- Email: [email protected]
The exponential method of angled derivatives
Applied Mathematics and Computation 177 (2006) 665–685 www.elsevier.com/locate/amc
The exponential method of angled derivatives Brian J. McCartin Applied Mathematics, Kettering University, 1700 West Third Avenue, Flint, MI 48504-4898, United States
Abstract A hybrid numerical method is presented for a linear, first order, hyperbolic partial differential equation (PDE): the nonhomogeneous one-way advection-reaction equation. This PDE is first reduced to an ordinary differential equation (ODE) along a characteristic emanating backward in time from each mesh point. These ODEs are then integrated numerically using exponential fitting to advance in time thus permitting wide disparities in size among the coefficients. Such temporal integration requires a spatial interpolation procedure dependent upon the slope of the local characteristic. If the local Courant number is less than one then the interpolation is afforded by the exponential angled derivative method while if greater than one then the exponential reflected angled derivative approximation is employed. Otherwise, if the Courant number is transitioning through one then the exponential box scheme is employed. At the boundary, parabolic tracing of the characteristics provides a comparable level of accuracy. The net result is a three-level explicit method which is second-order accurate and unconditionally stable. The efficacy of this exponential method of angled derivatives (EMAD) scheme for singular perturbation problems is demonstrated numerically. Ó 2005 Elsevier Inc. All rights reserved.
1. Introduction In what follows, attention is focused on the numerical solution of the nonhomogeneous one-way advectionreaction equation aðx; tÞ
ou ou þ bðx; tÞ þ cðx; tÞ u ¼ dðx; tÞ; ot ox
0 < x 6 1; t > 0; aðx; tÞ > 0; bðx; tÞ P 0; cðx; tÞ > 0; ð1Þ
subject to the initial-boundary conditions uðx; 0Þ ¼ /ðxÞ;
uð0; tÞ ¼ wðtÞ.
ð2Þ
For example, this problem arises in chemical engineering in the form of the plug flow reactor and heat exchanger equations [1]. Also, it plays a prominent role in the mathematical modelling of automotive catalytic converters [2]. In the latter application, u(x, t) is either a concentration or a temperature, a(x, t) is a measure of
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.09.097
666
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
porosity, b(x, t) is the convective velocity, c(x, t) is either a mass or heat transfer coefficient, and d(x, t) accounts for sources and sinks [3]. In the ensuing pages, we develop an approximate method for this problem which is nondissipative, secondorder accurate, explicit, and unconditionally stable that allows for wide disparities in the magnitude of the coefficients. We begin by rewriting Eq. (1) in the normalized form: ou ou þ bðx; tÞ þ cðx; tÞ u ¼ f ðx; tÞ; ot ox
ð3Þ
0 < x 6 1; t > 0;
where bðx; tÞ :¼ bðx; tÞ=aðx; tÞ P 0;
cðx; tÞ :¼ cðx; tÞ=aðx; tÞ > 0;
f ðx; tÞ :¼ dðx; tÞ=aðx; tÞ.
ð4Þ
The special case of pure advection, c(x, t) 0, has previously been treated by the method of angled derivatives (MAD) [4]. This MAD scheme is an amalgam of the angled derivative [5], box [6], and reflected angled derivative [4] schemes combined in such a way as to accentuate the best features of each of these component schemes. It is the express intent of the present paper to extend the MAD scheme to permit reaction processes. We begin by showing that straightforward methods of approximating the reaction term lead to undesirable numerical phenomena even in the homogeneous, constant coefficient case. A solution to this dilemma is provided by the application of exponential fitting [7]. It is important to realize that the numerical issues herein raised and resolved are not confined to the three schemes cited above but are in fact germane to any consistent approximation to Eq. (1) such as the Lax–Wendroff, upwind, and leap-frog schemes [8]. Once the need for exponential fitting has been firmly established, we develop the complete exponential method of angled derivatives (EMAD) for nonhomogeneous, variable coefficient problems. This new EMAD scheme is then exercised on a pair of interesting and challenging test problems which arise in the numerical simulation of automotive catalytic converters [3]. 2. The need for exponential fitting In order to firmly establish the need for exponential fitting for advection-reaction problems, we confine our attention, in this section only, to the homogeneous, constant coefficient problem, ou ou þ b þ c u ¼ 0; ot ox
0 < x 6 1; t > 0; b P 0; c > 0;
ð5Þ
subject to the initial-boundary conditions, Eq. (2). For pure advection, c = 0, a comprehensive analysis of the angled derivative (AD), reflected angled derivative (RAD), and box schemes was provided in [4]. The first obstacle which we must herein surmount is the existence of several seemingly natural generalizations of these schemes for the case of advection-reaction, c 5 0. We next consider these alternatives and argue that those based upon exponential fitting [9] possess superior numerical properties. 2.1. Exponential angled derivative (EAD) scheme With reference to Fig. 1, the angled derivative scheme may be derived as follows: 1 huiþ1;nþ1 þ ui;n uiþ1;n þ ui;n1 i uiþ1;n ui;n þ c ^uiþ1=2;n ¼ 0. þb Dt 2 2 Dx
ð6Þ
If we use the average value ^ uiþ1=2;n ðuiþ1;n þ ui;n Þ=2 then the following scheme results: uiþ1;nþ1 ¼ ui;n1 þ ð1 2l mÞ uiþ1;n ð1 2l þ mÞ ui;n ;
ð7Þ
where l :¼ bDt/Dx is the Courant number and m :¼ cDt. For future reference, we also introduce k :¼ cDx/b with the consequent identity l = m/k. We make the further assumption, pffiffiffiffiffiffiffi in this subsection only, that 0 6 l 6 1. Upon insertion of the Fourier mode [10] ui,n = gnejk(iDx) (j ¼ 1) into Eq. (7), we find that the amplification factor, g, satisfies g2 þ ð1 2l þ mÞ ejh ð1 2l mÞ g ejh ¼ 0; ð8Þ
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
667
Δ Δ
Fig. 1. Angled derivative scheme (0 6 l 6 1).
pffiffiffiffiffiffiffiffiffiffiffiffiffi where h :¼ kDx. Setting h = 0 in Eq. (8) and solving, we arrive at g ¼ m 1 þ m2 , so that one of its roots satisfies jgj > 1 for any m 5 0. Thus, Eq. (7) is unconditionally unstable for any m 5 0! A stable approximation to Eq. (5) may be derived by instead using the average value ^ uiþ1=2;n ðuiþ1;nþ1 þ ui;n1 Þ=2 resulting in the scheme: uiþ1;nþ1 ¼
1 ½ð1 mÞ ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ. 1þm
The amplification factor for this scheme satisfies 1 2l 1 m jh g2 þ ðejh 1Þ g e ¼ 0. 1þm 1þm
ð9Þ
ð10Þ
Application of the Schur criterion [11] establishes stability of this scheme for m > 0 and 0 6 l 6 1. However, observe that if m = 1 then one of the roots of Eq. (10) is zero. Furthermore, if, in addition, l = 1/2 then both roots vanish. In light of the fact that the exact amplification factor of a Fourier mode of Eq. (5) satisfies jgexactj = em, Eq. (9) is seen to be excessively dissipative. This difficulty can be alleviated by the use of the method of characteristics [8]. With reference to Fig. 2, dashed characteristics are traced forward from (xi, tn1) to (xT, tn) and backward from (xi+1, tn+1) to (xB, tn). Along a characteristic, Eq. (5) reduces to the ODE
Fig. 2. Characteristic interpretation of AD scheme (0 6 l 6 1).
668
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
du þ c u ¼ 0. dt
ð11Þ
This is numerically integrated from tn1 to tn by the trapezoidal rule to produce uT ¼
1 m=2 ui;n1 . 1 þ m=2
ð12Þ
Parabolic interpolation is then employed at time tn thereby yielding uB ¼ uT þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð13Þ
Finally, Eq. (11) is numerically integrated from tn to tn+1 by the trapezoidal rule to produce uiþ1;nþ1 ¼
1 m=2 uB . 1 þ m=2
Combining Eqs. (12)–(14) yields the angled derivative (AD) scheme 1 m=2 1 m=2 ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ . uiþ1;nþ1 ¼ 1 þ m=2 1 þ m=2
ð14Þ
ð15Þ
The amplification factor for this scheme satisfies 2 1 m=2 1 m=2 ð1 2lÞ ðejh 1Þ g g2 þ ejh ¼ 0; ð16Þ 1 þ m=2 1 þ m=2 thus implying that jgj ¼ 1m=2 Observe that this is the [1/1] Pade´ approximant [12] to jgexactj = em. Conse1þm=2. quently, the AD scheme is essentially nondissipative for small m but is highly dissipative for even moderate values of m. For example, m = 2 implies that the amplification factor vanishes. Yet, we may modify the AD scheme via exponential fitting [9] in such a way that the modulus of the amplification factor is exact. Returning to Eq. (11), we multiply by the integrating factor ect thereby producing d ct ðe uÞ ¼ 0. dt
ð17Þ
Thus, Eq. (12) becomes uT ¼ em ui;n1 ;
ð18Þ
while Eq. (14) becomes uiþ1;nþ1 ¼ em uB .
ð19Þ
Combining Eqs. (18), (13) and (19) yields the exponential angled derivative (EAD) scheme uiþ1;nþ1 ¼ em ½em ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð20Þ
The amplification factor for this scheme satisfies g2 þ ½em ð1 2lÞ ðejh 1Þ g e2m ejh ¼ 0; thus implying that jgj = em = jgexactj. For future reference, 8 9 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < = jh jh 2 1 e 1 e þ ð1 2lÞ gEAD ¼ em ð1 2lÞ þ ejh . : ; 2 2
ð21Þ
ð22Þ
Alternatively, this may be interpreted as a leap-frog approximation centered at (xi+1/2, tn) on an appropriately angled grid. The necessity of exponential fitting is underscored by Fig. 3 where Eqs. (15) and (20) are compared for l = 1, m = 4. Specifically, u(x, t) = ecÆt + ecÆx/b with b = 1, c = 40, Dx = Dt = .1 is displayed at t = 2. For subsequent comparison to the exponential box scheme in the range 0 6 l 6 1, we record here its local truncation error [6,10,13]:
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
669
Steady State: a=.025, b =.025, c = 1.0, d = 0.0 1
x = angled derivative scheme o = exponentially fitted scheme ___ = exact solution
0.5
0
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x Fig. 3. Need for exponential fitting: angled derivative scheme. 3
iþ1=2;n sEAD ¼
ðDxÞ uiþ1=2;n lð2l 1Þð1 lÞ. 12Dt xxx
ð23Þ
A further comment on the dispersive properties of the EAD scheme is appropriate at this junction. Since it is a three-level scheme, it inevitably possesses a ‘‘spurious mode’’ [6] travelling in the wrong direction (westward). However, as long as time-level t = Dt receives propitious treatment (as it will in the EMAD scheme developed below), this ‘‘parasitic solution’’ [13] will not be stimulated. 2.2. Exponential reflected angled derivative (ERAD) scheme Turning now to the range l P 1, perusal of Fig. 4 permits us to derive the reflected angled derivative scheme as follows: ui;nþ1 ui;n b huiþ1;nþ1 þ ui;n ui;nþ1 þ ui1;n i þ ð24Þ þ c ^ui;nþ1=2 ¼ 0. Dx Dt 2 2 If we use the average value ^ ui;nþ1=2 ðui;nþ1 þ ui;n Þ=2 then the following scheme results: uiþ1;nþ1 ¼ ui1;n þ ð1 2=l kÞ ui;nþ1 ð1 2=l þ kÞ ui;n . n jk(iDx)
Upon insertion of the Fourier mode [10] ui,n = g e given by
ð25Þ
into Eq. (25), we find that the amplification factor, g, is
Δ Δ
Fig. 4. Reflected angled derivative scheme (1 6 l 6 1).
670
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
g¼
2 þ lðejh 1 kÞ . 2 þ lðejh 1 þ kÞ
ð26Þ
, so that jgj > 1 for any m > 0 provided that l > 1. Thus, Eq. (25) is Setting h = p in Eq. (26) yields g ¼ 2ð1lÞm 2ð1lÞþm unconditionally unstable for any m > 0! A stable approximation to Eq. (5) may be derived by instead using the average value ^ui;nþ1=2 ðuiþ1;nþ1 þ ui1;n Þ=2 resulting in the scheme: uiþ1;nþ1 ¼
1 ½ð1 kÞ ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ. 1þk
ð27Þ
The amplification factor for this scheme is given by g¼
2 þ l½ð1 kÞejh 1 ; 2 þ l½ð1 þ kÞejh 1
ð28Þ
so that jgj 6 1 for any m > 0 provided that l P 1. The stability of this scheme is thus established for m > 0 and l P 1. However, observe that if k = 1 and l = 2 then g vanishes. In light of the fact that the exact amplification factor of a Fourier mode of Eq. (5) satisfies jgexactj = em, Eq. (27) is seen to be excessively dissipative. This difficulty can be alleviated by the use of the method of characteristics [8]. With reference to Fig. 5, dashed characteristics are traced forward from (xi1, tn) to (xi, tR) and backward from (xi+1, tn+1) to (xi, tL). Along a characteristic, Eq. (5) reduces to the ODE du c þ u ¼ 0. dx b
ð29Þ
This is numerically integrated from xi1 to xi by the trapezoidal rule to produce uR ¼
1 k=2 ui1;n . 1 þ k=2
ð30Þ
Parabolic interpolation is then employed at xi thereby yielding uL ¼ uR þ ð1 2=lÞ ðui;nþ1 ui;n Þ.
ð31Þ
Finally, Eq. (29) is numerically integrated from tn to tn+1 by the trapezoidal rule to produce uiþ1;nþ1 ¼
1 k=2 uL . 1 þ k=2
ð32Þ
Combining Eqs. (30)–(32) yields the scheme 1 k=2 1 k=2 ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ . uiþ1;nþ1 ¼ 1 þ k=2 1 þ k=2
ð33Þ
The amplification factor for this scheme is given by g¼
h2 ejh hð1 2=lÞ ; ejh hð1 2=lÞ
Fig. 5. Characteristic interpretation of RAD scheme (1 6 l 6 1).
ð34Þ
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
671
1k=2 where h ¼ 1þk=2 . Setting h = 0, observe that if l ! 1 with m fixed then g ! 1. This is in sharp contrast to the behavior of jgexactj = em and thus Eq. (33) is also of little general utility. Let us attempt to modify Eq. (33) via exponential fitting [9]. Returning to Eq. (29), we multiply by the integrating factor ecx/b thereby producing
d cx=b ðe uÞ ¼ 0. dx
ð35Þ
Thus, Eq. (30) becomes uR ¼ ek ui1;n ;
ð36Þ
while Eq. (32) becomes uiþ1;nþ1 ¼ ek uL .
ð37Þ
Combining Eqs. (36), (31) and (37) yields the exponentially fitted scheme uiþ1;nþ1 ¼ ek ek ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ .
ð38Þ
The amplification factor for this scheme is given by g¼
e2k ejh ek ð1 2=lÞ ; ejh ek ð1 2=lÞ
ð39Þ
which is identical to Eq. (34) if we replace ek by h which is its [1/1] Pade´ approximant [12]. Hence, Eq. (38) is equally as useless as Eq. (33). The dire nature of our plight is underscored by Fig. 6, which displays jgj/jgexactj for Eq. (38) with l = 1.5, m = 3.0 and l = 4.0, m = 2.0. It is now painfully obvious that exponential fitting in the case of l P 1 is a more delicate matter than it might first appear. Consequently, we take the following more circuitous route. Prior to the application of the method of characteristics, we multiply Eq. (5) by the integrating factor ect and introduce the auxiliary variable v(x, t) :¼ ect Æ u(x, t). This results in the pure advection problem ov ov þ b ¼ 0; ot ox
ð40Þ
which may then be approximated by the reflected angled derivative (RAD) scheme [4] viþ1;nþ1 ¼ vi1;n þ ð1 2=lÞ ðvi;nþ1 vi;n Þ.
ð41Þ
μ = 1.5, ν = 3.0 90
μ = 4.0, ν = 2.0 90
2
120
60
1
150
30
180
330
300
240 270
60
2
150
0
210
4
120
30
180
0
210
330
300
240 270
Fig. 6. Relative magnitude of Eq. (38) amplification factor.
672
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
In terms of the original variable u(x, t), this produces the exponential reflected angled derivative (ERAD) scheme uiþ1;nþ1 ¼ em ui1;n þ ð1 2=lÞ ðui;nþ1 em ui;n Þ;
ð42Þ
so named because its finite difference stencil is topologically equivalent to that of the EAD scheme reflected about the northeasterly diagonal. Alternatively, this may be interpreted as a leap-frog approximation centered at (xi, tn+1/2) on an appropriately angled grid. The amplification factor of this scheme is 2 þ l ðejh 1Þ gERAD ¼ em ; ð43Þ 2 þ l ðejh 1Þ thus implying that jgj = em = jgexactj. Observe that an analogous derivation of Eq. (20) could be given in the case of 0 6 l 6 1. Thus, in order to produce a nondissipative scheme, l P 1 requires multiplication by the integrating factor prior to the application of the method of characteristics whereas for 0 6 l 6 1 the order of these two operations is immaterial. For subsequent comparison to the exponential box scheme in the range l P 1, we record here its local truncation error [6,10,13]: i;nþ1=2
sERAD ¼
ðDxÞ3 i;nþ1=2 u lð2l 1Þð1 lÞ. 12Dt xxx
ð44Þ
As is clear on the grounds of symmetry, a prudent treatment of the grid line x = Dx (as provided by the EMAD scheme below) is necessary when employing the ERAD scheme. 2.3. Exponential box (EBOX) scheme With reference to Fig. 7, the box scheme may be derived as follows: 1 huiþ1;nþ1 uiþ1;n ui;nþ1 ui;n i b huiþ1;nþ1 ui;nþ1 uiþ1;n ui;n i þ þ þ þ c ^uiþ1=2;nþ1=2 ¼ 0. 2 2 Dt Dt Dx Dx
ð45Þ
If we use the average value ^ uiþ1=2;nþ1=2 ðuiþ1;nþ1 þ ui;nþ1 þ uiþ1;n þ ui;n Þ=4 then the following scheme results: uiþ1;nþ1 ¼
1 ½ð1 þ l m=2Þ ui;n þ ð1 l m=2Þ uiþ1;n ð1 l þ m=2Þ ui;nþ1 . 1 þ l þ m=2
ð46Þ
Upon insertion of the Fourier mode [10] ui,n = gnejk(iDx) into Eq. (46), we find that the amplification factor, g, is given by g¼
ð1 þ l m=2Þ þ ð1 l m=2Þejh ; ð1 l þ m=2Þ þ ð1 þ l þ m=2Þejh
ð47Þ
thereby implying that jgj 6 1 and thus that the box scheme, Eq. (46), is stable for any 0 6 l 6 1 provided that m P 0. However, its unsuitability is evident from Fig. 8, which displays jgj/jgexactj for l = 0.5, m = 2.0 and l = 2.0, m = 2.0.
Δ
Δ
Fig. 7. Box scheme (0 6 l 6 1).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 μ = 0.5, ν = 2.0 90
μ = 2.0, ν = 2.0 90
10
120
60
5
150
30
60
5
150
0
210
10
120
180
330
30
180
0
210
300
240
673
330
300
240
270
270
Fig. 8. Relative magnitude of box scheme amplification factor.
μ <−
μ −>
Fig. 9. Characteristic interpretation of box scheme (0 6 l 6 1).
With reference to Fig. 9, we may appeal to exponential fitting in order to resolve this dilemma. We first multiply Eq. (5) by the integrating factor ect and then introduce the auxiliary variable v(x, t) :¼ ect Æ u(x, t), resulting in the pure advection problem ov ov þ b ¼ 0; ot ox
ð48Þ
which may then be approximated by the method of characteristics (i.e., the box scheme) viþ1;nþ1 ¼ vi;n þ
1l ðviþ1;n vi;nþ1 Þ. 1þl
ð49Þ
In terms of the original variable u(x, t), this produces the exponential box (EBOX) scheme uiþ1;nþ1 ¼ em ui;n þ
1l ðem uiþ1;n ui;nþ1 Þ. 1þl
The amplification factor of this scheme is ð1 þ lÞ þ ð1 lÞejh gEBOX ¼ em ; ð1 lÞ þ ð1 þ lÞejh
ð50Þ
ð51Þ
thus implying that jgj = em = jgexactj. In order to compare the dispersive properties of the EAD, EBOX, and ERAD schemes, we consider their relative phase errors [10]. The relative phase errors are plotted in Fig. 10 for various values of 0 6 l 6 1. In this figure, the top curve corresponds to the EBOX scheme and the bottom curve to the EAD scheme. The adequacy of the EBOX scheme in the neighborhood of l = 1 is readily apparent as is the overall superiority of the EAD scheme. Observe that the EAD scheme is exact for l = 1/2 while both schemes are exact for l = 1. The relative phase errors are plotted in Fig. 11 for various values of
674
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 μ=.25
μ=.50 90 2
90 4 120
60 2
150
120 30
180
330 240
0
210
330 240
μ=1.0
90 2
90 1 60
1
150
0
330
30
180
0
210
330 240
300
270
0.5
150
30
210
60
120
180
240
300
270
μ=.75 120
30
180
300
270
1
150
0
210
60
300
270
Fig. 10. Comparison of relative phase errors (l 6 1) (top curve = EBOX scheme, bottom curve = EAD scheme).
μ=1.0 120
μ=1.5
90 1
120
60 0.5
150
30
180
0
330
210 240
1
0
330
210 240
μ=4.0 120
60 0.5
0
210
330 300 270
90 1
60
0.5
150
30
180
240
300
270
90 1
150
30
180
μ=2.0 120
60
150
300
270
90 2
30
0
180
210
330 240
270
300
Fig. 11. Comparison of relative phase errors (l P 1) (top curve = ERAD scheme, bottom curve = EBOX scheme).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
675
1 6 l 6 1. In this figure, the top curve corresponds to the ERAD scheme and the bottom curve to the EBOX scheme. The adequacy of the EBOX scheme in the neighborhood of l = 1 is once again readily apparent as is the overall superiority of the ERAD scheme. Observe that the ERAD scheme is exact for l = 2 while both schemes are exact for l = 1. As a vehicle for further assessing the relative merits of the EAD, EBOX, and ERAD schemes, consider their respective local truncation errors [6,10,13]: 3
iþ1=2;nþ1=2
sEBOX
¼
ðDxÞ uiþ1=2;nþ1=2 lðl 1Þðl þ 1Þ; 12Dt xxx
ð52Þ
together with Eqs. (23) and (44). As they all possess a common factor, we introduce the scaled local truncation errors t :¼
12Dt ðDxÞ3 uxxx
s.
ð53Þ
The scaled local truncation errors are plotted in Fig. 12 for the range 0 6 l 6 1. Observe the suitability of the EBOX scheme in the vicinity of l = 1 as well as the overall superiority of the EAD scheme in this regime. The scaled local truncation errors are plotted in Fig. 13 for the range 1 6 l 6 1. Observe once again the suitability of the EBOX scheme in the vicinity of l = 1 as well as the overall superiority of the ERAD scheme in this regime. Based upon these preliminary observations, we now turn to the construction of a new scheme for the nonhomogeneous, variable coefficient Eq. (3), which is a hybrid of the EAD, EBOX, and ERAD schemes that is nondissipative with good dispersive properties, second-order accurate, explicit, and unconditionally stable. 3. Exponential method of angled derivatives (EMAD) In light of the foregoing introductory observations, it would seem that the best course of action would be to use the EAD scheme, Eq. (20), if l 6 1 and the ERAD scheme, Eq. (42), if l P 1. This is indeed the case for
SCALED LOCAL TRUNCATION ERROR
0.05 EAD Scheme 0
−0.05
t
−0.1 −0.15 −0.2 EBOX Scheme −0.25 −0.3 −0.35 0
0.1
0.2
0.3
0.4
0.5
μ
0.6
0.7
0.8
Fig. 12. Comparison of scaled local truncation errors (l 6 1).
0.9
1
676
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 SCALED LOCAL TRUNCATION ERROR 25
20 EBOX Scheme
15
t
ERAD Scheme 10
5
0
1
1.5
2
2.5
μ
3
3.5
4
Fig. 13. Comparison of scaled local truncation errors (l P 1).
constant b in Eq. (3). However, if b is not constant then difficulties arise. With reference to Fig. 2, if neither li+1/2,n nor li,n+1/2 exceed one then the EAD scheme is stable and should be used. Likewise, with reference to Fig. 5, if both li+1/2,n and li,n+1/2 exceed one then the ERAD scheme is stable and should be used. Yet, consider the scenario pictured in Fig. 14. There, both the EAD and the ERAD schemes are stable and a seemingly arbitrary choice must be made between them. While this perplexity is venial at worst, that displayed in Fig. 15 is decidedly mortal. In this case, no matter which of the two schemes is utilized, the choice is condemned to instability. The resolution of these conundrums is provided by the employment of the EBOX scheme in these two instances. Since in both cases l passes through unity somewhere in the shaded cell, the EBOX scheme is an excellent choice because of its exemplary behaviour in the neighborhood of l = 1. We now develop these ideas into a complete algorithm for the initial-boundary value problem, Eqs. (2) and (3). In order to distinguish it from the EAD and ERAD schemes, which are merely components of this new scheme, it is dubbed the exponential method of angled derivatives (EMAD). An important feature of this development will be the preservation of the desirable properties described in the previous section, which are strictly applicable only to the homogeneous, constant coefficient problem, as we extend the component EAD, EBOX, and ERAD schemes to the nonhomogeneous, variable coefficient problem. Particular attention will be paid to approximation at the grid points adjacent to the spatiotemporal boundary.
μ
μ< Fig. 14. Indecisive scenario.
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
677
μ<
μ Fig. 15. Checkmate!
In the ensuing subsections, we extensively utilize the Bernoulli function z ; BðzÞ ¼ z e 1
ð54Þ
and the associated function 1 BðzÞ . ð55Þ z Their essential role in exponential fitting as well their stable numerical evaluation for small z are discussed in ^ and ^c, which will be specified [9]. Furthermore, we will freeze the coefficients of Eq. (3) at constant values, b ^ct below. In terms of the auxiliary variable vðx; tÞ :¼ e uðx; tÞ, Eq. (3) then becomes MðzÞ ¼
ov ^ ov þ b ¼ e^ct f ðx; tÞ. ot ox
ð56Þ
Subsequent application of the MAD scheme [4] to Eq. (56) with special attention to the required exponential quadrature on its right hand side will yield the EMAD scheme. 3.1. Nonhomogeneous, variable coefficient EAD scheme Let us begin with the case where neither li+1/2,n nor li,n+1/2 exceed one so that the EAD scheme is the ^ :¼ b ^ appropriate choice (see Fig. 2). Defining b c :¼ ciþ1=2;n and l :¼ bDt=Dx, m :¼ ^cDt, the dashed lines iþ1=2;n , ^ ^ In actuality, the characteristics are curved but, due to the symare approximate characteristics with slope 1=b. metry of the finite difference stencil about the point (xi+1/2, tn), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv ¼ e^ct f . dt
ð57Þ
Integrating this along the characteristics produces Z ðxT ;tn Þ m uT ¼ e ui;n1 þ e^cðtn tÞ f dt;
ð58Þ
ðxi ;tn1 Þ
uiþ1;nþ1 ¼ em uB þ
Z
ðxiþ1 ;tnþ1 Þ
e^cðtnþ1 tÞ f dt;
ð59Þ
ðxB ;tn Þ
where xT = xi + lDx and xB = xi + (1 l)Dx. Linear interpolation of the forcing term produces fT ¼ ð1 lÞfi;n þ lfiþ1;n ;
f B ¼ lfi;n þ ð1 lÞfiþ1;n .
ð60Þ
Employing the exponential trapezoidal rule [9] to perform the quadratures in Eqs. (58) and (59) now produces
678
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
Dt ½MðmÞ fi;n1 þ MðmÞ fT ; BðmÞ Dt ½MðmÞ fB þ MðmÞ fiþ1;nþ1 . ¼ em uB þ BðmÞ
uT ¼ em ui;n1 þ
ð61Þ
uiþ1;nþ1
ð62Þ
Parabolic interpolation to the values of v at (xi, tn), (xT, tn), (xi+1, tn) and subsequent evaluation at (xB, tn) produces uB ¼ uT þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð63Þ
Combination of Eqs. (60)–(63) now yields the EAD scheme: uiþ1;nþ1 ¼ em ½em ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ Dt fem MðmÞ fi;n1 þ ½em MðmÞ ð1 lÞ þ MðmÞ l fi;n þ BðmÞ þ ½em MðmÞ l þ MðmÞ ð1 lÞ fiþ1;n þ MðmÞ fiþ1;nþ1 g.
ð64Þ
3.2. Nonhomogeneous, variable coefficient ERAD scheme Next, we turn to the case where both li+1/2,n and li,n+1/2 exceed one so that the ERAD scheme is the appro^ :¼ b ^ ^ the dashed priate choice (see Fig. 5). Defining b c :¼ ci;nþ1=2 and l :¼ bDt=Dx, m :¼ ^cDt, k :¼ ^cDx=b, i;nþ1=2 , ^ ^ In actuality, the characteristics are curved but, due to the lines are approximate characteristics with slope 1=b. symmetry of the finite difference stencil about the point (xi, tn+1/2), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv 1 ^ct ¼ e f. ð65Þ ^ dx b Integrating this along the characteristics produces Z ðxi ;tR Þ uR ¼ ek ui1;n þ e^cðtR tÞ f dt;
ð66Þ
ðxi1 ;tn Þ
uiþ1;nþ1 ¼ ek uL þ
Z
ðxiþ1 ;tnþ1 Þ
e^cðtnþ1 tÞ f dt;
ð67Þ
ðxi ;tL Þ
where tR ¼ tn þ l1 Dt and tL ¼ tn þ ð1 l1ÞDt. Linear interpolation of the forcing term produces 1 1 1 1 fR ¼ 1 f L ¼ fi;n þ 1 fi;n þ fi;nþ1 ; fi;nþ1 . l l l l
ð68Þ
Employing the exponential trapezoidal rule [9] to perform the quadratures in Eqs. (66) and (67) now produces Dt=l ½MðkÞ fi1;n þ MðkÞ fR ; BðkÞ Dt=l ½MðkÞ fL þ MðkÞ fiþ1;nþ1 . ¼ ek uL þ BðkÞ
uR ¼ ek ui1;n þ
ð69Þ
uiþ1;nþ1
ð70Þ
Parabolic interpolation to the values of v at (xi, tn), (xi, tR), (xi, tn+1) and subsequent evaluation at (xi, tL) produces 2 mþ2k uL ¼ e uR þ 1 ð71Þ ðek uiþ1;n emþk ui;n Þ. l Combination of Eqs. (68)–(71) now yields the scheme:
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
uiþ1;nþ1
2 ðui;nþ1 em ui;n Þ ¼ e ui1;n þ 1 l
Dt=l 1 1 km km e MðkÞ fi1;n þ e MðkÞ 1 þ þ MðkÞ fi;n BðkÞ l l 1 1 þ ekm MðkÞ þ MðkÞ 1 fi;nþ1 þ MðkÞ fiþ1;nþ1 . l l
679
m
ð72Þ
Observe that Eq. (72) does not reproduce the constant solution u fc for constant c and f. Instead, it produces BðmÞ ekm þ1 f u BðkÞ 2 c. Thus, we modify it as follows to yield the ERAD scheme: 2 m uiþ1;nþ1 ¼ e ui1;n þ 1 ðui;nþ1 em ui;n Þ l
Dt=l 2 1 1 km þ ekm MðkÞ fi1;n þ ekm MðkÞ 1 þ MðkÞ fi;n BðmÞ e þ 1 l l 1 1 þ ekm MðkÞ þ MðkÞ 1 fi;nþ1 þ MðkÞ fiþ1;nþ1 . ð73Þ l l Eq. (73) has the same formal order of accuracy and dissipative properties as Eq. (72) yet possesses the correct constant solution. 3.3. Nonhomogeneous, variable coefficient EBOX scheme Finally, suppose that one of li+1/2,n and li,n+1/2 exceeds one while the other does not. In this case, the ^ :¼ b EBOX scheme is the appropriate choice (see Fig. 9). Defining b c :¼ ciþ1=2;nþ1=2 and iþ1=2;nþ1=2 , ^ ^ ^ ^ In actuall :¼ bDt=Dx, m :¼ ^cDt, k :¼ ^cDx=b, the dashed lines are approximate characteristics with slope 1=b. ity, the characteristics are curved but, due to the symmetry of the finite difference stencil about the point (xi+1/2, tn+1/2), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. 3.3.1. EBOX scheme (0 6 l 6 1) First consider the case l 6 1. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv ¼ e^ct f . dt Integrating this along the characteristic produces Z ðx2 ;tnþ1 Þ e^ct f dt; v2 ¼ v1 þ
ð74Þ
ð75Þ
ðx1 ;tn Þ
where x1 = xi + (1 l)Dx/2 and x2 = xi + (1 + l)Dx/2. Linear interpolation of the forcing term produces 1þl 1l 1l 1þl fi;n þ fiþ1;n ; fi;nþ1 þ fiþ1;nþ1 ; f2 ¼ ð76Þ f1 ¼ 2 2 2 2 while linear interpolation of the auxiliary variable produces 1þl 1l 1l 1þl vi;n þ viþ1;n ; vi;nþ1 þ viþ1;nþ1 . v2 ¼ ð77Þ v1 ¼ 2 2 2 2 Employing the exponential trapezoidal rule [9] to perform the quadrature in Eq. (75) now yields the EBOX scheme: 1l Dt ðem uiþ1;n ui;nþ1 Þ þ uiþ1;nþ1 ¼ em ui;n þ 1þl BðmÞ 1l fMðmÞ fiþ1;n þ MðmÞ fi;nþ1 g . MðmÞ fi;n þ MðmÞ fiþ1;nþ1 þ 1þl
ð78Þ
680
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
3.3.2. EBOX scheme (1 6 l 6 1) Next consider the case l P 1. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv 1 ^ct ¼ e f. ^ dx b
ð79Þ
Integrating this along the characteristic produces Z ðxiþ1 ;t2 Þ v2 ¼ v1 þ e^ct f dt;
ð80Þ
ðxi ;t1 Þ
where t1 = tn + (1 1/l)Dt/2 and t2 = tn + (1 + 1/l)D t/2. Linear interpolation of the forcing term produces 1 þ 1=l 1 1=l 1 1=l 1 þ 1=l fi;n þ fi;nþ1 ; fiþ1;n þ fiþ1;nþ1 ; f2 ¼ 2 2 2 2 while linear interpolation of the auxiliary variable produces f1 ¼
v1 ¼
1 þ 1=l 1 1=l vi;n þ vi;nþ1 ; 2 2
v2 ¼
1 1=l 1 þ 1=l viþ1;n þ viþ1;nþ1 . 2 2
ð81Þ
ð82Þ
Employing the exponential trapezoidal rule [9] to perform the quadrature in Eq. (80) now yields the scheme: 1l Dt eðkmÞ=2 ðem uiþ1;n ui;nþ1 Þ þ uiþ1;nþ1 ¼ em ui;n þ 1þl l BðkÞ 1l fMðkÞ fi;nþ1 þ MðkÞ fiþ1;n g . MðkÞ fi;n þ MðkÞ fiþ1;nþ1 1þl
ð83Þ
Observe that Eq. (83) does not reproduce the constant solution u fc for constant c and f. Instead, it produces BðmÞ u BðkÞ eðkmÞ=2 fc . Thus, we modify it as follows to yield the EBOX scheme: 1l Dt=l uiþ1;nþ1 ¼ em ui;n þ ðem uiþ1;n ui;nþ1 Þ þ 1þl BðmÞ 1l fMðkÞ fi;nþ1 þ MðkÞ fiþ1;n g . MðkÞ fi;n þ MðkÞ fiþ1;nþ1 1þl
ð84Þ
Eq. (84) has the same formal order of accuracy and dissipative properties as Eq. (83) yet possesses the correct constant solution. 3.4. Boundary treatment Depending on the numerical value of l, the EMAD approximation at grid points adjacent to the spatiotemporal boundary may require modification (see Fig. 16). Along t = 0 the EAD scheme is inapplicable while the same is true for the ERAD scheme along x = 0. Even though we know the exact solution there (see Eq.
μ
μ μ<
Fig. 16. Boundary treatment.
μ<
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
681
(2)), tracing a linear approximation to a characteristic from either t = Dt or x = Dx backwards to the spatiotemporal boundary would yield only first-order accuracy unless the coefficient b is constant in which case the characteristics are in fact straight. However, a parabolic approximation to the characteristics yields secondorder accuracy. 3.4.1. Southwestern corner (0 6 l 6 1) Let us begin with the southwestern-most grid cell (i = 1, n = 1). If l < 1 throughout this cell then, starting at (xi+1, tn+1), we trace backwards to t = 0 along a parabolic trajectory x(t) whose first two derivatives coincide with those of the actual characteristic through this point: 2
xB ¼ xiþ1 biþ1;nþ1 Dt þ ðbt þ bbx Þiþ1;nþ1
ðDtÞ . 2
ð85Þ
Then, as above, we integrate Eq. (74) along this trajectory with ^c :¼ cððxB þ xiþ1 Þ=2; Dt=2Þ and m :¼ ^cDt thereby producing the approximation uiþ1;nþ1 ¼ em /ðxB Þ þ
Dt ½MðmÞ fB þ MðmÞ fiþ1;nþ1 . BðmÞ
ð86Þ
In place of fB, we may instead employ the linear interpolant xiþ1 xB xB xi fi;n þ fiþ1;n ; f^ B ¼ Dx Dx
ð87Þ
without sacrifice of order of accuracy. In order to insure that the approximate characteristic slope, s, is correspondingly bounded (0 6 1/s 6 Dx/Dt), we must either impose the restriction 8 b db > > if > 0; > < db=dt dt iþ1;nþ1 iþ1;nþ1 Dt 6 > b db > > ð1=liþ1;nþ1 1Þ if < 0; : db=dt iþ1;nþ1 dt iþ1;nþ1 or utilize additional terms in the truncated Taylor series, Eq. (85). 3.4.2. Southwestern corner (1 6 l 6 1) On the other hand, if l > 1 throughout this cell then, starting at (xi+1, tn+1), we trace backwards to x = 0 along a parabolic trajectory t(x) whose first two derivatives coincide with those of the actual characteristic through this point: tL ¼ tnþ1
1 biþ1;nþ1
Dx
ðbt þ bbx Þiþ1;nþ1 ðDxÞ2 . 2 b3iþ1;nþ1
ð88Þ
Then, as above, we integrate Eq. (79) along this trajectory with ^c :¼ cðDx=2; ðtL þ tnþ1 Þ=2Þ, D :¼ tn+1 tL and ^m :¼ ^cD thereby producing the approximation uiþ1;nþ1 ¼ e^m wðtL Þ þ
D ½Mð^mÞ fB þ Mð^mÞ fiþ1;nþ1 . Bð^mÞ
In place of fL, we may instead employ the linear interpolant tnþ1 tL tL tn fi;n þ fi;nþ1 ; f^L ¼ Dt Dt
ð89Þ
ð90Þ
without sacrifice of order of accuracy. In order to insure that the approximate characteristic slope, s, is correspondingly bounded (0 6 s 6 Dt/Dx), we must either impose the restriction
682
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
8 b > > > < db=dx iþ1;nþ1 Dx 6 > b > > ðliþ1;nþ1 1Þ : db=dx iþ1;nþ1
db if < 0; dx iþ1;nþ1 db if > 0; dx iþ1;nþ1
or utilize additional terms in the truncated Taylor series, Eq. (88). 3.4.3. Remaining portion of the boundary In grid cells elsewhere along the southern boundary (i 5 1, n = 1) of the mesh, we use Eq. (86) if l < 1 throughout the cell, the ERAD scheme if l > 1 throughout the cell, and the EBOX scheme otherwise. Likewise, in grid cells along the western boundary (i = 1, n 5 1) of the mesh, we use Eq. (89) if l > 1 throughout the cell, the EAD scheme if l < 1 throughout the cell, and the EBOX scheme otherwise. This completes the specification of the EMAD algorithm. It is an explicit, three-level, unconditionally stable scheme that is nondissipative with good dispersive properties while maintaining second-order accuracy even at the boundary and in spite of the presence of variable coefficients and/or forcing function. We now present numerical evidence of the versatility of this new scheme. 4. Numerical examples 4.1. Spatiotemporally varying velocity with source Our first example is nonhomogeneous with forcing term f ðx; tÞ ¼ 2p ½2tð1 þ xÞ cosð2pxÞ cosð2ptÞ sinð2pxÞ sinð2ptÞ þ sinð2pxÞ cosð2ptÞ.
t=0.5 1
u
0.5
0
−0.5 −1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.6
0.7
0.8
0.9
1
x t=1.0 1
u
0.5
0
−0.5 −1
0
0.1
0.2
0.3
0.4
0.5
x Fig. 17. Spatiotemporally varying velocity with source (Dt = .05, Dx = .05).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
683
Also, there is a coefficient b(x, t) = 2t(1 + x) depending upon both space and time with c = 2p. The initial condition is /(x) = sin(2px) with boundary condition w(t) = 0. The exact solution is then u(x, t) = sin(2px) Æ cos(2pt). For this example, we have chosen Dt = .05 and Dx = .05. The numerical results for t = 0.5, halfway through a cycle, and t = 1.0, at the completion of a cycle, are displayed as dots in Fig. 17 where the solid line represents the exact solution. 4.2. Countercurrent wave propagation Our next example is of importance in the mathematical modelling of the warmup process for automotive catalytic converters [2,3]. In spite of the fact that the flow of air is down the length of the converter, chemical reactions along its surface give rise to a backward propagating wave known as the light-off point. Herewith, we model such countercurrent wave propagation through the forcing term 2 12 f ðx; tÞ ¼ 50 erfc½6ðx þ t 2Þ pffiffiffi e½6ðxþt2Þ ; p with b = 1 and c = 100. The initial condition is /ðxÞ ¼ 12 erfc½6ðx 2Þ with boundary condition wðtÞ ¼ 12 erfc½6ðt 2Þ. The exact solution is now uðx; tÞ ¼ 12 erfc½6ðx þ t 2Þ. For this example, we have chosen Dt = .025 and Dx = .025. The numerical results for t = 0.0, 1.0, 1.5, 2.0 are displayed as dots in Fig. 18 where the solid line represents the exact solution. The upstream propagation of the complementary error function profile is notable. 4.3. Oscillatory flow Our final example also arises in the numerical simulation of automotive catalytic converters [3]. If the converter is mounted in proximity to the exhaust manifold then the driving fluid flow decidedly pulsates. We
t = 0.0
t = 1.0
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
0.6
0.8
0
1
0
0.2
0.4
x
0.6
0.8
1
0.8
1
x
t = 1.5
t = 2.0
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
0.6
x
0.8
1
0
0
0.2
0.4
0.6
x
Fig. 18. Countercurrent wave propagation (Dt = .025, Dx = .025).
684
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 t = .06
t = .12
t = .18
1
1
1
0.5
0.5
0.5
0
0
0.5
1
0
0
t = .24
0.5
1
0
1
1
0.5
0.5
0.5
0
0.5
1
0
0
t = .42
0.5
1
0
1
0.5
0.5
0.5
0.5
1
0
0
0.5
1
0.5
1
t = .54
1
0
0
t = .48
1
0
0.5 t = .36
1
0
0
t = .30
1
0
0
0.5
1
Fig. 19. Oscillatory flow.
model this through the oscillating inlet velocity bðtÞ ¼ 1 þ 12 cos 5p t with c = 10 and f = 10. The initial con3 dition is /(x) = 0 with inlet boundary condition w(t) = 0. The exact solution is then 8 3 5p > 10ðtnÞ > sin t ; 1 e ; x 6 t þ < 10p 3 uðx; tÞ ¼ > 3 5p > : 1 e10t ; sin t ; xPtþ 10p 3 where n satisfies the transcendental equation 3 5p 3 5p nþ sin n þx tþ sin t ¼ 0. 10p 3 10p 3 For this example, we have chosen Dt = .01 and Dx = .02. Thus, every 12 time steps is equivalent to a complete pulsation cycle. The numerical results for t = .06, .12, .18, .24, .30, .36, .42, .48, .54 are displayed as dots in Fig. 19 where the solid line represents the exact solution, so that each successive frame is equivalent to half a cycle. Observe the periodic response, which develops once the transients have decayed. 5. Conclusion In this paper, we have presented a hybrid scheme for the numerical approximation of the nonhomogeneous one-way advection-reaction equation. This exponential method of angled derivatives (EMAD) scheme provides an explicit, three-level, nondissipative, second-order accurate, unconditionally stable numerical approximation to this important initial-boundary value problem. This new scheme provides the natural generalization of the method of angled derivatives (MAD) scheme [4] for pure advection problems.
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
685
We began by developing the EAD scheme for 0 6 l 6 1, where its interpretation in terms of characteristics, the explicit computation of its local truncation error, and the analysis of its dispersive properties took center stage. In addition, we have introduced the ERAD scheme for 1 6 l 6 1 and provided a detailed analysis of its numerical properties. The net result of this preliminary analysis was that, for the homogeneous, constant coefficient problem, the EAD and ERAD schemes provide superlative approximations within their respective parameter regimes. An ancillary benefit of the preceding analysis is a unified approach via exponential fitting and characteristics to the seemingly disparate box, angled derivative, and reflected angled derivative schemes. Likewise, the Lax–Wendroff, second-order upwind difference, Roberts–Weiss, and leap-frog schemes [5,6,8,10,11,13], as well as any other consistent scheme for Eq. (1), may also be so treated. For the nonhomogeneous, variable coefficient problem, the suitably extended EAD and ERAD schemes were then melded together using the EBOX scheme (whose characteristic features were also herein developed) as adhesive in order to provide a complete numerical algorithm unrestricted in its value of l. Particular attention has been paid to maintaining the order of accuracy and stability of the approximation in the vicinity of the spatiotemporal boundary. This new EMAD scheme was then exercised on a pair of important test problems in order to demonstrate its versatility. Acknowledgement The author would like to express his sincere gratitude to Mrs. Barbara A. McCartin for her dedicated assistance in the production of this paper. References [1] H. Rhee, R. Aris, N.R. Amundson, First Order Partial Differential Equations, Theory and Application of Single Equations, vol. I, Dover, New York, 2001. [2] B.A. Finlayson, Numerical Methods for Problems with Moving Fronts, Ravenna Park, Seattle, WA, 1992. [3] B.J. McCartin, P.D. Young, Numerical investigation of the validity of the quasi-static approximation in the modelling of catalytic converters, in: D. Ferguson, T. Peters (Eds.), Mathematics and Industry, SIAM, Philadelphia, PA, 2005. [4] B.J. McCartin, The method of angled derivatives, Applied Mathematics and Computation 170 (2005) 440–461. [5] A. Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Cambridge, 1996. [6] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, 1994. [7] B.J. McCartin, Seven deadly sins of numerical computation, American Mathematical Monthly 105 (10) (1998) 929–941. [8] J.D. Hoffman, Numerical Methods for Engineers and Scientists, McGraw-Hill, New York, 1992. [9] B.J. McCartin, Exponential fitting of the delayed recruitment/renewal equation, Journal of Computational and Applied Mathematics 136 (2001) 343–356. [10] J.C. Tannehill, D.A. Anderson, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, second ed., Taylor & Francis, New York, 1997. [11] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth, Belmont, CA, 1989. [12] G.A. Baker, P. Graves-Morris, Pade´ Approximants, second ed., Cambridge University Press, Cambridge, 1996. [13] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, Wiley-Interscience, New York, 1995.
The exponential method of angled derivatives Brian J. McCartin Applied Mathematics, Kettering University, 1700 West Third Avenue, Flint, MI 48504-4898, United States
Abstract A hybrid numerical method is presented for a linear, first order, hyperbolic partial differential equation (PDE): the nonhomogeneous one-way advection-reaction equation. This PDE is first reduced to an ordinary differential equation (ODE) along a characteristic emanating backward in time from each mesh point. These ODEs are then integrated numerically using exponential fitting to advance in time thus permitting wide disparities in size among the coefficients. Such temporal integration requires a spatial interpolation procedure dependent upon the slope of the local characteristic. If the local Courant number is less than one then the interpolation is afforded by the exponential angled derivative method while if greater than one then the exponential reflected angled derivative approximation is employed. Otherwise, if the Courant number is transitioning through one then the exponential box scheme is employed. At the boundary, parabolic tracing of the characteristics provides a comparable level of accuracy. The net result is a three-level explicit method which is second-order accurate and unconditionally stable. The efficacy of this exponential method of angled derivatives (EMAD) scheme for singular perturbation problems is demonstrated numerically. Ó 2005 Elsevier Inc. All rights reserved.
1. Introduction In what follows, attention is focused on the numerical solution of the nonhomogeneous one-way advectionreaction equation aðx; tÞ
ou ou þ bðx; tÞ þ cðx; tÞ u ¼ dðx; tÞ; ot ox
0 < x 6 1; t > 0; aðx; tÞ > 0; bðx; tÞ P 0; cðx; tÞ > 0; ð1Þ
subject to the initial-boundary conditions uðx; 0Þ ¼ /ðxÞ;
uð0; tÞ ¼ wðtÞ.
ð2Þ
For example, this problem arises in chemical engineering in the form of the plug flow reactor and heat exchanger equations [1]. Also, it plays a prominent role in the mathematical modelling of automotive catalytic converters [2]. In the latter application, u(x, t) is either a concentration or a temperature, a(x, t) is a measure of
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.09.097
666
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
porosity, b(x, t) is the convective velocity, c(x, t) is either a mass or heat transfer coefficient, and d(x, t) accounts for sources and sinks [3]. In the ensuing pages, we develop an approximate method for this problem which is nondissipative, secondorder accurate, explicit, and unconditionally stable that allows for wide disparities in the magnitude of the coefficients. We begin by rewriting Eq. (1) in the normalized form: ou ou þ bðx; tÞ þ cðx; tÞ u ¼ f ðx; tÞ; ot ox
ð3Þ
0 < x 6 1; t > 0;
where bðx; tÞ :¼ bðx; tÞ=aðx; tÞ P 0;
cðx; tÞ :¼ cðx; tÞ=aðx; tÞ > 0;
f ðx; tÞ :¼ dðx; tÞ=aðx; tÞ.
ð4Þ
The special case of pure advection, c(x, t) 0, has previously been treated by the method of angled derivatives (MAD) [4]. This MAD scheme is an amalgam of the angled derivative [5], box [6], and reflected angled derivative [4] schemes combined in such a way as to accentuate the best features of each of these component schemes. It is the express intent of the present paper to extend the MAD scheme to permit reaction processes. We begin by showing that straightforward methods of approximating the reaction term lead to undesirable numerical phenomena even in the homogeneous, constant coefficient case. A solution to this dilemma is provided by the application of exponential fitting [7]. It is important to realize that the numerical issues herein raised and resolved are not confined to the three schemes cited above but are in fact germane to any consistent approximation to Eq. (1) such as the Lax–Wendroff, upwind, and leap-frog schemes [8]. Once the need for exponential fitting has been firmly established, we develop the complete exponential method of angled derivatives (EMAD) for nonhomogeneous, variable coefficient problems. This new EMAD scheme is then exercised on a pair of interesting and challenging test problems which arise in the numerical simulation of automotive catalytic converters [3]. 2. The need for exponential fitting In order to firmly establish the need for exponential fitting for advection-reaction problems, we confine our attention, in this section only, to the homogeneous, constant coefficient problem, ou ou þ b þ c u ¼ 0; ot ox
0 < x 6 1; t > 0; b P 0; c > 0;
ð5Þ
subject to the initial-boundary conditions, Eq. (2). For pure advection, c = 0, a comprehensive analysis of the angled derivative (AD), reflected angled derivative (RAD), and box schemes was provided in [4]. The first obstacle which we must herein surmount is the existence of several seemingly natural generalizations of these schemes for the case of advection-reaction, c 5 0. We next consider these alternatives and argue that those based upon exponential fitting [9] possess superior numerical properties. 2.1. Exponential angled derivative (EAD) scheme With reference to Fig. 1, the angled derivative scheme may be derived as follows: 1 huiþ1;nþ1 þ ui;n uiþ1;n þ ui;n1 i uiþ1;n ui;n þ c ^uiþ1=2;n ¼ 0. þb Dt 2 2 Dx
ð6Þ
If we use the average value ^ uiþ1=2;n ðuiþ1;n þ ui;n Þ=2 then the following scheme results: uiþ1;nþ1 ¼ ui;n1 þ ð1 2l mÞ uiþ1;n ð1 2l þ mÞ ui;n ;
ð7Þ
where l :¼ bDt/Dx is the Courant number and m :¼ cDt. For future reference, we also introduce k :¼ cDx/b with the consequent identity l = m/k. We make the further assumption, pffiffiffiffiffiffiffi in this subsection only, that 0 6 l 6 1. Upon insertion of the Fourier mode [10] ui,n = gnejk(iDx) (j ¼ 1) into Eq. (7), we find that the amplification factor, g, satisfies g2 þ ð1 2l þ mÞ ejh ð1 2l mÞ g ejh ¼ 0; ð8Þ
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
667
Δ Δ
Fig. 1. Angled derivative scheme (0 6 l 6 1).
pffiffiffiffiffiffiffiffiffiffiffiffiffi where h :¼ kDx. Setting h = 0 in Eq. (8) and solving, we arrive at g ¼ m 1 þ m2 , so that one of its roots satisfies jgj > 1 for any m 5 0. Thus, Eq. (7) is unconditionally unstable for any m 5 0! A stable approximation to Eq. (5) may be derived by instead using the average value ^ uiþ1=2;n ðuiþ1;nþ1 þ ui;n1 Þ=2 resulting in the scheme: uiþ1;nþ1 ¼
1 ½ð1 mÞ ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ. 1þm
The amplification factor for this scheme satisfies 1 2l 1 m jh g2 þ ðejh 1Þ g e ¼ 0. 1þm 1þm
ð9Þ
ð10Þ
Application of the Schur criterion [11] establishes stability of this scheme for m > 0 and 0 6 l 6 1. However, observe that if m = 1 then one of the roots of Eq. (10) is zero. Furthermore, if, in addition, l = 1/2 then both roots vanish. In light of the fact that the exact amplification factor of a Fourier mode of Eq. (5) satisfies jgexactj = em, Eq. (9) is seen to be excessively dissipative. This difficulty can be alleviated by the use of the method of characteristics [8]. With reference to Fig. 2, dashed characteristics are traced forward from (xi, tn1) to (xT, tn) and backward from (xi+1, tn+1) to (xB, tn). Along a characteristic, Eq. (5) reduces to the ODE
Fig. 2. Characteristic interpretation of AD scheme (0 6 l 6 1).
668
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
du þ c u ¼ 0. dt
ð11Þ
This is numerically integrated from tn1 to tn by the trapezoidal rule to produce uT ¼
1 m=2 ui;n1 . 1 þ m=2
ð12Þ
Parabolic interpolation is then employed at time tn thereby yielding uB ¼ uT þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð13Þ
Finally, Eq. (11) is numerically integrated from tn to tn+1 by the trapezoidal rule to produce uiþ1;nþ1 ¼
1 m=2 uB . 1 þ m=2
Combining Eqs. (12)–(14) yields the angled derivative (AD) scheme 1 m=2 1 m=2 ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ . uiþ1;nþ1 ¼ 1 þ m=2 1 þ m=2
ð14Þ
ð15Þ
The amplification factor for this scheme satisfies 2 1 m=2 1 m=2 ð1 2lÞ ðejh 1Þ g g2 þ ejh ¼ 0; ð16Þ 1 þ m=2 1 þ m=2 thus implying that jgj ¼ 1m=2 Observe that this is the [1/1] Pade´ approximant [12] to jgexactj = em. Conse1þm=2. quently, the AD scheme is essentially nondissipative for small m but is highly dissipative for even moderate values of m. For example, m = 2 implies that the amplification factor vanishes. Yet, we may modify the AD scheme via exponential fitting [9] in such a way that the modulus of the amplification factor is exact. Returning to Eq. (11), we multiply by the integrating factor ect thereby producing d ct ðe uÞ ¼ 0. dt
ð17Þ
Thus, Eq. (12) becomes uT ¼ em ui;n1 ;
ð18Þ
while Eq. (14) becomes uiþ1;nþ1 ¼ em uB .
ð19Þ
Combining Eqs. (18), (13) and (19) yields the exponential angled derivative (EAD) scheme uiþ1;nþ1 ¼ em ½em ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð20Þ
The amplification factor for this scheme satisfies g2 þ ½em ð1 2lÞ ðejh 1Þ g e2m ejh ¼ 0; thus implying that jgj = em = jgexactj. For future reference, 8 9 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < = jh jh 2 1 e 1 e þ ð1 2lÞ gEAD ¼ em ð1 2lÞ þ ejh . : ; 2 2
ð21Þ
ð22Þ
Alternatively, this may be interpreted as a leap-frog approximation centered at (xi+1/2, tn) on an appropriately angled grid. The necessity of exponential fitting is underscored by Fig. 3 where Eqs. (15) and (20) are compared for l = 1, m = 4. Specifically, u(x, t) = ecÆt + ecÆx/b with b = 1, c = 40, Dx = Dt = .1 is displayed at t = 2. For subsequent comparison to the exponential box scheme in the range 0 6 l 6 1, we record here its local truncation error [6,10,13]:
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
669
Steady State: a=.025, b =.025, c = 1.0, d = 0.0 1
x = angled derivative scheme o = exponentially fitted scheme ___ = exact solution
0.5
0
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x Fig. 3. Need for exponential fitting: angled derivative scheme. 3
iþ1=2;n sEAD ¼
ðDxÞ uiþ1=2;n lð2l 1Þð1 lÞ. 12Dt xxx
ð23Þ
A further comment on the dispersive properties of the EAD scheme is appropriate at this junction. Since it is a three-level scheme, it inevitably possesses a ‘‘spurious mode’’ [6] travelling in the wrong direction (westward). However, as long as time-level t = Dt receives propitious treatment (as it will in the EMAD scheme developed below), this ‘‘parasitic solution’’ [13] will not be stimulated. 2.2. Exponential reflected angled derivative (ERAD) scheme Turning now to the range l P 1, perusal of Fig. 4 permits us to derive the reflected angled derivative scheme as follows: ui;nþ1 ui;n b huiþ1;nþ1 þ ui;n ui;nþ1 þ ui1;n i þ ð24Þ þ c ^ui;nþ1=2 ¼ 0. Dx Dt 2 2 If we use the average value ^ ui;nþ1=2 ðui;nþ1 þ ui;n Þ=2 then the following scheme results: uiþ1;nþ1 ¼ ui1;n þ ð1 2=l kÞ ui;nþ1 ð1 2=l þ kÞ ui;n . n jk(iDx)
Upon insertion of the Fourier mode [10] ui,n = g e given by
ð25Þ
into Eq. (25), we find that the amplification factor, g, is
Δ Δ
Fig. 4. Reflected angled derivative scheme (1 6 l 6 1).
670
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
g¼
2 þ lðejh 1 kÞ . 2 þ lðejh 1 þ kÞ
ð26Þ
, so that jgj > 1 for any m > 0 provided that l > 1. Thus, Eq. (25) is Setting h = p in Eq. (26) yields g ¼ 2ð1lÞm 2ð1lÞþm unconditionally unstable for any m > 0! A stable approximation to Eq. (5) may be derived by instead using the average value ^ui;nþ1=2 ðuiþ1;nþ1 þ ui1;n Þ=2 resulting in the scheme: uiþ1;nþ1 ¼
1 ½ð1 kÞ ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ. 1þk
ð27Þ
The amplification factor for this scheme is given by g¼
2 þ l½ð1 kÞejh 1 ; 2 þ l½ð1 þ kÞejh 1
ð28Þ
so that jgj 6 1 for any m > 0 provided that l P 1. The stability of this scheme is thus established for m > 0 and l P 1. However, observe that if k = 1 and l = 2 then g vanishes. In light of the fact that the exact amplification factor of a Fourier mode of Eq. (5) satisfies jgexactj = em, Eq. (27) is seen to be excessively dissipative. This difficulty can be alleviated by the use of the method of characteristics [8]. With reference to Fig. 5, dashed characteristics are traced forward from (xi1, tn) to (xi, tR) and backward from (xi+1, tn+1) to (xi, tL). Along a characteristic, Eq. (5) reduces to the ODE du c þ u ¼ 0. dx b
ð29Þ
This is numerically integrated from xi1 to xi by the trapezoidal rule to produce uR ¼
1 k=2 ui1;n . 1 þ k=2
ð30Þ
Parabolic interpolation is then employed at xi thereby yielding uL ¼ uR þ ð1 2=lÞ ðui;nþ1 ui;n Þ.
ð31Þ
Finally, Eq. (29) is numerically integrated from tn to tn+1 by the trapezoidal rule to produce uiþ1;nþ1 ¼
1 k=2 uL . 1 þ k=2
ð32Þ
Combining Eqs. (30)–(32) yields the scheme 1 k=2 1 k=2 ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ . uiþ1;nþ1 ¼ 1 þ k=2 1 þ k=2
ð33Þ
The amplification factor for this scheme is given by g¼
h2 ejh hð1 2=lÞ ; ejh hð1 2=lÞ
Fig. 5. Characteristic interpretation of RAD scheme (1 6 l 6 1).
ð34Þ
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
671
1k=2 where h ¼ 1þk=2 . Setting h = 0, observe that if l ! 1 with m fixed then g ! 1. This is in sharp contrast to the behavior of jgexactj = em and thus Eq. (33) is also of little general utility. Let us attempt to modify Eq. (33) via exponential fitting [9]. Returning to Eq. (29), we multiply by the integrating factor ecx/b thereby producing
d cx=b ðe uÞ ¼ 0. dx
ð35Þ
Thus, Eq. (30) becomes uR ¼ ek ui1;n ;
ð36Þ
while Eq. (32) becomes uiþ1;nþ1 ¼ ek uL .
ð37Þ
Combining Eqs. (36), (31) and (37) yields the exponentially fitted scheme uiþ1;nþ1 ¼ ek ek ui1;n þ ð1 2=lÞ ðui;nþ1 ui;n Þ .
ð38Þ
The amplification factor for this scheme is given by g¼
e2k ejh ek ð1 2=lÞ ; ejh ek ð1 2=lÞ
ð39Þ
which is identical to Eq. (34) if we replace ek by h which is its [1/1] Pade´ approximant [12]. Hence, Eq. (38) is equally as useless as Eq. (33). The dire nature of our plight is underscored by Fig. 6, which displays jgj/jgexactj for Eq. (38) with l = 1.5, m = 3.0 and l = 4.0, m = 2.0. It is now painfully obvious that exponential fitting in the case of l P 1 is a more delicate matter than it might first appear. Consequently, we take the following more circuitous route. Prior to the application of the method of characteristics, we multiply Eq. (5) by the integrating factor ect and introduce the auxiliary variable v(x, t) :¼ ect Æ u(x, t). This results in the pure advection problem ov ov þ b ¼ 0; ot ox
ð40Þ
which may then be approximated by the reflected angled derivative (RAD) scheme [4] viþ1;nþ1 ¼ vi1;n þ ð1 2=lÞ ðvi;nþ1 vi;n Þ.
ð41Þ
μ = 1.5, ν = 3.0 90
μ = 4.0, ν = 2.0 90
2
120
60
1
150
30
180
330
300
240 270
60
2
150
0
210
4
120
30
180
0
210
330
300
240 270
Fig. 6. Relative magnitude of Eq. (38) amplification factor.
672
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
In terms of the original variable u(x, t), this produces the exponential reflected angled derivative (ERAD) scheme uiþ1;nþ1 ¼ em ui1;n þ ð1 2=lÞ ðui;nþ1 em ui;n Þ;
ð42Þ
so named because its finite difference stencil is topologically equivalent to that of the EAD scheme reflected about the northeasterly diagonal. Alternatively, this may be interpreted as a leap-frog approximation centered at (xi, tn+1/2) on an appropriately angled grid. The amplification factor of this scheme is 2 þ l ðejh 1Þ gERAD ¼ em ; ð43Þ 2 þ l ðejh 1Þ thus implying that jgj = em = jgexactj. Observe that an analogous derivation of Eq. (20) could be given in the case of 0 6 l 6 1. Thus, in order to produce a nondissipative scheme, l P 1 requires multiplication by the integrating factor prior to the application of the method of characteristics whereas for 0 6 l 6 1 the order of these two operations is immaterial. For subsequent comparison to the exponential box scheme in the range l P 1, we record here its local truncation error [6,10,13]: i;nþ1=2
sERAD ¼
ðDxÞ3 i;nþ1=2 u lð2l 1Þð1 lÞ. 12Dt xxx
ð44Þ
As is clear on the grounds of symmetry, a prudent treatment of the grid line x = Dx (as provided by the EMAD scheme below) is necessary when employing the ERAD scheme. 2.3. Exponential box (EBOX) scheme With reference to Fig. 7, the box scheme may be derived as follows: 1 huiþ1;nþ1 uiþ1;n ui;nþ1 ui;n i b huiþ1;nþ1 ui;nþ1 uiþ1;n ui;n i þ þ þ þ c ^uiþ1=2;nþ1=2 ¼ 0. 2 2 Dt Dt Dx Dx
ð45Þ
If we use the average value ^ uiþ1=2;nþ1=2 ðuiþ1;nþ1 þ ui;nþ1 þ uiþ1;n þ ui;n Þ=4 then the following scheme results: uiþ1;nþ1 ¼
1 ½ð1 þ l m=2Þ ui;n þ ð1 l m=2Þ uiþ1;n ð1 l þ m=2Þ ui;nþ1 . 1 þ l þ m=2
ð46Þ
Upon insertion of the Fourier mode [10] ui,n = gnejk(iDx) into Eq. (46), we find that the amplification factor, g, is given by g¼
ð1 þ l m=2Þ þ ð1 l m=2Þejh ; ð1 l þ m=2Þ þ ð1 þ l þ m=2Þejh
ð47Þ
thereby implying that jgj 6 1 and thus that the box scheme, Eq. (46), is stable for any 0 6 l 6 1 provided that m P 0. However, its unsuitability is evident from Fig. 8, which displays jgj/jgexactj for l = 0.5, m = 2.0 and l = 2.0, m = 2.0.
Δ
Δ
Fig. 7. Box scheme (0 6 l 6 1).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 μ = 0.5, ν = 2.0 90
μ = 2.0, ν = 2.0 90
10
120
60
5
150
30
60
5
150
0
210
10
120
180
330
30
180
0
210
300
240
673
330
300
240
270
270
Fig. 8. Relative magnitude of box scheme amplification factor.
μ <−
μ −>
Fig. 9. Characteristic interpretation of box scheme (0 6 l 6 1).
With reference to Fig. 9, we may appeal to exponential fitting in order to resolve this dilemma. We first multiply Eq. (5) by the integrating factor ect and then introduce the auxiliary variable v(x, t) :¼ ect Æ u(x, t), resulting in the pure advection problem ov ov þ b ¼ 0; ot ox
ð48Þ
which may then be approximated by the method of characteristics (i.e., the box scheme) viþ1;nþ1 ¼ vi;n þ
1l ðviþ1;n vi;nþ1 Þ. 1þl
ð49Þ
In terms of the original variable u(x, t), this produces the exponential box (EBOX) scheme uiþ1;nþ1 ¼ em ui;n þ
1l ðem uiþ1;n ui;nþ1 Þ. 1þl
The amplification factor of this scheme is ð1 þ lÞ þ ð1 lÞejh gEBOX ¼ em ; ð1 lÞ þ ð1 þ lÞejh
ð50Þ
ð51Þ
thus implying that jgj = em = jgexactj. In order to compare the dispersive properties of the EAD, EBOX, and ERAD schemes, we consider their relative phase errors [10]. The relative phase errors are plotted in Fig. 10 for various values of 0 6 l 6 1. In this figure, the top curve corresponds to the EBOX scheme and the bottom curve to the EAD scheme. The adequacy of the EBOX scheme in the neighborhood of l = 1 is readily apparent as is the overall superiority of the EAD scheme. Observe that the EAD scheme is exact for l = 1/2 while both schemes are exact for l = 1. The relative phase errors are plotted in Fig. 11 for various values of
674
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 μ=.25
μ=.50 90 2
90 4 120
60 2
150
120 30
180
330 240
0
210
330 240
μ=1.0
90 2
90 1 60
1
150
0
330
30
180
0
210
330 240
300
270
0.5
150
30
210
60
120
180
240
300
270
μ=.75 120
30
180
300
270
1
150
0
210
60
300
270
Fig. 10. Comparison of relative phase errors (l 6 1) (top curve = EBOX scheme, bottom curve = EAD scheme).
μ=1.0 120
μ=1.5
90 1
120
60 0.5
150
30
180
0
330
210 240
1
0
330
210 240
μ=4.0 120
60 0.5
0
210
330 300 270
90 1
60
0.5
150
30
180
240
300
270
90 1
150
30
180
μ=2.0 120
60
150
300
270
90 2
30
0
180
210
330 240
270
300
Fig. 11. Comparison of relative phase errors (l P 1) (top curve = ERAD scheme, bottom curve = EBOX scheme).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
675
1 6 l 6 1. In this figure, the top curve corresponds to the ERAD scheme and the bottom curve to the EBOX scheme. The adequacy of the EBOX scheme in the neighborhood of l = 1 is once again readily apparent as is the overall superiority of the ERAD scheme. Observe that the ERAD scheme is exact for l = 2 while both schemes are exact for l = 1. As a vehicle for further assessing the relative merits of the EAD, EBOX, and ERAD schemes, consider their respective local truncation errors [6,10,13]: 3
iþ1=2;nþ1=2
sEBOX
¼
ðDxÞ uiþ1=2;nþ1=2 lðl 1Þðl þ 1Þ; 12Dt xxx
ð52Þ
together with Eqs. (23) and (44). As they all possess a common factor, we introduce the scaled local truncation errors t :¼
12Dt ðDxÞ3 uxxx
s.
ð53Þ
The scaled local truncation errors are plotted in Fig. 12 for the range 0 6 l 6 1. Observe the suitability of the EBOX scheme in the vicinity of l = 1 as well as the overall superiority of the EAD scheme in this regime. The scaled local truncation errors are plotted in Fig. 13 for the range 1 6 l 6 1. Observe once again the suitability of the EBOX scheme in the vicinity of l = 1 as well as the overall superiority of the ERAD scheme in this regime. Based upon these preliminary observations, we now turn to the construction of a new scheme for the nonhomogeneous, variable coefficient Eq. (3), which is a hybrid of the EAD, EBOX, and ERAD schemes that is nondissipative with good dispersive properties, second-order accurate, explicit, and unconditionally stable. 3. Exponential method of angled derivatives (EMAD) In light of the foregoing introductory observations, it would seem that the best course of action would be to use the EAD scheme, Eq. (20), if l 6 1 and the ERAD scheme, Eq. (42), if l P 1. This is indeed the case for
SCALED LOCAL TRUNCATION ERROR
0.05 EAD Scheme 0
−0.05
t
−0.1 −0.15 −0.2 EBOX Scheme −0.25 −0.3 −0.35 0
0.1
0.2
0.3
0.4
0.5
μ
0.6
0.7
0.8
Fig. 12. Comparison of scaled local truncation errors (l 6 1).
0.9
1
676
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 SCALED LOCAL TRUNCATION ERROR 25
20 EBOX Scheme
15
t
ERAD Scheme 10
5
0
1
1.5
2
2.5
μ
3
3.5
4
Fig. 13. Comparison of scaled local truncation errors (l P 1).
constant b in Eq. (3). However, if b is not constant then difficulties arise. With reference to Fig. 2, if neither li+1/2,n nor li,n+1/2 exceed one then the EAD scheme is stable and should be used. Likewise, with reference to Fig. 5, if both li+1/2,n and li,n+1/2 exceed one then the ERAD scheme is stable and should be used. Yet, consider the scenario pictured in Fig. 14. There, both the EAD and the ERAD schemes are stable and a seemingly arbitrary choice must be made between them. While this perplexity is venial at worst, that displayed in Fig. 15 is decidedly mortal. In this case, no matter which of the two schemes is utilized, the choice is condemned to instability. The resolution of these conundrums is provided by the employment of the EBOX scheme in these two instances. Since in both cases l passes through unity somewhere in the shaded cell, the EBOX scheme is an excellent choice because of its exemplary behaviour in the neighborhood of l = 1. We now develop these ideas into a complete algorithm for the initial-boundary value problem, Eqs. (2) and (3). In order to distinguish it from the EAD and ERAD schemes, which are merely components of this new scheme, it is dubbed the exponential method of angled derivatives (EMAD). An important feature of this development will be the preservation of the desirable properties described in the previous section, which are strictly applicable only to the homogeneous, constant coefficient problem, as we extend the component EAD, EBOX, and ERAD schemes to the nonhomogeneous, variable coefficient problem. Particular attention will be paid to approximation at the grid points adjacent to the spatiotemporal boundary.
μ
μ< Fig. 14. Indecisive scenario.
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
677
μ<
μ Fig. 15. Checkmate!
In the ensuing subsections, we extensively utilize the Bernoulli function z ; BðzÞ ¼ z e 1
ð54Þ
and the associated function 1 BðzÞ . ð55Þ z Their essential role in exponential fitting as well their stable numerical evaluation for small z are discussed in ^ and ^c, which will be specified [9]. Furthermore, we will freeze the coefficients of Eq. (3) at constant values, b ^ct below. In terms of the auxiliary variable vðx; tÞ :¼ e uðx; tÞ, Eq. (3) then becomes MðzÞ ¼
ov ^ ov þ b ¼ e^ct f ðx; tÞ. ot ox
ð56Þ
Subsequent application of the MAD scheme [4] to Eq. (56) with special attention to the required exponential quadrature on its right hand side will yield the EMAD scheme. 3.1. Nonhomogeneous, variable coefficient EAD scheme Let us begin with the case where neither li+1/2,n nor li,n+1/2 exceed one so that the EAD scheme is the ^ :¼ b ^ appropriate choice (see Fig. 2). Defining b c :¼ ciþ1=2;n and l :¼ bDt=Dx, m :¼ ^cDt, the dashed lines iþ1=2;n , ^ ^ In actuality, the characteristics are curved but, due to the symare approximate characteristics with slope 1=b. metry of the finite difference stencil about the point (xi+1/2, tn), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv ¼ e^ct f . dt
ð57Þ
Integrating this along the characteristics produces Z ðxT ;tn Þ m uT ¼ e ui;n1 þ e^cðtn tÞ f dt;
ð58Þ
ðxi ;tn1 Þ
uiþ1;nþ1 ¼ em uB þ
Z
ðxiþ1 ;tnþ1 Þ
e^cðtnþ1 tÞ f dt;
ð59Þ
ðxB ;tn Þ
where xT = xi + lDx and xB = xi + (1 l)Dx. Linear interpolation of the forcing term produces fT ¼ ð1 lÞfi;n þ lfiþ1;n ;
f B ¼ lfi;n þ ð1 lÞfiþ1;n .
ð60Þ
Employing the exponential trapezoidal rule [9] to perform the quadratures in Eqs. (58) and (59) now produces
678
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
Dt ½MðmÞ fi;n1 þ MðmÞ fT ; BðmÞ Dt ½MðmÞ fB þ MðmÞ fiþ1;nþ1 . ¼ em uB þ BðmÞ
uT ¼ em ui;n1 þ
ð61Þ
uiþ1;nþ1
ð62Þ
Parabolic interpolation to the values of v at (xi, tn), (xT, tn), (xi+1, tn) and subsequent evaluation at (xB, tn) produces uB ¼ uT þ ð1 2lÞ ðuiþ1;n ui;n Þ.
ð63Þ
Combination of Eqs. (60)–(63) now yields the EAD scheme: uiþ1;nþ1 ¼ em ½em ui;n1 þ ð1 2lÞ ðuiþ1;n ui;n Þ Dt fem MðmÞ fi;n1 þ ½em MðmÞ ð1 lÞ þ MðmÞ l fi;n þ BðmÞ þ ½em MðmÞ l þ MðmÞ ð1 lÞ fiþ1;n þ MðmÞ fiþ1;nþ1 g.
ð64Þ
3.2. Nonhomogeneous, variable coefficient ERAD scheme Next, we turn to the case where both li+1/2,n and li,n+1/2 exceed one so that the ERAD scheme is the appro^ :¼ b ^ ^ the dashed priate choice (see Fig. 5). Defining b c :¼ ci;nþ1=2 and l :¼ bDt=Dx, m :¼ ^cDt, k :¼ ^cDx=b, i;nþ1=2 , ^ ^ In actuality, the characteristics are curved but, due to the lines are approximate characteristics with slope 1=b. symmetry of the finite difference stencil about the point (xi, tn+1/2), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv 1 ^ct ¼ e f. ð65Þ ^ dx b Integrating this along the characteristics produces Z ðxi ;tR Þ uR ¼ ek ui1;n þ e^cðtR tÞ f dt;
ð66Þ
ðxi1 ;tn Þ
uiþ1;nþ1 ¼ ek uL þ
Z
ðxiþ1 ;tnþ1 Þ
e^cðtnþ1 tÞ f dt;
ð67Þ
ðxi ;tL Þ
where tR ¼ tn þ l1 Dt and tL ¼ tn þ ð1 l1ÞDt. Linear interpolation of the forcing term produces 1 1 1 1 fR ¼ 1 f L ¼ fi;n þ 1 fi;n þ fi;nþ1 ; fi;nþ1 . l l l l
ð68Þ
Employing the exponential trapezoidal rule [9] to perform the quadratures in Eqs. (66) and (67) now produces Dt=l ½MðkÞ fi1;n þ MðkÞ fR ; BðkÞ Dt=l ½MðkÞ fL þ MðkÞ fiþ1;nþ1 . ¼ ek uL þ BðkÞ
uR ¼ ek ui1;n þ
ð69Þ
uiþ1;nþ1
ð70Þ
Parabolic interpolation to the values of v at (xi, tn), (xi, tR), (xi, tn+1) and subsequent evaluation at (xi, tL) produces 2 mþ2k uL ¼ e uR þ 1 ð71Þ ðek uiþ1;n emþk ui;n Þ. l Combination of Eqs. (68)–(71) now yields the scheme:
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
uiþ1;nþ1
2 ðui;nþ1 em ui;n Þ ¼ e ui1;n þ 1 l
Dt=l 1 1 km km e MðkÞ fi1;n þ e MðkÞ 1 þ þ MðkÞ fi;n BðkÞ l l 1 1 þ ekm MðkÞ þ MðkÞ 1 fi;nþ1 þ MðkÞ fiþ1;nþ1 . l l
679
m
ð72Þ
Observe that Eq. (72) does not reproduce the constant solution u fc for constant c and f. Instead, it produces BðmÞ ekm þ1 f u BðkÞ 2 c. Thus, we modify it as follows to yield the ERAD scheme: 2 m uiþ1;nþ1 ¼ e ui1;n þ 1 ðui;nþ1 em ui;n Þ l
Dt=l 2 1 1 km þ ekm MðkÞ fi1;n þ ekm MðkÞ 1 þ MðkÞ fi;n BðmÞ e þ 1 l l 1 1 þ ekm MðkÞ þ MðkÞ 1 fi;nþ1 þ MðkÞ fiþ1;nþ1 . ð73Þ l l Eq. (73) has the same formal order of accuracy and dissipative properties as Eq. (72) yet possesses the correct constant solution. 3.3. Nonhomogeneous, variable coefficient EBOX scheme Finally, suppose that one of li+1/2,n and li,n+1/2 exceeds one while the other does not. In this case, the ^ :¼ b EBOX scheme is the appropriate choice (see Fig. 9). Defining b c :¼ ciþ1=2;nþ1=2 and iþ1=2;nþ1=2 , ^ ^ ^ ^ In actuall :¼ bDt=Dx, m :¼ ^cDt, k :¼ ^cDx=b, the dashed lines are approximate characteristics with slope 1=b. ity, the characteristics are curved but, due to the symmetry of the finite difference stencil about the point (xi+1/2, tn+1/2), first-order errors cancel and these linear approximations to the true characteristics are second-order accurate. 3.3.1. EBOX scheme (0 6 l 6 1) First consider the case l 6 1. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv ¼ e^ct f . dt Integrating this along the characteristic produces Z ðx2 ;tnþ1 Þ e^ct f dt; v2 ¼ v1 þ
ð74Þ
ð75Þ
ðx1 ;tn Þ
where x1 = xi + (1 l)Dx/2 and x2 = xi + (1 + l)Dx/2. Linear interpolation of the forcing term produces 1þl 1l 1l 1þl fi;n þ fiþ1;n ; fi;nþ1 þ fiþ1;nþ1 ; f2 ¼ ð76Þ f1 ¼ 2 2 2 2 while linear interpolation of the auxiliary variable produces 1þl 1l 1l 1þl vi;n þ viþ1;n ; vi;nþ1 þ viþ1;nþ1 . v2 ¼ ð77Þ v1 ¼ 2 2 2 2 Employing the exponential trapezoidal rule [9] to perform the quadrature in Eq. (75) now yields the EBOX scheme: 1l Dt ðem uiþ1;n ui;nþ1 Þ þ uiþ1;nþ1 ¼ em ui;n þ 1þl BðmÞ 1l fMðmÞ fiþ1;n þ MðmÞ fi;nþ1 g . MðmÞ fi;n þ MðmÞ fiþ1;nþ1 þ 1þl
ð78Þ
680
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
3.3.2. EBOX scheme (1 6 l 6 1) Next consider the case l P 1. Along a characteristic, the PDE, Eq. (56), reduces to the ODE dv 1 ^ct ¼ e f. ^ dx b
ð79Þ
Integrating this along the characteristic produces Z ðxiþ1 ;t2 Þ v2 ¼ v1 þ e^ct f dt;
ð80Þ
ðxi ;t1 Þ
where t1 = tn + (1 1/l)Dt/2 and t2 = tn + (1 + 1/l)D t/2. Linear interpolation of the forcing term produces 1 þ 1=l 1 1=l 1 1=l 1 þ 1=l fi;n þ fi;nþ1 ; fiþ1;n þ fiþ1;nþ1 ; f2 ¼ 2 2 2 2 while linear interpolation of the auxiliary variable produces f1 ¼
v1 ¼
1 þ 1=l 1 1=l vi;n þ vi;nþ1 ; 2 2
v2 ¼
1 1=l 1 þ 1=l viþ1;n þ viþ1;nþ1 . 2 2
ð81Þ
ð82Þ
Employing the exponential trapezoidal rule [9] to perform the quadrature in Eq. (80) now yields the scheme: 1l Dt eðkmÞ=2 ðem uiþ1;n ui;nþ1 Þ þ uiþ1;nþ1 ¼ em ui;n þ 1þl l BðkÞ 1l fMðkÞ fi;nþ1 þ MðkÞ fiþ1;n g . MðkÞ fi;n þ MðkÞ fiþ1;nþ1 1þl
ð83Þ
Observe that Eq. (83) does not reproduce the constant solution u fc for constant c and f. Instead, it produces BðmÞ u BðkÞ eðkmÞ=2 fc . Thus, we modify it as follows to yield the EBOX scheme: 1l Dt=l uiþ1;nþ1 ¼ em ui;n þ ðem uiþ1;n ui;nþ1 Þ þ 1þl BðmÞ 1l fMðkÞ fi;nþ1 þ MðkÞ fiþ1;n g . MðkÞ fi;n þ MðkÞ fiþ1;nþ1 1þl
ð84Þ
Eq. (84) has the same formal order of accuracy and dissipative properties as Eq. (83) yet possesses the correct constant solution. 3.4. Boundary treatment Depending on the numerical value of l, the EMAD approximation at grid points adjacent to the spatiotemporal boundary may require modification (see Fig. 16). Along t = 0 the EAD scheme is inapplicable while the same is true for the ERAD scheme along x = 0. Even though we know the exact solution there (see Eq.
μ
μ μ<
Fig. 16. Boundary treatment.
μ<
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
681
(2)), tracing a linear approximation to a characteristic from either t = Dt or x = Dx backwards to the spatiotemporal boundary would yield only first-order accuracy unless the coefficient b is constant in which case the characteristics are in fact straight. However, a parabolic approximation to the characteristics yields secondorder accuracy. 3.4.1. Southwestern corner (0 6 l 6 1) Let us begin with the southwestern-most grid cell (i = 1, n = 1). If l < 1 throughout this cell then, starting at (xi+1, tn+1), we trace backwards to t = 0 along a parabolic trajectory x(t) whose first two derivatives coincide with those of the actual characteristic through this point: 2
xB ¼ xiþ1 biþ1;nþ1 Dt þ ðbt þ bbx Þiþ1;nþ1
ðDtÞ . 2
ð85Þ
Then, as above, we integrate Eq. (74) along this trajectory with ^c :¼ cððxB þ xiþ1 Þ=2; Dt=2Þ and m :¼ ^cDt thereby producing the approximation uiþ1;nþ1 ¼ em /ðxB Þ þ
Dt ½MðmÞ fB þ MðmÞ fiþ1;nþ1 . BðmÞ
ð86Þ
In place of fB, we may instead employ the linear interpolant xiþ1 xB xB xi fi;n þ fiþ1;n ; f^ B ¼ Dx Dx
ð87Þ
without sacrifice of order of accuracy. In order to insure that the approximate characteristic slope, s, is correspondingly bounded (0 6 1/s 6 Dx/Dt), we must either impose the restriction 8 b db > > if > 0; > < db=dt dt iþ1;nþ1 iþ1;nþ1 Dt 6 > b db > > ð1=liþ1;nþ1 1Þ if < 0; : db=dt iþ1;nþ1 dt iþ1;nþ1 or utilize additional terms in the truncated Taylor series, Eq. (85). 3.4.2. Southwestern corner (1 6 l 6 1) On the other hand, if l > 1 throughout this cell then, starting at (xi+1, tn+1), we trace backwards to x = 0 along a parabolic trajectory t(x) whose first two derivatives coincide with those of the actual characteristic through this point: tL ¼ tnþ1
1 biþ1;nþ1
Dx
ðbt þ bbx Þiþ1;nþ1 ðDxÞ2 . 2 b3iþ1;nþ1
ð88Þ
Then, as above, we integrate Eq. (79) along this trajectory with ^c :¼ cðDx=2; ðtL þ tnþ1 Þ=2Þ, D :¼ tn+1 tL and ^m :¼ ^cD thereby producing the approximation uiþ1;nþ1 ¼ e^m wðtL Þ þ
D ½Mð^mÞ fB þ Mð^mÞ fiþ1;nþ1 . Bð^mÞ
In place of fL, we may instead employ the linear interpolant tnþ1 tL tL tn fi;n þ fi;nþ1 ; f^L ¼ Dt Dt
ð89Þ
ð90Þ
without sacrifice of order of accuracy. In order to insure that the approximate characteristic slope, s, is correspondingly bounded (0 6 s 6 Dt/Dx), we must either impose the restriction
682
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
8 b > > > < db=dx iþ1;nþ1 Dx 6 > b > > ðliþ1;nþ1 1Þ : db=dx iþ1;nþ1
db if < 0; dx iþ1;nþ1 db if > 0; dx iþ1;nþ1
or utilize additional terms in the truncated Taylor series, Eq. (88). 3.4.3. Remaining portion of the boundary In grid cells elsewhere along the southern boundary (i 5 1, n = 1) of the mesh, we use Eq. (86) if l < 1 throughout the cell, the ERAD scheme if l > 1 throughout the cell, and the EBOX scheme otherwise. Likewise, in grid cells along the western boundary (i = 1, n 5 1) of the mesh, we use Eq. (89) if l > 1 throughout the cell, the EAD scheme if l < 1 throughout the cell, and the EBOX scheme otherwise. This completes the specification of the EMAD algorithm. It is an explicit, three-level, unconditionally stable scheme that is nondissipative with good dispersive properties while maintaining second-order accuracy even at the boundary and in spite of the presence of variable coefficients and/or forcing function. We now present numerical evidence of the versatility of this new scheme. 4. Numerical examples 4.1. Spatiotemporally varying velocity with source Our first example is nonhomogeneous with forcing term f ðx; tÞ ¼ 2p ½2tð1 þ xÞ cosð2pxÞ cosð2ptÞ sinð2pxÞ sinð2ptÞ þ sinð2pxÞ cosð2ptÞ.
t=0.5 1
u
0.5
0
−0.5 −1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.6
0.7
0.8
0.9
1
x t=1.0 1
u
0.5
0
−0.5 −1
0
0.1
0.2
0.3
0.4
0.5
x Fig. 17. Spatiotemporally varying velocity with source (Dt = .05, Dx = .05).
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
683
Also, there is a coefficient b(x, t) = 2t(1 + x) depending upon both space and time with c = 2p. The initial condition is /(x) = sin(2px) with boundary condition w(t) = 0. The exact solution is then u(x, t) = sin(2px) Æ cos(2pt). For this example, we have chosen Dt = .05 and Dx = .05. The numerical results for t = 0.5, halfway through a cycle, and t = 1.0, at the completion of a cycle, are displayed as dots in Fig. 17 where the solid line represents the exact solution. 4.2. Countercurrent wave propagation Our next example is of importance in the mathematical modelling of the warmup process for automotive catalytic converters [2,3]. In spite of the fact that the flow of air is down the length of the converter, chemical reactions along its surface give rise to a backward propagating wave known as the light-off point. Herewith, we model such countercurrent wave propagation through the forcing term 2 12 f ðx; tÞ ¼ 50 erfc½6ðx þ t 2Þ pffiffiffi e½6ðxþt2Þ ; p with b = 1 and c = 100. The initial condition is /ðxÞ ¼ 12 erfc½6ðx 2Þ with boundary condition wðtÞ ¼ 12 erfc½6ðt 2Þ. The exact solution is now uðx; tÞ ¼ 12 erfc½6ðx þ t 2Þ. For this example, we have chosen Dt = .025 and Dx = .025. The numerical results for t = 0.0, 1.0, 1.5, 2.0 are displayed as dots in Fig. 18 where the solid line represents the exact solution. The upstream propagation of the complementary error function profile is notable. 4.3. Oscillatory flow Our final example also arises in the numerical simulation of automotive catalytic converters [3]. If the converter is mounted in proximity to the exhaust manifold then the driving fluid flow decidedly pulsates. We
t = 0.0
t = 1.0
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
0.6
0.8
0
1
0
0.2
0.4
x
0.6
0.8
1
0.8
1
x
t = 1.5
t = 2.0
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
0.6
x
0.8
1
0
0
0.2
0.4
0.6
x
Fig. 18. Countercurrent wave propagation (Dt = .025, Dx = .025).
684
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685 t = .06
t = .12
t = .18
1
1
1
0.5
0.5
0.5
0
0
0.5
1
0
0
t = .24
0.5
1
0
1
1
0.5
0.5
0.5
0
0.5
1
0
0
t = .42
0.5
1
0
1
0.5
0.5
0.5
0.5
1
0
0
0.5
1
0.5
1
t = .54
1
0
0
t = .48
1
0
0.5 t = .36
1
0
0
t = .30
1
0
0
0.5
1
Fig. 19. Oscillatory flow.
model this through the oscillating inlet velocity bðtÞ ¼ 1 þ 12 cos 5p t with c = 10 and f = 10. The initial con3 dition is /(x) = 0 with inlet boundary condition w(t) = 0. The exact solution is then 8 3 5p > 10ðtnÞ > sin t ; 1 e ; x 6 t þ < 10p 3 uðx; tÞ ¼ > 3 5p > : 1 e10t ; sin t ; xPtþ 10p 3 where n satisfies the transcendental equation 3 5p 3 5p nþ sin n þx tþ sin t ¼ 0. 10p 3 10p 3 For this example, we have chosen Dt = .01 and Dx = .02. Thus, every 12 time steps is equivalent to a complete pulsation cycle. The numerical results for t = .06, .12, .18, .24, .30, .36, .42, .48, .54 are displayed as dots in Fig. 19 where the solid line represents the exact solution, so that each successive frame is equivalent to half a cycle. Observe the periodic response, which develops once the transients have decayed. 5. Conclusion In this paper, we have presented a hybrid scheme for the numerical approximation of the nonhomogeneous one-way advection-reaction equation. This exponential method of angled derivatives (EMAD) scheme provides an explicit, three-level, nondissipative, second-order accurate, unconditionally stable numerical approximation to this important initial-boundary value problem. This new scheme provides the natural generalization of the method of angled derivatives (MAD) scheme [4] for pure advection problems.
B.J. McCartin / Applied Mathematics and Computation 177 (2006) 665–685
685
We began by developing the EAD scheme for 0 6 l 6 1, where its interpretation in terms of characteristics, the explicit computation of its local truncation error, and the analysis of its dispersive properties took center stage. In addition, we have introduced the ERAD scheme for 1 6 l 6 1 and provided a detailed analysis of its numerical properties. The net result of this preliminary analysis was that, for the homogeneous, constant coefficient problem, the EAD and ERAD schemes provide superlative approximations within their respective parameter regimes. An ancillary benefit of the preceding analysis is a unified approach via exponential fitting and characteristics to the seemingly disparate box, angled derivative, and reflected angled derivative schemes. Likewise, the Lax–Wendroff, second-order upwind difference, Roberts–Weiss, and leap-frog schemes [5,6,8,10,11,13], as well as any other consistent scheme for Eq. (1), may also be so treated. For the nonhomogeneous, variable coefficient problem, the suitably extended EAD and ERAD schemes were then melded together using the EBOX scheme (whose characteristic features were also herein developed) as adhesive in order to provide a complete numerical algorithm unrestricted in its value of l. Particular attention has been paid to maintaining the order of accuracy and stability of the approximation in the vicinity of the spatiotemporal boundary. This new EMAD scheme was then exercised on a pair of important test problems in order to demonstrate its versatility. Acknowledgement The author would like to express his sincere gratitude to Mrs. Barbara A. McCartin for her dedicated assistance in the production of this paper. References [1] H. Rhee, R. Aris, N.R. Amundson, First Order Partial Differential Equations, Theory and Application of Single Equations, vol. I, Dover, New York, 2001. [2] B.A. Finlayson, Numerical Methods for Problems with Moving Fronts, Ravenna Park, Seattle, WA, 1992. [3] B.J. McCartin, P.D. Young, Numerical investigation of the validity of the quasi-static approximation in the modelling of catalytic converters, in: D. Ferguson, T. Peters (Eds.), Mathematics and Industry, SIAM, Philadelphia, PA, 2005. [4] B.J. McCartin, The method of angled derivatives, Applied Mathematics and Computation 170 (2005) 440–461. [5] A. Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Cambridge, 1996. [6] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, 1994. [7] B.J. McCartin, Seven deadly sins of numerical computation, American Mathematical Monthly 105 (10) (1998) 929–941. [8] J.D. Hoffman, Numerical Methods for Engineers and Scientists, McGraw-Hill, New York, 1992. [9] B.J. McCartin, Exponential fitting of the delayed recruitment/renewal equation, Journal of Computational and Applied Mathematics 136 (2001) 343–356. [10] J.C. Tannehill, D.A. Anderson, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, second ed., Taylor & Francis, New York, 1997. [11] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth, Belmont, CA, 1989. [12] G.A. Baker, P. Graves-Morris, Pade´ Approximants, second ed., Cambridge University Press, Cambridge, 1996. [13] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, Wiley-Interscience, New York, 1995.