Numerical Solution of Partial Differential Equations
Numerical Solution of Partial Differential Equations MA 584
Popular in Course
verified elite notetaker
Popular in Mathematics (M)
This 22 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 584 at North Carolina State University taught by Staff in Fall. Since its upload, it has received 11 views. For similar materials see /class/223724/ma-584-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for Numerical Solution of Partial Differential Equations
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/15/15
Finite Difference Alethod 33 6 Finite difference methods for two dimensional second order elliptic problems There are many important applications of elliptic partial differential equations Below are some examples Laplace equation um uw 0 61 a Poisson equation um uw f 62 The solution to a Laplace or Poisson equation sometimes is called a potential equation This is because for a vector eld v that is conservative meaning that curlv 0 its potential function u satis es um uw V v where V 3 IT is the gradient operator Helmholtz equation um uw u f If lt 0 it is called generalized Helmholtz equation and it is easy to solve If gt U is big then the problem is hard to solve a General selfadjoint elliptic PDEs V am 2Vucv y MicaM I x 2 64 or mark 112 Macv yu 27135331v 65 We should assume that Many does not change the sign in the solution domain a General elliptic PDEs diffusion and advection equations avyuxx 2bauxy Cxzyluyy dc ylux diallug 9557 3 37 Z Z 3 It can be rewritten as V 12 yVuacv y Wacv 3 Vu Cacv yu ing 66 after some transformation where wcy is a vector Below are some examples of nonlinear elliptic PDEs Diffusion and reaction equation V 39 no 2Vuc 21 u 67 The V ac yVuac y is called the diffusion term the nonlinear term u is called the reaction term If 0101531 E 1 it is also called a nonlinear Poisson equation p Laplacian equation div mayo 2m 0 p 2 2 68 0 Minimal surface equation div U 69 W Note that an elliptic PDE P au U can be regarded as the steady state solution of a corresponding parabolic equation m P u If a linear PDE is de ned on a rectangle domain then we can use a nite difference approximation dimension by dimension both for the equation and the boundary condition The most difficult part of a finite difference method is how to solve the resulting linear system of equation efficiently 61 Boundary conditions and compatibility condition Assume a second order elliptic PDE is defined on a domain S2 Let 8S2 be the boundary of S2 and n is the unit normal direction positive according to the right side rule see Fig 8 11 882 Figure 8 A diagram of a two dimensional domain S2 its boundary 8S2 and its unit normal direction Some commonly used boundary conditions are listed below Finite Difference Allethod 35 Dirichlet boundary condition the solution ucy3g u0cy is given along the boundary Neumann or ux boundary condition the normal derivative Vu n an uxnx uyny gcy is given along the boundary where n mmy is the unit normal direction 7377 y is 35 Robin or mixed boundary condition crcyucy xy given along the boundary 8S2 For some domains for example we can de ne a periodic boundary condition For example if S2 11 b X 0 the periodic boundary condition in the c direction is uay uby the periodic boundary condition in the y direction is ucc uc d We can have different boundary condition on different parts of the boundary For example for a channel ow in a domain 11 x 051 where b a gtgt at c it is reasonable to assume that the ux boundary condition at a a and a b but nonslip u 0 along the boundaries 3 c and y d For a Poisson equation with pure Neumann boundary condition the solution does not exist unless the compatibility condition is satis ed If u is a solution to the following problem 3 Auf39r7y7 6S2 Integrate and use the Green s theorem we can get ff Audacdy fcydcdy 6937315 fltcydmdy 39995st L amdm If the compatibility condition is satis ed and 8S2 is smooth then the solution do exist but it is not unique If uc y is a solution then uc y C is also a solution with arbitrary constant C In this case we can specify the solution at a particular point for example uco yo U to make the solution well de ned 62 The central nite difference method with ve point stencil for Poisson equation Consider the Poisson equation umx uyy f7yv E Z a X Cad u0cy Dirichlet BC 611 36 Z Li If f E L2SZ then the solution exists and it is unique Analytic solution is rarely available Now we discuss how to use the nite difference equation to solve the Poisson equation a Step 1 Generate a grid A uniform Cartesian grid can be used I xaih i012m hm ma 612 d yjzcfjhy jU12n 1197 613 We want to nd an approximate solution U15 to the exact solution at all the grid points 12213 where Mingj is unknown So there are m 1n 1 unknown for Dirichlet boundary condition Step 2 Substitute the partial derivatives with a nite difference formula in terms of the function values at grid points to get Min 1213 2uiayj 37139le 323 Mindy I 2uiayj Mini31 2 2 hr by fijnja i17quot39m 17 jZIvH39n lv where fig an 21339 The local truncation error satis es hm2 3421 hy2 3421 T 7 7 14 1 12 854 12 834 6 De ne h max hr by 615 The nite difference discretization is consistent if lim 0 616 h W Therefore the discretization is consistent and second order accurate If we remove the error term in the equation above and replace the exact solution Man 313 with the approximate solution U15 which is the solution of the linear system of equations Ui m Uma h 2 U 7 1U 7 2 2 U at mm w U 2 f 63917 The nite difference scheme at a grid point mmj involves ve grid points east 6 north west south and the center The center is called the master grid point Solve the linear system of equations to get an approximate solution at grid points how 0 Error analysis implementation visualization etc Finite Difference Allethod 37 63 Matrixvector form of the nite difference equations Generally if one wants to use a direct method such as Gaussian elimination method or sparse matrix techniques then one needs to nd out the matrix structure If one use an iterative method such as Jacobi Gauss Seidel SORw methods then it may be not necessarily to have the matrix and vector form In the matrix vector form AU F the unknown is a one dimensional array For the two dimensional Poisson equations the unknowns U15 are a two dimensional array Therefore we need to order it to get a one dimensional array We also need to order the nite difference equations It is common practice that we use the same ordering for the equations and for the unknowns There are two commonly used ordering One is called the natural ordering that ts sequential computers The other one is called the red and black ordering that ts parallel computers Figure 9 The natural ordering left and the redblack ordering right The natural row ordering In the natural row ordering we order the unknowns equations rowwise therefore the kth equation corresponding to j with the following relation kim 1j 1 i12 m l j12 n l 618 We use the following example to verify the matrixvector form of the nite difference equations Assume that hm hy h m n 4 so we will have nine equations and nine unknowns The coefficient matrix is 9 by 9 To write down the matrixvector form we use 38 Z Li a onedimensional array x to express the unknown U15 351 2 U11 352 2 U21 353 2 U31 354 I U12 355 I U227 619 6 I U32 7 I U13 its 2 U23 359 I U33 If we order the equations the same way as we order the unknowns then the nine equations from the standard central nite difference scheme using the ve point stencil are 1 HO 1110 4ml24 I fH T7 1 U20 h7l 423 55 1 21 1 u30u4l h72 4 E3 E6 2 f3 T 1 no h71 4 E4 E5 E7 2 132 722 1 72 52 54 4355 56 58 13922 1 7 U42 h735 4 E6 E9 Ina hi2 1 u03ul4 h74 4 E7 E8 I fl3 T 1 u 4 h757 4 E8 E9 I f23 hi2 1 U34 2143 h7lt 6 s 49 I f33 T Now we can write down the coef cient matrix easily It is block tridiagonal and has the following form B I U 1 A 72 I B I 620 U I B where I is a 3 X 3 identity matrix 4 1 U B 1 4 1 Finite Difference Alethod 39 For a general n by n grid we will have B I 4 1 I B I 1 4 1 IB 1 4 Note that A is a symmetric positive de nite matrix and it is weakly diagonally dominant Therefore A is nonsingular and there is a unique solution The matrixvector form is useful to understand the structure of the linear system of equations and it may be necessary if a direct method such as Gaussian elimination or sparse matrix techniques are used for solving the system However it is more convenient sometimes to use the two parameters system j especially if an iterative method is used to solve the system It is more intuitive and useful to visualize the data using two index system The eigenvalues and eigenvectors of A can be indexed by two parameters p and k corre sponding to wave numbers in the a and y directions The p kth eigenvector um has n2 elements for a n by n matrix of the form above ugh sinp7rih sink7rjh j 1 2 n 621 for p k 1 2 n The corresponding eigenvalues are N cosp7rh 1 cosk7rh 1 622 The least dominant eigenvalue the smallest in the magnitude is A 27r 0012 623 The dominant eigenvalue the largest in the magnitude is 4 Win2 624 Therefore we have the following estimates 1 4 i min vk 27139v hi2 comm1 HAH2HAquot H2 2 0022 HAH2 max WM 2 HAquot H2 2 625 Since the condition number is considered to be large we should use double precision to reduce the effect of round off errors 40 Z Li 7 Maximum Principle and Error Analysis Given an elliptic operator 2 2 2 21 8 L xi 7 Ci 8152 3533 832v b2 aclt0 for 3231 ESZ agtUCgtU The maximum principle is the following theorem Theorem 71 Ifuc y 6 C362 satis es Luc y 2 U is a bounded domain SZ than uc n has its maximum on Me boundary Proof If the theorem is not true then there is an interior point 0210 6 S2 such that ucgy0 2 ucu for all 3531 6 S2 Therefore form the necessary condition of the local extrema we have an an 75 0 75 0 as 010 32 010 Since 350310 is not on the boundary and uc y is continuous there is a neighborhood of 0310 within the domain S2 that we can use the Taylor expansion 1 um Am Ag WEN0 5 mm 2AmAyu2 AgFug cam Ag3 where the superscript of 0 indicating that the functions are evaluated at 50 go for example mgr gang0 Since we have ucg Argo Ag 3 ucgy0 for all A515 and Ag that are small enough We have 1 i AmVuEx 2AmAgugy Ag2ugy lt 0 71 On the other hand from the given condition we know that Lu0 aougx 2192129 cougy 2 U 72 where a0 acgy0 and so forth In order to match the Taylor expansion to get an contradiction we rewrite the inequality above as 2 0 0 0 0 2 0 02 a 0 a b 0 b 0 w 0b VM MHVMxMuwbm W M C an 20 7393 where M gt U is a constant The role of M is to make some choices of A3 and Ag that are small enough Next we take A 1 0 A 0 3 7 I 7 M39 l vaDM Finite Difference Alethod 41 From 71 we know that 0 0 0 a 0 2b 0 b 0 Mum YUM Mum S U 74 Now we take A3 0 Ag co 1W M a0 From 71 again we know that 1 AZFug 2 b0 2 M c0 a0 ugy s U 75 From 74 and 711 we know that the left hand side of 73 should not be positive which contradicts the the condition Lu0 aougx 2192129 cougy 2 0 This completes the proof D Similarly if Lu 3 0 then n has its minimum value on the boundary of 2 For general elliptic equations with no cross derivative term um1 the maximum principle is the following Let Lu an 2buw mw dlux dguy an 0 521 6 S2 2 aclt0v agt03cgt0 6S0 where S2 is a bounded domain Then uc y can not have a positive local maximum or a negative local minimum in the interior of 2 This can be easily proved using Theorem 71 71 Discrete maximum principle Theorem 72 Given a grid function U15 239 U1 m j U12 71 If the discrete Laplacian operator using the central vepoint stencil satis es 7 Ui lj Uilj Uij l Uijl 4023 AhUij I12 2 07 76 i12m 1 j12n 1 Then U15 attains its maximum on the boundary IfAhUij s 0 then U15 attains its minimum on the boundary Proof Assume the theorem is not true then U15 has its maximum at an interior grid point 2390 jg We have inu gt UM for all 239 and j s Therefore we have 1 Uiogjo 2 1U1 o ljo UileJo Uiogjo l Uimjo 39 42 Z Li On the other hand from the condition AMIquot 2 0 we know that 1 Uiogjo g 1Uio UU UileJo Uimjo l Uimjo 39 This contradicts the inequality above except that all U15 at the four neighbors of 2390 jg have the same value as Ui0 jg So Uir 1731 is also a maximum If we apply the same arguments several times then we will reach the boundary and know that UDJU is also a maximum From the proof we can see that if U13 has its maximum in interior then U13 is a constant If AU g 0 then we can consider Uij which gives us the conclusion D 72 Error estimates of the FD method for Poisson equations Lemma 71 Let U15 be a grid function that satis es AhUij fijv Z 07 an Ui Ij Umu Uiu I Uio H 4 h2 with homogeneous boundary condition Then 1 1 Z lt 7 Z 7 39 HUHoc Ogig nllfm 8 lgquot lAnUm 8 1173 an 7 8 Proof De ne a grid function wij I 2 CM a2 le g 79 where 2 mi 1zh yj 1h 170717quot39n7hgv corresponding to the continuous function wc 4i Q2 y 2 We have 2 I 4 hi Q 34quot 1 710 l to A Z hwy may WM 12 354 3214 where 33 is some point near iyj Therefore we have AMUu HfHoc39wij AhUij Hfo S 0 AMUij HfHoc39wij I AhUij Hfo 2 0 711 From the discrete maximum principle we know that U15 Hf waij has its maximum on the boundary while U15 Hf xwij has its minimum on the boundary That is Hfo HwinaQ S U13 HfHonij S U13 and Uij S Uij HfHoc39wij S Hfo wallow Finite Difference Allethod 43 since U15 is zero on the boundary and Hfowij 2 0 It is easy to check that 1 llwijllao g therefore we get 1 1 Fume s Us guinea 712 That completes the proof of the lemma D Theorem 73 Let U15 be the solution of the nite di erence equations using the standard central ne point stencil for Poisson equation with a Dirichlet boundary condition Then the global error is second order accurate and satis es h2 HU 11HQC niaxUij uciyj g max max uww 713 Proof We know that AhUij I fij T1 AhEij 2 T1 where T15 is the local truncation error at an 313 and satis es h2 T1 j g E max max WWW Therefore from Lemma we have 1 h2 HEHOC s 4 er s maxlumyl 73 Finite difference discretization for general second order elliptic PDE using dimension by dimension If the domain of the interest is a rectangle a b X c d and there is no cross derivative term uw in the equation then we can discretize the PDE dimension by dimension We verify this for the example below Consider V 39 MimiWu 1675721 I x 2 or mm 2mm m I f with Dirichlet boundary condition at a b y c and y d but a Neumann boundary condition um gy along a a If we use a uniform Cartesian grid acizaH39hm i20313mv hm yjzcjhyv jZUvLH39Tl by 44 Z Li If we discretize the PDE dimension by dimension the nite difference equation then is pijUilj ngg pi jUiJ PijUi Ij w Fi Uijl Fi pi7j Uij magma I 2 lijUij I fij by 714 for 239 12 m 1 and j 12 n 1 where pii ij 21235 ilk23131 and so forth For the indices 239 U j 12 n 1 we use the ghost point method to deal with the Neumann boundary condition Using the central nite difference scheme for the ux boundary condition U1 U l W 2 go or UM 2 Um 2211mm J 212 n 21 5C Plugging this into the nite difference equation at U j we get 2L3 pjU1J 2 ng p jUm hr P073 U041 2104 P0jU03 p07jU0j l hm 2pjgyj h 39 1031103 I ij For general second order elliptic PDE with no cross derivative term aw V 39 MimiWu W 39 W 2 dis 2u 2 x 2 we can still use the central nite difference scheme of ltlt 1 h otherwise we should use the upwinding scheme to deal with the advection term Finite difference formula for cross derivative um If there is a cross derivative term um we can not simply use the method of dimension by dimension However it is easy to derive a centered nite difference scheme for aw ucy by uc y by 2119 uc hwy by uc hwyl by uc hmy by uc hwyl by 411x119 aw N S lyuz x 716 Using Taylor expansion at 3531 we can show the nite difference formula is consistent and it it second order accurate The nite difference formula for the PDE that involves a cross derivative term has nine grid points involved The resulting linear system of equations for PDEs with a cross derivative term is more difficult to solve because it is not symmetric and there is no diagonally dominance any more and there is no upwinding scheme to deal with the cross derivative term Finite Difference Alethod 45 74 Solving the resulting linear system of equations The linear system of equations resulted from a nite nite difference discretization is very large for two dimensional problems For example using an n X n grid for an elliptic PDE the linear system has about n2 equations The coef cient matrix has the size of n2 X n2 When we take n 100 a modest number the size of the matrix would be 10 1 x 104 which can not be stored in most of computers if double precision is used However the matrix is very sparse for a selfadjoint elliptic PDE the nonzero entries are about 5722 Therefore usually iterative methods or sparse matrix techniques are used Typically for an elliptic PDE de ned on a rectangle domain or a disk several methods listed below are used frequently Fast Poisson solvers such as the FFT cyclic reduction see Usually the implemen tation is not so easy It is recommended to use software packages for example the Fishpack in Fortran which is free on the Netlib Multigrid solver for example MGDQV using a structured grid for example a nine point stencil which includes the vepoint stencil as a special case AMG algebraic multigrid solver if the coef cient Matrix is an Mmatrix Sparse matrix techniques It was quite popular before the multigrid method was developed Simple iterative method such as Jacobi GaussSeidel SORw They are easy to implement and fewer restrictions on the coef cient matrix But they are generally slow in convergence Other iterative method such as CG or PCG conjugate gradient method with pre conditioning GMRES generalized minimized residual and BICG methods for non symmetric system of equations One of the most important advantage of an iterative method is that it only needs the matrixvector multiplication So zero entries play no role In the implementation we use the component form instead of manipulating the matrix and vector In other words the equations are used directly 75 Basic stationary iterative methods Give a linear system of equation Ax b detA 0 Assume that we can write A as A M N where M is an invertible matrix Then we have Nx b or Mx Nx b or x M INX M lb We can form an iteration starting from an initial guess 3 x m MquotNx Mquotb k12 717 46 Z Li Such iteration may converge or diverge depending the spectral radius of MMquot max i quot Since T M I N is a constant matrix such a method is called a stationary iterative method 751 The Jacobi iterative method Solve for the diagonals If we solve for as in the rst equation and 52 in the second equation we can get 1 it I fltbl al2m2 a133 53quot39 alnn ll 1 2 I i b2 a2licl a23ft3quot39 a2nn 0622 1 351 I bi auxI anmm amxn 0m 1 in I i bn a ilxl an2x2quot39 ann ln l arm Given an initial guess x the Jacobi iterative method is k l 1 k k k cl L b1 a2m2 a13m3 alnmn all k l 1 k k k 3521L 52 06mm a23 53quot39 a2n 5n 1122 k l 1 k k k 321 L 7 bi anac ai2c2 amacn an k l 1 k k k as 1L bu aux an2m2 an7nlmnl arm In a compact form it can be written as 1 quot as f I E 11535 z 12 n 718 an 7 3441 which provides the basis for easy programming For the nite difference equation Um 2UiUil 7 1 f1 with a Dirichlet boundary condition we have UM uaU h2f 39 2 2 Uk UC 2 Ugh 1 I2 1IL2f1 23n1 UM 7 U 2ubh2fnl39 n l 2 2 Finite Difference Alethod For a two dimensional Poisson equation it is Uk Uk Uk Uk h2f MI I 1 I l 7 Um va 2 12 11 752 The GaussSeidel iterative method Solve for the diagonals and use the most updated information In the Jacobi iterative method we update all the components of Xkl based on x In the GaussSeidel iterative method we use the most updated information as follows k l 1 k k k cl L b1 al2m2 a13m3 alnmn all k l 1 k l k k 322 b2 d2lml L a23m3 a2nmn 0622 kl 1 kl kl kl k k 35 bi WIIE ai2552 am 1551 aiilil quot39 amn k l 1 k l k l k 1 337i bu ail art23 521L 39 39 39 anm lmnil arm or in compact form 1 i l n kli4 k H k 2 H 51 a b E 113353 E away 2 12 n 719 31 PM Below is a pseudocode of the GaussSeidel iterative method for solving the Poisson equation Give u0ij and a tolerance tol err 1000 k 0 u u0 while err gt tol for i1n for j1n uijgt ui 1jui1juij 1uij1 h 2fij 4 end end err maxmaxabsuu0 110 u k k 1 7 Next iteration if err gt tol end 48 Z Li 753 The successive over relaxation SORw iterative method an extrapo lation technique Let x21 be the new update of the GaussSeidel method from x The new update using SORw is the linear combination of x 6 and 323 Xkl 1 wx c wxg gl 720 If the parameter w lt 1 the iteration above is called an interpolation if u gt 1 it is called an extrapolation or over relaxation For elliptic problems usually we should choose 1 g u lt 2 if u 1 we are back to the GaussSeidel method In the component form the SORw method can be written as w i l n as 1 bi Emir E 11535 i12 n 721 31 j1 l We just need to change one line in the pseudocode of the GaussSeidel method to get the SORw method uij l omegau0ij omega ui 1j ui1j uij 1 uij1 h 2fij gt4 The convergence of the SORw method depends on the choice of w For the linear system of equations obtained from the standard vepoint stencil applied to a Poisson equation with h hm hy 172 it has been shown that the optimal w is 2 1 sin7rn N 1 7rn39 73922 wept Z Note that the optimal w approaches number two as n approaches in nity For general elliptic PDEs we do not know the optimal u but we can use the optimal w for the Poisson equation as a trial value It is also recommended to use larger w than smaller ones If u is too large that the iterative method diverges we would nd it quickly because the solution will blowup Also the optimal w is independent of the right hand side 754 Convergence of the stationary iterative methods For a stationary iterative method the discussion of the convergence of the method is based on the following theorem suf cient and necessary condition Theorem 74 Given a stationary iteration Xkl TX C c 723 Finite Difference blethod 49 where T is a constant matrim and c is a constant vector The vector sequence xk converges for arbitrary x0 if and only if pT lt 1 where pT is the spectral radius ofT de 39ned as pT max 139Tv 724 that is the largest magnitude of all the eigenvalues of T Another commonly used suf cient condition to check the convergence of a stationary iterative method is the following theorem Theorem 75 If there is a matril39 norm such that lt 1 then the stationary iterative method converges for arbitrary initial guess x0 We often check HTHp norm for p 1 23 so Now let us discuss the convergence of Jacobi GaussSeidel and SORw methods Given a linear system Ax I Let D be the diagonal matrix formed from the diagonal elements of A L be the lower triangular part of A and U be the upper triangular part of A The iteration matrices for the three basic iteration methods are Jacobi method T Dquot L U c D lb GaussSeidel method T D Lquot U c D Lquot b I SORw method T I wLquot 1 wI wU c wI wLquotDquotb Theorem 76 IfA is strictly diagonally dominant that is n 1 gt Z lag 725 warn Then both Jacobi and GaussSeidel iterative methods converge The conclusion is also true when I A is weakly diagonally dominant 71 11111 2 E 10615 726 warn 2 the inequality holds for at least one row 3 A is irreducible From the theorem above we can conclude that both Jacobi and GaussSeidel iterative methods converge when they are applied to solve the linear system of equations obtained from the standard central nite difference method for Poisson equations 50 Z Li 76 Nine point discrete Laplacian for Poisson equations The purpose of the ninepoint discrete Laplacian is to get a compact 4th order accurate scheme Hu UH g Oh for a Poisson equation One of advantages of high order method is that we can use fewer grid points to get the same order accuracy compared with lower order method We will have smaller system of equations One of the disadvantages is that we will have a denser system of equations While other methods may be used we use a symbolic derivation from the second order central scheme for Recall that 2 7 7 if 4 6 352 12 8324 001 7 27 7 1 i2 0014 12 352 8322 u 39 Plugging the operator relation 82 Six 0012 into equation 727 we get 2 12 2 2 82 4 5mm 1 001 852 u 001 h22 32 4 1 mu0h We solver 61 from the equation above to get 82 h 39 h2 I g 2 1 E6350 53w 1 E63 0014 Note that 112 quot 112 1 35 1 E63 004 if h is suf cient small We get the symbolic relation 2 2 8 ML 5 xuoh4 or W 2 12 m a2 112 2 1 air 33sz 0014 With a Cartesian grid we can approximate the Poisson equation Au f to get I I12 I h2 1 166 lirak 1 659 35 fcy 0014 Finite Difference Alethod where h maxhx by We multiply 2 hb 1 1 E699 to the expression above and use the fact that Ax A202 A202 Ax 1H12 33 1H12 359 1lt12 39 1lt12 33 That is the are commutative we can get hi 2 2 hi 2 2 hz2r 2 hi 2 4 1 E6 3mm 1 T26 Swu 1 T26 1 w fcy 001 h 2 2 I 1172396m 6W fcy0h Expand the expression above we get the ninepoint scheme for the Poisson equation for example h 112 U 2UU I J 2 2 7 y 2 1 173 13 Hm G M G m Ui Ij 2UijUiIj 1 hwy 12hy2 Ui lg lz 2Ui lg U1 ljl 2024 401339 2014 Uilj l 2Uilj Uiljl For the special case when hm by h the nite difference coefficients and the linear combination of f are expressed in Fig 10 MW ing 948 O H o quotbe 0 Ho Figure 10 The coefficients of the finite difference scheme using the ninepoint stencil Z Li Possible advantages and disadvantages of ninepoint finite difference methods for Poisson equations include o It is fourth order accurate and it is still compact The coefficient matrix is still block tridiagonal Less grid orientation effects compared with the standard five point nite difference scheme Note that if we apply a2 112 2 1 ES to the Poisson equation we will get another ninepoint finite difference scheme which is not gt53 0014 compact and stronger grid orientation effects 77 A nite difference method for Poisson equation using polar coordi nates If the domain of the interest is a circle or an annulus or a fan etc see Fig 11 It is much easier to use polar coordinates a 7quot cos 0 y 7quot sin 0 728 Under the polar coordinates the Poisson equation is 1 8 an 1 32a 77 T7 if 7 0 T3quot 3quot 7 2 8492 39 32u13u182u 8T2 7 8quot 7 2 8492 If the origin is not in the domain R g r g R2 9 g H g HT we can use a uniform grid in polar coordinates to discretize the Poisson equation 2 f0quot OI anIH39AT i01m ArziRZ Rl m sz jA6L jU1N AH The discretized equation is using the conservative form 1thi m Ti g 72ng m ma 7 M2 729 l UiJ I 2 UAJH I I T12 A602 fT1793 where U15 is an approximation to the solution um 0339 Finite Difference Alethod 53 BC BC 2 T BC T BC pe iodic p0 iOdic B BC No BC 27 0 27 0 0 Figure 11 A diagram of domains and boundary conditions that may be better solved in polar coordinates 78 Treating the boundary conditions If the origin is within the domain then 0 g H g 27139 we should use the pmiodic boundary condition in 0 direction that is 22036 22030 27c In the 2quot direction R 0 needs special attention There are different methods in the literature in dealing with the pole singularity Some methods will lead to an undesirable structure of the coefficient matrix of the nite difference equation One clever approach discussed here is to use the staggered grid n 1A2quot Aquot R2 2 m 23912 m 730 2 Notice that 2quot Ar 2 and rm R2 We can use the conservative form of the discretization at 239 2 m 1 except for 239 1 At 239 1 we use the nonconservative form to take care of the pole singularity Uoj 2UjU2j 1U23 2gU0j 1 Ulj l 2UjUljl fhhgj 2quot ATP 2quot r12 A602 Note that 7 0 Ar2 and 2quot Ar2 We can see that the coefficient of U03 which is the approximation at the ghost point To is zero The above finite difference equation is 54 Z Li simpli ed to 2UljU2j 1 U2j 1Uj 2UljUljl i 7 i 7 7 I my 7 2m T A9 3 We still have a diagonally dominant system of linear equations 79 Use FFT to solve Poisson equation in polar coordinates Since the solution uquot H is periodic in H is the origin is an interior point we can approximate it by the truncated Fourier series as N2 l uquot 49 Z unrem9 731 n N2 where 239 l l and unr is the complex Fourier coefficient given by N l 1 unquot E uquot Elk mke 732 k70 Plugging 731 in to the Poisson equation we get 1 8 1 an n2 E 7un fnquot n N2 N2 1 733 where fnr is the nth coefficient of the Fourier series of fr 0 de ned in 732 For each n we can discretize in the 7quot direction using the staggered grid to get a tridiagonal system of equations which can be solved easily Assuming a Dirichlet boundary condition rm 0 uBC9 at r rm we can use the Fourier transform 2 1 5C Tmazc W URL 96 quotmy 734 2 l l o to nd af rm which is the boundary condition for the ordinary differential equation Once we have the Fourier coefficient an we can use the inverse Fourier transform 731 to get an approximate solution to the Poisson equation
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'