More numerical elliptic PDEs March 30, 2009ME 501B – Engineering Analysis1Additional Topics in Numerical Additional Topics in Numerical Solutions of Elliptic EquationsSolutions of Elliptic EquationsLarry CarettoMechanical Engineering 501BSeminar in Engineering AnalysisMarch 30, 20092Outline• Review last class• Treatment of boundary conditions– Second kind (Neumann)– Third kind a∂u/∂s+ bu= c• s is coordinate normal to boundary (x or y)• Treatment using different gradient expressions• Compact difference expressions for higher order accuracy3Review Finite Differences• Second-order second derivatives])[()(2221122xOxuuuxuijjijiijΔ+Δ−+=∂∂−+ε+Δ−+=∂∂−+21122)(2yuuuyuijijijij0)(2)(22112112222=⎟⎟⎠⎞⎜⎜⎝⎛+Δ−++Δ−+≈⎥⎥⎦⎤⎢⎢⎣⎡+∂∂+∂∂−+−+ijijijijijjijiijkQyuuuxuuukQyuxu&&• Poisson-Laplace equation O[(Δx)2,(Δy)2]()()()0122211211=⎟⎟⎠⎞⎜⎜⎝⎛Δ+β+−+β++−+−+ijijijijjijikQxuuuuu&• Multiply by (Δx)2and define β = Δx/Δy4Review Small Grid (N = 6, M = 5)i = 0 i = 1 i = 2 i = 3 i = 4 i = 5 i = 6j = 5j = 4j = 3j = 2j = 1j = 0u43u01u11u10u21u12Boundarynodesu51u44u53u42u33Compu-tationalMoleculeu44+ u33–4u43+ u53+ u44= 0Review N = 6, M = 5 Matrix⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−−−−−−−−−−−−−−−54443424145343332313524232221251413121114111411141114111411411114111141111411114114111141111411114111141141114111411141114uuuuuuuuuuuuuuuuuuuuZero coefficients because of x boundaryu44+ u33+ -4u43+ u53+ u44= 06More numerical elliptic PDEs March 30, 2009ME 501B – Engineering Analysis27Review Iterative Solutions• 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()[]')(1')(1')1(1')1(1')()1(1ijnijNijnjiEijnjiWijnijSijnijnijbuAuAuAuAuu−+++ω−ω−=+++−+−+8Effect of Relaxation Factor on Execution Time0.010.11101.4 1.5 1.6 1.7 1.8 1.9 2Relaxation FactorExecution Time (seconds16x16 grid32x32 grid64x64 grid128x128 gridOther codeSquare ( L= H = 1)64 by 64 GridZero boundary on left, right and bottomTop boundary has u(x,H) = sin(πx/L)"Other" is different code with u(x,H) = 1⎟⎟⎠⎞⎜⎜⎝⎛−−=ξξω112opt()()2221coscos⎥⎥⎦⎤⎢⎢⎣⎡β+πβ+π=ξMN9Review Errors())1()()1()1(+++−=nijnijnijnijuuuChangeRelative()')1(1')1(1')1(1')1(1')1()1(ijnijNijnjiEijnjiWijnijSijnijnijbuAuAuAuAu−++++−=+++++−+−++Residual())1()()1(+∞+−=nijijnijuuErrorIteration()())1()1(,++−=nijjiPDEnijuyxuErrorExactComputable10Effects of Iterations on Laplace Equation Errors1.E-151.E-141.E-131.E-121.E-111.E-101.E-091.E-081.E-071.E-061.E-051.E-041.E-031.E-021.E-011.E+000 100 200 300 400 500 600 700IterationsErrorsDifferenceResidualIteration errorExact ErrorSquare ( L= H = 1)64 by 64 GridZero boundary on left, right and bottomTop boundary has u(x,H) = sin(πx/L)11Other Boundary Conditions• General condition a∂u/∂s+ bu= c– a = 0, b = 1 for Dirichlet (value given)– a = 1, b = 0 for Neumann (gradient given)– Mixed has both a and b nonzero• Write general boundary condition using a finite difference expression for ∂u/∂s– Two approaches• Using second order forward or backward difference for boundary node• Add fictitious node outside boundary and use central differences at boundary12Notation• Label the boundaries at x = x0and x = xNas the (W)est and (E)ast boundaries– a, b, and c can be different at each yj– Notation: ajW, bjW, cjW, ajE, bjE, and cjE• Boundaries at y = y0and y = yMare the (S)outh and (N)orth boundaries– a, b, and c can be different at each xi– Notation: aiS, biS, ciS, aiN, biN, and ciNMore numerical elliptic PDEs March 30, 2009ME 501B – Engineering Analysis313Central Difference Approach• General condition a∂u/∂s + bu = c by central differences– Write general equation for boundary nodes that contains fictitious node outside region– Write central difference equation for boundary condition– Solve this equation for potential at fictitious node and use result to replace this value in general equation at boundary– Include boundary nodes in iterations14(N)orth Boundary Example• Finite-difference equation at y = yMiMiMNiMMiEiMiMPiMMiWiMiMSiMQuAuAuAuAuA =++++++−− 1111• Fictitious uiM+1from central-difference boundary condition equation at y = yMNiiMNiiMiMNicubyuua =+Δ−−+211()iMNiNiNiiMiMubcayuu −Δ+=−+21115(N)orth Boundary Example II• Combine equations to eliminate uiM+1()iMiMNiNiNiiMNiMMiEiMiMPiMMiWiMiMSiMQubcayuAuAuAuAuA=⎥⎦⎤⎢⎣⎡−Δ+++++−+−−21111• Modified coefficients at (iterated) boundaries: new AS, AP, and Q; AN= 0; AWand AEmultiplied by aN()()NiNiMiMNiMiEiMNiiMNiNiMPiMNiMiWiMNiiMNiMSiMNicyAQauAaubyAAauAauAAaΔ−=+Δ−++++−−2211116CD General Boundary• Modified east boundary coefficients()022←←←Δ−←+←Δ−←ENjSNjEjSNjNNjEjNNjWNjEjNjEjNjWNjENjEjWNjEjENjPNjEjPNjAAaAAaAAxcQaQAAaAxbAAaA• Modified north boundary coefficients()022←←←Δ−←+←Δ−←NiMWiMNiWiMEiMNiEiMNiMNiiMNiiMSiMNiMNiSiMNiNiMPiMNiPiMAAaAAaAAycQaQAAaAybAAaA17CD General Boundary II• Modified west boundary coefficients()02200000000000000←←←Δ+←+←Δ+←WjSjWjSjNjWjNjWjWjjWjjWjEjWjEjWjWjPjWjPjAAaAAaAAxcQaQAAaAxbAAaA• Modified south boundary coefficients()02200000000000000←←←Δ+←+←Δ+←SiWiSiWiEiSiEiSiSiiSiiSiNiSiNiSiSiPiSiPiAAaAAaAAycQaQAAaAybAAaA18Forward/Backward Difference• General condition, a∂u/∂s+ bu= c– Use forward differences at i = 0 or j = 0 and backward differences at i = N or j = M– Obtain equation for boundary potential in terms of two nodes in from boundary– Combine this equation with general equation for first node in from the boundary to eliminate unknown boundary potential– No iteration on boundary values, which are found at end of iterations– Derivation at end of this presentationMore numerical elliptic PDEs March 30, 2009ME 501B – Engineering