The properties of an implicit three-level finite difference scheme are   investigated.  The modified equivalent partial differential equation   is used to determine the speed and amplitude characteristics of the   scheme, and these are used to determine the optimal value of the   weighting between a central difference and an upwind biased   difference.  Results of numerical experiments that confirm that the   predicted value of the weight does minimise the error are presented

May, Robert L.

Polanco, Frank G.

English

Australian Mathematical Society (AustMS): E-Journals

ANZIAM J. 42 (E) ppC1179–C1198, 2000 C1179An implicit finite difference approximation tothe one-dimensional transport equationFrank G. Polanco∗ Robert L. May†(Received 7 August 2000)AbstractThe properties of an implicit three-level finite difference schemeare investigated. The modified equivalent partial differential equa-tion is used to determine the speed and amplitude characteristics ofthe scheme, and these are used to determine the optimal value of the∗ Aeronautical and Maritime Research Laboratory, GPO Box 4331, Melbourne,Victoria 3000, Australia. mailto:Frank.Polanco@dsto.defence.gov.au† Department of Mathematics, RMIT University, GPO Box 2476V, Melbourne,Victoria 3000, Australia. mailto:rob@rmit.edu.au0See http://anziamj.austms.org.au/V42/CTAC99/Pola for this article and ancillaryservices, c© Austral. Mathematical Soc. 2000. Published 27 Nov 2000.Contents C1180weighting between a central difference and an upwind biased differ-ence. Results of numerical experiments that confirm that the pre-dicted value of the weight does minimise the error are presented.Contents1 Introduction C11802 Finite Difference Approximation C11813 Modified Equivalent PDE C11834 Numerical Experiments C11875 Conclusion C1197References C11981 IntroductionConvection and diffusion occur simultaneously in many physical situations,so it is not surprising that the accurate modelling of these processes hasbeen the subject of much research over the last few decades [1]. We consider2 Finite Difference Approximation C1181a finite difference (fd) approximation to the one-dimensional transport (orconvection-diffusion) equation∂T∂t+ C ∂T∂x= D∂2T∂x2, (1)where t is time, x is space, C and D are the convective velocity and thediffusive coefficient respectively, and T is some scalar quantity (temperaturefor example). We restrict the convective velocity to be non-negative (C ≥ 0),and obviously the diffusion coefficient is also non-negative (D ≥ 0).2 Finite Difference ApproximationThe fd scheme is implicit and second order in space and time; the diffusiveterm is approximated by a central difference and the convective term bya linear combination between a central and an upwind biased difference.A three-level approximation employing a variable time step is used for thetemporal derivative. The fd approximation is[(1 + 2κ)T(n)j − (1 + κ)2T (n−1)j + κ2T (n−2)j(1 + κ)∆t(n)]+ C[(1 − ω)(T(n)j+1 − T (n)j−12∆x)+ω(3T(n)j − 4T (n)j−1 + T (n)j−22∆x)]= D(T(n)j+1 − 2T (n)j + T (n)j−1(∆x)2), (2)2 Finite Difference Approximation C1182where T(n)j is the approximation of T at time tn and x = xj , tn = tn−1+∆t(n),xj = xj−1 + ∆x and ∆t(n) = κ(n)∆t(n−1). The superscript on κ is not shown;we understand that κ = κ(n). Note that κ < 1 and κ > 1 respectively yielda decrease or an increase in the time step. The weight ω, which controls theamount of upwind bias, is assumed to be non-negative. Choosing ω = 0 orω = 1 approximates the convective term by a central difference or an upwindbiased difference respectively.The fd scheme blends an upwind biased scheme with a central scheme.Thus it may be of advantage in problems with time-dependent convective (ordiffusive) terms, where the flow varies from being convection dominated todiffusive dominated. A reduction in the number of nodes needed to evaluatethe fd scheme was the motivation behind deriving Equation (2). The aboveimplicit scheme requires only 3n + 3 nodes in n-dimensions compared witha Crank-Nicolson type of scheme (using the same spatial discretisation asEquation (2)) which would require 6n + 2 nodes. The variable time stepovercomes the limitation of some three-level schemes with fixed time steps.Setting κ = 1 (constant time steps) gives the backward differentiationformula (which is second order in time), while setting κ = 0 gives a schemethat is first order in time. As the finite difference approximation involvesthree time levels if κ 6= 0, it can not be used for the first time step. Anytwo-level scheme could be used, but all results that follow were computedusing the difference equation with κ = 0 for the first time step.3 Modified Equivalent PDE C11833 Modified Equivalent PDEWe use the modified equivalent partial differential equation (mepde) to in-vestigate the numerical diffusive and convective characteristics of the fdscheme. The mepde technique was first outlined by Warming and Hyett [5],and later improved Noye and Hayman [3]. The mepde technique involvesre-writing the Taylor series expansion of the fd scheme (or the equivalentpde) with all temporal derivatives eliminated using the equivalent pde. Theequivalent pde of the fd scheme given by Equation (2) is∂T∂t+ C ∂T∂x−D∂2T∂x2+ TE = 0,where TE is the truncation error of the fd scheme. The truncation errorfrom the equivalent pde is given byTE =∆x212[2C(1 − 3ω)∂3T∂x3−D∂4T∂x4]− ∆t26(1 + κ)κ∂3T∂t3+ O(∆x)3 + O(∆t)3.Note that the above expression is only valid for κ 6= 0, since κ = 0 results in afirst order approximation in time. It would appear that taking ω = 13wouldgive the highest accuracy as this makes the approximation to the convectiveterm third order. However, we shall see that this is not the best choice.Using the above equivalent pde to eliminate temporal derivatives yields3 Modified Equivalent PDE C1184the truncation error from the mepdeTE =∆t26(1 + κ)κ[C3 ∂3T∂x3− 3C2D∂4T∂x4+ 3CD2∂5T∂x5−D3∂6T∂x6]+∆x212[2C(1 − 3ω)∂3T∂x3−D∂4T∂x4]+ O(∆x)3 + O(∆t)3 + O(∆x)2O(∆t)2. (3)The above expression may be written in the more compact formTE ={∆t26(1 + κ)κ(C − D ∂∂x)3+∆x212[2(1 − 3ω)C − D ∂∂x]}∂3T∂x3+ O(∆x)3 + O(∆t)3 + O(∆x)2O(∆t)2.Warming and Hyett [5] were the first authors to note the relation betweenthe mepde’s truncation error and the artificial damping and phase shiftingproperties of the numerical solution. Noye [2, pp. 242–247] formalised therelation and later improved the presentation [4, pp. 236–237] that we nowfollow.Let the truncation of the mepde be given byTE =∞∑k=3C(∆x)k−1 ηkk!∂kT∂xk. (4)3 Modified Equivalent PDE C1185The above series begins at k = 3 since the scheme we have developed issecond order in time and space (unless κ = 0).Consider the amplitude and speed of an infinite travelling wave of unitamplitude propagating according to both the transport equation (1), whichwe term the true solution, and fd equation (2), which we term the numericalsolution. The ratio of the numerical to true wave amplitudes and the ratio ofthe numerical to true wave speed after one period are collectively known asthe wave propagation characteristics, and denoted by ζ and µ respectively.Noye [4] provides the following two expressionsζ = exp2π ∑k=4(2)∞(−1)(k/2)+1(2πNλ)k−1ηkk!andµ = 1 +∑k=3(2)∞(−1)(k−1)/2(2πNλ)k−1ηkk!,which relate the truncation error coefficients to the wave amplitude and speedpropagation characteristics respectively. The parameter Nλ = λ/∆x definesthe number of grid spacings in one wavelength λ.Using the truncation error given by Equation (3), the truncation error3 Modified Equivalent PDE C1186coefficients η (as defined by Equation (4)) areη33!= [1 − 3ω + c2(1 + κ)/κ] /6, η44!= − [s/c + 6cs(1 + κ)/κ] /12,η55!= s2(1 + κ)/(2κ), andη66!= −s3(1 + κ)/(6cκ),where the Courant number, diffusion number, and Peclet number are definedasc = C ∆t∆x, s = D ∆t∆x2, and Pe =cs=CD∆xrespectively.Then the wave amplitude characteristic isζ = exp{π6(2πNλ)3sc[1 + 6(1 + κκ)c2 − 2(2πNλ)2(1 + κκ)s2 + · · ·]},while the wave speed characteristic isµ = 1 − 16(2πNλ)2 [1 − 3ω +(1 + κκ)c2 − 3(2πNλ)2(1 + κκ)s2 + · · ·].As the time step ratio κ increases, both the wave amplitude and speedcharacteristics improve marginally. However, as κ increases so does the re-4 Numerical Experiments C1187sulting time step and associated truncation error. The amplitude character-istic does not depend on ω, but choosingω = ωopt =13[1 +(1 + κκ)c2], (5)results in cancellation of the first 3 terms in the square brackets in the expres-sion for the wave speed characteristic. For small c the optimal weight ωoptis approximately 13. As expected, for convection dominated flows (c large)the wave speed characteristics improve as ω increases, since more weightingis given to the upwind biased difference than the central difference.4 Numerical ExperimentsIt is easily verified thatT (x, t) =0.1√0.01 + 4Dt exp[−(x −12− Ct)20.01 + 4Dt](6)is an exact solution of the transport equation. At t = 0 it gives a narrowGaussian pulse of unit height centred on x = 12(see Figure 1). This solutionis used to examine the accuracy of the fd approximation (2).As time passes the Gaussian pulse moves to the right with velocity C, atthe same time becoming shorter and broader due to diffusion. Without loss4 Numerical Experiments C1188of generality we take C = 1, so the pulse will be centred on x = 3.5 at timet = 3. Thus an interval of at least [0, 5] would be required for large valuesof D. The fd approximation (2) is implicit, and requires that a system oflinear equations is solved for each time step. To avoid the large interval andconsequent large number of linear equations, a unit interval is used with thecyclic boundary conditionT (1, t) = T (0, t).Thus the pulse enters the left boundary as it leaves the right boundary, as ifthe unit interval was bent into a circle. An alternative interpretation is thatwe start with an infinite train of pulses spaced a unit distance apart, and wewatch through the “window” [0, 1] as they move to the right.The use of the cyclic boundary condition does reduce the size of thesystem of linear equations to be solved at each time step, but unfortunatelythe coefficient matrix is no longer banded. In addition to the band of width4, it has 4 extra non-zero elements; 1 in the lower left corner and 3 in thetop right corner (in the uppermost 2 diagonals). An LU factorisation of thecoefficient matrix is used to solve the linear equations, but instead of therebeing non-zero elements in only 3 diagonals of L and 2 diagonals of U , the“extra” elements require that there are also non-zero elements in the lastrow of L and last two columns of U . An algorithm has been constructed sothat only the non-zero elements of L and U are calculated. The factorisationmust be found at each step unless κ = 1 for consecutive steps.The initial values T(0)j = T (xj , 0) are obtained from Equation (6). As4 Numerical Experiments C1189mentioned earlier, the first step is carried out using the fd approximationwith κ = 0. Subsequent steps are calculated using a fixed value of κ (usually1), and at specified times the root mean square error is calculated.A grid with ∆x = 140and a convective velocity C = 1 was used forall calculations. Unless stated otherwise, the time step was kept constant(κ = 1). The Courant number c is determined by the choice of the time step∆t, and the diffusion number s (and hence the Peclet number) determinedby the value of D.Figure 1 shows a comparison of the exact solution and the numericalsolution. The values ∆t = 0.01 and D = 0.008 were used, resulting in aCourant number of c = 0.4 and a Peclet number of Pe = 3.125. The value ofω given by Equation (5) (that is ωopt) was used for the calculation.Initially the Gaussian pulse is at the centre of the interval. At t = 0.5 thepulse has moved a distance of 0.5, so half of the pulse has disappeared outthe right side of the interval and re-entered the left side of the interval. Analternative, and perhaps simpler interpretation, is that we can see the lefthalf of the pulse at the right end of the interval and can see the right half ofthe pulse which follows at the left. Due to the diffusive effects, their maximaare about 0.6 and they are correspondingly broader (the integral of T overthe interval should be conserved). At t = 1 the pulses have diffused furtherand moved a distance of 1, so the following pulse is centred on x = 12. Thefinal plot is at t = 3; the pulse that was the third behind the initial pulse iscentred on x = 12, and due to the large amount of diffusion the pulses have4 Numerical Experiments C1190t = 0 t = 0.5T0.2 0.4 0.6 0.8 10.20.40.60.810.2 0.4 0.6 0.8 10.20.40.60.81x xt = 1 t = 3T0.2 0.4 0.6 0.8 10.20.40.60.810.2 0.4 0.6 0.8 10.20.40.60.81x xFigure 1: The exact solution (solid line) and numerical solution (dots) forc = 0.4 and Pe = 3.125 at the indicated times. The numerical solution wascalculated using ω = ωopt = 0.44 .4 Numerical Experiments C1191Error0.1 0.2 0.3 0.4 0.50.020.040.060.08ωFigure 2: Error dependence on ω for c = 0.1 using Pe = 25 and Pe = 2.5(dots and line-plus-dots respectively). The value of ωopt is indicated by thedashed line.4 Numerical Experiments C1192“coalesced”. The numerical solution is clearly in excellent agreement withthe exact solution at all times.To determine the best value of the weight ω the numerical solution andthe norm of its error were calculated with ∆t = 0.0025 (c = 0.1) for ω =0(0.025)0.5. Two values of D, corresponding to Pe = 2.5 and Pe = 25, wereused. The results are shown in Figure 2. For both values of the Pecletnumber, the minimum error occurs for a value of ω extremely close to thevalue of ωopt defined by Equation (5). The same outcome was produced whentime steps of ∆t = 0.01 and ∆t = 0.03 (c = 0.4 and c = 1.2) were used (seeFigures 3 and 4). Note that in the latter case the optimal weight is greaterthan 1.It is apparent from Figures 2, 3 and 4 that the error is much smallerfor the smaller Peclet number for any value of ω. This is not surprisingas increasing the rate of diffusion lowers the Peclet number. ComparingFigures 2, 3 and 4 shows that the errors corresponding to the optimal valueof ω increase with the Courant number. However, the results are still quiteacceptable for c = 1.2, particularly for the smaller Peclet number. It shouldbe noted that the best way to improve the accuracy for given values of Cand D is to reduce ∆x. Although this increases the Courant number, itreduces the value of the Peclet number and the Peclet number has a strongerinfluence on the error.It was observed that the numerical wave speed was too low if ω < ωopt,and too high if ω > ωopt. This can be seen in Figure 5, where the solution at4 Numerical Experiments C1193Error0.2 0.4 0.6 0.8 10.020.040.060.080.1ωFigure 3: Error dependence on ω for c = 0.4 using Pe = 25 and Pe = 2.5(dots and line-plus-dots respectively). The value of ωopt is indicated by thedashed line.4 Numerical Experiments C1194Error0.5 1 1.5 20.0250.050.0750.10.1250.15ωFigure 4: Error dependence on ω for c = 1.2 using Pe = 25 and Pe = 2.5(dots and line-plus-dots respectively). The value of ωopt is indicated by thedashed line.4 Numerical Experiments C1195ω = 0.2 ω = ωopt = 0.44T0.2 0.4 0.6 0.8 1-0.20.20.40.60.80.2 0.4 0.6 0.8 1-0.20.20.40.60.8x xFigure 5: The exact solution (solid line) and numerical solution (dots) attime t = 3 for c = 0.4 and Pe = 25. The numerical solution was calculatedusing two values of ω.4 Numerical Experiments C1196T0.2 0.4 0.6 0.8 10.20.40.60.81xFigure 6: The numerical solution for for T at Pe = 3.125 calculated withan increasing time step (κ = 1.01). The optimal value of ω was used at eachstep.t = 3 has been calculated using two values of ω. A time step of ∆t = 0.01and diffusion coefficient D = 0.001 were used, giving the Courant numberof c = 0.4 and Peclet number Pe = 25. The numerical solution calculatedwith ω = 0.2 clearly lags the exact solution; the numerical wave speed is toolow. When the optimal value ω = 0.44 is used, the numerical wave speed isreasonably accurate. However, in both cases the peak value is too low; thenumerical solution has suffered from slightly too much diffusion.5 Conclusion C1197The numerical solution shown in Figure 6 was calculated using a variabletime step. The initial time step was ∆t = 0.01, and at each time step itwas increased by the factor κ = 1.01 resulting in t = 3 being reached in 139steps with a final time step of 0.0395. The Peclet number was 3.125, and theCourant number varied from 0.4 to about 1.58. The optimal value of ω wascalculated at each time step. As expected, the solution is not as accurateas the one obtained with a constant time step of 0.01 (see Figure 1); theerror is approximately 3 times greater. It is also more accurate to take thesame number of steps of constant length, but in practice one would only usea value of κ > 1 for a few steps.5 ConclusionWe have investigated the accuracy of an implicit fd approximation of thetransport equation and found that the highest accuracy is obtained not bychoosing the weight that maximises the order of the approximation to theconvective term, but by choosing the value that makes the numerical wavespeed as close as possible to the true wave speed. Although we have only con-sidered one spatial dimension and a constant convective velocity, the methodcould be adapted to two or more dimensions and a variable velocity by cal-culating the optimal value of the weight at each grid point.References C1198References[1] K. W. Morton. Numerical Solution of Convection-Diffusion Problems.Chapman & Hall, 1996. C1180[2] J. Noye. Finite Difference Techniques for Partial Differential Equations.Computational Techniques for Differential Equations, pages 95–354.North-Holland, 1984. C1184[3] J. Noye and K. Hayman. Accurate Finite Difference Methods forSolving the Advection-Diffusion Equation. In J. Noye and R. May,editors, Computational Techniques and Applications: CTAC-85, pages137–158. Elsevier Science, 1986. C1183[4] J. Noye. Finite Difference Methods for Solving the One-DimensionalTransport Equation. Numerical Modelling: Applications to MarineSystems, pages 231–256. North-Holland, 1987. C1184, C1185[5] R. Warming and B. Hyett. The Modified Equation Approach to theStability and Accuracy of Finite-Difference Methods. J. ComputationalPhysics, 14:159–179, 1974. C1183, C1184

An implicit finite difference approximation to the one-dimensional transport equation

https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/download/641/509

An implicit finite difference approximation to the one-dimensional transport equation

Abstract

Similar works

Full text

Available Versions

Australian Mathematical Society (AustMS): E-Journals