Using High Order Finite Differences/Definitions and Basics

From testwiki
Jump to navigation Jump to search

Template:Status

Definitions and Basics

Vector Norms

Definition of a Vector Norm

The most ordinary kind of vectors are those consisting of ordered n-tuples of real or complex numbers. They may be written in row      <xyz>,      [xyz],       (x,y,z),     

or column      [xyz]

forms. Commas or other seperators of components or coordinates may or may not be used. When a vector has many elements, notations like

[x1x2xn]   or   [x1x2xn]T   are often used.

A most popular notation to indicate a vector is  v=<v1v2vn>.

Vectors are usually added component-wise, for

 v=<v1v2vn>  and   w=<w1w2wn>,

 v+w=<(v1+w1)(v2+w2)(vn+wn)>.

Scalar multiplication is defined by   αv=<αv1αv2αvn>.

A vector norm is a generalization of ordinary absolute value  |x|  of a real or complex number.

For  u  and  vvectors, and  α,  a scalar, a vector norm is a real value   associated with a vector for which the following properties hold.

(i)v0(ii)v=0v=0(iii)αv=|α|v(iv)v+wv+w.

Common Vector Norms

The most commonly used norms are:

v2=|v1|2+|v2|2++|vn|2v1=|v1|+|v2|++|vn|v=max|vi|1invp=(|v1|p+|v2|p++|vn|p)1p.

Any two norms on n dimensional vectors of complex numbers are topologically equivalent in the sense that, if  a  and  b  are two different norms, then there exist positive constants  c1  and  c2  such that  c1abc2a.

Inner Product of Vectors

The inner product (or dot product), of two vectors

 v=<v1v2vn>  and   w=<w1w2wn>,

is defined by

 vw=v1w1+v2w2++vnwn,

or when  v  and  w  are complex valued, by

Template:RefQuote


It is often indicated by any one of several notations:

vw,vTw,<v,w>,    or    (v,w).

Besides the dot product, other inner products are defined to be a rule that sssigns to each pair of vectors  v,w,  a complex number with the following properties.

<αv,w>=α<v,w>

<v,w>=<w,v>

<v1+v2,w>=<v1,w>+<v2,w>

for  v0<v,v>  is real valued and positive

<v,v>>0

and

|<v,w>|2<v,v><w,w>

An inner product defines a norm by

Template:RefQuote

Inequalities Involving Norms

The Cauchy Schwarz and Holder's inequalities are commonly employed.

|vw|v2w2

|vw|vpwq    for  1p+1q=1

Matrix Norms

Definition of a Matrix Norm

The most ordinary kind of matices are those consisting of rectangular arrays of real or complex numbers. They may be written in element form      A=(aij),1m,1n      and be considered as collections of column  or  row  vectors.

Matrices are usually added element-wise, for

A=(aij),B=(bij),A+B=(aij+bij)

Scalar multiplication is defined by   αA=(αaij).

The notation  A=0  means that all elements of  A  are identically zero.

A matrix norm is a generalization of ordinary absolute value  |x|  of a real or complex number, and can be considered a type of a vector norm.

For  A  and  Bmatrices, and  α,  a scalar, a matrix norm is a real value   associated with a matrix for which the following properties hold.

(i)A0(ii)A=0A=0(iii)αA=|α|A(iv)A+BA+B.

Common Matrix Norms

The most commonly used norms are:

AF=i=1mj=1n|aij|2A=max1im(|ai1|+|ai2|++|ain|)A1=max1jn(|a1j|+|a2j|++|amj|)Amax=max{|aij|}1im,1jnAp=(i=1mj=1n|aij|p)1pA2=maxx2=1Ax2.

Like vector norms any two matrix norms on m×n matrices of complex numbers are topologically equivalent in the sense that, if  a  and  b  are two different norms, then there exist positive constants  c1  and  c2  such that  c1abc2a.

The norm  A2  on matrices  A  is an example of an induced norm. An induced norm is for a vector norm xb  defined by

Ab=maxxb=1Axb.

This could cause the same subscript notation to be used fot two different norms, sometimes.

Positive definite matrices

n×n  matrix  A  is said to be positive definite if for any vector  x

xTAxαxTx

for some positive constant  α , not depending on  x .

The property of being positive definite insures the numerical stability of a variety of common numerical techniques used to solve the equation  Ax=y .

Taking  x=A1y  then

(A1y)TA(A1y)α(A1y)TA1y .

so that

(A1y)TyαA1y22     and     αA1y22A1y2y2  .

This gives

αA1y2y2

and

A12α1 .

Consistency of Norms

A matrix norm  M  and a vector norm  V  are said to be consistent when

AxMAMxV .

When  M  is the matrix norm induced by the vector norm  V  then the two norms will be consistent.

When the two norms are not consistent there will still be a positive constant  α  such that

AxMαAMxV .

Finite Difference Operators

approximation of a first derivative

The defition of the derivative of a function gives the first and simplest finite difference.

f(x0)=limx1x0f(x1)f(x0)x1x0orf(x0)=limh0f(x0+h)f(x0)h.

So the finite difference 

dhdhxf(x)=f(x+h)f(x)h=h1(f(x+h)f(x))

can be defined. It is an approximation to  f(x)  when  h is near 0.

The finite difference approximation  dhndhxnf(x)  for  f(n)(x)  is said to be of order  k,  if there exists  M>0  such that

|dhndhxnf(x)f(n)(x)|Mhk,

when  h is near 0.

For practical reasons the order of a finite difference will be described under the assumption that  f(x)  is sufficiently smooth so that it's Taylor's expansion up to some order exists. For example, if

f(x+h)=f(x)+f(x)h+12f(zh)h2

then

f(x+h)f(x)h=f(x)+12f(zh)h

so that

dhdhxf(x)f(x)=12f(zh)h,

meaning that the order of the approximation of  f(x)  by  dhdhxf(x)  is  1.

The finite difference so far defined is a 2-point operator, since it requires 2 evaluations of f(x). If

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(zh)h3

then another 2-point operator

dhdhxf(x)=f(x+h)f(xh)2h=(2h)1(f(x+h)f(xh))

can be defined. Since

f(x+h)f(xh)2h=f(x)+12(16f(zh)+16f(zh))h2,

this  dhdhxf(x)  is of order 2, and is referred to as a centered difference operator. Centered difference operators are usually one order of accuracy higher than un-centered operators for the same number of points.


More generally, for  m  points  x+α1h,x+α2h,,x+αmh  a finite difference operator

dhdhxf(x)=h1(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

is usually defined by choosing the coefficients  a1,a2,,am  so that  dhdhxf(x)  has as high of an order of accuracy as possible. Considering

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(x)h3++1(m1)!f(m1)(x)hm1+1m!f(m)(zh)hm.

Then

dhdhxf(x)=h1(c1f(x)+c2f(x)h+12c3f(x)h2+16c4f(x)h3++1(m1)!cmf(m1)(x)hm1+1m!Rmhm).

where the  c1,c2,,cm  are the right hand side of the Vandermonde system

[1111α1α2αm1αmα12α22αm12αm2α1m2α2m2αm1m2αmm2α1m1α2m1αm1m1αmm1][a1a2a3am1am]=[c1c2c3cm1cm]

and

Rm=a1α1mf(m)(zα1h)+a2α2mf(m)(zα2h)+a3α3mf(m)(zα3h)++am1αm1mf(m)(zαm1h)+amαmmf(m)(zαmh).

When the  ai's  are chosen so that 

c1=0,c2=1,c3==cm=0

then

dhdhxf(x)=f(x)+1m!Rmhm1.

so that the operator is of order  m1.

At the end of the section a table for the first several values of  m,  the number of points, will be provided. The discussion will move on to the approximation of the second derivative.

approximation of a second derivative

The definition of the second derivative of a function

f(x)=limh0f(x+h)f(x)h.

used together with the finite difference approximation for the first derivative

dhdhxf(x)=f(x+h)f(x)h=h1(f(x+h)f(x))

gives the finite difference 

dh2dhx2f(x)=h1(h1(f(x+2h)f(x+h))h1(f(x+h)f(x)))=h2(f(x+2h)2f(x+h)+f(x)))

In view of

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(zh)h3

for the operator just defined

dh2dhx2f(x)=f(x)+(43f(z2h)13f(zh))h.

If instead, the difference operator

dhdhxf(x)=(2h)1(f(x+h)f(xh))

is used

dh2dhx2f(x)=h1((2h)1(f(x+2h)f(x))(2h)1(f(x+h)f(xh)))=12h2(f(x+2h)f(x+h)f(x)+f(xh))=f(x)+(23f(z2h)112f(zh)112f(zh))h

If the other obvious possibility is tried

dh2dhx2f(x)=h1(h1(f(x+h)f(x))h1(f(x)f(xh)))=h2(f(x+h)2f(x)+f(xh))

In view of

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(x)h3+124f(iv)(zh)h4,

dh2dhx2f(x)=f(x)+(124f(iv)(zh)+124f(iv)(zh))h2.

So   dh2dhx2f(x)=h2(f(x+h)2f(x)+f(xh))

is a second order centered finite difference approximation for  f(x).


The reasoning applied to the approximation of a first derivative can be used for the second derivative with only a few modifications.

For  m  points  x+α1h,x+α2h,,x+αmh  a finite difference operator

dh2dhx2f(x)=h2(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

is usually defined by choosing the coefficients  a1,a2,,am  so that  dh2dhx2f(x)  has as high of an order of accuracy as possible. Considering

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(x)h3++1(m1)!f(m1)(x)hm1+1m!f(m)(zh)hm.

Then

dh2dhx2f(x)=h2(c1f(x)+c2f(x)h+12c3f(x)h2+16c4f(x)h3++1(m1)!cmf(m1)(x)hm1+1m!Rmhm).

where the  c1,c2,,cm  are the right hand side of the Vandermonde system

[1111α1α2αm1αmα12α22αm12αm2α1m2α2m2αm1m2αmm2α1m1α2m1αm1m1αmm1][a1a2a3am1am]=[c1c2c3cm1cm]

and

Rm=a1α1mf(m)(zα1h)+a2α2mf(m)(zα2h)+a3α3mf(m)(zα3h)++am1αm1mf(m)(zαm1h)+amαmmf(m)(zαmh).

When the  ai's  are chosen so that 

c1=0,c2=0,c3=2,c4==cm=0

then

dh2dhx2f(x)=f(x)+1m!Rmhm2.

so that the operator is of order  m2.

The effect of centering the points will be covered somewhere below.

approximation of higher derivatives

Although approximations to higher derivatives can be defined recursively from those for derivatives of lower order, the end result is the same finite difference operators. The Vandermonde type system will be used again for this purpose.

f(n)(x)=limh0f(n1)(x+h)f(n1)(x)h.

The number of points needed to approximate  f(n)(x)  by finite differences is at least  n+1.


For  m  points  x+α1h,x+α2h,,x+αmh  a finite difference operator

dhndhxnf(x)=hn(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

is usually defined by choosing the coefficients  a1,a2,,am  so that  dhndhxnf(x) approximates  f(n)(x)  to as high of an order of accuracy as possible. Considering

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(x)h3++1(m1)!f(m1)(x)hm1+1m!f(m)(zh)hm.

Then

dhndhxnf(x)=hn(c1f(x)+c2f(x)h+12c3f(x)h2+16c4f(x)h3++1n!cn+1f(n)(x)hn++1(m1)!cmf(m1)(x)hm1+1m!Rmhm).

where the  c1,c2,,cm  are the right hand side of the Vandermonde system

[1111α1α2αm1αmα12α22αm12αm2α1m2α2m2αm1m2αmm2α1m1α2m1αm1m1αmm1][a1a2a3am1am]=[c1c2c3cm1cm]

and

Rm=a1α1mf(m)(zα1h)+a2α2mf(m)(zα2h)+a3α3mf(m)(zα3h)++am1αm1mf(m)(zαm1h)+amαmmf(m)(zαmh).

When the  ai's  are chosen so that 

cn+1=n!,andci=0,forin+1

then

dhndhxnf(x)=f(n)(x)+1m!Rmhmn.

so that the operator is of order  mn.

An alternative analysis is to require that the finite difference operator differentiates powers of  x  exactly, up to the highest power possible.

effect of the placement of points

Usually the  αi's  are taken to be integer valued, since the points are intended to coincide with those of some division of an interval or 2 or 3 dimensional domain. If these points and hence  αi's  are chosen with only accuracy in mind, then a higher accuracy of only one order can be achieved.

So start by seeing how high is the accuracy that  f(x)  can be approximated with three points.

dhdhxf(x)=h1(a1f(x+α1h)+a2f(x+α2h)+a3f(x+α3h))

Then accuracy of order 4 can not be achieved, because it would require the solution of

[111α1α2α3α12α22α32α13α23α33α14α24α34][a1a2a3]=[01000]

which can not be solved since the matrix

[α12α22α32α13α23α33α14α24α34]

is non-singular. The possibility of an  αi  being 0  can be ruled out otherwise.

For accuracy of order 3

[111α1α2α3α12α22α32α13α23α33][a1a2a3]=[0100]

So the matrix

[111α12α22α32α13α23α33]

is singular and  α1,α2,α3  are the roots of some polynomial 

p(α)=α3+b1α2+b0 .

Two examples are next.

dhdhxf(x)=h1((33+12)f(x+(331)h)+233f(x+33h)(3312)f(x+(33+1)h))

dhdhxf(x)=h1(92f(x12h)+163f(x+34h)56f(x+32h))


To see what the accuracy that  f(x)  can be approximated to with three points.

dh2dhx2f(x)=h2(a1f(x+α1h)+a2f(x+α2h)+a3f(x+α3h))

Then accuracy of order 3 can not be achieved, because it would require the solution of

[111α1α2α3α12α22α32α13α23α33α14α24α34][a1a2a3]=[00200]

which can not be solved since the matrices

[111α1α2α3α13α23α33]      and      [111α1α2α3α14α24α34]

would both need to be singular.

If the matrix

[111α1α2α3α13α23α33]

is singular, then  α1,α2,α3  are the roots of some polynomial 

p(α)=α3+b1α+b0 ,

implying

[α14α24α34]=b1[α12α22α32]b0[α1α2α3]

meaning that elementary row operations can transform

[111α1α2α3α14α24α34]

to

[111α1α2α3α12α22α32]

which is non-singular.

Conversely, if  α1,α2,α3  are the roots of some polynomial 

p(α)=α3+b1α+b0 , then

[111α1α2α3α12α22α32α13α23α33][a1a2a3]=[0020]

can be solved and  f(x)  approximated to an order of 2 accuracy.

See how high is the accuracy that  f(x)  can be approximated with  m  points.

dhdhxf(x)=h1(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

Then accuracy of order  m+1  can not be achieved, because it would require the solution of

[1111α1α2α3αmα12α22α32αm2α13α23α33αm3α1mα2mα3mαmmα1m+1α2m+1α3m+1αmm+1][a1a2a3am]=[0100]

which can not be solved since the matrix

[α12α22α32αm2α13α23α33αm3α1mα2mα3mαmmα1m+1α2m+1α3m+1αmm+1]

is non-singular. The possibility of an  αi  being 0  can be ruled out otherwise, because, for example, if  α1=0,  then the non-singularity of the block

[α22α32αm2α23α33αm3α2mα3mαmm]

would force  a2=a3=am=0.

For accuracy of order  m

[1111α1α2α3αmα12α22α32αm2α13α23α33αm3α1mα2mα3mαmm][a1a2a3am]=[0100]

So the matrix

[1111α12α22α32αm2α13α23α33αm3α1mα2mα3mαmm]

is singular and  α1,α2,α3,,αm  are the roots of some polynomial 

p(α)=αm+bm2αm1++b2α3+b1α2+b0 .

The progression for the second, third, ... derivatives goes as follows.

If  α1,α2,α3,,αm  are the roots of some polynomial 

p(α)=αm+bm2αm1++b2α3+b1α+b0 

then the system

[1111α1α2α3αmα12α22α32αm2α13α23α33αm3α1mα2mα3mαmm][a1a2a3a4am]=[00200]

can be solved, and

dh2dhx2f(x)=h2(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

approximates  f(x)  to an order of accuracy of  m1.

If  α1,α2,α3,,αm  are the roots of some polynomial 

p(α)=αm+bm2αm1++b3α4+b2α2+b1α+b0 

then the system

[1111α1α2α3αmα12α22α32αm2α13α23α33αm3α14α24α34αm4α1mα2mα3mαmm][a1a2a3a4a5am]=[000600]

can be solved, and

dh3dhx3f(x)=h3(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

approximates  f(x)  to an order of accuracy of  m2.

Now, the analysis is not quite done yet. Returning to the approximation of  f(x).  If for the polynomial

p(α)=αm+bm2αm1++b2α3+b1α+b0

it were that  b1=0,  then the system can be solved for one more order of accuracy. So the question arises as to whether polynomials of the form

p(α)=αm+bm2αm1++b2α3+b0

exist that have  m  distinct real roots. When  m=3  there is not. So consider  m=4.

p(α)=α4+b2α3+b0

If  p(α)  has 4 distinct real roots, then

p(α)=4α3+3b2α2

has 3 distinct real roots, which it does not. So the order of approximation can not be improved. This is generally the case.

Returning to the approximation of  f(x).  If for the polynomial

p(α)=αm+bm2αm1++b3α4+b2α2+b1α+b0

it were that  b2=0,  then the system can be solved for one more order of accuracy. So the question arises as to whether polynomials of the form

p(α)=αm+bm2αm1++b3α4+b1α+b0

exist that have  m  distinct real roots.

If  p(α)  has  m  distinct real roots, then

p(α)=m(m1)αm2+(m1)(m2)bm2αm3++12b3α2

has  m2  distinct real roots, which it does not. So the order of approximation can not be improved.

For functions of a complex variable using roots of unity, for example, may obtain higher orders of approximation, since complex roots are allowed.

centered difference operators

For  m  points  x+α1h,x+α2h,,x+αmh  a finite difference operator

dhndhxnf(x)=hn(a1f(x+α1h)+a2f(x+α2h)++amf(x+αmh))

is said to be centered when the points are symmetrically placed about  x.

αi=αmi+1fori=1,2,,[m/2]

When  m  is odd  α[m/2]+1=0.

To find the centered difference operators, consider

f(x+h)=f(x)+f(x)h+12f(x)h2+16f(x)h3++1m!f(m)(x)hm+1(m+1)!f(m+1)(zh)hm+1.

Then

dhndhxnf(x)=hn(c1f(x)+c2f(x)h+12c3f(x)h2+16c4f(x)h3++1(m1)!cmf(m1)(x)hm1+1m!cm+1f(m)(x)hm+1(m+1)!Rm+1hm+1).

where the  c1,c2,,cm,cm+1  are the right hand side of the over-determined system

[1111α1α2αm1αmα12α22αm12αm2α1m2α2m2αm1m2αmm2α1m1α2m1αm1m1αmm1α1mα2mαm1mαmm][a1a2a3am1am]=[c1c2c3cm1cmcm+1]

and

Rm+1=a1α1m+1f(m+1)(zα1h)+a2α2m+1f(m+1)(zα2h)+a3α3m+1f(m+1)(zα3h)++am1αm1m+1f(m+1)(zαm1h)+amαmm+1f(m+1)(zαmh).

When the  ai's  are chosen so that 

cn+1=n!,ci=0,forin+1

then

dhndhxnf(x)=f(n)(x)+1(m+1)!Rm+1hmn+1.

so that the operator is of order  mn+1.

Since, for the centered case, the system is over-determined, some restriction is needed for the system to have a solution. A solution occurs when the  α1,α2,,αm  are the roots of a polynomial

p(α)=αm+bm1αm1++b2α2+b1α+b0

with

bn=0.

Observing that when  m  is even

p(α)=i=1m2(α2αi2)


and when  m  is odd

p(α)=αi=1[m2](α2αi2).

So when  m  is even  p(α)  has  bn=0  for all odd  n and when  m  is odd  p(α)  has  bn=0  for all even  n.

So a centered difference operator will achieve the one order extra of accuracy when the number of points  m  is even and the order of the derivative  n  is odd or when the number of points  m  is odd and the order of the derivative  n  is even.


Shift Operators

Shift of a Function

Trigometric polynomials

Let 

Template:RefQuote


be trigometric polynomial defined on  1x1.

Define the inner product on such trigometric polynomials by

Template:RefQuote

In light of the orthogonalities

<sin(kπx),sin(rπx)>=0,<cos(kπx),cos(rπx)>=0,

and  <sin(kπx),cos(rπx)>=0,      when  kr,

inner products can be calculated easily.

Template:RefQuote

and for

p1(x)=a0,1+k=1m(ak,1cos(kπx)+bk,1sin(kπx))

p2(x)=a0,2+k=1m(ak,2cos(kπx)+bk,2sin(kπx))

the inner product is given by

Template:RefQuote

Definition of a shift operator

Define the shift operator  (sft)h  on  p(x)  by


Template:RefQuote


Since

sin(kπ(xh))=cos(kπh)sin(kπx)sin(kπh)cos(kπx)

and

cos(kπ(xh))=cos(kπh)cos(kπx)+sin(kπh)sin(kπx),

so that

Template:RefQuote

Approximation by trigometric polynomials

Let  f(x)  be a function defined on and periodic with respect to the interval  1x1.  That is  f(x+2)=f(x).

The  mth  degree trigometric polynomial approximation to  f(x)  is given by

p(x)=a0+k=1m(akcos(kπx)+bksin(kπx))

where

Template:RefQuote

p(x)  approximates  f(x)  in the sense that

11|f(x)p(x)|2dx

is minimized over all trigometric polynomials, of degree  m  or less, by  p(x).

In fact

11|f(x)p(x)|2dx=11(f(x)p(x))(f(x)p(x))dx

=11(|f(x)|22(f(x)p(x))+|p(x)|2)dx.

The term in the center

11f(x)p(x)dx =a011f(x)dx+k=1m(ak11f(x)cos(kπx)dx+bk11f(x)sin(kπx))dx

=|a0|2+k=1m(|ak|2+|bk|2)=11|p(x)|2dx,

so that

Template:RefQuote

Template:RefQuote

If  p(x)  and  q(x)  are the  mth  degree trigometric polynomial approximations to  f(x)  and  g(x),  then the  mth  degree trigometric polynomial approximation to  f(x)+αg(x)  is given by  p(x)+αq(x).

This follows immediately from (3.4.0) since

11(f(x)+αg(x))cos(kπx)dx=11f(x)cos(kπx)dx+α11g(x)cos(kπx)dx      and      11(f(x)+αg(x))sin(kπx)dx=11f(x)sin(kπx)dx+α11g(x)sin(kπx)dx.

Fundamental Property

Template:RefQuote

Generally if  p(x)  is the  mth  degree trigometric polynomial approximation to a function  f(x)periodic on  1x1,  then  (sft)hp(x)  is the  mth  degree trigometric polynomial approximation to  f(xh).

To see this calculate the trigometric polynomial approximations for  f(xh).

11f(xh)cos(kπx)dx=1h1hf(x)cos(kπ(x+h))dx

=1h1hf(x)(cos(kπh)cos(kπx)sin(kπh)sin(kπx))dx

=cos(kπh)1h1hf(x)cos(kπx)dxsin(kπh)1h1hf(x)sin(kπx)dx

=cos(kπh)11f(x)cos(kπx)dxsin(kπh)11f(x)sin(kπx)dx

=akcos(kπh)bksin(kπh),

where

ak=11f(x)cos(kπx)dx      and      bk=11f(x)sin(kπx)dx.

11f(xh)sin(kπx)dx=1h1hf(x)sin(kπ(x+h))dx

=1h1hf(x)(cos(kπh)sin(kπx)+sin(kπh)cos(kπx))dx

=cos(kπh)1h1hf(x)sin(kπx)dx+sin(kπh)1h1hf(x)cos(kπx)dx

=cos(kπh)11f(x)sin(kπx)dx+sin(kπh)11f(x)cos(kπx)dx

=aksin(kπh)+bkcos(kπh).

Comparing the results with (3.3.1) finishes the observation.

Another detail of use is

11|f(xh)(sft)hp(x)|2dx

=1h1h|f(x)p(x)|2dx=11|f(x)p(x)|2dx

which is

Template:RefQuote

Error Estimation

A result used in error estimation is

Template:RefQuote

When  s(x)  is a sine polynomial

s(x)=k=1mbksin(kπx)

then

Template:RefQuote

and

Template:RefQuote

Simple Functions

Let  0=x0<x1<x2<<xn1<xn=1

so that

{[x0,x1],(x2,x2],,(xn2,xn1],(xn1,xn]}

is a partition of  0x1.

The function  f(x)  is said to be simple on  0x1,  if

Template:RefQuote

Of particular interest is when the points have equal spacing  xixi1=1n.

The intent is to make estimates of  01|f(x)f(xh)|2dx.

Begin by making an odd extension of  f(x)  to  1x1 by setting  f(x)=f(x)  and continue the definition by extending  f(x)  periodically.

Then approximate  f(x)  with a sine polynomial

s(x)=k=1mbksin(kπx)

where

bk=201f(x)sin(kπx)dx.

When  m  is large enough so that some  k  are divided by  2n then, for  k=2nr

bk=2i=1nxi1xif(x)sin(kπx)dx=2i=1nfi1(i1)1ni1nsin(kπx)dx,

and letting  x=1nru,

(i1)1ni1nsin(kπx)dx=(i1)rirsin(2πu)du=0 ,

so that

bk=b2nr=0 .

Now, return to the sum

s(x)(sft)hs(x)2=k=1m2|bk|2(1cos(kπh))

with  b2nr=0  for  r=1,2,3,.

If  gcd(j,2n)=1  and  (2n)k,  then  (2n)kj,  and in this case

cos(kπj1n)0  and  (1cos(kπj1n))(1cos(π1n)).

So for  h=j1n   with  gcd(j,2n)=1 ,

s(x)(sft)hs(x)2k=1m2|bk|2(1cos(π1n))

Template:RefQuote

Next observe that if  s(x)  is the  mth  degree sine polynomial approximation to  f(x)  then  (sft)hs(x)  is the  mth  degree trigometric polynomial approximation to  f(xh).  The assumption that  f(x)  is odd and periodic is still in effect.

Finally the intended results follow.

s(x)(sft)hs(x)

=(f(x)f(xh))+(s(x)f(x))((sft)hs(x)f(xh))

f(x)f(xh)+f(x)s(x)+f(xh)(sft)hs(x)

=f(x)f(xh)+2f(x)s(x)

so that

f(x)f(xh)s(x)(sft)hs(x)2f(x)s(x).

Making use of (3.7.1)

f(x)f(xh)2(1cos(π1n))s(x)2f(x)s(x).

Being that simple functions can be approximated by trigometric polynomials to arbitrary accuracy,

Template:RefQuote

Now, for  h=1n,  and using the definition of the simple function  f(x)

f(x)2=2ni=0n1|fi|2

To find the sum for  f(x)f(xh)2  list tthe values of  f(x)  over the values of  f(xh)  on the whole interval  1x1.

fn1fn2f1f0f0f1f2fn2fn1fn1fn1f2f1f0f0f1fn3fn2

This gives

f(x)f(xh)2=2n(2|f0|2+i=1n1|fifi1|2+2|fn1|2).

So the inequality follows

Template:RefQuote


Specialized Norms

Energy and Heat Norms

There are a number of less well known, but important norms. These norms are important in the analysis of many physical problems and are used in error estimation for finite difference and finite element methods. Examples are the energy and heat norms.

These norms are usually expressed in an integral form.

yE=ab(12ky2(t)+12m(y)2(t))dt

yH=ab(y2(t)+(y)2(t)+(y)2(t)++(y(k))2(t))dt

When   y(a)=y(b)=0   the inequality below holds.

ab(y)2(t)dt43(ba)aby2(t)dt

This follows from a completely elementary, but lengthy calculation, as shown in appendix a). When additional assumptions on  y  are made this inequality can be improved somewhat.

ab(y)2(t)dtπ(ba)aby2(t)dt

See appendix b) for an explanation.

Discrete Heat Norms

In the analysis of finite difference methods for partial differential equations it is useful to have discrete analogs of norms like the energy and heat norms.

In an attempt not to have notation too cumbersome some indices are suppressed when from the discussion it is clear what they refer to. When  v=<v1v2vn>  is a  n  dimensional vector of complex numbers, finite difference operators are defined as discrete approximations or analogs of derivatives.

The most appropriate definition of a discrete energy or heat norm may vary due to differences in the handling of initial or boudary conditions. So for this reason the reader should make appropriate adjustments, when needed, to apply the same reasoning to another problem.

Before the discrete versions of energy or heat norms can be defined,  finite differece operations need to be defined and explained. This was done in the section on finite difference operators.


The next discrete version of the preceding inequality has important applications to the estimation of the error when approximating a second derivative with a finite difference operator.

If  v0=vn+1=0  then

i=1n+1(vivi1)2βnn+1i=1nvi2

with

β1=22,β2=3,

and generally the following under-estimate holds.

βn2n+1n+3

See appendix c) for a proof.

The inequality can be improved for general  n  by using (3.7.3) with  n  increased by 2.

2|f0|2+i=1n1|fifi1|2+2|fn1|22(1cos(π1n))i=0n1|fi|2

If  v0=vn+1=0  then

i=1n+1(vivi1)22(1cos(π1n+2))i=1nvi2

and

2(1cos(π1n+2))πn+2.


appendix a)

When   y(a)=y(b)=0   the inequality below holds.

ab(y)2(t)dt43(ba)aby2(t)dt

First apply the Cauchy Schwarz inequality.

12y2(x)=axy(t)y(t)dtaxy2(t)dtax(y)2(t)dt,

Next observe the integrals on the right are increasing with  x.

12y2(u)axy2(t)dtax(y)2(t)dt,foraux

Integrate, make cancellations, and reapply the first inequality.

12axy2(t)dt(xa)axy2(t)dtax(y)2(t)dt .

axy2(t)dt2(xa)ax(y)2(t)dt .

y2(x)4(xa)ax(y)2(t)dt .

After integrating again the inequality is improved.

axy2(t)dt4ax(ta)dtax(y)2(t)dt=2(xa)2ax(y)2(t)dt .

axy2(t)dt2(xa)ax(y)2(t)dt .

Now, assume the inequality immediately below holds for some  α.

axy2(t)dtα(xa)ax(y)2(t)dt

After substituting the inequality above into the next

12axy2(t)dt(xa)axy2(t)dtax(y)2(t)dt .

the following observations are made.

axy2(t)dt2α(xa)2ax(y)2(t)dt .

axy2(t)dt2α(xa)ax(y)2(t)dt .

y2(x)22α(xa)ax(y)2(t)dt .

axy2(t)dt22αax(ta)dtax(y)2(t)dt=2α(xa)2ax(y)2(t)dt .

axy2(t)dt2α4(xa)ax(y)2(t)dt .

Repeating this iteration leads to a sequence  αn+1=2αn4

which converges to  α=23,  the solution of  α=2α4.

So:

axy2(t)dt23(xa)ax(y)2(t)dt

and

y2(x)2223(xa)ax(y)2(t)dt .

For  c=(a+b)/2 ,

acy2(t)dt2223ac(ta)dtac(y)2(t)dt=223(ca)2ac(y)2(t)dt=2234(ba)2ac(y)2(t)dt .

acy2(t)dt434(ba)2ac(y)2(t)dt .

To handle the part  cby2(t)dt

12y2(x)=xby(t)y(t)dtxby2(t)dtxb(y)2(t)dt,

12y2(u)xby2(t)dtxb(y)2(t)dt,forxub

12xby2(t)dt(bx)xby2(t)dtxb(y)2(t)dt .

xby2(t)dt2(bx)xb(y)2(t)dt .

y2(x)4(bx)xb(y)2(t)dt .

Reasoning as before this inequality can be strengthend to

y2(x)2223(bx)xb(y)2(t)dt .

cby2(t)dt2223cb(bt)dtcb(y)2(t)dt=223(bc)2cb(y)2(t)dt=2234(ba)2cb(y)2(t)dt .

cby2(t)dt434(ba)2cb(y)2(t)dt .

Adding the results of the two main calculations

aby2(t)dt=acy2(t)dt+cby2(t)dt434(ba)2ac(y)2(t)dt+434(ba)2cb(y)2(t)dt=434(ba)2ab(y)2(t)dt .

appendix b)

When   y(a)=y(b)=0   the inequality below holds.

ab(y)2(t)dtπ(ba)aby2(t)dt


Assume the conditions that  y(x)  can be approximated by a trigometric polynomial

s(x)=a1sin(π(ba)x)+a2sin(π(ba)2x)+a3sin(π(ba)3x)++ansin(π(ba)nx).

such that  y(x)  is also approximated by the polynomials derivative

s(x)=a1π(ba)cos(π(ba)x)+2a2π(ba)cos(π(ba)2x)+3a3π(ba)cos(π(ba)3x)++nanπ(ba)cos(π(ba)nx)

That is to say, given  ϵ>0,  there exixts  a trigometric sine polynomial  s(x),  such that

ab|f(x)s(x)|2dx<ϵ   and   ab|f(x)s(x)|2dx<ϵ.

Now,

2baab|s(x)|2dx=|a1|2+|a2|2+|a3|2++|an|2

  and  

2baab|s(x)|2dx=π2(ba)2|a1|2+22π2(ba)2|a2|2+32π2(ba)2|a3|2++n2π2(ba)2|an|2.

So it is easy to see that

ab|s(x)|2dxπ2(ba)2ab|s(x)|2dx.

In this case the inequality is sharp, since it holds with equality for  f(x)=sin(π(ba)x).

Since

ab|s(x)|2dxϵ<ab|f(x)|2dx<ab|s(x)|2dx+ϵ,

and  

ab|s(x)|2dxϵ<ab|f(x)|2dx<ab|s(x)|2dx+ϵ,

the inequality holds for  y(x)  and  y(x).

appendix c)

If  v0=vn+1=0  then

i=1n+1(vivi1)2βnn+1i=1nvi2

with

β1=22,β2=3,

and generally the following under-estimate holds.

βn2n+1n+3

The cases for 1 and 2 are simple.

(v1v0)2+(v2v1)2=2v12

and

(v1v0)2+(v2v1)2+(v2v1)2=v12+(v2v1)2+v22v12+v22

with equality, when  v1=v2.

To prove the general under-estimate do as follows.

Use  v0=0  and apply the Cauchy Schwartz inequality together with inequality ().

vk2=i=1k(vi2vi12)=i=1k(vi+vi1)(vivi1)

i=1k(vi+vi1)2i=1k(vivi1)2

2i=1kvi2i=1k(vivi1)2

(1)

So for  ik,

vi22j=1ivj2j=1i(vjvj1)2

and

i=1kvi22ki=1k(vivi1)2.

After inserting the inequality above into (1)

vk24ki=1k(vivi1)2.

So for  ik,

vi24ij=1i(vjvj1)2

i=1kvi24(i=1ki)j=1k(vjvj1)2

and after using formula ()

i=1kvi22k(k+1)i=1k(vivi1)2.

Using  vn+1=0  the right end of the sum is estimated using nearly the same procedure.

vk2=i=kn(vi2vi+12)=i=kn(vi+vi+1)(vivi+1)

i=k+1n+1(vi+vi1)2i=k+1n+1(vivi1)2

2i=knvi2i=k+1n+1(vivi1)2

So for  ik,

vi22j=invj2j=i+1n+1(vjvj1)2,

i=knvi22(nk+1)i=k+1n+1(vivi1)2,

and

vk24(nk+1)i=k+1n+1(vivi1)2.

So for  ik,

vi24(ni+1)j=i+1n+1(vjvj1)2,

i=knvi24((n+1)(nk+1)(i=kni))i=k+1n+1(vivi1)2

=2(2(n+1)(nk+1)(n(n+1)(k1)k))i=k+1n+1(vivi1)2

=2((n2k+3)n+(k1)(k2))i=k+1n+1(vivi1)2.

Combining the two results:

i=1nvi2=i=1kvi2+i=k+1nvi2

2k(k+1)i=1k(vivi1)2

+2((n2k+1)n+k(k1))i=k+2n+1(vivi1)2.

When  n  is odd, use  k=(n+1)/2  so that the inequality becomes

2n+12n+32i=1k(vivi1)2

+2(n+12n12)i=k+2n+1(vivi1)2

(n+1)(n+3)2i=1n+1(vivi1)2.

When  n  is even, use  k=n/2  so that the inequality becomes

2n2n+22i=1k(vivi1)2

+2(n2n+22)i=k+2n+1(vivi1)2

(n)(n+2)2i=1n+1(vivi1)2.

Special Formulas

The rule named L´Hospital's Rule named by a French math teacher who published it without proof, because he said you need to be taken to the hospital if you can't prove it yourself, states.

If  limxaf(x)=0  and  limxag(x)=0  then,

limxaf(x)g(x)=limxaf(x)g(x).

The rule named Leibniz's Rule states how to differentiate a product of two functions  n  times.

(uv)(n)=k=0n(nk)u(k)v(nk)

There are a number of ways to find the formulas for the sums of the powers of the first  n  integers. One of the most convenient is as follows. Begin with the geometric formula.

1xn+11x=(1xn+1)(1x)1=k=0nxk

By convention  x0=1.

Apply Leibniz's Rule rule to the right hand side of the equivalent identity

(1xn+1)=(1x)k=0nxk

to differentiate both sides  m  times.

(n+1)n(nm+2)xnm+1=j=0m(mj)(1x)(j)(k=0nxk)(mj)

Upon noticing that only the first two terms of the sum on the right are non-zero

(n+1)n(nm+2)xnm+1=(1x)(k=0nxk)(m)m(k=0nxk)(m1).

Now take the limit as  x1  to establish the formula

(n+1)n(nm+2)=limx1m(k=0nxk)(m1).

From term by term differentiation.

(k=0nxk)=k=1nkxk1

and

(n+1)n=2k=1nk.

More generally

(k=0nxk)(j)=k=jnk(k1)(kj+1)xkj.

and thus

(n+1)n(nm+2)=mk=m1nk(k1)(km+2).

A progression of formulas follow by setting  m=3,4,.

k=1nk=n(n+1)2

(n+1)n(n1)=3k=2nk(k1)=3k=1n1(k+1)k.

=3k=1n1k2+3k=1n1k.

and

(n+2)(n+1)n=3k=1nk2+3k=1nk.

k=1nk2=n(n+1)(n+2)3n(n+1)2

After adjusting the index of summation and increasing  n  by  m2  the progression has the form

(n+m1)(n+m2)(n+1)n=mk=1n(k+m2)(k+m3)(k+1)k.

Letting  m=4

(n+3)(n+2)(n+1)n=4(k=1nk3+3k=1nk2+2k=1nk).

k=1nk3=n(n+1)(n+2)(n+3)4n(n+1)(n+2)+n(n+1)2


The formula called  summation by parts  is a discrete analog of  integration by parts . It is verified simply by inspection.

i=1mwi(vivi1)=wmvmw0v0i=0m1vi(wi+1wi)

One form of the statement of the Fundamental Theorem of Integral Calculus is

For  F(x)=axf(t)dt  it holds  F(x)=f(x).

Another rule which slightly generalizes the Fundamental Theorem of Integral Calculus says

ddxaxf(x,t)dt=f(x,x)+axxf(x,t)dt.


The following simple inequality can be used to bound sums of products.

(ab)2=a22ab+b202aba2+b2ab12(a2+b2)

and also

ab=(αa)(b/α)12(α2a2+b2/α2),

so that

i=1naibi12i=1n(αi2ai2+bi2/αi2).

Letting  bi=b  and  αi2=n 

i=1naib12n((i=1nai2)+b2).

The particular case with  n=4  gives

(c1+c2+c3+c4)c0(c02+c12+c22+c32+c42).

save scratch



i=1kvi22i=1kj=1ivj2j=1i(vjvj1)2

2i=1kj=1ivj2i=1kj=1i(vjvj1)2

=2i=1k(ki+1)vi2i=1k(ki+1)(vivi1)2

i=1kvi22ki=1kvi2i=1k(vivi1)2

Now, suppose for some  α

i=1kvi2αki=1k(vivi1)2

vk22αki=1k(vivi1)2

So for  ik

vi22αij=1i(vjvj1)2

i=1kvi22α(i=1ki)j=1k(vjvj1)2

i=1kvi2αk(k+1)i=1k(vivi1)2

i=1kvi2α(1+1k)ki=1k(vivi1)2


i=1kvi22ki=1k(ki+1)(vivi1)2 zzzzzzzzzz

So for  ik

j=1ivj22j=1i(ij+1)vj2j=1i(ij+1)(vjvj1)2

i=1kj=1ivj22i=1kj=1i(ij+1)vj2j=1i(ij+1)(vjvj1)2

i=1k(ki+1)vi22ki=1k(ki+1)vi2i=1k(ki+1)(vivi1)2

i=1k(ki+1)vi22ki=1k(ki+1)(vivi1)2

i=1kvi24ki=1k(ki+1)(vivi1)2

i=1kvi22ki=1k(ki+1)(vivi1)2

xxxxxxxxxxxxxxxx i=1k(ki+1)vi22i=1kj=1i(ij+1)vj2i=1kj=1i(ij+1)(vjvj1)2

i=1k(ki+1)vi22i=1kj=1i(ij+1)vj2i=1kj=1i(ij+1)(vjvj1)2


2ki=1k(ki+1)vi2i=1k(vivi1)2

i=1kvi22ki=1kvi2i=1k(vivi1)2

i=1kvi22ki=1k(vivi1)2

Template:BookCat