Numeric solutions of elliptic PDEs March 25, 2009ME 501B – Engineering Analysis 1Numerical Solutions of Elliptic Numerical Solutions of Elliptic EquationsEquationsLarry CarettoMechanical Engineering 501BSeminar in Engineering AnalysisMarch 25, 20092Outline• Review last class• Numerical solution of elliptic equations• Laplace’s equation as an example• Finite difference forms• Iterative solution of algebraic equations• Sample results3Review Properties of Solutions• 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. 4Review Two Space Dimensions• Extension of one space dimension• Have grid in yjas well as xiand 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)2and fy= αΔt/(Δy)2)1,2/1,(2222++=⎥⎦⎤⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂=∂∂nnnsteptimeijyTxTtTα5Finite-difference Grid• Node related to four nearest neighbors at both time stepsunij-1Tnijuni-1junij+1uni+1jΔx ΔxΔyΔyun+1ij-1Tn+1ijun+1i-1jun+1ij+1un+1i+1jΔxΔxΔyΔyOld timeNew time• CN uses all T values at old time step6Review Algorithms• Explicit method()()()niyxnijnijynjinjixnijTffTTfTTfT22111111−−++++=−+−++()()()()nijnijyxnijnijynjinjixnijynjixnijyxnjixnijyRTffTTfTTfTfTfTffTfTf=−−++++=−−+++−−−+−++−+−+++++12121111111111111• Crank-Nicholson• Fully implicitnijnijynjixnijyxnjixnijyTTfTfTffTfTf =−−+++−−++++++−+−111111111)221(Numeric solutions of elliptic PDEs March 25, 2009ME 501B – Engineering Analysis 27Review Explicit Stability()()()niyxnijnijynjinjixnijTffTTfTTfT 22111111−−++++=−+−++•1 –2fx–2fymust be positive for stability so fx+ fymust be less than 1/221)()(22≤ΔΔ+ΔΔ=+ytxtffyxαα•To cut Δxand Δy by a factor of 2 we must cut Δt by a factor of 4• Work increases by a factor of 168Alternating Direction Implicit • ADI is technique for allowing use of Thomas algorithm for simple solution– At each time step, treat one direction as explicit and one as implicit– Even number of time steps required⎥⎥⎦⎤⎢⎢⎣⎡Δ−++Δ−+α=Δ−−+−+2112**1*1*)(2)(2yuuuxuuutuunijnijnijijjijinijij⎥⎥⎦⎤⎢⎢⎣⎡Δ−++Δ−+=Δ−++−++−++2221212**1*1*2)(2)(2yuuuxuuutuunijnijnijijjijiijnijα9Elliptic Equations• 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)(022spaceFuu =∇=∇• Laplace equation in two-dimensional Cartesian coordinates02222=∂∂+∂∂yuxu102D Cartesian Laplace• Use finite-difference grid where xi= x0+ iΔx and yj= y0+ jΔy–x0≤ xi≤ xNand 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• South j = 0, 0 ≤ i ≤ N• North j = M, 0 ≤ i ≤ N• West i = 0, 0 ≤ j ≤ M• East i = N, 0 ≤ j ≤ M),(jiijyxuu =• Dependent variable, u11Finite Differences• Use second-order expressions for second derivatives])[()(2221122xOxuuuxuijjijiijΔ+Δ−+=∂∂−+])[()(2221122yOyuuuyuijijijijΔ+Δ−+=∂∂−+0])(,)[()(2)(2222112112222=ΔΔ+Δ−++Δ−+=⎥⎦⎤⎢⎣⎡∂∂+∂∂−+−+yxOyuuuxuuuyuxuijijijijjijiijyxΔΔ=β• Multiply by (Δx)212Finite-difference Equation• Links central node to four nearest neighbors()()012211211=+−+++−+−+ ijijijjijiuuuuuββuij-1uijui-1juij+1ui+1jΔx ΔxΔyΔy•Only parameter: β = Δx/Δy•For β = 1, 2(1 + β2) = 4Numeric solutions of elliptic PDEs March 25, 2009ME 501B – Engineering Analysis 313Finite-difference Equation• For 2D Cartesian Laplace equation with β = Δx/Δy= 1, uijis the average of its four nearest neighbors41111 −+−++++=ijijjijiijuuuuu• Consider Dirichlet boundary conditions–uijknown at all boundary nodes– Need to find (N – 1)(M – 1) unknown values of uijat central nodes14General Equation• Most general case has five diagonals, but can have different terms– Occurs with variable properties• Notation Aijpointis general coefficient– ij refers to a particular node– Point = N(orth), S(outh), E(ast), W(est) refers to neighboring nodes by direction• General equation shown belowijijNijjiEijijPijjiWijijSijbuAuAuAuAuA =++++++−− 11110=+∂∂∂∂+∂∂∂∂QyTkyxTkx&15Solving the Equations• Typically have large number of equations forming sparse matrix–For Δx= Δy = .01 have 992equations so matrix has 96x106potential 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 solutions16Iterative Solutions• Simplest examples are Jacobi, Gauss-Seidel, and Successive Over Relaxation• Move from iteration n to iteration n+1• Iteration 0 is initial guess (often all zero)• Solve equation for uijand use this as basis for iterationPijijNijjiEijjiWijijSijPijijijAuAuAuAuAAbu1111 ++−−+++−=1'1'1'1''++−−−−−−=ijNijjiEijjiWijijSijijijuAuAuAuAbu17Iterative Solutions II• Use superscript (n) for iteration number• Jacobi iteration uses all old values)(1')(1')(1')(1'')1( nijNijnjiEijnjiWijnijSijijnijuAuAuAuAbu++−−+−−−−=• Gauss–Seidel uses most-recent values)(1')(1')1(1')1(1'')1( nijNijnjiEijnjiWijnijSijijnijuAuAuAuAbu+++−+−+−−−−=• Relaxation basis: Gauss–Seidel provides a correction that can be adjusted[])(),1()()1( nijGSnijnijnijuuuu −+=++ωRelaxation Factor18Relaxation Methods• 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 well-behaved problems[]()()[]')(1')(1')1(1')1(1')(),1()()(),1()()1(11ijnijNijnjiEijnjiWijnijSijnijGSnijnijnijGSnijnijnijbuAuAuAuAuuuuuuu−+++−−=+−=−+=+++−+−+++ωωωωωNumeric solutions of elliptic PDEs March 25, 2009ME 501B – Engineering Analysis 419Excel Relaxation Code I• Gauss-Seidel solution to Poisson equation with constant source