- Email: [email protected]
Special methods for problems whose oscillatory solution is damped
161
SPECIAL METHODS FOR PROBLEW WRXiE OSCILLATORYSOLUTION IS DAMPED Beny Neta Naval Postgradua+eSchool Departmentof Mathematics Code 53Nd Monterey, CA 93943 ABSTRACT This paper introducesmethods tailored especially for prcblems whose solution behaves like eXX, where h is complex. The shallow water equations with topographyadmit such solution. This paper complementsthe results of Pratt and others on exponentialfitted methods and those of Gautschi,Neta, van der Wouwen and others on trigonometrically-fitted methods. 1.
Introduction In this paper we consider linear multistep methods h
1, n 2 k-l f bff(xn+l_f,~~+l_~),k _?. !2,=0
(1)
for integratingthe initial value problem (2)
Y'(X) = f(x,y(x)), Y(Xo) = Y. * This linear multistep method is characterizedby the polynomials P(S) = krO aetk-' ,
~(0
=
F beck-” 11=0
.
(3)
The main assumptionof this paper is that it is a priori known that the solution is approximatelyof the form m iXjt 1 cje j=l.
Y(X) - no +
(4)
where hj = Wj + iJlj,and the frequencieswj are in a given interval [wc,wu]. The special case where . = jwo with WC given was considered first by Gautschi [B]. His approach ias the following. Let: 9(Z) = p(e') - za(e'>
(5)
then the local truncationerror of (1) is given by Lambert [ill Tn+k = $(h &)
y(tn) .
(6)
Inserting (4) in (6) yields Tn+k - cP(0)cO+
I
iXjt , ij = jw, ci4(ihXj)e
(7)
j=l The coefficientsbi are chosen in such a way that $(ihjw,) = 0 ,
j =
O,l,...,q ,
@Elsevier Scrence Publmhlng Co.,lnc.,1999 655Avenue of ths Americas, New York, NY 10010
(8)
for the largest value of q possible. q is then called the trigonometric order of the method. Gautschi has chosen ai such that the methods are those of Adams and Stormer type. However, these methods are sensitive to changes in the frequency w0. Neta and Ford [13] developed Nystrom and generalized Milne-Simpson type methods. These methods showed less sensitivity to perturbation in w0 but require the eigenvalues of the Jacobian to be purely imaginary. Neta [14] has developed families of backward differentiation methods that overcome the above-mentioned restrrctfon. Salzer [17] has developed predictor-corrector methods based on trigonometric polynomials. See also Steifel and Bettis [18] and Bettis [3]. Van der Houwen and Sommeijer [lo] have developed an alternative approach. The conditions (8) were replaced by
wx = 0
I
(9) $(ihX(j)) = 0 ,
j = l,i,...,q ,
where the A(j) are appropriately chosen points in the fnterval [wa,wU]. An advantage of this so-called minimax approach over the fitting approach is the increased accuracy in cases where no accurate estimate of wD is available or when the frequency is varying fn time. The other special case considered in the literature is where Xj = i$j. Probably the first article on the subject is due to Brock and Murray [5]. They discuss the use of exponential sums in the integration of a system of first order ordinary differential equations. Dennis [7] also suggested special methods for problems whose solution Is exponential. He suggested a transformation of variables. More recently, Carroll 163 has developed exponentfally fitted one-step methods for the scalar Riccati equation. For the general first order system of equations, Pratt [Hi] suggests methods based on the three parameter exponential function f(x)
= A + BeZX .
(10)
The parameters A, B are given in terms of values of y and f. $everal possibilities for z are given based on results of Brandon [2] and Babcock et al. [I]. Lyche [12] analyzes multistep methods which exactly integrate the set wx {xme n ], where wn is real or imaginary. In this article we developed various methods fitting exponentials and methods obtained via the minimax approach. 2.
Construction of Methods
2.1 Fitting Methods In this subsection we discuss various fitting methods. To this end, we separate $(ihjA) = 0, j = l,Z,...,q, into real and imaginary parts. This yields the following equations relating the coefficients aE,bt, cos
jv(k-ll)-
i R=O
b .ju(k-a) [-JV . cos jv(k-ll) ,c - jv sin jv(k-!Z)]= 0 ,
[ju sin jV(k-a) sin jV(k-ll)- t aeejp(k-k) bllejU(k-J1) !2=0
ll=o
f jv cos jV(k-r?)]= 0 , where X =w+i$,
u = -h$, V = hw, j = l,Z,...,q.
163
For explicitmethods, b0 = 0. For Adams type methods a0 = 1, al = -1, ai = 0 for i = Z,...,k. For Nystrom or generalizedMilne-Simpsonmethods a0 = 1, a2 = -1 and other ai = 0. k = 1 Implicit Adams b. 5
we-’ +
sin V - v cos v
1-I
9
sin V
b2+v2)
bl = ve' - p sin v - v cos v ’
(p2+v2) sin v
1 - cos v =b 1, which agree with sin v Gautschi [S] if the coefficientsare expanded in Taylor series with respect to v).
For $ = 0, the coefficientsbecome b0 =
k=
2 Explicit
Adams bl =
(n sin 2v - v cos 2v)ep + v cos v - p sin V (p2+v') sin v (13)
= (v cos v - u sin V)e2' - we'
b 2
.
(p2*V2) sin V
Nystrom (n sin 2v - v cos 2")s' + Ve-' bl =
I
(v2+V2) sin V (14)
b2 =
A (v
cos
2lJ
v -
(p’+V2)
. sin
V
sin V For $ = 0, the coefficientsbecome bl = 2 7, Neta [14]. k=
b2 = 0 which agree with
2 Implicit
In this case, one obtains a one-parameterfamily of (Adams, generalized Milne-Simpson)methods of trigonometricorder I. The free parameter can be used to increase the algebraic order of the method as in 1131. Backward Differentiation e2’
cos
2v + a e' co6 v + a - b e2'(p cos 2v 1 2 0 - bOe2'(p sin
ezp sin 2v + ale' sin V 1
+ al
+ a2
2~
+
=Z 0
,
V
sin
2V)
v
COS
?V) = 0 3
(15)
=o.
This system can be solved by MACSYMA (ProjectNAG's SYmbolic MAnipulation system written in LISP and used for performing symbolic as well as numerical mathematicalmanipulation [4]) or by REDUCE 191. The solution is
ve*n - Ilsin 2v - v cos 2v al= -eu(v cos V + u sin V) + u sin 2V + V cos 2V ' (16)
a2 = -1 - al , _e*u sin
b0 =
V + e’
sin
Zv -
sin V
-e2qv cos v + D sin L) + eu(u sin 2v + v coo 2v) * For 9 = 0, the coefficientsagree with those given in 1141. k = 3 Explicit Again here, one obtaiqs a one-parameterfamily of methods of trigonometric order 1. In ordc.,to get methods of trigonometricorder 2, one has to constructa 3 step -mpiicitmethod of Adams or generalizedMilne-Simpson q-se. In order to in:rease the trigonometricorder without going to a higher step number, 3ne can constructliusar multistepmethods for which the coefficientsallare also functionsof J1,w. Some examples are given in the next subsection. 2.2 GeneralizedFitting Methods In this section,we constructsome linear multistepmethods of the form f
all(h)yel_ll
= h
!fba(A)fn+l_E
R=O
L=O
WI
k>l,nlk-1. -
3
Since a% are functions of X one has more free parametersfor his disposal which can be used to obtain higher trigonometricorder methods with relatively lower step number. k = 2 Implicit In this case, one has to solve the following linear system of five equations for the parameters al, a2, b0, bl, b2 to obtain a method of trigonometric order 2. = 0, l+ al + a2 e5 cos
e*D sin
e4' cos
2~ +
2V +
4~ +
ale' cos v + a - boeLr(p cos 2v - v sin 2
aleu sin
v
-
bleu(u cos I -
-
bOe2u(p sin
-
ble'(u sin
e41.1 sin 4v + a e2% sin 2V 1
ble2’(2u
v +
cos
sin V) - nb2 = 0 ,
2V + V cos
ale*' cos 2V+a2 - boe4u(2p cos -
v
4~
2V)
v
cos
-
2~
2V)
V)
-
sin
vb2 = 0 ,
(18)
@J)
2v - 2v sin 2v) - 2ub2 = 0 ,
sin 4V + 20 cos 4V) - ooe4p(21.1 - ble*u(2u sin 2v +
2J
cos
2V)
-
2vb2 = 0 .
The system was solved by REDJCE [9]. The expressionsfor :he coefficients are complicatedbut SEDUCE produces an output in the form of Fortran Statements that can be incorporatedin a computer program for numerical experiments with such a method.
165
2.3 Minimax Methods In this section we discuss minimax methods, i.e., methods obtained by satisfyingconditions (9). These conditions can be wrftten in terms of all, bllas follows:
cos(k-f)"(j)
=
e(k~~)~(j’~_“(j)s~n(k_~)~(~)+~(~)co8(k_~)~(j f, (19) j = 1,2,.*.,q s
j
= 1,2,...rq ,
ih~(j) are the zeros of the function (p(ih.hx) such that it has a small maximum norm in the rectanglewI12 w ~wu, $'a5 JI( qu. To obtain the best approximationin this case is certainly not easy; but we will assume that one can write +(ihX(j)) = $(!!(j),w(j)) as a product of 2 one variable functions. Thus $(j) and w(j) can be taken as Chebyshev's points on the correspondinginterval,i.e.
+Cii = $JJ, 1 + q,j +
$(& - qpos y
w(j) = $(Wk + wuj + +wu - Wk)COS F
ll
,
(21) 71
)
j = 1,2,...,q . For this choice of points, one can evaluate the coefficientsaf, bf by solving Q(O) = 0 , Q(&),w(j))
(22) = 0 , ; = 1,2,...,q .
We call such methods product min;max (pM2). The number of free narametersfor i&licit methods 2k+l and the number of equations is 2q+l, thus, the trigonometricorder q is equal to the step number k. k = 1 Implicit In this case
(23)
and the system of equations can be solved for bC, bl. This yields the coefficientsgiven by (12) where
166
P = _,p
)
”
=
hw(l) .
(24)
Thus, the product minimax method would suggest using the center of the rectangle I~~,&.Jx~~~,wulas ho* To obtain a product minimax method of trigonometric order 2, one has to solve a system of 5 equationssimilar to equation (18) with the unknowns b0, bl, b2, al, a2. The difference is that in the last 2 equations one should replace 2u by D(2) and 2~ by ~(~1. In the second and third equations of (18), the 1-1, V, should be replaced by u(l), V(l) respectively. The resultingsystem can be solved by MACSMA [4] or REDUCE 191. In the next section we implementtwo methods cf trigonometricorder 1 and 2 and see how the product minimax methods compare with fitting methods. 3. Numerical Example In this sectionwe compare Adams fitting method of trigonometricorder 1, the generalizedfitting method of trigonometricorder 2 obtained when solving system (18), the product minimax method of trigonometricorder 1 given by (23) and of trigonometricorder 2 obtained when solving the system (19)-(21). Both systems (18) and (19)-(Lllwere solved by REDUCE which produced a F3RTRAN subroutinefor the evaluationof the coefficients. This subroutine is called only once during the integration. The methods were compared for the solution of the initial value problem . 5 i(1 + i)x = 0 W O
(26)
thus
x =;(l+i),
$--Jj,
w=-.71 2
(27)
Iirorder to avoid complex arithmetic,we rewrite the differentialequation as a system of equations for the real and imaginarypart of z = u + iv. I;+;(u + v) = 0 , . v- ;(u - v) = 0 ,
oit<4 (28)
u(0) = 1 , v(0) = 0 . The system is solved by fitting methods of trigonometricorder 1 and 2 witn h = .Ol and various values of II, and w. In Table 1 we list the Euclidean norm of the error at t = 4. It is clear that the method is not sensitive to perturbationsin the values of 9 and w.
167
error
0 0 .1 .1 0 .2
hw
first order
second order
0 .1 0 .l .2 0
.3678 .4482 .4346 .6342 .9241 .8706
.4433 -3508 .2544 .4421 .8126 .5095
(-12) (-7) (-7) (-7) (-7) (-7)
(-12) (-11) (-11) (-11) (-11) (-11)
Table 1 Using the product minimax methods of trigonometric order 1 and 2 with h = .Ol and various squares centered at $ = -r/2, w = */2, the error is much larger but again is insensitive to small perturbations in the length of the sides of the squares. In Table 2 we list the Euclidean norm of the error at t = 4.
length of side .4 .8 1.2 1.6 2
I]
el m x
first order .3678 .3678 .3678 -3678 .3678
Table 2
(-12) (-12) (-12) (-12) (-12)
second order
1 .1469 .2348 .1186 .3728 .9001
(-6) (-5) (-4) (-4) (-4)
Note that the perturbations in the product minimax methods are lnrger than those allowed in the fitting methods. It is possible that the larger errors in the product minimax methods are due to the assumption that $ can be written as a product of 2 one variable functions. Also note that for the first order method, one always gets a good result since PM2 always uses the centei of the square. Acknowledgment The author would like to thank the NPS Foundation Research Program for its support of this research. REFERENCES 1.
Babcock, P.D., Stutzman, L. F., Brandon, D. M., Imprwements in a single-step integration algorithm, Simulation 28, 1 (1979), l-10.
2.
BrandJn, D. M., A new single-step integration algorithm wfth Astability and improved accuracy, Simulation 23, 1 (1974), 17.
3.
Bettis, D. G., Numerical integration of products of Fourier and ordinary polynomials, Numer. Math., g (1970), 421-434.
4.
Bogen, R. et al., MACSYMA Reference Manual, MIT Press, Cambridge, MA, 1975.
5.
Brock, P., Murray, F. J., The use of exponential sums in step by step integration, Math. Tables Aids Comput. 5 (1952), 63-78.
6.
Carroll, J., Exponentially fitted one-step methods for the numerical solution of the scalar Riccatl Equation, IMA J. Numer. Anal., to appear.
168
7_
Dennis, S. C. R., The numerical integration of ordinary differential equations possessing exponential type solutions, Proc. Cambridge Phil. Sot. -65 (1960), 240-246.
8.
Gautschi, W., Numerical integration of crdinary differential equations based on trigonometric polynomials, Numer. Math., 2 (1961), 381-397.
9.
Hearn, A. C. (ea.), REDUCE user's manual, Rand Corporation, Santa Monica, CA, 1983.
10:
van Eer Houwen, P. J., Sommeijer, B. P., Linear multistep methods with minimized truncation error for periodic initial value problem, IMA J. Numer. -Anal.,5 (1984), 479-489.
11.
Lambert, J. D., Computational Methods in Ordinary Differential Equations, Wiley, London, 1977.
12.
Lyche, T., Chebyshevian multisiep methods for ordinary differential equations, Numer. Math. -19 (1972), 65-75.
13.
Neta, B., Ford, C. H., Families of methods for ordinary differential equations based on trigonometric polynomials, J. Comp. Appl. Math., -10 (1984), 33-38.
14.
Neta, B., Families of backward differentiation methods based on trigonometric polynomials, Intern. J. Computer Math., 20(1986). 67-75.
15.
Neta, B., Williams, R. T., Hinsman, D. E., Studies in a shallow water fluid model with topography, Proc. Eleventh IMXS World Congress, Oslo, Norway, August 5-9. 1985 (R. Vichnevetsky, ed.).
16.
Pratt, D. T., Exponential-fitted methods for solving stiff systems of ordinary differential equations, Proc. Eleventh IMACS World Congress, Oslo, Norway, August 5-9, 1985 (R. Vichnevetsky, ed.).
17.
Salzer, B. E., Trigonometric interpolation and predictor-corrector formulas for numerical integration, ZAMM, -42 11962). 403-412.
18.
Steifel, E., Bettis, D. G., Stabilization of Cowell's method, Numer. Math., -13 (1969), 154-175.
APPENDIX Here we show that the shallow water equations with topography have a solution cf the form eh*, where A is complex. This system of equations consists of three equations with three forecast var;ables, u, v and @. The equations are:
$+,?E+,au-
fv+$U, ay
$+u%+v g
(A. 2)
av,f,,g=,, ay
+ $po$-lbBN
+ $p(rt-Q,)1
= 0 ,
where Q = gh is the geopotential height (h = height of free surface), $B is the bottom topography (assumed to be independent of time), u, v are the components of the -wind velocity in the x, y direction, respectively, and f is the Coriolis parameter. Linearizing the equations by letting u=U+u',
v=V+v',
$I=@+$'.
where U, V are the constant mean flow and @ is independent of time. ing that U, V are related to @ via the geostrophic relations U
se--- 1 a@ f ay 3
v=55
i a@
Assum-
one obtains the linear system (after dropping the primes):
g+Ug++_
fv+%=O,
aY g+IJ
*+,u++o, ay
g+V
B.+W+x+"~+v~, ay ax
2?L $$+l_L&+v
(A.5) (~4
(A.71
where y = Q,- $b. If the flow is assumed to be along the topographyas in 1153, then the right hand side of (A.7) is zero. In such a case, one can wrfte the solution in the form " = u eS%+rly
-ut)
0 v
_
v
.i(Sx
+rly-0t)
eiKx
+Ily-bt)
0 cp ^- 0 0
,
,
(A.81
,
where S=P-io,
(A.91 n=V-i6. In order for (A.8) to be a solution for (A-5)-(A.7),one must have
x = -0 + CJJ+ nv
(A.lO)
satisfying iX[(iX)2- in(iny + $1
- f[-ihf - i<(iny 4-z)]
+ (i
(A.ll)
The real and imaginaryparts of xr = -5 + pll+ vv ) (A.12) xi = -pu - ev , satisfy the following system of equations (after dropping nonlinear terms iu X)
In general, X is complex and, thus, the shallow water equations have a solution in the class of problems to be discussed here.
SPECIAL METHODS FOR PROBLEW WRXiE OSCILLATORYSOLUTION IS DAMPED Beny Neta Naval Postgradua+eSchool Departmentof Mathematics Code 53Nd Monterey, CA 93943 ABSTRACT This paper introducesmethods tailored especially for prcblems whose solution behaves like eXX, where h is complex. The shallow water equations with topographyadmit such solution. This paper complementsthe results of Pratt and others on exponentialfitted methods and those of Gautschi,Neta, van der Wouwen and others on trigonometrically-fitted methods. 1.
Introduction In this paper we consider linear multistep methods h
1, n 2 k-l f bff(xn+l_f,~~+l_~),k _?. !2,=0
(1)
for integratingthe initial value problem (2)
Y'(X) = f(x,y(x)), Y(Xo) = Y. * This linear multistep method is characterizedby the polynomials P(S) = krO aetk-' ,
~(0
=
F beck-” 11=0
.
(3)
The main assumptionof this paper is that it is a priori known that the solution is approximatelyof the form m iXjt 1 cje j=l.
Y(X) - no +
(4)
where hj = Wj + iJlj,and the frequencieswj are in a given interval [wc,wu]. The special case where . = jwo with WC given was considered first by Gautschi [B]. His approach ias the following. Let: 9(Z) = p(e') - za(e'>
(5)
then the local truncationerror of (1) is given by Lambert [ill Tn+k = $(h &)
y(tn) .
(6)
Inserting (4) in (6) yields Tn+k - cP(0)cO+
I
iXjt , ij = jw, ci4(ihXj)e
(7)
j=l The coefficientsbi are chosen in such a way that $(ihjw,) = 0 ,
j =
O,l,...,q ,
@Elsevier Scrence Publmhlng Co.,lnc.,1999 655Avenue of ths Americas, New York, NY 10010
(8)
for the largest value of q possible. q is then called the trigonometric order of the method. Gautschi has chosen ai such that the methods are those of Adams and Stormer type. However, these methods are sensitive to changes in the frequency w0. Neta and Ford [13] developed Nystrom and generalized Milne-Simpson type methods. These methods showed less sensitivity to perturbation in w0 but require the eigenvalues of the Jacobian to be purely imaginary. Neta [14] has developed families of backward differentiation methods that overcome the above-mentioned restrrctfon. Salzer [17] has developed predictor-corrector methods based on trigonometric polynomials. See also Steifel and Bettis [18] and Bettis [3]. Van der Houwen and Sommeijer [lo] have developed an alternative approach. The conditions (8) were replaced by
wx = 0
I
(9) $(ihX(j)) = 0 ,
j = l,i,...,q ,
where the A(j) are appropriately chosen points in the fnterval [wa,wU]. An advantage of this so-called minimax approach over the fitting approach is the increased accuracy in cases where no accurate estimate of wD is available or when the frequency is varying fn time. The other special case considered in the literature is where Xj = i$j. Probably the first article on the subject is due to Brock and Murray [5]. They discuss the use of exponential sums in the integration of a system of first order ordinary differential equations. Dennis [7] also suggested special methods for problems whose solution Is exponential. He suggested a transformation of variables. More recently, Carroll 163 has developed exponentfally fitted one-step methods for the scalar Riccati equation. For the general first order system of equations, Pratt [Hi] suggests methods based on the three parameter exponential function f(x)
= A + BeZX .
(10)
The parameters A, B are given in terms of values of y and f. $everal possibilities for z are given based on results of Brandon [2] and Babcock et al. [I]. Lyche [12] analyzes multistep methods which exactly integrate the set wx {xme n ], where wn is real or imaginary. In this article we developed various methods fitting exponentials and methods obtained via the minimax approach. 2.
Construction of Methods
2.1 Fitting Methods In this subsection we discuss various fitting methods. To this end, we separate $(ihjA) = 0, j = l,Z,...,q, into real and imaginary parts. This yields the following equations relating the coefficients aE,bt, cos
jv(k-ll)-
i R=O
b .ju(k-a) [-JV . cos jv(k-ll) ,c - jv sin jv(k-!Z)]= 0 ,
[ju sin jV(k-a) sin jV(k-ll)- t aeejp(k-k) bllejU(k-J1) !2=0
ll=o
f jv cos jV(k-r?)]= 0 , where X =w+i$,
u = -h$, V = hw, j = l,Z,...,q.
163
For explicitmethods, b0 = 0. For Adams type methods a0 = 1, al = -1, ai = 0 for i = Z,...,k. For Nystrom or generalizedMilne-Simpsonmethods a0 = 1, a2 = -1 and other ai = 0. k = 1 Implicit Adams b. 5
we-’ +
sin V - v cos v
1-I
9
sin V
b2+v2)
bl = ve' - p sin v - v cos v ’
(p2+v2) sin v
1 - cos v =b 1, which agree with sin v Gautschi [S] if the coefficientsare expanded in Taylor series with respect to v).
For $ = 0, the coefficientsbecome b0 =
k=
2 Explicit
Adams bl =
(n sin 2v - v cos 2v)ep + v cos v - p sin V (p2+v') sin v (13)
= (v cos v - u sin V)e2' - we'
b 2
.
(p2*V2) sin V
Nystrom (n sin 2v - v cos 2")s' + Ve-' bl =
I
(v2+V2) sin V (14)
b2 =
A (v
cos
2lJ
v -
(p’+V2)
. sin
V
sin V For $ = 0, the coefficientsbecome bl = 2 7, Neta [14]. k=
b2 = 0 which agree with
2 Implicit
In this case, one obtains a one-parameterfamily of (Adams, generalized Milne-Simpson)methods of trigonometricorder I. The free parameter can be used to increase the algebraic order of the method as in 1131. Backward Differentiation e2’
cos
2v + a e' co6 v + a - b e2'(p cos 2v 1 2 0 - bOe2'(p sin
ezp sin 2v + ale' sin V 1
+ al
+ a2
2~
+
=Z 0
,
V
sin
2V)
v
COS
?V) = 0 3
(15)
=o.
This system can be solved by MACSYMA (ProjectNAG's SYmbolic MAnipulation system written in LISP and used for performing symbolic as well as numerical mathematicalmanipulation [4]) or by REDUCE 191. The solution is
ve*n - Ilsin 2v - v cos 2v al= -eu(v cos V + u sin V) + u sin 2V + V cos 2V ' (16)
a2 = -1 - al , _e*u sin
b0 =
V + e’
sin
Zv -
sin V
-e2qv cos v + D sin L) + eu(u sin 2v + v coo 2v) * For 9 = 0, the coefficientsagree with those given in 1141. k = 3 Explicit Again here, one obtaiqs a one-parameterfamily of methods of trigonometric order 1. In ordc.,to get methods of trigonometricorder 2, one has to constructa 3 step -mpiicitmethod of Adams or generalizedMilne-Simpson q-se. In order to in:rease the trigonometricorder without going to a higher step number, 3ne can constructliusar multistepmethods for which the coefficientsallare also functionsof J1,w. Some examples are given in the next subsection. 2.2 GeneralizedFitting Methods In this section,we constructsome linear multistepmethods of the form f
all(h)yel_ll
= h
!fba(A)fn+l_E
R=O
L=O
WI
k>l,nlk-1. -
3
Since a% are functions of X one has more free parametersfor his disposal which can be used to obtain higher trigonometricorder methods with relatively lower step number. k = 2 Implicit In this case, one has to solve the following linear system of five equations for the parameters al, a2, b0, bl, b2 to obtain a method of trigonometric order 2. = 0, l+ al + a2 e5 cos
e*D sin
e4' cos
2~ +
2V +
4~ +
ale' cos v + a - boeLr(p cos 2v - v sin 2
aleu sin
v
-
bleu(u cos I -
-
bOe2u(p sin
-
ble'(u sin
e41.1 sin 4v + a e2% sin 2V 1
ble2’(2u
v +
cos
sin V) - nb2 = 0 ,
2V + V cos
ale*' cos 2V+a2 - boe4u(2p cos -
v
4~
2V)
v
cos
-
2~
2V)
V)
-
sin
vb2 = 0 ,
(18)
@J)
2v - 2v sin 2v) - 2ub2 = 0 ,
sin 4V + 20 cos 4V) - ooe4p(21.1 - ble*u(2u sin 2v +
2J
cos
2V)
-
2vb2 = 0 .
The system was solved by REDJCE [9]. The expressionsfor :he coefficients are complicatedbut SEDUCE produces an output in the form of Fortran Statements that can be incorporatedin a computer program for numerical experiments with such a method.
165
2.3 Minimax Methods In this section we discuss minimax methods, i.e., methods obtained by satisfyingconditions (9). These conditions can be wrftten in terms of all, bllas follows:
cos(k-f)"(j)
=
e(k~~)~(j’~_“(j)s~n(k_~)~(~)+~(~)co8(k_~)~(j f, (19) j = 1,2,.*.,q s
j
= 1,2,...rq ,
ih~(j) are the zeros of the function (p(ih.hx) such that it has a small maximum norm in the rectanglewI12 w ~wu, $'a5 JI( qu. To obtain the best approximationin this case is certainly not easy; but we will assume that one can write +(ihX(j)) = $(!!(j),w(j)) as a product of 2 one variable functions. Thus $(j) and w(j) can be taken as Chebyshev's points on the correspondinginterval,i.e.
+Cii = $JJ, 1 + q,j +
$(& - qpos y
w(j) = $(Wk + wuj + +wu - Wk)COS F
ll
,
(21) 71
)
j = 1,2,...,q . For this choice of points, one can evaluate the coefficientsaf, bf by solving Q(O) = 0 , Q(&),w(j))
(22) = 0 , ; = 1,2,...,q .
We call such methods product min;max (pM2). The number of free narametersfor i&licit methods 2k+l and the number of equations is 2q+l, thus, the trigonometricorder q is equal to the step number k. k = 1 Implicit In this case
(23)
and the system of equations can be solved for bC, bl. This yields the coefficientsgiven by (12) where
166
P = _,p
)
”
=
hw(l) .
(24)
Thus, the product minimax method would suggest using the center of the rectangle I~~,&.Jx~~~,wulas ho* To obtain a product minimax method of trigonometric order 2, one has to solve a system of 5 equationssimilar to equation (18) with the unknowns b0, bl, b2, al, a2. The difference is that in the last 2 equations one should replace 2u by D(2) and 2~ by ~(~1. In the second and third equations of (18), the 1-1, V, should be replaced by u(l), V(l) respectively. The resultingsystem can be solved by MACSMA [4] or REDUCE 191. In the next section we implementtwo methods cf trigonometricorder 1 and 2 and see how the product minimax methods compare with fitting methods. 3. Numerical Example In this sectionwe compare Adams fitting method of trigonometricorder 1, the generalizedfitting method of trigonometricorder 2 obtained when solving system (18), the product minimax method of trigonometricorder 1 given by (23) and of trigonometricorder 2 obtained when solving the system (19)-(21). Both systems (18) and (19)-(Lllwere solved by REDUCE which produced a F3RTRAN subroutinefor the evaluationof the coefficients. This subroutine is called only once during the integration. The methods were compared for the solution of the initial value problem . 5 i(1 + i)x = 0 W O
(26)
thus
x =;(l+i),
$--Jj,
w=-.71 2
(27)
Iirorder to avoid complex arithmetic,we rewrite the differentialequation as a system of equations for the real and imaginarypart of z = u + iv. I;+;(u + v) = 0 , . v- ;(u - v) = 0 ,
oit<4 (28)
u(0) = 1 , v(0) = 0 . The system is solved by fitting methods of trigonometricorder 1 and 2 witn h = .Ol and various values of II, and w. In Table 1 we list the Euclidean norm of the error at t = 4. It is clear that the method is not sensitive to perturbationsin the values of 9 and w.
167
error
0 0 .1 .1 0 .2
hw
first order
second order
0 .1 0 .l .2 0
.3678 .4482 .4346 .6342 .9241 .8706
.4433 -3508 .2544 .4421 .8126 .5095
(-12) (-7) (-7) (-7) (-7) (-7)
(-12) (-11) (-11) (-11) (-11) (-11)
Table 1 Using the product minimax methods of trigonometric order 1 and 2 with h = .Ol and various squares centered at $ = -r/2, w = */2, the error is much larger but again is insensitive to small perturbations in the length of the sides of the squares. In Table 2 we list the Euclidean norm of the error at t = 4.
length of side .4 .8 1.2 1.6 2
I]
el m x
first order .3678 .3678 .3678 -3678 .3678
Table 2
(-12) (-12) (-12) (-12) (-12)
second order
1 .1469 .2348 .1186 .3728 .9001
(-6) (-5) (-4) (-4) (-4)
Note that the perturbations in the product minimax methods are lnrger than those allowed in the fitting methods. It is possible that the larger errors in the product minimax methods are due to the assumption that $ can be written as a product of 2 one variable functions. Also note that for the first order method, one always gets a good result since PM2 always uses the centei of the square. Acknowledgment The author would like to thank the NPS Foundation Research Program for its support of this research. REFERENCES 1.
Babcock, P.D., Stutzman, L. F., Brandon, D. M., Imprwements in a single-step integration algorithm, Simulation 28, 1 (1979), l-10.
2.
BrandJn, D. M., A new single-step integration algorithm wfth Astability and improved accuracy, Simulation 23, 1 (1974), 17.
3.
Bettis, D. G., Numerical integration of products of Fourier and ordinary polynomials, Numer. Math., g (1970), 421-434.
4.
Bogen, R. et al., MACSYMA Reference Manual, MIT Press, Cambridge, MA, 1975.
5.
Brock, P., Murray, F. J., The use of exponential sums in step by step integration, Math. Tables Aids Comput. 5 (1952), 63-78.
6.
Carroll, J., Exponentially fitted one-step methods for the numerical solution of the scalar Riccatl Equation, IMA J. Numer. Anal., to appear.
168
7_
Dennis, S. C. R., The numerical integration of ordinary differential equations possessing exponential type solutions, Proc. Cambridge Phil. Sot. -65 (1960), 240-246.
8.
Gautschi, W., Numerical integration of crdinary differential equations based on trigonometric polynomials, Numer. Math., 2 (1961), 381-397.
9.
Hearn, A. C. (ea.), REDUCE user's manual, Rand Corporation, Santa Monica, CA, 1983.
10:
van Eer Houwen, P. J., Sommeijer, B. P., Linear multistep methods with minimized truncation error for periodic initial value problem, IMA J. Numer. -Anal.,5 (1984), 479-489.
11.
Lambert, J. D., Computational Methods in Ordinary Differential Equations, Wiley, London, 1977.
12.
Lyche, T., Chebyshevian multisiep methods for ordinary differential equations, Numer. Math. -19 (1972), 65-75.
13.
Neta, B., Ford, C. H., Families of methods for ordinary differential equations based on trigonometric polynomials, J. Comp. Appl. Math., -10 (1984), 33-38.
14.
Neta, B., Families of backward differentiation methods based on trigonometric polynomials, Intern. J. Computer Math., 20(1986). 67-75.
15.
Neta, B., Williams, R. T., Hinsman, D. E., Studies in a shallow water fluid model with topography, Proc. Eleventh IMXS World Congress, Oslo, Norway, August 5-9. 1985 (R. Vichnevetsky, ed.).
16.
Pratt, D. T., Exponential-fitted methods for solving stiff systems of ordinary differential equations, Proc. Eleventh IMACS World Congress, Oslo, Norway, August 5-9, 1985 (R. Vichnevetsky, ed.).
17.
Salzer, B. E., Trigonometric interpolation and predictor-corrector formulas for numerical integration, ZAMM, -42 11962). 403-412.
18.
Steifel, E., Bettis, D. G., Stabilization of Cowell's method, Numer. Math., -13 (1969), 154-175.
APPENDIX Here we show that the shallow water equations with topography have a solution cf the form eh*, where A is complex. This system of equations consists of three equations with three forecast var;ables, u, v and @. The equations are:
$+,?E+,au-
fv+$U, ay
$+u%+v g
(A. 2)
av,f,,g=,, ay
+ $po$-lbBN
+ $p(rt-Q,)1
= 0 ,
where Q = gh is the geopotential height (h = height of free surface), $B is the bottom topography (assumed to be independent of time), u, v are the components of the -wind velocity in the x, y direction, respectively, and f is the Coriolis parameter. Linearizing the equations by letting u=U+u',
v=V+v',
$I=@+$'.
where U, V are the constant mean flow and @ is independent of time. ing that U, V are related to @ via the geostrophic relations U
se--- 1 a@ f ay 3
v=55
i a@
Assum-
one obtains the linear system (after dropping the primes):
g+Ug++_
fv+%=O,
aY g+IJ
*+,u++o, ay
g+V
B.+W+x+"~+v~, ay ax
2?L $$+l_L&+v
(A.5) (~4
(A.71
where y = Q,- $b. If the flow is assumed to be along the topographyas in 1153, then the right hand side of (A.7) is zero. In such a case, one can wrfte the solution in the form " = u eS%+rly
-ut)
0 v
_
v
.i(Sx
+rly-0t)
eiKx
+Ily-bt)
0 cp ^- 0 0
,
,
(A.81
,
where S=P-io,
(A.91 n=V-i6. In order for (A.8) to be a solution for (A-5)-(A.7),one must have
x = -0 + CJJ+ nv
(A.lO)
satisfying iX[(iX)2- in(iny + $1
- f[-ihf - i<(iny 4-z)]
+ (i
(A.ll)
The real and imaginaryparts of xr = -5 + pll+ vv ) (A.12) xi = -pu - ev , satisfy the following system of equations (after dropping nonlinear terms iu X)
In general, X is complex and, thus, the shallow water equations have a solution in the class of problems to be discussed here.