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()

2 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.

3 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.

4 ,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.

5 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!

6 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 .

7 () 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!

8 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.

9 P(x)u (x) p (x)u (x) +q(x)u=f(x)LetUi u(xi) so thatU0= ,U1= .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).

10 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.