Numerical Diffusion Equation March 18, 2009ME 501B – Engineering Analysis 1Numerical Solutions of the Numerical Solutions of the Diffusion EquationDiffusion EquationLarry CarettoMechanical Engineering 501ABSeminar in Engineering AnalysisMarch 18, 20092Outline• Review last class• Numerical solutions of the diffusion equation in one space dimension– Explicit algorithm– Stability of algorithms– Crank-Nicholson algorithm– (Fully) implicit algorithm– DuFort-Frankel algorithm3Review Numerical Analysis• Transform differential equation into a system of algebraic equations• Obtain solution for discrete points in domain• Two basic approaches: finite differences and finite elements• Start with finite elements• Get expressions for derivatives and measure of error with their use4Review Finite Difference Grids• Grid notation for four independent variables: x, y, z and tx0= xminxN= xmaxxi–xi-1= Δxiy0= ymjnyM= ymaxyj–yj-1= Δyjz0= zmknzK= zmaxzk–zk-1= Δzkt0= tmintL= tmaxtn–tn-1= Δtn• Dependent variable u(x,y,z,t) at discrete points u(xi, yj, zk, tn)• Use notation below for this value of u),,,(nkjinijktzyxuu =5Review Truncation Error• If we truncate series after m terms∑∑∞+====+=10)-(!1)-(!1)(mnnaxnnmnnaxnnaxdxfdnaxdxfdnxfTerms used Truncation error, εm• Can write truncation error as single term at unknown location (derivation based on the theorem of the mean)1111)-(!)1(1)-(!1+=++∞+==+==∑mxmmmnnaxnnmaxdxfdmaxdxfdnξε6Review Order of the Error• Derivative expressions have error that is proportional to hn• This power, n, is called the order of the error• Use notation O(hn) to indicate this error• Reducing step size by a factor of a reduces nthorder error by annhh⎟⎟⎠⎞⎜⎜⎝⎛≈1212εεNumerical Diffusion Equation March 18, 2009ME 501B – Engineering Analysis 27Review Derivative Expressions• First-order error, first derivatives)(2211'hOhfffiii+−=−+)(1'hOhfffiii+−=+)(1'hOhfffiii+−=−• Second-order error, first derivatives)(234212'hOhffffiiii+−+−=++)(234212'hOhffffiiii++−=−−• Second derivative()2211''2hOhffffiiii+−+=−+8Find f’ and f’’ for sin at x = 1)(2211'hOhfffiii+−=−+Second order central()2211''2hOhffffiiii+−+=−+)1(.2)1.1sin()1.1sin('−−+=if2)1(.)1sin(2)1.1sin()1.1sin(''−−++=if)01(.2)01.1sin()01.1sin('−−+=if)001(.2)001.1sin()001.1sin('−−+=if2)01(.)1sin(2)01.1sin()01.1sin(''−−++=if2)001(.)1sin(2)001.1sin()001.1sin(''−−++=if9Review Roundoff Error• Possible in derivative expressions from subtracting close differences• Example f(x) = ex: f’(x) ≈ (ex+h–ex-h)/(2h) and error at x = 1 is (e1+h–e1-h)/(2h) - e3105.4718282.2)1.0(2722815.2004166.3−=−−= xE9105.4597182818284.2)0001.0(27180100139.27185536702.2−=−−= xE9109.5718281828.2)0000001.0(203887182815566.287247182821002.2−=−−= xESecond order errorFigure 2-1. Effect of Step Size on Error1.E-111.E-101.E-091.E-081.E-071.E-061.E-051.E-041.E-031.E-021.E-011.E+001.E+011.E-17 1.E-15 1.E-13 1.E-11 1.E-09 1.E-07 1.E-05 1.E-03 1.E-01Step SizeErrorSlope of a log error – log h plot is order nΔ(log ε) = 10; Δ(log h) = 5Slope = order = 211Numerical PDE Solutions• Define a finite-difference grid in the independent variables (x, y, z, t)• Place grid points on region boundary whose values are found from boundary conditions for the problem• At some grid location convert differential equation into a finite difference equation– Observe truncation error in process– Neglect truncation error to get set of algebraic equations to solve12Diffusion Equation• Apply difference formulas derived for ordinary derivatives to partial derivatives• Use notation to consider different coordinate directions• Apply to diffusion equation•Grids xi= x0+ iΔx and tn= t0+ nΔt• Try finite difference expressions below to get simple finite-difference equationnixutu⎥⎦⎤⎢⎣⎡∂∂=∂∂22α])[()(2)(2211221xOxuuuxuandtOtuutunininininininiΔ+Δ−+=∂∂Δ+Δ−=∂∂−++Numerical Diffusion Equation March 18, 2009ME 501B – Engineering Analysis 313Diffusion Equation II• Substitute finite difference expressions into differential equation])(,[)(222111xtOxuuutuunininininiΔΔ+Δ−+=Δ−−++α• Ignore truncation error, solve for uin+1()ninininiuxtuuxtu⎟⎟⎠⎞⎜⎜⎝⎛ΔΔ−++ΔΔ=−++21121)(21)(αα• Obtain potential at x = xiand t = tn+1in terms of u values at old time step14Explicit (FTCS) Method() ()()nininininininiufuufuxtuuxtu 21)(21)(1121121−++=⎟⎟⎠⎞⎜⎜⎝⎛ΔΔ−++ΔΔ=−+−++αα• Method just derived is called explicit method; can solve one equation at a timeuni-1 uinuni+1●----------●----------●● uin+1•uin+1 does not depend on other u values at the new time step (n+1)2)( xtfΔΔ≡αxt15Explicit Method Example• Pick α = 1, Δx = 0.25, Nx= 4, Δt = 0.01•f = αΔt/(Δx)2= 1(.01)/(.25)2= 0.16• Pick initial ui0= 1000 and boundaries, un0= un4= 0 for time > 0 (n ≥ 0)()()ninininiufuufuApply 21111−++=−++[][][]840]1000[68.0]01000[16.0)21(1000]1000[68.0]10001000[16.0)21(840]1000[68.0]10000[16.0)21(030402130203011201020011=++=−++==++=−++==++=−++=ufuufuufuufuufuufu• Repeat for subsequent time steps 16Explicit Method Results f = 0.160388.7548.9388.70t = 0.08n = 80429.2605.3429.20t = 0.07n = 70474.2667474.20t = 0.06n = 60524.6734524.60t = 0.05n = 50582805.55820t = 0.04n = 40649879.26490t = 0.03n = 30731.2948.8731.20t = 0.02n = 2084010008400t = 0.01n = 101000100010000t = 0+n = 010001000100010001000t = 0x = 1.00x = 0.75x = 0.50x = 0.25x = 0.00i = 4i = 3i = 2i = 1i = 017Explicit Method Results f = 0.1605.88.25.80t = 0.20Error0125.1176.9125.10t = 0.20Exact0119.2168.6119.20t = 0.20n = 200131.6186.1131.60t = 0.19n = 190145.2205.3145.20t = 0.18n = 180160.2226.5160.20t = 0.17n = 170176.8250176.80t = 0.16n = 160195275.81950t = 0.15n = 150215.2304.4215.20t = 0.14n = 140237.5335.8237.50t = 0.13n = 130262370.52620t = 0.12n = 12x = 1.00x = 0.75x = 0.50x = 0.25x = 0.00i = 4i = 3i = 2i = 1i = 018Explicit Results f = 0.3201825.4180t = 0.20Error0125.1176.9125.10t = 0.20Exact0107.1151.4107.10t = 0.20n = 100131.8186.4131.80t = 0.18n = 90162.2229.4162.20t = 0.16n = 80457.9647.7457.90t = 0.06n = 30564.8795.2564.80t = 0.04n = 2068010006800t = 0.02n = 101000100010000t = 0+n = 010001000100010001000t = 0x = 1.00x = 0.75x = 0.50x = 0.25x = 0.00i = 4i = 3i = 2i = 1i = 0Numerical Diffusion Equation March 18, 2009ME 501B – Engineering Analysis