- Email: [email protected]
A second-order accurate difference method for systems of hyperbolic partial differential equations
COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 1 (1972) 173- 188. 0 NORTH-HOLLAND PUBLISHING COMPANY
A SECOND-ORDER ACCURATE DIFFERENCE METHOD FOR SYSTEMS OF HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS 1 S. RANGANATH 2 and R.J. CLIFTON 3
Received 16 November 1971 A second order accurate difference method is presented for systems of first order hyperbolic differential equations. The method is analogous to the Courant, Isaacson, Rees (CIR) method, except that the error introduced in one time step is O(Ats) instead of O(At’) as is the case for the CIR method. Convergence of the proposed method is established. The cumulative error at a fixed time is shown to be O(APt*).The proposed method is compared with several other second order accurate methods by considering in detail the special case of the system of equations governing flexural wave propagation in elastic beams. These comparisons indicate that the proposed method has substantial advantages over the other methods considered in that the method is computationally stable for larger mesh sizes. As a result, less computing time is required to obtain solutions with a given accuracy.
1. Introduction In 1953, Courant, Isaacson and Rees [ 11 presented a difference method for the solution of quasi-linear hyperbolic partial differential equations of first order. The method makes use of a fixed mesh and integration along characteristics to determine the solution at an arbitrary mesh point from data at neighboring mesh points at the previous time step. The error introduced in going from the exact solution at time t to the difference solution at time t + At was shown to be O(At2) which indicates that the method corresponds to what is generally called a first order accurate method. In order to save computing time for a given accuracy it is frequently advantageous to make use of higher order accurate difference schemes for which the same accuracy can be obtained with a larger mesh size. A well-known example of such higher order accurate schemes is the second order accurate difference method developed by Lax and Wendroff [2] . The LaxWendroff method is based on expanding the solution vector in a Taylor’s series and retaining terms up to second order. Richtmeyer [ 3 I has presented a two step Lax-Wendroff scheme which is also second order accurate but slightly simpler to use. When the coefficients are constant it reduces to the original Lax-Wendroff scheme [2]. The generalized Heun’s method [4] is also a two-step difference scheme which is second order accurate. Stetter [51 has constructed higher order accurate methods based on characteristic finite difference schemes for systems of equations
‘. ‘IhIs research was supported by the Air Force Materials Laboratory, Wright-Patterson Airforce Base under Contract No. F 33615-70-C-1101 with Brown University. 2. Research Associate, Brown University. Now at Advanced Technology Center, Inc., Dallas, Texas. 3. Professor of Engineering, Brown University.
174
S. Rangmath, R.J. Oifton, A second-order accurate difference method
which have only two independent characteristics. Using the Adam’s procedure for ordinary differential equations Ansorge 161 has presented higher order accurate methods for semilinear equations. We present here a second order accurate method which is a natural extension of the first order method proposed by Courant, Isaacson and Rees [ I] for general quasilinear hyperbolic systems in two independent variables. This method (referred to hereafter as the second order accurate CIR method) is based on integration along characteristics, preserving second order accuracy throughout. The purpose of this paper is to describe the second order accurate CIR method and establish the convergence of the difference scheme. This method is then compared with other second order accurate methods by considering in detail the special case of the propagation of flexural waves in elastic beam. It will be shown that, at least for this case, the new method has definite advantages over other second order accurate methods.
2. Description of the difference method In this section we develop a difference method for quasi-linear systems of hyperbolic differential equations in two independent variables. Convergence of the method will be established in the next section. Quasi-linear systems of first order partial differential equations in two independent variables have the form * *$ {A%;
+ B’h) + di} = 0
j=l,---,y1
(2.1)
where A’j, Bii and di are in general, functions of the (n+2) variables (x, t, u’). If the system is hyperbolic, n distinct linear combinations of equations (2.1) may be formed to produce an equivalent system of equations in the normal form n
C
aji(uf+
,iu;
) =:
bi
j= 1,2, ---,n
i=l
(2.2)
where aif, ci and bj are functions of x, t and ui. (We assume here that the direction t = constant is not a characteristic.) In the j-th equation of the normal form (2.2) the terms (uf+cj,L) are simply the derivatives, du’fdt, of ui along the curve x = x(t) for which dx/dt = ci (the j-th characteristic direction). Thus, (2.2) can now be written in the form ze a’idu’ = bidt
(2.3)
where du’ represents the change in ui along the j-th characteristic for a change dt in t. We assume that the coefficients in (2.1) and the initial data are sufficiently smooth for a unique solution, u’(t,x) E C3 to exist in some finite region t 2 0. In particular we shall assume det aii # 0 * Subscripts x and 2 denote differentiation with respect to n and t respectively.
S. Ranganath, R.J. CIifton, A second-order accurate difference method
175
Fig. 1. Grid and Characteristics.
and &, cj and bj are Lipschitz continuous functions of all their arguments. Our objective is to develop a finite difference method for computing functions w’(t,x) which can be made arbitrarily close to u’(t,x) in some finite region t > 0 by taking the mesh size sufficiently small. We select a mesh width Ax and consider a rectangular network of points in the x-t plane. From the known difference solution at mesh points on the line IV, we obtain the solution at mesh points on the line IV+i by integration of (2.3) along characteristics. (Fig. 1). The time step between the lines IV and IV+i is chosen such that, if P is a net point on IV+i and Q its nearest neighboring point on IV, then the line segments drawn backwards through P with slopes dx/dt = cj intersect IV at points which lie between Q and its nearest neighboring points Q’ and Q” on 1”. In this way the Courant, Friedrichs, Lewy [ 71 necessary condition for stability is automatically satisfied. The difference solution at a mesh point P is computed by means of a two-step procedure. First, a provisional solution G’(P) is obtained using the CIR first order accurate method. Then, this solution is used in evaluating coefficients for a second order accurate system of difference equations for the difference solution w’(P). The provisional solution G(P) is given by
dj(Q, w(Q))[G,‘(P) - w’@)j,l = b’(Q, w(Q)) At
(2.4)
where ~~~=Q+%(Q)A_x ij
=g
c'(Q,w(Q))
and wi(Qj) = w’(Q) + j;,(w’(Q”) - w’(Q))
O
Wi(Qj) = w’(Q) + ‘j (w’(Q) - w’(Q’))
_l
1
The final difference solution U(P) is calculated from the following system of equations:
( f ){a”&,, MC,)) + a”(P, G(P)) [W’(P) - w’(Oj)l = {bjG’,WV + bj(Qj, d$H}
$
(2.5)
S. Ranganath, R.J. Clifton, A second-order
176
accurate difference
method
where
c’(P,ti(P)) ii =!!! Ax L
+ C’(ijl,W(Q.)) 2
I
1
and d(&)
i. = w’(Q) +; {d(Q”)
XT - w’(Q’)} +$ {w’(Q”) - 2w’(Q) + w’(Q’))
.
3. Proof of convergence From (2.3) the exact solution
ui satisfies
{aii(Qj, U*(Qj)) + a’i(P, U(P))} 2
=
[u’(P) -
hj(Qj, u*(Qj)) + bjU’,0)) 2
~*~(a,,1 At + 0(At3 )
(3.1)
where Qj = Q+ XjAx Qi Xi = (l/Ax)
s
cjdt
P
and U* is the newly defined
function
given by
u*‘(Q) = u’(Q) for mesh points Q and uxi(Qj) = u’(Q) + (Aj/2)[ui(Q”)
- u’(Q’)] + (~;/2)[u’(Q”)
for intermediate points Qj 7 Q + hj Ax. If only first order terms m At are retained, a’j(Q, u(Q))[u’(P)
+x.
-
the exact solution
- 2u’(Q) + u’(Q’)]
ui satisfies
ui(aj)l = b’(Q, u(Q)) At + O(At’)
where 6. = Q Ax and “x. = (At/Ax) c.(Q, Our aim is to dake the dilfference of (i. 1) ui = wi - Use. Before doing so we first derive (i) g(P) - wita,, = {a”(Q, w(Q))}-’
u(Q)) . and (2.5) in order to obtain a difference some order of magnitude estimates. b’(Q, w(Q))At
(3.2)
equation
for
(3.3)
S. Ranganath, R.J. Clifton, A second-order accurate difference method
u’(P) - u*I(Gj) = (&Q,
177
u(Q)>}-’ b’(Q, u(Q)) At + 0(At2)
(3.4)
Subtracting (3.4) from (3.3) gives G,‘(P) - u’(P) = u’((ii) + O(Atu’(Q)) + 0(At2)
(3.5)
where, as before, we define ui(Qj) by Ui(Qj) = Wi(Qj) - u*i(Qj) . (ii) From (3.5) we can write &P,
G(P)) = a”(P, u(P)) + O(ui(Qj)) + O(Atd(Q)) + O(At’ )
(3.6)
@(P,&(P)) = bj(P, u(P)) + O(u’(Q/)) + O(Atd(Q)) + O(At’)
(3.7)
u*‘(e)
= Ups
(3.8)
Qj -Q
= (‘t/2)
(iii) + O(Q, - sj, + O(At3).
Now
{c’(P, U(P)) + C’(Qj> I)
- c’(P, ‘(PI> - Ci(6j, ~(Gj)>I
(3.9)
or Qj - $ = O(Atui(Qj)) + 0(At2ui(Q)) + 0(At3) + (At/2) 1O - W’(Oj)) 1.
(3.10)
Also,
Q, -6, = At&Q,
u(Q)) + 0W2 ) - At&Q,
w(Q))
(3.11)
or Qj -Qj = O(Atr+(Q)) + O(At’) .
(3.12)
Using (3.12) in (3.10) we have Q, -Q/ = O(Atu’(Qj)) + O(At’u’(Q)) + 0(At3) .
(3.13)
Equation (3.8) can now be written as ~*‘(a,) = u*I(Qj) + O(Atr+(Qj)> + 0(At2u’(Q)) + 0(At3) . (iv)
ufw3,9 4&N =u"(Qr
u*(Q,))
+O(dt~'t~~))
+O(At'd(Q))
(3.14)
+Otdt@> +0W3)
.
Having obtained the order of magnitude estimates, we re-write equation (2.5) in the following form :
178
S. Ranganath, R.J. Clifton, A second-order
uri(Qj, u*(Qj)) + a”(P, u(P))
method
{w’(P) - wita,),
2
=
accurate difference
“(P, u(P)>+ “(Qj, u*(Qj)) At + O(v’(P)ui(Qj)) + O(u’(P)u’(Qj)) + O(Atd(P)u’(Q))
2
+ O(Atu’(P)b’(Qj)) + 0(At2~‘(P)ui(Q))
+ O(At” u’(P)) + 0(ui2 (6j) + O(At”‘(Gj)“‘(Q))
+ Q(t~‘(t$~~)d(e~)) + O(Atu’(Gj)u’(C&)) + 0(At2 “‘(Gj)u’(Q)) + Q(At2”‘(Gj))
(3.15)
+ O(Atvi(Qj)) + O(Atui(Qj)) + O(At2 u’(Q)) + O(At’ ~it6~>, + O(At3 u’(Q)) + ON3
).
Similarly, eq. (3.1) can be rewritten as follows:
a'j(Qj, u*(Qj)> + a”(P, U(P)) 2
=
b’(Qj, ‘“(Qj))
{U’(P) - U*i(Qj)}
+ “(P, u(P)) At + O(Atu’(ii,)) + 0(At2 u’(Q)) + 0(At3). 2
(3.16)
Subtracting (3.16) from (3.15) gives
a'j(Qj, u*(Qj)) + a’j(P, u(P)> {u’(P) - ui(ej)} = O(Atu’(oj)) + O(Atu’(Gj))
2
+ O(vi(P)vi(Qj)) + O(v’(P)u’(Qj)) + 0(ui2 (Qj)) + 0(Atui(P)ui(6j)) + O(Atu’(P)u’(Q)) + O(U’(~~~)U’(~~)) + O(Atui(ijj)u’(Q)) + 0(Atui(6j)ui(6j)) + 0(At2u’(ijj)v’(Q)) + 0(At2u’(dj))
+ 0(At2u’(P>u’(Q)) + O(At2 u’(P))
+ O(At’u’(Q)> + O(At’ vi(iij)> + 0(At3 u’(Q)>+ O(At3 ) .
(3.17)
It is convenient to introduce a new measure of the error, Vi, defined by Vi(Q) = a”(Q, u*(Q)) u’(Q) at mesh points and vi(Qj) = vi(Q) + 2 at intermediate
points.
(V’(Q”) - Vj(Q’)) + 3
(Vj(Q”) - 2Vj(Q) + Vj(Q’))
(3.18)
S. Ranganath, R.J. Clifton, A second-order accurate difference method
In terms of Vi, the following order of magnitude estimates are valid: ; a@(P, u(P)) v’(P) = $ F(P)
J a”(P, U(P)) u@, = y’ V@)
+
O(AWj(Q)) + O(A#(Q’))
+ O(ArVj(@‘))
3 a’i(Qi, u*(Qj))uz(P) = + V’(P) + O(AbVi(P)) + a”i(Qj, ~*(Qj))~~(~j) = + F’j(Qj) + Q(A~~~(Q)) + Q(A~~~(Q’)) + O(A~~~(~“)).
(3.19)
Then eq. (3.17) becomes V(P) - V($)
= O(A~~(~~)) + O(Atvf(P)) + O(AtV?Q)) + O(AtYi(Q’H
+ O(Aa Vi(Q”)) + O(~’ (~~>>+ 0( Vi(P) vita,,> f Q(A~~(2li)) + O( y’(P) YiCa,>) + O( ‘(
+ O(A”(P)‘t~jt)
+ 0(At2 vi(P) vi(Q)) + O(At’vi(a,)
+ O(At’(P)In’(Q)) + O(A~~(#j)
‘t(i))
vi(Q)) + O(Ar’ vi($,)
+ 0(A.t2Vi(P)) + 0(At2 Vj(Q)) + 0(At2 Vj(#$) + O(At’) + O(Aht3Vi(Q)) .
(3.20)
The right-hand side of eq. (3.20) can be shortened by-eliminating redundant terms, such as O(At Vi(Q.) Vj(Qj)) which can be included in 0( V(Qj) Vj(Qj)), etc. In this way, eq. (3.20) can be reduce d to Tri(p) = vj(Qj) + O(A~~(~j))
+ O(ArYi(P)) + O(A~~~(Q)) + O(Afvj(Q’))
+ O(ArV’(Q”)) + 0(Vi2 (b,>> f Q(F”(P)“(Qj))
+ O(“(P)“(Qj))
+ o( V(Qj) V(Qj)) + O(AtVj(P)Vj(Q)> + O(LW~(~~V~(QN + O(A~~j(#j)) + 0(At3 1.
(3.21)
As a norm for Vi and IV we take the maximum absolute value of Vi at any point along the line t = vAt. That is, we take as our norm, My, defined by
180
S. Ranganath,
R.J. Clifton, A second-order
accurate difference
method
From (3.22) and (3.21) we obtain A,!+1 < MV+ a {AtMV+
At&+ 1+M;+MyMV+,+AtMpMV+,+AtM;+At3}
(3.23)
where (Yis a positive constant. The terms AtMVMV+Iand AtME can be omitted since for large (Ythey are dominated by the terms M,,MV+l and MS. Thus the inequality to be satisfied reduces to MV+l < MV+ Q:{AtM” + AtMV+, + M;5+ MJV~+~ + At3>.
(3.24)
We now show that for fixed t = vAt the bound on the error MY is O(At*). For this, consider the function HV which satisfies the equation H v+l = HV + a {AtH” + AtH”+, + q
+H,,HV+l + At? .
(3.25)
If we set Ho = To, it is clear that Hy is an upper bound on MY. From (3.25) it is seen that HV+l > H, since (Yis positive. Thus we can replace q in (3.25) by the term HVHV+Ito obtain the following inequality Hv+, < H, + a { AtH, + AtHV+l + 2H,H,,+l+ At3 1 .
(3.26)
The initial value for M, , and consequently H,, , can be written as Ho = M. = pAt3, where fi is a positive constant, since interpolation of the initial data at intermediate points is second order accurate. For all v, and sufficiently small At, H,, is less than the function G
(v+P14~At3
=
”
I-
(v+P/a) aAt
which satisfies the difference equation G“+, = Gv+a {AtG”+ AtGV+l + k
GVG,,+I + At”)
(3.27)
and the initial condition G - flat3 0 l-pat
’
Thus,
M/C HV< GV
(3.28)
and at any time t = vAt G,, =
atAt* + pAt 3 = O(At*) 1 -cd-/3At
.
(3.29)
S. Ranganath, R.J. CZifton,A second-order accurate difference method
181
From (3.28) and (3.29) we conclude that (3.30)
Mv = 0(At2> .
This shows that the difference method is convergent and that for fixed t the error is 0(At2) as At+ 0.
4. Comparison of various second order accurate schemes In this section we will describe two second order accurate difference schemes, the Lax-Wendroff method and the generalized Heun’s method [4], and then compare them with the second order CIR method for a problem in flexural wave propagation. For convenience of presentation we will consider a linear system of first order partial differential equations with constant coefficients. The difference method can, however, be used for solving quasilinear equations. The general linear system of first order partial differential equations with constant coefficients can be written in the form j= 1,2,---,n where Aii, Bii and di are constant coefficients. Assuming that the matrix [B”] is non-singular, n U; = C DiiU: + Eiiui j= l,---,y1 i=l
(4.1)
(4.2)
where [D’j] and [J!?‘] are constant matrices. The Lax-Wendroff scheme involves the expansion of u in a Taylor’s series in t and retaining terms up to second order in At At2 u(x, t + At) = u(x, t) + u&x, t)At + u,,(x, t)z .
(4.3)
From (4.3) j=l,---,It.
(4.4)
Writing the matrix of coefficients [Dii] and [Eii I as D and E, (4.4) can be expressed in vector notation as follows U tt
=D{Duxx+Eux}
+E{Dux+Eu)
.
(4.5)
The Lax-Wendroff scheme is obtained by replacing partial derivatives in (4.5) by the correspond-
182
S. Ranganath,
R.J. Clifton, A second-order
accurate difference
method
ing centered difference expressions
(4.6) where the superscript y1and subscript j indicate the value of u at a point P(jAx, nAt). For quasi-linear systems of equations, the two-step Lax-Wendroff scheme [ 3 ] is more convenient. First, provisional values are calculated at the centers of the rectangular meshes of the net in the x-t plane:
(T+l-“I) +
n
in
n-+1/2 ‘j+
Ax
(4.7)
112
where q+1,2
= $;+,
‘ui”)
and D”j+l,2 = D((i++)Ax,
nAt, $+1/Z )
Enj+l,z = E((j+$)Ax,
nAt, uY+~,~).
The final values are obtained from the equation u:+’ = ui” + At {Di"+"'(u;~;;;
- u”+~/~) + E;+112u;+1~2}
(4.8)
J--1/2
I
The generalized Heun’s method resembles the two-step Lax-Wendroff scheme closely. If we write (4.2) as ut=Dux+Eu=f(x,t,u,uX)
(4.9)
then Heun’s method has the form (see ref. [4] ) u(x, t+At) = u(x, t) + At
x, t +$,
u +$
f, uX + $
fX 11
(4.10)
1
(4.11)
The difference scheme would then be
“+’ zci
= zci” + At
n+ l/2 (uj+l/Z
Ax
f4”+1/2 j-112 >
+ E”+l12 J
u”+l/2 J
S. Rangenuth, R.J. Clifton, A second-order accurate difference method
183
where
The convergence of this method has been proved in 141 for the semi-linear case, for mesh ratios satisfying condition (4.14). For linear systems, the Heun’s method and the Lax-Wendroff scheme yield the same difference equations. In order to investigate the stability of the various difference schemes, we consider the solution of an initial value problem with periodic initial data. In this case U” can be represented by a Fourier series in x. Substitut~g the Fourier series for u” and zP1 in the systems of difference equations for the various schemes (e.g. (4.6) for the Lax-Wendroff scheme) we find that in evaluating the solution at mesh points on x = jAx, the typical term is $k@@ax where 0 = + 1 or 0. Cancelling the common part eijkAx and denoting the Fourier coefficients by y(k) we obtain (see e.g. 131 p. 67)
yn+l (k) = G(At, k,
k)y” (k)
(4.12)
where the matrix G(At, AX, k) is the amplification matrix for the associated difference method. If Ri(i= 1, ---, n) are the eigenvalues of G, then the spectral radius R is the maximum of IR,{ (i = 1, ---, n). For stability, the Von Neumann necessary condition
RG
1
+ O(A0
(4.13)
must be satisfied. For the Lax-Wendroff scheme, the necessary condition (4.13) reduces to IXI &
(4.14)
where lhl = Max]+] (i= 1, ---, n) and X, (i= 1, ---, n) are the eigenvalues of the matrix D in (4.2). We now consider the application of the various difference schemes to the problem of flexural waves in an elastic beam. The equation governing the system (see ref. [8]) can be written as Aut+Bux+Cu=o
where A, B and C are constant matrices given by
(4.15)
184
S. Ranganatk, R.J. Clifton, A second-order
accurate difference
method
and
in which K = k’ A b G where k’ is a dimensionless shear coefficient A, = cross sectional area of the beam G = shear modulus = time t x = distance along beam M = bending moment Q = shear force = transverse velocity V 0 = angular velocity F = mass density of the mate~al of the beam I = moment of inertia of the cross section of the beam P = mass per unit length of the beam E = modulus of elasticity We choose the numerical constants in the matrix A to correspond to a steel beam. The ratio of the mesh sizes ~/A~ is taken to be equal to the longitudin~ wave velocity ci = m, which is the largest value of this ratio for which the necessary condition for stability (71 is satisfied. This choice of the mesh ratio also satisfies the sufficient conditions for all the schemes considered here (see p. 184 ref. 14 I). Thus all the schemes considered are convergent and stable as AX becomes small. However, for finite Ax the various schemes may or may not exhibit computational instability, Also, time required for computing the solution for the same range oft and to the same accuracy is in general different for the different schemes. In what follows we shall compare the schemes on the basis of their stability for finite Ax and on the basis of the computing time required to obtain a given accuracy. The difference equations for the Lax-Wendroff method and the generalized Heun’s method are identical and are given below.
185
S. Ranganath, R.J. Clifton, A second-order accumte difference method
M”+’ I
-M;
n+l n = ui -vi
=
JW;+~ - w;_l )
+
2c 1
(q+l -Q;_,)
+-
2Wl
--
c2 2
()
(u;+1-2u~+u;_,)
(“i”+r-“i”_r)
2
c1
-4
cp 2 ( ch 1 1
(4.16)
K Qi” $-KayAt
PI
where
c2 =m.
The difference equations for the higher order CIR method are similar to those given in (4.16) but a few additional terms are introduced. We show here only the additional terms while the parenthesis is understood to include all the terms on the right hand side of the corresponding equation in (4.16) cd.n+l-q
=
I
[
( ) + ((:r +1)G
(a” /+1
M
i”+‘-“;
= 0
Qj’+l-Qj
= [
( ) -y
-2”;
+ q-1
(l +(32)(w;tl-i?w;+q=l)
KAt2
-s,lte;,l-2Q;+Q;_l)
(4.17)
Comparison of eqs. (4.16) and (4.17) shows that the two schemes differ at most by terms of order A.x3. This is to be expected as both schemes considered are second order accurate.
186
-.
S. Ranganath, R.J. Clifton, A second-order accurate difference
0.1
1I
---
-
I
5-
4
Second Method 1
f
Generalized Heun’s Method and Lax - Wendroff Method
,
--2s 3
I
sx 6
Order
*
method
04 J
CIR
t
L
a
,
f
+-
f36
---
Lax - Wendroff Method
-
Second Order CIR Method
I
I
t
--2A
5%
n
c
01
Fig. 3. Comparison of spectral radii for different values of At/Ax keeping Ax fixed.
Fig. 2. Comparison of spectral radii for the case of flexural waves in beams.
In order to compare the two difference schemes the spectral radius of the amplification matrix was plotted for different values of the mesh size, keeping Ax/At = c1 for both difference schemes. Figure 2 shows a plot of the spectral radius for values of (Y= kAx varying from 0 to r (only the range O-n need be considered because of symmetry about the line 01= n). For Heun’s method and the Lax-Wendroff scheme, the spectral radius increases with increasing mesh size and the maximum value over all Q, is always greater than unity. While the difference scheme is stable in the mathematical sense (i.e. R + 1 as Ax -+ 0), from the computational point of view, the scheme is unstable except for very small mesh sizes. Thus, computational stability can be achieved only by using a small mesh size, which in turn makes the computing time large. For the second order CIR method the spectral radius is always strictly less than unity for all mesh sizes considered. This means that larger mesh sizes can be used with substantial savings in computing time while still maintaining stability, Increasing the mesh size does not lead to numerical instability, but has the effect of decreasing the accuracy, Figure 3 shows the spectral radius of the Lax-Wendroff method for different mesh ratios At/Ax keeping Ax fixed. Although there is a reduction in the spectral radius as a result of the smaller time step, the spectral radius remains greater than unity. It is easy to show that for quite general systems the spectral radius for the Lax-Wendroff method is greater than unity for all non zero At. To show this, consider the amplification matrix G co~espond~g to eq. (4.2): G=
Z+E&+E2
where Z is the unit matrix.
$f
At
ED+(DE+ED)rh
At2
aP( +b2
cosa!-1)P
(4.18)
187
S. Ranganath, R.J. Clifton,A second-order accurate difference method
For E = 0, G is a normal matrix and the Von Neumann necessary condition (4.14) is also a sufficient condition for stability. For E # 0, we can obtain a lower bound on the spectral radius R, by taking (Y= 0 in (4.18). ThUs’,
R>
Max k=l,---,n
1 +p,At+Er,
Zti
2
(4.19)
where p, (k= 1, ---, n) are the eigenvalues of E. For the particular case when E is skew-sym& metric (for symmetric D, this corresponds to non-dissipative systems) its eigenvalues are purely imaginary. Writing p = Max ICC,1, (4.19) becomes k= l,m-,n
(4.20) Thus it is clear that R is always greater than unity and a lower bound on R increases with At. The advantages of the second order CIR method over the other methods considered here were confirmed by a numerical example in which eqs. (4.15) were solved for the case of normal impact of a beam by a semi-infinite elastic rod [8]. With a mesh size A.x = 0.01” corresponding to At = 0.05 microsecond the solution from the Lax-Wendroff scheme was stable only up to 30 /.Is. The computing time for getting the solution up to the first 30 ps was about nine minutes on an IBM 360 model 67 computer. However, for the second order accurate CIR method, the solution was stable for the full range of 50 ~.ts.For a mesh-size five times as large, Ax = O.OS”, the difference scheme was stable throughout the computation and the computing time required was less than one minute. The accuracy of computation was checked by an energy balance (see ref. [ 91) and throughout the computation, the error in the energy balance was less than 1%. Comparisons of numerical solutions obtained by the various methods have also been made for an example in which eqs. (4.15) are quasi-linear. This situation arises in the case of plastic wave propagation for which the material behavior is nonlinear [9]. The problem was solved for plastic wave propagation in an aluminum beam using both the Lax-Wendroff scheme and the second order CIR method. Here again the superiority of the second order CIR method over the LaxWendroff method was evident. For a mesh size of Ax = 0.01” the Lax-Wendroff scheme was stable only up to 10 ps, while for a mesh size five times as large, the extended CIR method was stable for the full range of 100 ps. Thus, we can conclude that the second order accurate CIR method is a useful method having some advantages over other high order accurate difference schemes.
References [l] R. Courant, E. Isaacson and M. Rees, “On the solution of nonlinear hyperbolic differential equations by finite differences,” Comm. Pure Appl. Math. 5.(1952) 243-255. [ 21 P.D. Lax and P.D. Wendroff, “Difference Schemes with higher order accuracy for solving hyperbolic equations,” Comm. Pure Appl. Math. 17 (1964) 381-398. [ 31 R.D. Richtmeyer and K.W. Morton, “Difference methods for initial value problems,” second edition (Interscience publishers, (1967)pp. 302-303.
188
S. Ranganath. R.J. Clifton, A second-order accurate difference method
[4] B. Wendroff, “Theoretical Numerical Analysis” (Academic Press, 1966) pp. 183-193. [5] H.J. Stetter, “On the convergence of characteristic finite-difference methods of high accuracy for quasi-linear hyperbolic equations,” Numer. Math., Vol. 3 (1961) 321. [6] R. Ansorge, “Die Adams-Verfahren als Charakteristikenverfahren hoherer Ordnung der Losung von hyperbolischen Systemen halblinellrer Differentialgleichungen,” Numer. Math., Vol. 5 (1963) 443. [ 71 R. Courant, K.O. Friedrichs and H. Lewy, “Ueber die partlellen Differenzengleichungen der mathematischen Physik,” Math. Ann. 100 (1928) pp. 32-74. [$I S. Ranganath, “Normal impact of an infinite elastic beam by a semi-infinite elastic rod,” Jour. App. Mech., Vol. 38, No. 2 (1971) 455-460. [9] S. Ranganath and R.J. Clifton, “Normal impact of an infinite elastic-plastic beam by a semi-infinite elastic rod,” Int’l. J. Solids and Structs. 8 (1972) 41-67.
A SECOND-ORDER ACCURATE DIFFERENCE METHOD FOR SYSTEMS OF HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS 1 S. RANGANATH 2 and R.J. CLIFTON 3
Received 16 November 1971 A second order accurate difference method is presented for systems of first order hyperbolic differential equations. The method is analogous to the Courant, Isaacson, Rees (CIR) method, except that the error introduced in one time step is O(Ats) instead of O(At’) as is the case for the CIR method. Convergence of the proposed method is established. The cumulative error at a fixed time is shown to be O(APt*).The proposed method is compared with several other second order accurate methods by considering in detail the special case of the system of equations governing flexural wave propagation in elastic beams. These comparisons indicate that the proposed method has substantial advantages over the other methods considered in that the method is computationally stable for larger mesh sizes. As a result, less computing time is required to obtain solutions with a given accuracy.
1. Introduction In 1953, Courant, Isaacson and Rees [ 11 presented a difference method for the solution of quasi-linear hyperbolic partial differential equations of first order. The method makes use of a fixed mesh and integration along characteristics to determine the solution at an arbitrary mesh point from data at neighboring mesh points at the previous time step. The error introduced in going from the exact solution at time t to the difference solution at time t + At was shown to be O(At2) which indicates that the method corresponds to what is generally called a first order accurate method. In order to save computing time for a given accuracy it is frequently advantageous to make use of higher order accurate difference schemes for which the same accuracy can be obtained with a larger mesh size. A well-known example of such higher order accurate schemes is the second order accurate difference method developed by Lax and Wendroff [2] . The LaxWendroff method is based on expanding the solution vector in a Taylor’s series and retaining terms up to second order. Richtmeyer [ 3 I has presented a two step Lax-Wendroff scheme which is also second order accurate but slightly simpler to use. When the coefficients are constant it reduces to the original Lax-Wendroff scheme [2]. The generalized Heun’s method [4] is also a two-step difference scheme which is second order accurate. Stetter [51 has constructed higher order accurate methods based on characteristic finite difference schemes for systems of equations
‘. ‘IhIs research was supported by the Air Force Materials Laboratory, Wright-Patterson Airforce Base under Contract No. F 33615-70-C-1101 with Brown University. 2. Research Associate, Brown University. Now at Advanced Technology Center, Inc., Dallas, Texas. 3. Professor of Engineering, Brown University.
174
S. Rangmath, R.J. Oifton, A second-order accurate difference method
which have only two independent characteristics. Using the Adam’s procedure for ordinary differential equations Ansorge 161 has presented higher order accurate methods for semilinear equations. We present here a second order accurate method which is a natural extension of the first order method proposed by Courant, Isaacson and Rees [ I] for general quasilinear hyperbolic systems in two independent variables. This method (referred to hereafter as the second order accurate CIR method) is based on integration along characteristics, preserving second order accuracy throughout. The purpose of this paper is to describe the second order accurate CIR method and establish the convergence of the difference scheme. This method is then compared with other second order accurate methods by considering in detail the special case of the propagation of flexural waves in elastic beam. It will be shown that, at least for this case, the new method has definite advantages over other second order accurate methods.
2. Description of the difference method In this section we develop a difference method for quasi-linear systems of hyperbolic differential equations in two independent variables. Convergence of the method will be established in the next section. Quasi-linear systems of first order partial differential equations in two independent variables have the form * *$ {A%;
+ B’h) + di} = 0
j=l,---,y1
(2.1)
where A’j, Bii and di are in general, functions of the (n+2) variables (x, t, u’). If the system is hyperbolic, n distinct linear combinations of equations (2.1) may be formed to produce an equivalent system of equations in the normal form n
C
aji(uf+
,iu;
) =:
bi
j= 1,2, ---,n
i=l
(2.2)
where aif, ci and bj are functions of x, t and ui. (We assume here that the direction t = constant is not a characteristic.) In the j-th equation of the normal form (2.2) the terms (uf+cj,L) are simply the derivatives, du’fdt, of ui along the curve x = x(t) for which dx/dt = ci (the j-th characteristic direction). Thus, (2.2) can now be written in the form ze a’idu’ = bidt
(2.3)
where du’ represents the change in ui along the j-th characteristic for a change dt in t. We assume that the coefficients in (2.1) and the initial data are sufficiently smooth for a unique solution, u’(t,x) E C3 to exist in some finite region t 2 0. In particular we shall assume det aii # 0 * Subscripts x and 2 denote differentiation with respect to n and t respectively.
S. Ranganath, R.J. CIifton, A second-order accurate difference method
175
Fig. 1. Grid and Characteristics.
and &, cj and bj are Lipschitz continuous functions of all their arguments. Our objective is to develop a finite difference method for computing functions w’(t,x) which can be made arbitrarily close to u’(t,x) in some finite region t > 0 by taking the mesh size sufficiently small. We select a mesh width Ax and consider a rectangular network of points in the x-t plane. From the known difference solution at mesh points on the line IV, we obtain the solution at mesh points on the line IV+i by integration of (2.3) along characteristics. (Fig. 1). The time step between the lines IV and IV+i is chosen such that, if P is a net point on IV+i and Q its nearest neighboring point on IV, then the line segments drawn backwards through P with slopes dx/dt = cj intersect IV at points which lie between Q and its nearest neighboring points Q’ and Q” on 1”. In this way the Courant, Friedrichs, Lewy [ 71 necessary condition for stability is automatically satisfied. The difference solution at a mesh point P is computed by means of a two-step procedure. First, a provisional solution G’(P) is obtained using the CIR first order accurate method. Then, this solution is used in evaluating coefficients for a second order accurate system of difference equations for the difference solution w’(P). The provisional solution G(P) is given by
dj(Q, w(Q))[G,‘(P) - w’@)j,l = b’(Q, w(Q)) At
(2.4)
where ~~~=Q+%(Q)A_x ij
=g
c'(Q,w(Q))
and wi(Qj) = w’(Q) + j;,(w’(Q”) - w’(Q))
O
Wi(Qj) = w’(Q) + ‘j (w’(Q) - w’(Q’))
_l
1
The final difference solution U(P) is calculated from the following system of equations:
( f ){a”&,, MC,)) + a”(P, G(P)) [W’(P) - w’(Oj)l = {bjG’,WV + bj(Qj, d$H}
$
(2.5)
S. Ranganath, R.J. Clifton, A second-order
176
accurate difference
method
where
c’(P,ti(P)) ii =!!! Ax L
+ C’(ijl,W(Q.)) 2
I
1
and d(&)
i. = w’(Q) +; {d(Q”)
XT - w’(Q’)} +$ {w’(Q”) - 2w’(Q) + w’(Q’))
.
3. Proof of convergence From (2.3) the exact solution
ui satisfies
{aii(Qj, U*(Qj)) + a’i(P, U(P))} 2
=
[u’(P) -
hj(Qj, u*(Qj)) + bjU’,0)) 2
~*~(a,,1 At + 0(At3 )
(3.1)
where Qj = Q+ XjAx Qi Xi = (l/Ax)
s
cjdt
P
and U* is the newly defined
function
given by
u*‘(Q) = u’(Q) for mesh points Q and uxi(Qj) = u’(Q) + (Aj/2)[ui(Q”)
- u’(Q’)] + (~;/2)[u’(Q”)
for intermediate points Qj 7 Q + hj Ax. If only first order terms m At are retained, a’j(Q, u(Q))[u’(P)
+x.
-
the exact solution
- 2u’(Q) + u’(Q’)]
ui satisfies
ui(aj)l = b’(Q, u(Q)) At + O(At’)
where 6. = Q Ax and “x. = (At/Ax) c.(Q, Our aim is to dake the dilfference of (i. 1) ui = wi - Use. Before doing so we first derive (i) g(P) - wita,, = {a”(Q, w(Q))}-’
u(Q)) . and (2.5) in order to obtain a difference some order of magnitude estimates. b’(Q, w(Q))At
(3.2)
equation
for
(3.3)
S. Ranganath, R.J. Clifton, A second-order accurate difference method
u’(P) - u*I(Gj) = (&Q,
177
u(Q)>}-’ b’(Q, u(Q)) At + 0(At2)
(3.4)
Subtracting (3.4) from (3.3) gives G,‘(P) - u’(P) = u’((ii) + O(Atu’(Q)) + 0(At2)
(3.5)
where, as before, we define ui(Qj) by Ui(Qj) = Wi(Qj) - u*i(Qj) . (ii) From (3.5) we can write &P,
G(P)) = a”(P, u(P)) + O(ui(Qj)) + O(Atd(Q)) + O(At’ )
(3.6)
@(P,&(P)) = bj(P, u(P)) + O(u’(Q/)) + O(Atd(Q)) + O(At’)
(3.7)
u*‘(e)
= Ups
(3.8)
Qj -Q
= (‘t/2)
(iii) + O(Q, - sj, + O(At3).
Now
{c’(P, U(P)) + C’(Qj> I)
- c’(P, ‘(PI> - Ci(6j, ~(Gj)>I
(3.9)
or Qj - $ = O(Atui(Qj)) + 0(At2ui(Q)) + 0(At3) + (At/2) 1O
(3.10)
Also,
Q, -6, = At&Q,
u(Q)) + 0W2 ) - At&Q,
w(Q))
(3.11)
or Qj -Qj = O(Atr+(Q)) + O(At’) .
(3.12)
Using (3.12) in (3.10) we have Q, -Q/ = O(Atu’(Qj)) + O(At’u’(Q)) + 0(At3) .
(3.13)
Equation (3.8) can now be written as ~*‘(a,) = u*I(Qj) + O(Atr+(Qj)> + 0(At2u’(Q)) + 0(At3) . (iv)
ufw3,9 4&N =u"(Qr
u*(Q,))
+O(dt~'t~~))
+O(At'd(Q))
(3.14)
+Otdt@> +0W3)
.
Having obtained the order of magnitude estimates, we re-write equation (2.5) in the following form :
178
S. Ranganath, R.J. Clifton, A second-order
uri(Qj, u*(Qj)) + a”(P, u(P))
method
{w’(P) - wita,),
2
=
accurate difference
“(P, u(P)>+ “(Qj, u*(Qj)) At + O(v’(P)ui(Qj)) + O(u’(P)u’(Qj)) + O(Atd(P)u’(Q))
2
+ O(Atu’(P)b’(Qj)) + 0(At2~‘(P)ui(Q))
+ O(At” u’(P)) + 0(ui2 (6j) + O(At”‘(Gj)“‘(Q))
+ Q(t~‘(t$~~)d(e~)) + O(Atu’(Gj)u’(C&)) + 0(At2 “‘(Gj)u’(Q)) + Q(At2”‘(Gj))
(3.15)
+ O(Atvi(Qj)) + O(Atui(Qj)) + O(At2 u’(Q)) + O(At’ ~it6~>, + O(At3 u’(Q)) + ON3
).
Similarly, eq. (3.1) can be rewritten as follows:
a'j(Qj, u*(Qj)> + a”(P, U(P)) 2
=
b’(Qj, ‘“(Qj))
{U’(P) - U*i(Qj)}
+ “(P, u(P)) At + O(Atu’(ii,)) + 0(At2 u’(Q)) + 0(At3). 2
(3.16)
Subtracting (3.16) from (3.15) gives
a'j(Qj, u*(Qj)) + a’j(P, u(P)> {u’(P) - ui(ej)} = O(Atu’(oj)) + O(Atu’(Gj))
2
+ O(vi(P)vi(Qj)) + O(v’(P)u’(Qj)) + 0(ui2 (Qj)) + 0(Atui(P)ui(6j)) + O(Atu’(P)u’(Q)) + O(U’(~~~)U’(~~)) + O(Atui(ijj)u’(Q)) + 0(Atui(6j)ui(6j)) + 0(At2u’(ijj)v’(Q)) + 0(At2u’(dj))
+ 0(At2u’(P>u’(Q)) + O(At2 u’(P))
+ O(At’u’(Q)> + O(At’ vi(iij)> + 0(At3 u’(Q)>+ O(At3 ) .
(3.17)
It is convenient to introduce a new measure of the error, Vi, defined by Vi(Q) = a”(Q, u*(Q)) u’(Q) at mesh points and vi(Qj) = vi(Q) + 2 at intermediate
points.
(V’(Q”) - Vj(Q’)) + 3
(Vj(Q”) - 2Vj(Q) + Vj(Q’))
(3.18)
S. Ranganath, R.J. Clifton, A second-order accurate difference method
In terms of Vi, the following order of magnitude estimates are valid: ; a@(P, u(P)) v’(P) = $ F(P)
J a”(P, U(P)) u@, = y’ V@)
+
O(AWj(Q)) + O(A#(Q’))
+ O(ArVj(@‘))
3 a’i(Qi, u*(Qj))uz(P) = + V’(P) + O(AbVi(P)) + a”i(Qj, ~*(Qj))~~(~j) = + F’j(Qj) + Q(A~~~(Q)) + Q(A~~~(Q’)) + O(A~~~(~“)).
(3.19)
Then eq. (3.17) becomes V(P) - V($)
= O(A~~(~~)) + O(Atvf(P)) + O(AtV?Q)) + O(AtYi(Q’H
+ O(Aa Vi(Q”)) + O(~’ (~~>>+ 0( Vi(P) vita,,> f Q(A~~(2li)) + O( y’(P) YiCa,>) + O( ‘(
+ O(A”(P)‘t~jt)
+ 0(At2 vi(P) vi(Q)) + O(At’vi(a,)
+ O(At’(P)In’(Q)) + O(A~~(#j)
‘t(i))
vi(Q)) + O(Ar’ vi($,)
+ 0(A.t2Vi(P)) + 0(At2 Vj(Q)) + 0(At2 Vj(#$) + O(At’) + O(Aht3Vi(Q)) .
(3.20)
The right-hand side of eq. (3.20) can be shortened by-eliminating redundant terms, such as O(At Vi(Q.) Vj(Qj)) which can be included in 0( V(Qj) Vj(Qj)), etc. In this way, eq. (3.20) can be reduce d to Tri(p) = vj(Qj) + O(A~~(~j))
+ O(ArYi(P)) + O(A~~~(Q)) + O(Afvj(Q’))
+ O(ArV’(Q”)) + 0(Vi2 (b,>> f Q(F”(P)“(Qj))
+ O(“(P)“(Qj))
+ o( V(Qj) V(Qj)) + O(AtVj(P)Vj(Q)> + O(LW~(~~V~(QN + O(A~~j(#j)) + 0(At3 1.
(3.21)
As a norm for Vi and IV we take the maximum absolute value of Vi at any point along the line t = vAt. That is, we take as our norm, My, defined by
180
S. Ranganath,
R.J. Clifton, A second-order
accurate difference
method
From (3.22) and (3.21) we obtain A,!+1 < MV+ a {AtMV+
At&+ 1+M;+MyMV+,+AtMpMV+,+AtM;+At3}
(3.23)
where (Yis a positive constant. The terms AtMVMV+Iand AtME can be omitted since for large (Ythey are dominated by the terms M,,MV+l and MS. Thus the inequality to be satisfied reduces to MV+l < MV+ Q:{AtM” + AtMV+, + M;5+ MJV~+~ + At3>.
(3.24)
We now show that for fixed t = vAt the bound on the error MY is O(At*). For this, consider the function HV which satisfies the equation H v+l = HV + a {AtH” + AtH”+, + q
+H,,HV+l + At? .
(3.25)
If we set Ho = To, it is clear that Hy is an upper bound on MY. From (3.25) it is seen that HV+l > H, since (Yis positive. Thus we can replace q in (3.25) by the term HVHV+Ito obtain the following inequality Hv+, < H, + a { AtH, + AtHV+l + 2H,H,,+l+ At3 1 .
(3.26)
The initial value for M, , and consequently H,, , can be written as Ho = M. = pAt3, where fi is a positive constant, since interpolation of the initial data at intermediate points is second order accurate. For all v, and sufficiently small At, H,, is less than the function G
(v+P14~At3
=
”
I-
(v+P/a) aAt
which satisfies the difference equation G“+, = Gv+a {AtG”+ AtGV+l + k
GVG,,+I + At”)
(3.27)
and the initial condition G - flat3 0 l-pat
’
Thus,
M/C HV< GV
(3.28)
and at any time t = vAt G,, =
atAt* + pAt 3 = O(At*) 1 -cd-/3At
.
(3.29)
S. Ranganath, R.J. CZifton,A second-order accurate difference method
181
From (3.28) and (3.29) we conclude that (3.30)
Mv = 0(At2> .
This shows that the difference method is convergent and that for fixed t the error is 0(At2) as At+ 0.
4. Comparison of various second order accurate schemes In this section we will describe two second order accurate difference schemes, the Lax-Wendroff method and the generalized Heun’s method [4], and then compare them with the second order CIR method for a problem in flexural wave propagation. For convenience of presentation we will consider a linear system of first order partial differential equations with constant coefficients. The difference method can, however, be used for solving quasilinear equations. The general linear system of first order partial differential equations with constant coefficients can be written in the form j= 1,2,---,n where Aii, Bii and di are constant coefficients. Assuming that the matrix [B”] is non-singular, n U; = C DiiU: + Eiiui j= l,---,y1 i=l
(4.1)
(4.2)
where [D’j] and [J!?‘] are constant matrices. The Lax-Wendroff scheme involves the expansion of u in a Taylor’s series in t and retaining terms up to second order in At At2 u(x, t + At) = u(x, t) + u&x, t)At + u,,(x, t)z .
(4.3)
From (4.3) j=l,---,It.
(4.4)
Writing the matrix of coefficients [Dii] and [Eii I as D and E, (4.4) can be expressed in vector notation as follows U tt
=D{Duxx+Eux}
+E{Dux+Eu)
.
(4.5)
The Lax-Wendroff scheme is obtained by replacing partial derivatives in (4.5) by the correspond-
182
S. Ranganath,
R.J. Clifton, A second-order
accurate difference
method
ing centered difference expressions
(4.6) where the superscript y1and subscript j indicate the value of u at a point P(jAx, nAt). For quasi-linear systems of equations, the two-step Lax-Wendroff scheme [ 3 ] is more convenient. First, provisional values are calculated at the centers of the rectangular meshes of the net in the x-t plane:
(T+l-“I) +
n
in
n-+1/2 ‘j+
Ax
(4.7)
112
where q+1,2
= $;+,
‘ui”)
and D”j+l,2 = D((i++)Ax,
nAt, $+1/Z )
Enj+l,z = E((j+$)Ax,
nAt, uY+~,~).
The final values are obtained from the equation u:+’ = ui” + At {Di"+"'(u;~;;;
- u”+~/~) + E;+112u;+1~2}
(4.8)
J--1/2
I
The generalized Heun’s method resembles the two-step Lax-Wendroff scheme closely. If we write (4.2) as ut=Dux+Eu=f(x,t,u,uX)
(4.9)
then Heun’s method has the form (see ref. [4] ) u(x, t+At) = u(x, t) + At
x, t +$,
u +$
f, uX + $
fX 11
(4.10)
1
(4.11)
The difference scheme would then be
“+’ zci
= zci” + At
n+ l/2 (uj+l/Z
Ax
f4”+1/2 j-112 >
+ E”+l12 J
u”+l/2 J
S. Rangenuth, R.J. Clifton, A second-order accurate difference method
183
where
The convergence of this method has been proved in 141 for the semi-linear case, for mesh ratios satisfying condition (4.14). For linear systems, the Heun’s method and the Lax-Wendroff scheme yield the same difference equations. In order to investigate the stability of the various difference schemes, we consider the solution of an initial value problem with periodic initial data. In this case U” can be represented by a Fourier series in x. Substitut~g the Fourier series for u” and zP1 in the systems of difference equations for the various schemes (e.g. (4.6) for the Lax-Wendroff scheme) we find that in evaluating the solution at mesh points on x = jAx, the typical term is $k@@ax where 0 = + 1 or 0. Cancelling the common part eijkAx and denoting the Fourier coefficients by y(k) we obtain (see e.g. 131 p. 67)
yn+l (k) = G(At, k,
k)y” (k)
(4.12)
where the matrix G(At, AX, k) is the amplification matrix for the associated difference method. If Ri(i= 1, ---, n) are the eigenvalues of G, then the spectral radius R is the maximum of IR,{ (i = 1, ---, n). For stability, the Von Neumann necessary condition
RG
1
+ O(A0
(4.13)
must be satisfied. For the Lax-Wendroff scheme, the necessary condition (4.13) reduces to IXI &
(4.14)
where lhl = Max]+] (i= 1, ---, n) and X, (i= 1, ---, n) are the eigenvalues of the matrix D in (4.2). We now consider the application of the various difference schemes to the problem of flexural waves in an elastic beam. The equation governing the system (see ref. [8]) can be written as Aut+Bux+Cu=o
where A, B and C are constant matrices given by
(4.15)
184
S. Ranganatk, R.J. Clifton, A second-order
accurate difference
method
and
in which K = k’ A b G where k’ is a dimensionless shear coefficient A, = cross sectional area of the beam G = shear modulus = time t x = distance along beam M = bending moment Q = shear force = transverse velocity V 0 = angular velocity F = mass density of the mate~al of the beam I = moment of inertia of the cross section of the beam P = mass per unit length of the beam E = modulus of elasticity We choose the numerical constants in the matrix A to correspond to a steel beam. The ratio of the mesh sizes ~/A~ is taken to be equal to the longitudin~ wave velocity ci = m, which is the largest value of this ratio for which the necessary condition for stability (71 is satisfied. This choice of the mesh ratio also satisfies the sufficient conditions for all the schemes considered here (see p. 184 ref. 14 I). Thus all the schemes considered are convergent and stable as AX becomes small. However, for finite Ax the various schemes may or may not exhibit computational instability, Also, time required for computing the solution for the same range oft and to the same accuracy is in general different for the different schemes. In what follows we shall compare the schemes on the basis of their stability for finite Ax and on the basis of the computing time required to obtain a given accuracy. The difference equations for the Lax-Wendroff method and the generalized Heun’s method are identical and are given below.
185
S. Ranganath, R.J. Clifton, A second-order accumte difference method
M”+’ I
-M;
n+l n = ui -vi
=
JW;+~ - w;_l )
+
2c 1
(q+l -Q;_,)
+-
2Wl
--
c2 2
()
(u;+1-2u~+u;_,)
(“i”+r-“i”_r)
2
c1
-4
cp 2 ( ch 1 1
(4.16)
K Qi” $-KayAt
PI
where
c2 =m.
The difference equations for the higher order CIR method are similar to those given in (4.16) but a few additional terms are introduced. We show here only the additional terms while the parenthesis is understood to include all the terms on the right hand side of the corresponding equation in (4.16) cd.n+l-q
=
I
[
( ) + ((:r +1)G
(a” /+1
M
i”+‘-“;
= 0
Qj’+l-Qj
= [
( ) -y
-2”;
+ q-1
(l +(32)(w;tl-i?w;+q=l)
KAt2
-s,lte;,l-2Q;+Q;_l)
(4.17)
Comparison of eqs. (4.16) and (4.17) shows that the two schemes differ at most by terms of order A.x3. This is to be expected as both schemes considered are second order accurate.
186
-.
S. Ranganath, R.J. Clifton, A second-order accurate difference
0.1
1I
---
-
I
5-
4
Second Method 1
f
Generalized Heun’s Method and Lax - Wendroff Method
,
--2s 3
I
sx 6
Order
*
method
04 J
CIR
t
L
a
,
f
+-
f36
---
Lax - Wendroff Method
-
Second Order CIR Method
I
I
t
--2A
5%
n
c
01
Fig. 3. Comparison of spectral radii for different values of At/Ax keeping Ax fixed.
Fig. 2. Comparison of spectral radii for the case of flexural waves in beams.
In order to compare the two difference schemes the spectral radius of the amplification matrix was plotted for different values of the mesh size, keeping Ax/At = c1 for both difference schemes. Figure 2 shows a plot of the spectral radius for values of (Y= kAx varying from 0 to r (only the range O-n need be considered because of symmetry about the line 01= n). For Heun’s method and the Lax-Wendroff scheme, the spectral radius increases with increasing mesh size and the maximum value over all Q, is always greater than unity. While the difference scheme is stable in the mathematical sense (i.e. R + 1 as Ax -+ 0), from the computational point of view, the scheme is unstable except for very small mesh sizes. Thus, computational stability can be achieved only by using a small mesh size, which in turn makes the computing time large. For the second order CIR method the spectral radius is always strictly less than unity for all mesh sizes considered. This means that larger mesh sizes can be used with substantial savings in computing time while still maintaining stability, Increasing the mesh size does not lead to numerical instability, but has the effect of decreasing the accuracy, Figure 3 shows the spectral radius of the Lax-Wendroff method for different mesh ratios At/Ax keeping Ax fixed. Although there is a reduction in the spectral radius as a result of the smaller time step, the spectral radius remains greater than unity. It is easy to show that for quite general systems the spectral radius for the Lax-Wendroff method is greater than unity for all non zero At. To show this, consider the amplification matrix G co~espond~g to eq. (4.2): G=
Z+E&+E2
where Z is the unit matrix.
$f
At
ED+(DE+ED)rh
At2
aP( +b2
cosa!-1)P
(4.18)
187
S. Ranganath, R.J. Clifton,A second-order accurate difference method
For E = 0, G is a normal matrix and the Von Neumann necessary condition (4.14) is also a sufficient condition for stability. For E # 0, we can obtain a lower bound on the spectral radius R, by taking (Y= 0 in (4.18). ThUs’,
R>
Max k=l,---,n
1 +p,At+Er,
Zti
2
(4.19)
where p, (k= 1, ---, n) are the eigenvalues of E. For the particular case when E is skew-sym& metric (for symmetric D, this corresponds to non-dissipative systems) its eigenvalues are purely imaginary. Writing p = Max ICC,1, (4.19) becomes k= l,m-,n
(4.20) Thus it is clear that R is always greater than unity and a lower bound on R increases with At. The advantages of the second order CIR method over the other methods considered here were confirmed by a numerical example in which eqs. (4.15) were solved for the case of normal impact of a beam by a semi-infinite elastic rod [8]. With a mesh size A.x = 0.01” corresponding to At = 0.05 microsecond the solution from the Lax-Wendroff scheme was stable only up to 30 /.Is. The computing time for getting the solution up to the first 30 ps was about nine minutes on an IBM 360 model 67 computer. However, for the second order accurate CIR method, the solution was stable for the full range of 50 ~.ts.For a mesh-size five times as large, Ax = O.OS”, the difference scheme was stable throughout the computation and the computing time required was less than one minute. The accuracy of computation was checked by an energy balance (see ref. [ 91) and throughout the computation, the error in the energy balance was less than 1%. Comparisons of numerical solutions obtained by the various methods have also been made for an example in which eqs. (4.15) are quasi-linear. This situation arises in the case of plastic wave propagation for which the material behavior is nonlinear [9]. The problem was solved for plastic wave propagation in an aluminum beam using both the Lax-Wendroff scheme and the second order CIR method. Here again the superiority of the second order CIR method over the LaxWendroff method was evident. For a mesh size of Ax = 0.01” the Lax-Wendroff scheme was stable only up to 10 ps, while for a mesh size five times as large, the extended CIR method was stable for the full range of 100 ps. Thus, we can conclude that the second order accurate CIR method is a useful method having some advantages over other high order accurate difference schemes.
References [l] R. Courant, E. Isaacson and M. Rees, “On the solution of nonlinear hyperbolic differential equations by finite differences,” Comm. Pure Appl. Math. 5.(1952) 243-255. [ 21 P.D. Lax and P.D. Wendroff, “Difference Schemes with higher order accuracy for solving hyperbolic equations,” Comm. Pure Appl. Math. 17 (1964) 381-398. [ 31 R.D. Richtmeyer and K.W. Morton, “Difference methods for initial value problems,” second edition (Interscience publishers, (1967)pp. 302-303.
188
S. Ranganath. R.J. Clifton, A second-order accurate difference method
[4] B. Wendroff, “Theoretical Numerical Analysis” (Academic Press, 1966) pp. 183-193. [5] H.J. Stetter, “On the convergence of characteristic finite-difference methods of high accuracy for quasi-linear hyperbolic equations,” Numer. Math., Vol. 3 (1961) 321. [6] R. Ansorge, “Die Adams-Verfahren als Charakteristikenverfahren hoherer Ordnung der Losung von hyperbolischen Systemen halblinellrer Differentialgleichungen,” Numer. Math., Vol. 5 (1963) 443. [ 71 R. Courant, K.O. Friedrichs and H. Lewy, “Ueber die partlellen Differenzengleichungen der mathematischen Physik,” Math. Ann. 100 (1928) pp. 32-74. [$I S. Ranganath, “Normal impact of an infinite elastic beam by a semi-infinite elastic rod,” Jour. App. Mech., Vol. 38, No. 2 (1971) 455-460. [9] S. Ranganath and R.J. Clifton, “Normal impact of an infinite elastic-plastic beam by a semi-infinite elastic rod,” Int’l. J. Solids and Structs. 8 (1972) 41-67.