Finite Difference Methods for Boundary Value Problems

1 Finite Difference Methods for Boundary Value ProblemsOctober 2, 2013() Finite DifferencesOctober 2, 20131 / 52 GoalsLearn steps to approximate BVPs using the Finite Difference MethodStart with two-point BVP (1D)Investigate common FD approximations foru (x) andu (x) in 1 DUse FD quotients to write a system of Difference equations to solvetwo-point BVPH igher order accurate schemesSystems of first order BVPsUse what we learned from 1D and extend to Poisson s equation in 2D& 3 DLearn how to handle different Boundary conditions() Finite DifferencesOctober 2, 20132 / 52 Steps in the Finite Difference Approach to linear DirichletBVPsOverlay domain with gridChoose Difference quotients to approximate derivatives in DEWrite a Difference equation at each node where there is an unknownSet up resulting system of equations as a matrix problemSolving resulting linear system efficientlyCompute error when solution is known() Finite DifferencesOctober 2, 20133 / 52 Prototype Dirichlet BVP in 1D ddx(p(x)dudx)+q(x)u=f(x)a<x<b,(1)u(a) = u(b) = Whenp(x) =p, a constant, we have pu (x) +q(x)u=f(x)a<x<b,Whenp= 1 andq= 0 we have the Poisson equation u (x) =f(x)a<x<bin one dimension.

2 We will see that the minus sign is important.() Finite DifferencesOctober 2, 20134 / 52 Herep(x),q(x) are required to satisfy the bounds0<pmin p pmaxandqmin= 0 q(x) qmax.(2)For existence and uniqueness we also require thatfandqbe continuousfunctions ofxon the domain [a,b] and thatphas a continuous firstderivative there.() Finite DifferencesOctober 2, 20135 / 52 Step 1: Overlay domain with a gridSuppose that we subdivide our domain [a,b] inton+ 1 subintervals usingthe (n+ 2) uniformly spaced pointsxi,i= 0,1,..n+ 1 withx0=a,x1=x0+h,..,xi=xi 1+h,..,xn+1=xn+h=bwhereh=b an+ 1() Finite DifferencesOctober 2, 20136 / 52 The pointsxiare called the grid points or nodesx1,x2,..,xnare interior nodes (denoted by open circles inthe diagram below)The two nodesx0andxn+1are Boundary nodes (denoted by solidcircles in the diagram).e e e e e e e e euux0x1x1xi 1xixi+1xnxn+1abh- () Finite DifferencesOctober 2, 20137 / 52 Step 2: Choose Difference quotients to approximatederivatives in DEFor generalp(x) we have p(x)u (x) p (x)u (x) +q(x)u=f(x)a<x< we need Difference quotient approximations for both the first andsecond derivatives.

3 So far we have approximations for the first Difference :u (x) =u(x+h) u(x)h+O(h)Backward Difference :u (x) =u(x) u(x h)h+O(h)() Finite DifferencesOctober 2, 20138 / 52 Difference quotient foru (x)A Taylor series expansion foru(x+h) isu(x+h) =u(x) +h u (x) +h22!u (x) +h33!u (x) +O(h4).(3)Now we want an approximation foru (x) but if we solve for it we geth22!u (x) =u(x+h) u(x) hu (x) h33!u (x) +O(h4).in (3) then we still have theu (x) term. However if we consider the Taylorseries expansion foru(x h)u(x h) =u(x) h u (x) +h22!u (x) h33!u (x) +O(h4)(4)then we can eliminate theu (x) term by adding the two expansions; wehave() Finite DifferencesOctober 2, 20139 / 52u(x+h) =u(x) +h u (x) +h22!u (x) +h33!u (x) +O(h4).u(x h) =u(x) h u (x) +h22!u (x) h33!u (x) +O(h4)impliesu(x+h) +u(x h) 2u(x) =h2u (x) +O(h4)Note that the terms Difference quotient is called a second centered differencequotient or a second order central Difference approximation tou (x)It is second order centered Difference :u (x) =u(x+h) 2u(x) +u(x h)h2+O(h2)(5)() Finite DifferencesOctober 2, 201310 / 52 Another way to derive this approximation is to Difference the forward andbackward approximations to the first derivative, ,u (x) 1h[Forward Difference foru (x) Backward Difference foru (x)]which impliesu (x) 1h[u(x+h) u(x)h u(x) u(x h)h]u (x) u(x+h) 2u(x) +u(x h)hhence the name second Difference .

4 () Finite DifferencesOctober 2, 201311 / 52 Finite Difference StencilFinite Difference approximations are often described in a pictorial format bygiving a diagram indicating the points used in the approximation. Theseare called Finite Difference stencils and this second centered Difference iscalled a three point stencil for the second derivative in one kxi 1xixi+1-211() Finite DifferencesOctober 2, 201312 / 52 Finite Difference quotient foru (x)The forward or backward Difference quotients foru (x) are first orderThe second centered Difference foru (x) is second orderSo we need a second order approximation tou (x)If we subtract the expansionsu(x+h) =u(x) +h u (x) +h22!u (x) +h33!u (x) +O(h4).u(x h) =u(x) h u (x) +h22!u (x) h33!u (x) +O(h4)we getu(x+h) u(x h) = 2hu (x) +O(h3)() Finite DifferencesOctober 2, 201313 / 52which gives the (first) centered differenceFirst centered Difference :u (x) =u(x+h) u(x h)2h+O(h2)It is described by the stencilkkxi 1xixi+1-11() Finite DifferencesOctober 2, 201314 / 52 Step 3: Writing the Difference EquationWe have the ODE p(x)u (x) p (x)u (x) +q(x)u=f(x)with the approximationsu (xi) u(xi+1) 2u(xi) +u(xi 1)h2u (xi) u(xi+1) u(xi 1)2h() Finite DifferencesOctober 2, 201315 / 52 ODE: p(x)u (x) p (x)u (x) +q(x)u=f(x)LetUi u(xi) so thatU0= ,U1=.

5 Using our Difference quotients at each interior grid pointxi,i= 1,..,nwe havep(xi)( Ui+1+ 2Ui Ui 1h2) p (xi)(Ui+1 Ui 12h)+q(xi)Ui=f(xi).Ati= 1 we havep(x1)( U2+ 2U1 U0h2) p (x1)(U2 U02h)+q(x1)U1=f(x1),U0= is known so we take it to the right hand side of the equationp(x1)( U2+ 2U1h2) p (x1)(U22h)+q(x1)U1=f(x1)+p(x1) h2+p (x1) 2h,() Finite DifferencesOctober 2, 201316 / 52 Step 4: Write Difference equations as linear system ofequationsFirst consider the simple case whenp= 1 andq= 0 then we have theequations2U1 U2=h2f(x1) + U3+ 2U2 U1=h2f(x2) U4+ 2U3 U2=h2f(x3).. Un+ 2Un 1 Un 2=h2f(xn 1)2Un Un 1=h2f(xn) + () Finite DifferencesOctober 2, 201317 / 52 The corresponding matrix problem isA~U=~FwhereAis the matrix 2 100 0 12 10 00 12 10 0 12 10 0 12 (6)with the vector of unknowns~U= 1Un ~F= h2f(x1) + h2f(x2)..h2f(xn 1)h2f(xn) + () Finite DifferencesOctober 2, 201318 / 52 The linear system istridiagonalsymmetricpositive definiteO(n) operations to solveCholesky for tridiagonal system can be usedA=LLTthen forwardsolveL~Y=~Fand back solveLT~U=~Ystorage required is two vectors for matrix and one for~FNote that if we didn t have the minus sign in u (x) =f(x) then thematrix would not be positive definite.

6 () Finite DifferencesOctober 2, 201319 / 52 Example 1 - Homogeneous Dirichlet Boundary ConditionsWe want to use Finite differences to approximate the solution of the BVP u (x) = 2sin( x)0<x<1u(0) = 0,u(1) = 0usingh= 1/4. Our grid will contain five total grid pointsx0= 0,x1= 1/4,x2= 1/2,x3= 3/4,x4= 1and three interior pointsx1,x2,x3. Thus we have three unknownsU1,U2,U3. We will write the equation at each interior node todemonstrate that we get the tridiagonal system. We have() Finite DifferencesOctober 2, 201320 / 522U1 U2= 216sin( 4) U1+ 2U2 U3= 216sin( 2) U2+ 2U3= 216sin(3 4).Writing these three equations as a linear system gives 2 10 12 10 12 U1U2U3 = 216 sin( 4)sin( 2)sin(3 4) = .Solving this system givesU1= ,U2= andU3= ; theexact solution to this problem isu= sin( x) so at the interior nodes wehave the exact solution ( ,1, ).() Finite DifferencesOctober 2, 201321 / 52 Example 2 - Inhomogeneous Dirichlet BCsConsider the BVP u (x) = 2cos( x)0<x<1u(0) = 1,u(1) = 1whose exact solution isu(x) = cos( x).

7 Using the same grid (h= 1/4) asin the previous example, we still have three unknowns so we write theequations at the three interior nodes U0+ 2U1 U2= 216cos( 4) U1+ 2U2 U3= 216cos( 2)) U2+ 2U3 U4= 216cos(3 4))NowU0= 1 andU4= 1 so we simply substitute these values into theequations and move the terms to the right hand side to get() Finite DifferencesOctober 2, 201322 / 522U1 U2= 216cos( 4) + 1 U1+ 2U2 U3= 216cos( 2)) U2+ 2U3= 216cos(3 4) 1 Writing these three equations as a linear system gives 2 10 12 10 12 U1U2U3 = .Solving this system givesU1= ,U2= 0 andU3= ; theexact solution at these interior nodes is ( , , ).() Finite DifferencesOctober 2, 201323 / 52 Higher order accurate schemeTo get a higher order scheme we need to include more points in the stencilFive point stencil: xi 1xixi+1xi 2xi+2 524343 112 112u (x) =1h2[ 112u(x 2h) +43u(x h) 52u(x) +43u(x+h) 112u(x+ 2h)]+O(h4)where we derive this by combining Taylor series expansions foru(x 2h),u(x h),u(x+h),andu(x+ 2h).

8 () Finite DifferencesOctober 2, 201324 / 52If we have uniform pointsxi,i= 0,..,n+ 1, how many unknownsdo we have for a Dirichlet 2-point BVP?What do you think the structure of the resulting matrix is?Do we handle the boundaries in the same way as the three-pointstencil?() Finite DifferencesOctober 2, 201325 / 52 Summary of FD approximationsu (x)forwardu(x+h) u(x)hO(h)differenceu (x)backwardu(x) u(x h)hO(h)differenceu (x)centeredu(x+h) u(x h)2hO(h2)differenceu (x)2ndcenteredu(x+h) 2u(x) +u(x h)h2O(h2) Difference () Finite DifferencesOctober 2, 201326 / 52u (x)five-point112h2[ u(x 2h) + 16u(x+h)stencil (1D) 30u(x) + 16u(x h) u(x+ 2h)]O(h4)Recall that for IVPs we made a LTE at each step and then the total orglobal error was the accumulated error over all the time IVPs, for BVPs the total error (assuming no roundoff and using adirect solver) is just due to the error in replacing the DE with a differenceequation.

9 () Finite DifferencesOctober 2, 201327 / 52 Calculating the numerical rate of convergenceWe want to calculate the numerical rate of convergence for our simulationsas we did for IVPs. However, in this case our solution is a vector ratherthan a single number. To calculate the numerical rate using the formular=lnE1E2lnh1h2we need a single number which represents the error so we use a vectornorm. A commonly used norm is the standard Euclidean norm defined by ~x 2=[n i=1x2i]1/2or ~E 2=[1nn i=1E2i]1/2for a vector in~x Rn. Other choices include the maximum norm orthe one-norm 1() Finite DifferencesOctober 2, 201328 / 52 Structure of code for Dirichlet 1D BVPUser specifiesn, the number of interior grid points (alternately the grid spacingh);aandb, the right and left endpoints of interval;the Boundary Value atx=aand atx=b,a routine for the forcing functionf(x)a routine for the exact solution, if known.() Finite DifferencesOctober 2, 201329 / 52 Code could be structured as follows:computeh= (b a)/(n+ 1);compute grid pointsx(i),i= 0,1,2.

10 ,n+ 1;set up the coefficient matrix and store efficiently; for example, for thethree-point stencil the matrix can be stored as two vectors;set up the right hand side for all interior points;modify the first and last entries of the right hand side to account forinhomogeneous Dirichlet Boundary data;solve the resulting linear system using an appropriate solver;output solution to file for plotting, if desired;compute the error vector and output a norm of the error (normalized)if the exact solution is known.() Finite DifferencesOctober 2, 201330 / 52 Example 3 u (x) = 2sin( x)0<x<1u(0) = 0,u(1) = 0h E 2 u 2numerical 10 10 10 10 10 () Finite DifferencesOctober 2, 201331 / 52 Example 4 u (x) = 20<x<1u(0) = 0,u(1) = 0with exact solutionu(x) =x2 E 2 u 2numerical 10 10 16 Why are we getting essentially zero for the error?() Finite DifferencesOctober 2, 201332 / 52 From the Taylor series derivation of the second centered differenceu(x+h) +u(x h) 2u(x) =h2u (x) + 2h44!