Numeric solutions of elliptic PDEs March 25 2009 Outline Numerical Solutions of Elliptic Equations Larry Caretto Mechanical Engineering 501B Seminar in Engineering Analysis Review last class Numerical solution of elliptic equations Laplace s equation as an example Finite difference forms Iterative solution of algebraic equations Sample results March 25 2009 2 Review Properties of Solutions Review Two Space Dimensions time step n n 1 2 n 1 T 2T 2T 2 2 y ij x t Consistency truncation error becomes zero as step sizes approach zero Stability errors remain bounded Convergent tends to the exact solution of the differential equation as the grid size tends towards zero Physical reality Solutions produce physically realistic results Accuracy Many sources of error in numerical solutions 3 Extension of one space dimension Have grid in yj as well as xi and time Explicit method has almost no difference but has more stringent stability limit Cannot apply Thomas algorithm directly to implicit method Define fx t x 2 and fy t y 2 4 Finite difference Grid Review Algorithms Node related to four Old nearest neighbors at both y time time steps x x n CN uses all u i 1j uni 1j Tnij T values at y un 1ij 1 old time y unij 1 step un 1i 1j New time x y x n 1 Tn 1ij u i 1j un 1ij 1 ME 501B Engineering Analysis n 1 n n n n Explicit Tij f x Ti 1 j Ti 1 j f y Tij 1 Tij 1 method 1 2 f x 2 f y Ti n Crank Nicholson unij 1 f yTijn 11 f xTi n1 1j 2 1 f x f y Tijn 1 f xTi n1 1j f yTijn 11 f x Ti n1 j Ti n1 j f y Tijn 1 Tijn 1 2 1 f x f y Tijn Rijn Fully implicit f yTijn 11 f xTi n1 1j 1 2 f x 2 f y Tijn 1 f xTi n1 1j f yTijn 11 Tijn 5 6 1 Numeric solutions of elliptic PDEs March 25 2009 Review Explicit Stability Alternating Direction Implicit ADI is technique for allowing use of Thomas algorithm for simple solution 1 2fx 2fy must be positive for stability so fx fy must be less than 1 2 Tijn 1 f x Ti n1 j Ti n1 j f y Tijn 1 Tijn 1 1 2 f x 2 f y Ti n fx f y 1 t t x 2 y 2 2 At each time step treat one direction as explicit and one as implicit Even number of time steps required ui 1 j ui 1 j 2uij uijn 1 uijn 1 2uijn t x 2 y 2 n 2 n 2 n 2 n 2 ui 1 j ui 1 j 2uij uij 1 uij 1 2uij uij uij t x 2 y 2 uij uijn To cut x and y by a factor of 2 we must cut t by a factor of 4 Work increases by a factor of 16 7 8 Elliptic Equations 2D Cartesian Laplace Solutions at any point in the region depend on conditions at all boundaries No open boundaries like parabolic case Main examples are Laplace and Poisson Equations Use finite difference grid where xi x0 i x and yj y0 j y 2u 0 x0 xi xN and y0 yj yM 0 i N and 0 j M x xN x0 N and y yM y0 M Must specify boundary conditions on u at all boundaries 2u F space Laplace equation in two dimensional Cartesian coordinates u u 0 x 2 y 2 2 2 South j 0 0 i N North j M 0 i N West i 0 0 j M East i N 0 j M Dependent variable u uij u xi y j 9 Finite Differences Finite difference Equation Use second order expressions for second derivatives ui 1 j ui 1 j 2uij 2u O x 2 x 2 ij x 2 uij 1 uij 1 2uij u O y 2 y 2 ij y 2 2 x y ui 1 j ui 1 j 2uij uij 1 uij 1 2uij 2u 2u O x 2 y 2 0 2 2 x y x 2 y 2 ij Multiply by x 2 11 ME 501B Engineering Analysis 10 Links central node to four nearest neighbors u ij 1 Only y parameter x x x y ui 1j ui 1j uij For 1 y 2 1 2 4 uij 1 ui 1 j ui 1 j 2 uij 1 uij 1 2 1 2 uij 0 12 2 Numeric solutions of elliptic PDEs March 25 2009 Finite difference Equation General Equation For 2D Cartesian Laplace equation with x y 1 uij is the average of its four nearest neighbors uij ui 1 j ui 1 j uij 1 uij 1 Most general case has five diagonals but can have different terms k T k T Occurs with variable properties x x y Q 0 Notation Aijpoint is general coefficient 4 Consider Dirichlet boundary conditions uij known at all boundary nodes Need to find N 1 M 1 unknown values of uij at central nodes ij refers to a particular node Point N orth S outh E ast W est refers to neighboring nodes by direction General equation shown below AijS uij 1 AijW ui 1 j AijP uij AijE ui 1 j AijN uij 1 bij 13 Solving the Equations 14 Iterative Solutions Typically have large number of equations forming sparse matrix Simplest examples are Jacobi GaussSeidel and Successive Over Relaxation Move from iteration n to iteration n 1 Iteration 0 is initial guess often all zero Solve equation for uij and use this as basis for iteration bij AijS uij 1 AijW ui 1 j AijE ui 1 j AijN uij 1 uij P Aij AijP For x y 01 have 992 equations so matrix has 96x106 potential coefficients Only 48609 0 051 are nonzero Want data structure and algorithm for handling sparse matrices Gauss elimination uses storage for banded matrices Iterative methods used for solutions uij bij AijS uij 1 AijW ui 1 j AijE ui 1 j AijN uij 1 15 Iterative Solutions II uij n 1 bij AijS uij n 1 AijW ui n1 j AijE ui n1 j AijN uij n 1 Gauss Seidel uses most recent values u bij AijS uij n 11 AijW ui n1 j1 AijE ui n1 j AijN uij n 1 Relaxation basis Gauss Seidel provides a correction that can be adjusted uij n 1 uij n uij n 1 GS uij n Relaxation Factor ME 501B Engineering Analysis 16 Relaxation Methods Use superscript n for iteration number Jacobi iteration uses all old values n 1 ij y Relaxation factor greater than or less than 1 is over or underrelaxation Underrelaxation procures stability in problems that will not converge Overrelaxation procures speed in wellbehaved problems uij n 1 uij n uij n 1 GS uij n 1 uij n uij n 1 GS 1 uij n AijS uij n 11 AijW ui n1 j1 AijE ui n1 j AijN …