Using High Order Finite Differences/Preliminary Estimates

From testwiki
Jump to navigation Jump to search

Introduction

The One-dimensional Problem

Continuous Problem

This book is primarily concerned with finite difference approximations to the solutions of partial differential equations. However, it is both useful in itself and instructive to study the solution of some one-dimensional problems by finite differences.

Consider the meager problem to solve

y(x)=f(x)    on  axb

subject to the boundary conditions

y(a)=y(b)=0.

The following inequality derived next motivates the analysis of the higher-dimensional problems involving the Laplacian operator.

aby(x)y(x)dx=aby(x)y(x)+ab|y(x)|2dx=ab|y(x)|2dx

Making use of inequality () that

ab|y(x)|2dxπ2(ba)2ab|y(x)|2dx

and using the equation together with the Cauchy Schwartz inequality

aby(x)y(x)dx=aby(x)f(x)dxab|y(x)|2dxab|f(x)|2dx.

It follows

π2(ba)2ab|y(x)|2dxab|y(x)|2dxab|f(x)|2dx.

and that

ab|y(x)|2dx(ba)2π2ab|f(x)|2dx.

Now, consider when  z(x)  solves the approximate problem

z(x)=g(x)    on  axb

subject to the boundary conditions  z(a)=z(b)=0,  with  g(x)  near  f(x).
Since  (y(x)z(x))=(f(x)g(x)) ,


ab|y(x)z(x)|2dx(ba)2π2ab|f(x)g(x)|2dx.

This inequality just above will be generalized to both discrete and higher dimensional analogs. This will enable the analysis of the accuracy of finite difference methods for several kinds of equations.

Another point to make is, if non-zero boundary conditions

y(a)=A,y(b)=B,

are wanted, then

y(x)+(BA)(ba)(xa)+A

will solve the new problem. So there will not be a loss of generality by assuming zero boundary conditions, in most cases. In the higher dimensional case this will have the effect of separating the problem into two parts.

Solution by Finite Differences

Return to () the problem to solve from the previous section.

y(x)=f(x)    on  axb

subject to the boundary conditions

y(a)=y(b)=0.

This problem may be better by means other than finite differences. For example,

letting  g(x)=ax(xt)f(t)dt,  and applying the rule (),  g(x)=(xx)f(x)+axx(xt)f(t)dt=axf(t)dt  so that   g(x)=f(x).


Then   y(x)=g(x)g(b)ba(xa)   is a solution of the problem.

However, it will illustrate the finite difference method without being to difficult to follow, for a start. Partial differential equations can on the other hand be difficult to solve by direct analytic methods.

Begin with a partition of the interval given by

a=x0<x1<x2<<xn<xn+1=b.

For simplicities sake, assume a uniform mesh  h=(ba)/(n+1).

That is to say  xixi1=h    for   i=1,2,,n+1.

For the beginning approximate the second derivative by the second order accurate difference operator

dh2dhx2y(x)=h2(y(xh)2y(x)+y(x+h)).

Define  fi=f(xi).  It will be shown that there exist unique

y0,y1,y2,,,yn,yn+1   with   y0=yn+1=0

that solve the equations

h2(yi12yi+yi+1)=fi   for   i=1,2,,n.

In addition, and most importantly, the  yi  approximate  y(x)  in the sense

i=1n|y(xi)yi|2nBh2

for some bound  B  independent of  h.

At this point some clarification must be made and some notation introduced. The term finite difference operator is used for two different, but related operators. One is the difference quotient applied to the function  y(x),  namely

dh2dhx2y(x)=h2(y(xh)2y(x)+y(x+h))

and the other is the linear operator

h2(yi12yi+yi+1)

applied to the vector  y=<y0y1y2ynyn+1>.

The notation

(dh2dhx2)iy=h2(yi12yi+yi+1)

will be used to distinguish between the two. The equations () can be written as

(dh2dhx2)iy=fi.

The linear operator   (dh2dhx2)I  is defined by

(dh2dhx2)Iy=<(dh2dhx2)1y(dh2dhx2)2y(dh2dhx2)ny>.

Then   (dh2dhx2)I  can be thought of as

(dh2dhx2)I=<(dh2dhx2)1(dh2dhx2)2(dh2dhx2)n>.

This representation of a linear operator differs from the matrix notation usually used. This has the advantage of allowing estimates to be generalized more easily to higher order operators and two or three dimensional domains.

In fact this linear operator is well studied and has the matrix representation

[21000000121000000000012100000012],

when applied to the vector  <y1y2yn>T.

The eigenvalues, eigenvectors, and inverse are known for this matrix, and can be found in some references. If it were just for the sake of this one introductory example the details of this matrix would be used for the analysis. As is being pointed out a method of analysis that generalizes to higher order operators and domains is being developed instead.

Finally the equations () can be written as

(dh2dhx2)Iy=f,

where the vector  f=<f1f2fn>.

One last piece of notation, for the vector  y=<y0y1y2ynyn+1>

the interior points of  y,  denoted by  yo,  is the  n  dimensional vector

yo=<y1y2yn>.

Restating the problem, it will be shown that there exist unique

y=<y0y1y2ynyn+1>   with   y0=yn+1=0

that solve the equation ()

(dh2dhx2)Iy=f,

In addition, and most importantly, the  yi  approximate  y(x)  in the sense

i=1n|y(xi)yi|2nBh2

where the bound  B  is independent of  h  and if fact a good estimate of  B  is

B=M(αn)1h2

where   M=112maxaxb(f(iv)(x)),

αn=2(1cos(π1n+2))π2(n+2)2,

and   BMπ2.

The remainder of this section is a proof of the claim () immediately above.

The proof is done by first showing that the operator   (dh2dhx2)I   is positive definite.

In particular for  αn=2(1cos(π1n+2))π2(n+2)2

(yo)T(dh2dhx2)Iyαnh2yo22.

To prove what is said immediately above

h2(yo)T(dh2dhx2)Iy=i=1nyi(yi1+2yiyi+1)=i=1nyi((yi+1yi)(yiyi1)).

Rearrange  summation by parts  ()

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

as

i=0m1vi(wi+1wi)=wmvm+w0v0+i=1mwi(vivi1)

Now, let  m=n+1  to get

i=0nvi(wi+1wi)=wn+1vn+1+w0v0+i=1n+1wi(vivi1)

Set   vi=yi   for   i=0,1,,n+1.

Set   w0=y0   and   wi=yiyi1   for   i=1,2,,n+1.

Since   v0=vn+1=0,   the identity becomes

i=1nvi(wi+1wi)=i=1n+1wi(vivi1),

which is then in turn

i=1nyi((yi+1yi)(yiyi1))=i=1n+1(yiyi1)(yiyi1).

This has proven the equality

h2(yo)T(dh2dhx2)Iy=i=1n+1(yiyi1)2.

Making use of ()

If  v0=vn+1=0  then

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

and

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

We have the following

i=1n+1(yiyi1)22(1cos(π1n+2))i=1nyi2,

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

This finishes the proof of the claim ().

Since the operator is linear and positive definite,   y,   the solution to (), exists and is unique.

The next part is to observe

(yo)T(dh2dhx2)Iy=(yo)Tf(yo)Tyo(f)Tf.

and putting this together with ()     (yo)T(dh2dhx2)Iyαnh2yo2

αnh2(yo)Tyo(yo)Tyo(f)Tf

and

(yo)Tyo(αn)1h2(f)Tf.

The notation

y(x)=<y(x0)y(x1)y(x2)y(xn)y(xn+1)>

will be used for the exact values of  y(x),  as well as the notation

yo(x)=<y(x1)y(x2)y(xn)>

for the interior points of  y(x).

The vector  (y(x)y)  satisfies the problem

(dh2dhx2)I(y(x)y)=(dh2dhx2)Iy(x)f,

with  (y(a)y0)=(y(b)yn+1)=0.

So the estimate () can be applied to get

(yo(x)yo)T(yo(x)yo)(αn)1h2((dh2dhx2)Iy(x)f)T((dh2dhx2)Iy(x)f).

Recall ()

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

where  xh<zh<x  and  x<zh<x+h.

Now,

(dh2dhx2)Iy(x)=<dh2dhx2y(x1)dh2dhx2y(x2)dh2dhx2y(xn)>

and

f=<y(x1)y(x2)y(xn)>.

So each component of    ((dh2dhx2)Iy(x)f)    has the form

(124f(iv)(zh)+124f(iv)(zh))h2.

This leads to the estimate

((dh2dhx2)Iy(x)f)T((dh2dhx2)Iy(x)f)nMh2,

where

  M=112maxaxb(f(iv)(x)).

After combining the inequalities () and ()

yo(x)yo2nM(αn)1h4.


Generalities about Difference Solutions

Before moving on to describe the effect of using a higher order finite difference operator, it is useful to study some generalities about the solutions of linear equations. This provides the needed motivation to the design of methods and proofs.

Suppose  L  is a linear operator that is positive definite, which means to say there exists a constant  α>0,  not depending on  x  such that

xTLxαxTx.

Then as was explained for matrices at ()

L12(α)1.

Now, if the exact solution  y(x)  to some problem is given by

Ly(x)=f(x)

except that  f(x)  is only known up to some degree of approximation by  f,  then the solution  y  to the equation

Ly=f,

is an approximation to the desired exact solution  y(x).  Since

L(y(x)y)=(f(x)f),

the closeness of the approximation can be estimated with

(y(x)y)2L12(f(x)f)2(α)1(f(x)f)2.

With regards to finite difference methods the strategy will be to define a linear operator  Lh  such that

yTLhyαyTy,

with  α  independent of  h.  Then for the exact solution  y(x),  to whatever problem,

(y(x)y)2(α)1(Lhy(x)Lhy)2.

When  Lhy(x)  is known to some degree of approximation by  f, that is when

 Lhy(x)f2Bh

then for  y,  the solution to  Lhy=f,

(y(x)y)2(α)1(Lhy(x)f)2(α)1Bh.

Third Order Estimation

In this section the properties of the five-point finite difference operator for the approximation to the second derivative is studied. The notations introduced in the previous sections of this chapter will be reused. The notations   y   and   yo   will have the same meaning as in the section Solution by Finite Differences. The notations

(dh2dhx2)i   and   (dh2dhx2)I

are redefined to represent the five-point operator.

h2(h2hx2)1 y =.
112y413y312y2+53y11112y0


h2(h2hx2)i y =.
112yi+243yi+1+52yi43yi1+112yi2


for i=2, 3, n1and


h2(h2hx2)n y =.

1112yn+1+53yn12yn113yn2+112yn3

The five-point operator for the second derivative is third-order accurate at the ponts nearest the endpoints of the interval and being a centered difference operator, is fourth-order accurate for the more internal points.

The expressions for   (dh2dhx2)i  are rearranged to make the summation procedure to be performed easier to follow. These identities are verified simply by comparing coefficients.

h2(h2hx2)1 y=76(y2+2y1y0)


+(112(y4y3)13(y3y2)+512(y2y1)16(y1y0))

112(y2+2y1y0)

112(y3+2y2y1)


h2(h2hx2)i y=76(yi+1+2yiyi1)


112(yi+2+2yi+1yi)

112(yi+2yi1yi2)


for i=2, 3, n1and


h2(h2hx2)n y=76(yn+1+2ynyn1)


+(112(yn2yn3)+13(yn1yn2)512(ynyn1)+16(yn+1yn))

112(yn+2yn1yn2)

112(yn+1+2ynyn1)

The intent of this section is to establish the inequality.

h2(yo)T(h2hx2)Iy23i=1n+1(yiyi1)2.

This is the most technically difficult part of making estimates as to the accuracy of finite difference approximations. The remainder of the analysis follows by applying the reasoning described in the section Generalities about Difference Solutions. The estimates of the accuracy to which five-point finite differences estimate the second derivative of a function are easier and will be covered in a separate section.

Take into account that  y0=0.   So  y1=(y1y0)  and then

h2y1(h2hx2)1y=76y1(y2+2y1y0)

+y1(112(y4y3)13(y3y2)+512(y2y1)16(y1y0))

112y1(y2+2y1y0)112y1(y3+2y2y1)


=76y1(y2+2y1y0)

112y1(y2+2y1y0)112y2(y3+2y2y1)

+112(y2y1)2+16(y1y0)2

112(y3y2)(y2y1)+112(y2y1)(y1y0)

+(112(y4y3)13(y3y2)+13(y2y1)13(y1y0))(y1y0)



h2yi(h2hx2)iy=76yi(yi+1+2yiyi1)

112yi(yi+2+2yi+1yi)112yi(yi+2yi1yi2)


=76yi(yi+1+2yiyi1)

112yi+1(yi+2+2yi+1yi)112yi1(yi+2yi1yi2)

+112(yi+1yi)2+112(yiyi1)2

112(yi+2yi+1)(yi+1yi)112(yiyi1)(yi1yi2)


for i=2, 3, n1and

Take into account that  yn+1=0.   So  yn=(yn+1yn)  and then

h2yn(h2hx2)ny=76yn(yn+1+2ynyn1)

+yn(112(yn2yn3)+13(yn1yn2)512(ynyn1)+16(yn+1yn))

112yn(yn+2yn1yn2)112yn(yn+1+2ynyn1)


=76yn(yn+1+2ynyn1)

112yn1(yn+2yn1yn2)112yn(yn+1+2ynyn1)

+16(yn+1yn)2+112(ynyn1)2

112(ynyn1)(yn1yn2)+112(yn+1yn)(ynyn1)

+(yn+1yn)(112(yn2yn3)13(yn1yn2)+13(ynyn1)13(yn+1yn))

This organization of the terms leads to

h2(yo)T(h2hx2)Iy=h2i=1nyi(h2hx2)iy=

i=1nyi(yi+1+2yiyi1)+16i=1n+1(yiyi1)216i=2n1(yi+1yi)(yiyi1)

+(112(y4y3)13(y3y2)+13(y2y1)13(y1y0))(y1y0)

+(yn+1yn)(112(yn2yn3)13(yn1yn2)+13(ynyn1)13(yn+1yn))

It is known from () that the sum   i=1nyi(yi+1+2yiyi1)=i=1n+1(yiyi1)2.   Using this and moving the first and last terms of the second sum

h2(yo)T(h2hx2)Iy=

i=1n+1(yiyi1)2+16i=2n(yiyi1)216i=2n1(yi+1yi)(yiyi1)

+(112(y4y3)13(y3y2)+13(y2y1)16(y1y0))(y1y0)

+(yn+1yn)(112(yn2yn3)13(yn1yn2)+13(ynyn1)16(yn+1yn))

Using ()


|i=3n2(yi+1yi)(yiyi1)|12(y3y2)2+i=4n2(yiyi1)2+12(yn1yn2)2,

|(y3y2)(y2y1)|12α12(y2y1)2+12(y3y2)2/α12,

|(ynyn1)(yn1yn2)|12β12(ynyn1)2+12(yn1yn2)2/β12,

the inequality next follows.

h2(yo)T(h2hx2)Iyi=1n+1(yiyi1)216(y1y0)216(yn+1yn)2

+112(11/α12)(y3y2)2+16(112α12)(y2y1)2

+16(112β12)(ynyn1)2+112(11/β12)(yn1yn2)2

+(112(y4y3)13(y3y2)+13(y2y1))(y1y0)

+(yn+1yn)(112(yn2yn3)13(yn1yn2)+13(ynyn1))

Now, using ()

|(112(y4y3)13(y3y2)+13(y2y1))(y1y0)|

124(α22(y4y3)2+1α22(y1y0)2)+16(α32(y3y2)2+1α32(y1y0)2)

16(α42(y2y1)2+1α42(y1y0)2).

Taking   α12=1+52,α22=8,   and   α32=α42=2+354   yields

16(y1y0)2+112(11/α12)(y3y2)2+16(112α12)(y2y1)2

+(112(y4y3)13(y3y2)+13(y2y1))(y1y0)

13((y4y3)2+(y3y2)2+(y2y1)2+(y1y0)2).

In identical fashion

16(yn+1yn)2+16(112β12)(ynyn1)2+112(11/β12)(yn1yn2)2

+(yn+1yn)(112(yn2yn3)13(yn1yn2)+13(ynyn1))

13((yn+1yn)2+(ynyn1)2+(yn1yn2)2+(yn2yn3)2).

The sought inequality has been established.

h2(yo)T(h2hx2)Iy23i=1n+1(yiyi1)2.


Two-dimensional Domain

h2(h2hx2)1,j e =.

112e4,j13e3,j12e2,j+53e1,j1112e0,j



=76(11e2,j+21e1,j11e0,j)

+(112(e4,je3,j)14(e3,je2,j)+512(e2,je1,j)14(e1,je0,j))



h2(h2hx2)i,j e =.

112ei+2,j43ei+1,j+52ei,j43ei1,j+112ei2,j

=76(11ei+1,j+21ei,j11ei1,j)

+112(ei+2,jei+1,j)112(ei+1,jei,j)

+112(ei,jei1,j)112(ei1,jei2,j)



for i=2, 3, m1and

h2(h2hx2)m,j e =.

1112em+1,j+53em,j12em1,j13em2,j+112em3,j



=76(11em+1,j+21em,j11em1,j)

+(112(em2,jem3,j)+14(em1,jem2,j)512(em,jem1,j)+14(em+1,jem,j))



Grid Vectors

Discrete Laplacian

To approximate u(x, y) numerically, use the grid

(xi,yj), 0im+1, 0jn+1.
with
xi=ih,yi=j k
and
h=a/(m+1),k=b/(n+1)

The second partial derivatives

2x2 u(x,y)      and      2y2 u(x,y)

can be approximated on the grid by difference quotients

(h2hx2) u(xi,yj)      and      (k2ky2) u(xi,yj).

These difference quotients can be chosen by the number of points or the order of accuracy. In either case they will be the same as explained in the section on difference quotients in the chapter Definitions and Basics. The possibilities arising from choosing difference quotients with less than the maximum order of accuracy using some other criteria such as a minimization of the size of or differences in the coefficients, relative to order, is not analyzed at this point.

Since it is cumbersome to include many indices in notations for difference operators the same expression for a difference quotient that approximates a second derivative will be reused for different order operators. The matter as to which one will be made clear when needed. Certain generalities apply to any one of them and can be discussed in this light.

The Laplacian  Δu(x,y)  then can be approximated on the interior of the grid by

Δh,ku(xi,yj)=h2hx2u(xi,yj)+k2ky2u(xi,yj)

(h2hx2) u(x1,yj) =.

112u(x1+3h,yj)+13u(x1+2h,yj)+12u(x1+h,yj)53u(x1,yj)+1112u(x1h,yj) h2

(h2hx2) u(xi,yj) =.

112u(xi+2h,yj)+43u(xi+h,yj)52u(xi,yj)+43u(xih,yj)112u(xi2h,yj) h2

for i=2, 3, m1and

(h2hx2) u(xm,yj) =.

1112u(xm+h,yj)53u(xm,yj)+12u(xmh,yj)+13u(xm2h,yj)112u(xm3h,yj) h2


The second partial derivative
2y2 u(x,y)
can be approximated on the grid by difference quotients
(k2ky2) u(xi,yj).

These difference quotients are given by

(k2ky2) u(xi,y1) =.

112u(xi,y1+3k)+13u(xi,y1+2k)+12u(xi,y1+k)53u(xi,y1)+1112u(xi,y1k) k2

(k2ky2) u(xi,yj) =.

112u(xi,yj+2k)+43u(xi,yj+k)52u(xi,yj)+43u(xi,yjk)112u(xi,yj2k) k2

for j=2, 3, n1and

(k2ky2) u(xi,yn) =.

1112u(xi,yn+k)53u(xi,yn)+12u(xi,ynk)+13u(xi,yn2k)112u(xi,yn3k) k2

Template:BookCat