- Email: [email protected]
A method for finding the periodic solution of autonomous ordinary differential equations
A Method for Finding the Periodic Solution of Autonomous Ordinary Differential Equations Chang
Jin Jiang”
International Centre for Theoretical Trieste, Italy
Transmitted
Physics
by John C&i
ABSTRACT Using the least square method, the problem of computing periodic solutions of autonomous ordinary differential equations has been reduced to a problem of minimization on a normal plane. The conjugate gradient method is used to yield a successive initial value, and the RungeKutta method is used to solve the initial value problem in the ordinary differential equations. So a new method for this problem is formed. To illustrate, it is applied to the Brusselator model and Field-Noyes model of Belousov-Zhabotinskii reactions. The numerical results obtained show that this method is very efficient.
1.
INTRODUCTION
The existence of stable oscillation is well known in some ecological, chemical, and biochemical systems [l, 21. Modeling mathematically for these systems the autonomous ordinary differential equations are obtained and the oscillations are described by periodic solutions of the systems. By using a qualitative method one can demonstrate the existence and outline the periodic solution in some special region [3, 41, but only by using the *Permanent China.
APPLIED
address:
MATHEMATICS
Department
AND
of Mathematics,
COMPUTATION
@ Elsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010
USTC
66:161-173
Hefei, Anhui,
(1994)
230026,
161
0096-3003/94/$7.00
162
C. 3. JIANG
numerical numerical
method one can find the solution concretely and minutely. A method for finding the periodic solution of nonlinear autonomous
differential equations was described in [5], but it was also pointed out that the application of this method to more than two-dimensional problems is difficult. We know the right side of the differential equations defines a direction field. A point can be chosen as an approximate initial value and the direction at this point can be taken as a normal vector of the plane passing this point; here the plane is called the normal plane. In this paper the problem for computing the periodic solution of nonlinear autonomous differential equations in n-dimension is reduced to a minimization problem on the normal plane, which is (n - 1)-dimensional. The conjugate gradient method is used to yield successive
initial
values,
and the Runge-Kutta
method
is
used to solve the initial value problem in ordinary differential equations. By this method the computational efficiency can be greatly improved and the demand on accuracy of initial value is very little. This method will be successful numerical
only if the normal plane cuts the periodic solution practice in this paper shows that our approaches
curve. The are feasible
and effective.
2.
ALGORITHM A general
autonomous
ordinary
g here u = (~1, ~2, it has a periodic
, u,) T, 44 solution
U(t),
differential
equation
is given below:
(2.1)
= A(u); = VI(U),
h(u),
with period
T,
. . . , .f,(~))~. i.e.,
u(t)
We wwe
= U(t + T)
for
arbitrary t. If a point on the solution curve u(t) denoted by 6(O) is given, the problem for periodic solution could become a usual initial value problem in ordinary differential equations. However, such a point and the period T are unknown. The right side of (2.1) defines a direction field in a n-dimensional region. Take a point u(‘)(O), which is close to U(O), but not necessary on the periodic solut.ion curve. Let
A(d"'(0)) n =
be the unit vector
of A(u(‘)(O)).
(24
I14~‘“‘@))l12
The pl ane with the normal
vector
n and
Periodic
Solutions of ODES
passing
the point u(O) (0) is
163
nT(u which is called the normal
- u(O)(O)) = 0,
plane.
(2.3)
We take a disc region with radius R and
center u(‘)(O) on the normal plane. When t E [0, T] the periodic solution
is a closed curve in n-dimensional
space. We suppose the curve intersects with the normal plane but does not contact it. Generally speaking, the solution curve of (2.1), which goes from point u(O)(O), is not closed. According to continuity we can suppose it will intersect with the disc again when the t increases if the point u(O) (0) is sufficiently near to the periodic solution curve. We are now in the position to solve the initial problem differential
$ = A(u) u(0) = uqo)
i=o,1,2
(2.4)
,....
The solution curve goes from UC’)(0)) which is in the the normal plane. When it comes back to the disc again increase oft, the intersection point and the correspondent t T ci) can be regarded as an ,u,‘;hz )a;z,(i), respectively. According that
of the ordinary
equations:
to the law of the least square and taking
u(T(%)), u(O) lie on the normal
disc region on along with the are denoted by approximation
into consideration
plane, let
J(4))) = ;(u(O)- u(T), ~(0) - u(T)) nT(u(T)
- r&O’(O)) = 0
nT(u(0)
- u(O) (0)) = 0.
(2.5)
When we chose u(O) on the normal plane to minimize the inner product, namely the right side of first formula in (2.5), a minimization problem in (n - l)-dimension space is formed. If the solution curve has intersection point with the normal plane, the minimization problem at least has a solution u(t) and J(u(0)) = 0. By using the u(O) as initial value of (2.4), the solution u(t) of (2.4) and T(‘) are approximations of the periodic solution and the period, respectively. We use the conjugate gradient method to solve the minimization problem (2.5), here u(O) is an independent variable on the normal plane. To this end we have to compute 6’J(u(O))/&r(O) in every step. We use the perturbation method to compute it. Suppose u(O) has a small perturbation
C. J. JIANG
164
au(O) in th e d’rsc on the normal plane. When ult=a = u(O), u(0) + &u(O), the solutions and the approximate periods of (2.4) are denoted by u, ii and T, T, respectively. Let 6u(t) = ii(t) - u(t), ST = ?; - T, SJ(u(0)) = J(u(0) + 640)) - J(u(O)), p(t) 1s . a n-dimensional vector function and We have p(T) =u(T) -u(O), a = pT(T)A(u(T))/nTA(u(T)). THEOREM 1. If u(0) h as a small perturbation plane, then the minimization problem (2.5) becomes Wu(0))
= (-p(T),WO))
Su(0)
on the normal
+ (P(T) - an,WT)).
(2.6)
PROOF. SJ(u(0))= (u(O)- u(T), Su(0) + u(T) - ii(T)) +~@u(O)+u(T)-u(T)&(O)+u(T)-s(T))
are small, by omitting the higher-order
Since both &r(O) and u(T) -U(T) small terms we obtain SJ(u(0))
= (u(O) - u(T), Su(0) + u(T) -U(F)).
Expanding A@(T)) at the point u(T) terms, we obtain also u(T)
-ii(T)=u(T)
(2.7)
and omitting the higher-order small
-ii(T)
+ii(T)
= -Su(T)
- A(a(T))6T
=-&I(T)
- A(u(T))ST.
-ii(T) (2.8)
The second formula in (2.5) nT(U( T) - u(T))
= nT(6u(T)
results in ST = -
+ A(u(T))ST)
= 0
nT6u(T) nTA(u(T))
’
place it in (2.7), (2.8) and note that
-P(T)> n;y($(u('.)))
= (-
P~(T)A(~(T)), &(T) nTA(u(T)) ’
a=~T(T)A(u(T)) nTA(u(T))
Periodic
Solutions
of ODES
We have SJ(u(0))
=
- P(T),WO) +A(u(T))
nTSu(T) nTA(G”))
= (-P(T)>WO))
+
-WT)
>
+(P(T)>wT))
-P(W(~(T))~;~($;)) (
= (-p(T),WO)) + (P(T) -an,WT)). So (2.6)
holds.
To consider
The proof is complete. (2.4),
0
the following results
d&(t) dt
are used:
aA
= dubu(t)
(2.9) = h(O)
hI(t) Here afl -
au, af2
afl -
af2
aA(U>= -au1 -au2
-a&f1
...
-af2
au,
au, ... ... ... ... afn afn ... -afn au, au2 au,
t3U
Letting
. . .
au2
(2.10)
t E [0, T], we have
THEOREMS.
Let p(t)
be a solution of the following
adjoint
differential
equations:
T(p(t) -
an)
= 0
(2.11)
p(T) = u(T) - u(O); then
aJ(U(o)) = p(O)-p(T)-
WO)
an.
(2.12)
C. J. JIANG
166
PROOF.
From (2.9) we have
(
%,p(t)-un
Integrating
(hu,P
aA(u) 3u6u,p--an
)
(
on [0, T] by parts the following
- un)IT - l’
(au(t),
=O. )
$)dt
- lT
is yielded:
(%6u(t),p
- an)dt = CL
Since (
aA(u) _)UEu,p-an)
= ([F]T(p-an),Eu)
then
(WT),
By using (2.11)
P(T)
- un) - (Su(O), P(O) - un)
we have
(WV, P(T) - 4 and putting
P(O) - un);
it into (2.6) (-p(T),
640))
= (-p(T),
640))
~J(u@)) =
= (P(O) - P(T) we obtain
= (NO),
aJ(u(o)) au(o)
+ (P(T)
- un, au(T))
+ (P(O) - un, 640)) - un, 640))
=p(O)-p(T)--an. 0
This proves the theorem.
We will look for the minimum point u(0) only on the normal plane. Both u(0) and u(T) are on the normal plane. Whether p(0) - an is on it or not, taking into consideration errors in calculation we have to compute the project of p(0) - an on the normal plane and choose the gradient vector
g = P(O) - (n, p(O))n
- p(T).
(2.13)
Periodic
Solutions
of ODES
167
Particularly, if the problem is two-dimensional, the normal plane becomes a normal line (the neighboring region is line segment) and the gradient direction coincides with the normal line direction, the procedure for solving p(0) can be omitted and can let g = u(0) -u(T). The procedure for finding out optimal u(0) is simplified as solving the minimization problem on the normal line, which is a one-dimensional problem.
3.
THE COMPUTATION
PROCEDURE
According to the preceding analysis of the algorithm, the procedure for finding the period solution can be described as below. Let u(‘)(O) approximate to u(O), the normal vector n and the normal plane are achieved from the right side of (2.1). By using u(‘)(O) as an initial value we solve If necessary we can solve (2.4); then u(t), T(O), and u(T(‘)) are obtained. the adjoint linear differential equations (2.11)by backward pattern about t, i.e., to solve p(0) from p(T(‘)). g(O) is yielded from (2.13). Let w(O) = &O). When u(‘) (0)) gci), wci) have b een obtained. The (i + l)th step solving the one-dimensional minimization problem may begin. Let
l$+l)(o) = u(qq
J(ui’+”(0))
In order to minimize
- aLw(i).
then the procedure can be completed; solution and T(ifl) is an approximate can solve the adjoint linear ordinary out P(~+')(O), let [6]
yG+l) = w(i+l)
E
(3.2)
ujitl) (t) is an approximate periodic period. If (3.2) does not hold, we differential equations (2.11) to find
(p(i+l)(o),
n)n
_
p(“+l)(T(i+l))
(gG+l), g(“+l) _ g’“‘)
z-g
(g’“‘, g(i)) (i+l) +
yG+l)w(i)
i=o,1,2,...; then turn
(3.1)
we have to find out al and solve (2.4). If
J(Uj”f”(0))<
$$“+I)= p@+l)(o)-
for
to (i + 2)th step.
(3.3)
168
C. J. JIANG
We use second-order minimization
interpolation
problem.
Denote
polynomial
for the one-dimensional
briefly
J(tJi) (0)
-
cyw(q =
J(a).
(3.4)
Given CXO,al, CYZto compute J(ao), J(cti), J(IYz), the three points can define a quadratic order interpolation polynomial of J(o),
UC?+ ba + c;
$(a) =
(3.5)
here (3.6)
The
J(Q~,
W-1),
The minimum
J(w,
CYL--~, ~1-2)
point of q(o)
Qlfl
=
denote first and second difference
of ~(a!).
is
al + QL-1
2
J(w,
-
QL-1)
2J(W,QL-13%2)’
(3.7)
is The LIAR” ’ d e t ermined by ol+i and (3.1). We take ~~~~“(0) as an initial value, then solve (2.4), and compute J(or+i). The largest one among J(al+i), J(ol), J(cr_i), J(al_2) can be omitted. So the new three points and their values are obtained again. We can stop when (3.2) holds or do it some times. We take the largest time equal to five in this paper. If the quadrative interpolation polynomial $(a) does not have minimum point or just J(al+i) J(QI), J(w-l), its zero point
is the largest one, then we take two small values among polynomial and take J(w-2) t o f orm a linear interpolation as QL+~ in order to accelerate convergence. We take the
a0 = O,cyl = -0.8(lu (i)(O) - u(~)(T(~))~~/~I(w(~)II, a2 = 0.5~1 in this paper. The numerical method, for example, the implicit Runge-Kutta method, is used for solving the initial value problem in ordinary differential equations (2.4). The discrete step of t is denoted by h. When uj+i is yielded from uj, we have to check that uj+] is across the disc on the normal plane. When u.j+i is across it, the intersect and the time t, which are denoted by u(Tci)), T(‘), can be defined by using linear interpolation. The coefficient matrix of adjoint linear differential equations (2.11) is obtained by using interpolation too. The steps of t for solving numerical solution of (2.11) are taken to coincide with ones of (2.4) and the computational procedure is taken backward pattern about t.
Periodic
Solutions
169
of ODES
TABLE ZJ(U(O)(O))
u(O) (0)
4.
Time
1.
2J(d”)(O))
T
u(i) (0) 2.806571,2.403286
2.0,2.0
4.013856e-01
5
1.480771e-11
7.1565
2.0,4.0
9.026728e-02
3
3.184155e-10
7.1565
1.684179,3.715761
2.0,6.0
1.733149e+OO
3
8.380585e-10
7.1565
0.646447,4.721644
4.0,2.0
2.397104e-01
5
2.948539e-10
7.1566
3.472134,1.551305
4.0,4.0
1.887851e+OO
3
1.659828e-10
7.1565
2.585821,2.667408
4.0,6.0
5.283966e+OO
3
3.437961e-10
7.1565
1.659904,3.743479
NUMERICAL EXAMPLE
1.
PRACTICE
The famous dx dt
=
Brusslator A
-
model in dissipation
is
(B + 1)~ + 2’1~
dy - = Bx - x’y, dt
(4.1)
where A, B are positive coefficients. When B > 1 + A2 there is a limit cycle or periodic solution [3]. We take A = 1, B = 3. Since n = 2, the normal plane degenerates into a straight line, which is the normal line. The procedure for solving adjoint differential equations can be omitted. The change of sign in the inner product (n, u(t) - u(0)) can judge if the curve is across the normal segment. So the Tci) can be computed. By using the presented method and second-order implicit Runge-Kutta method for solving (2.4) we obtain the results in Table 1 for some different initial values. All real numbers are taken single length arithmetic here. In Table 1 the third column, namely time, indicates the number of times to solve the initial problem (2.4). The last column U(~)(O) is the final approximation of u(0). Table accuracy the limit inside of
1 shows that the convergence rate is very fast and the demand on of initial values is very little. According to computational result cycle is drawn in Figure 1. There is the singular point (1, 3) in the limit cycle,
EXAMPLE 2. B. P. Belousov and A. M. Zhabotinskii noted that different organic acids and metal ions could be used in reaction and still preserve oscillation. Under some assumptions R. J. Field and R. M. Noyes [l] have given a simple model which mechanism retains quantitative features.
170
C. J. JIANG
FIG. 1. The limit cycle.
By
using some transformation
Zhabotinskii
reaction
of variable
their
model
of the Belousov-
is dx - =
s(y -
&
i(fz - y-
dt
xy +
- = dt s dz - = w(x - z),
5
-
4x2) (4.2)
xy)
dt
where s = 77.27, q = 8.375e - 6, w = f = 1.0. It is an ordinary differential equation without initial conditions. S. P. Hastings and J. D. Murray showed that the model equations have an oscillatory or periodic solution [4]. We use the present method to solve the periodic solution of this problem. Since n = 3, we have to solve the adjoint differential equations (2.11). The stiffness of the Field-Noyes model of the Belousov-Zhabotinskii reaction should be conquered by using self-adaptive fourth-order implicit Rung+ Kutta method with variable step. The idea of a self-adaptive method is described
as below:
with variable
step mentioned
here
Periodic
Solutions
of ODES
171
TABLE
g
w
i
2
CY
0
8.333739e-03 -2.00036le-02 1.87661&-01
2J(u@)(O)) 3.566940e-02
T(i) 231.503
-8,333739e-03 2.000361e-02 -1.876618e-01
1.294955e-05
-1.300868e-05
4.923687e-04
-4.922267e-04
5.186635e-05
-5.319774e-05
1
-0.799808
1.426972e-03
231.830
2
-0.399904
1.284108e-02
231.667
3
-0.999756
2.448174e-07
231.911
4
-0.799235
9.793161e-09
231.902
Suppose the step of t is h at moment of t,, we use h to compute u,+r, from u,, denote u,+r by uk+r. We use h/2 to compute two times, namely from u, to compute u,+r/z and then to compute u,+r from u,+~/Z, denote h’2 The ]]ukyr -u~+r((~ can estimate the truncation error. We u,+r by u,+r. restrict h by h,,,, which equals 0.04 in this paper, and Eis the requirement of accuracy. So the procedure is > E, th’ 1s st e P 1s not successful; we omit it, replace (1) If rlu;/;“1 -4+1/12 h by h/2, and do it again.
(2) If llu%
- ui+r]]2 < E, this step is successful. Let u,+r = uiyr and
(a) if e/10 I llu$?r - u~+r]]2 < E, the next step of t is unchanged. (b) if ll$$‘r - di+J2 it is just taken h,,,
< e/10,
the next step oft is taken as 2h, but if 2h > h,,,.
We can use the extrapolation method to improve the accuracy of u,+r C. W. Gear points out [7] that all results about convergency, stability, and
error estimate of the implicit Runge-Kutta method with fixed step is also effective to the method with variable step if the largest step is restrict. This gives a theoretical basis of self-adaptive method with variable step. The computational procedures are recorded in Table 2 and Table 3 for u(‘)(O) = (2.5, 1.65, 2.59) and (10.0, 50.0, 50.0), the final approximate values of u(0) are (2.4916579, 1.6696053, 2.4023413) and (1.0030284, 331.20989, 7651.67908), respectively. Here all real numbers are doubled precision. According to the computational result the periodic solution curve of the model is drawn in Figure 2, where the y-axis is magnified 50 times.
172
C. J. JIANG
TABLE
g
w
i
3 a
0 8.989289e+OO -2.439098efOZ
2.439098e+OZ 7.601184e+03
5.682839e-03 1.195153e+oo
-2.156741e-04
5,78646848e+07
149.009
-8.989289e+OO
-7,601184e+03
-3.724665efOl
T(i)
ZJ(u’i’(O))
1
-0.800188
2.315698e+06
2
-0.400094
2.083184e+07
199.811
3
-1.000222
1,388678e+03
228.425
4
-0.799972
5.557380efOl
231.230
5
-0.399986
4.999374e+OZ
229.853
6
-0.999953
3.334272e-02
231.900
221.021
-5.898636e-03 3.725250e+Ol -l.O12680e+OO
2.156904e-04
4,859032e-03
-4.961698e-03
1.825308e-01
-1.825280e-01 7
-0.800019
1.333710e-03
231.900
8
-0.400009
1.200338e-02
231.900
9
-1.000024
l.OlZOOZe-11
231.899
0
x-120000 FIG. 2. The periodic
solution
curve.
Periodic
The
Solutions
author
of ODES
thanks
Professor
173
Abdus Salam,
Energy Agency, and UNESCO for hospitality for Theoretical Physics, Trieste.
the International at the International
Atomic
Centre
REFERENCES 1
2 3
4
5
6 7
R. Field and R. M. Noyes, Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction, J. Chem. Phys. 160:18771884 (1974). G. Nicolis and I. Prigogina, Self-organization in nonequilibrium systems, Wiley, New York, 1977. Q. Yuaxun and Z. Xianwu, Qualitative investigation of the differential equations of Brusslator in biochemistry, A Monthly Journal of Science 25 (4):273276 (1980). S. P. Hastings and J. D. Murray, The existence of oscillatory solutions in the Field-Noyes model for the Belousov-Zhabotinski reaction, SIAM. J. Appl. Math. 28 (3):678-688 (1975). Z. Suo-Chun and G. Auchmuty, Variational principles and numerical methods for periodic solutions of autonomous ordinary differential equations, Numer. Calculation and App. Cornput., 3:141-148 (1989). W. De Ren, Solution of Nonlinear Equations and Optimization Method, The People’s Education Publishing House, Beijing, 1979. C. W. Gear, Numerical Initial Value Problems in Ordinary Difjerential Equations, Prentice-Hall, 1971.
Jin Jiang”
International Centre for Theoretical Trieste, Italy
Transmitted
Physics
by John C&i
ABSTRACT Using the least square method, the problem of computing periodic solutions of autonomous ordinary differential equations has been reduced to a problem of minimization on a normal plane. The conjugate gradient method is used to yield a successive initial value, and the RungeKutta method is used to solve the initial value problem in the ordinary differential equations. So a new method for this problem is formed. To illustrate, it is applied to the Brusselator model and Field-Noyes model of Belousov-Zhabotinskii reactions. The numerical results obtained show that this method is very efficient.
1.
INTRODUCTION
The existence of stable oscillation is well known in some ecological, chemical, and biochemical systems [l, 21. Modeling mathematically for these systems the autonomous ordinary differential equations are obtained and the oscillations are described by periodic solutions of the systems. By using a qualitative method one can demonstrate the existence and outline the periodic solution in some special region [3, 41, but only by using the *Permanent China.
APPLIED
address:
MATHEMATICS
Department
AND
of Mathematics,
COMPUTATION
@ Elsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010
USTC
66:161-173
Hefei, Anhui,
(1994)
230026,
161
0096-3003/94/$7.00
162
C. 3. JIANG
numerical numerical
method one can find the solution concretely and minutely. A method for finding the periodic solution of nonlinear autonomous
differential equations was described in [5], but it was also pointed out that the application of this method to more than two-dimensional problems is difficult. We know the right side of the differential equations defines a direction field. A point can be chosen as an approximate initial value and the direction at this point can be taken as a normal vector of the plane passing this point; here the plane is called the normal plane. In this paper the problem for computing the periodic solution of nonlinear autonomous differential equations in n-dimension is reduced to a minimization problem on the normal plane, which is (n - 1)-dimensional. The conjugate gradient method is used to yield successive
initial
values,
and the Runge-Kutta
method
is
used to solve the initial value problem in ordinary differential equations. By this method the computational efficiency can be greatly improved and the demand on accuracy of initial value is very little. This method will be successful numerical
only if the normal plane cuts the periodic solution practice in this paper shows that our approaches
curve. The are feasible
and effective.
2.
ALGORITHM A general
autonomous
ordinary
g here u = (~1, ~2, it has a periodic
, u,) T, 44 solution
U(t),
differential
equation
is given below:
(2.1)
= A(u); = VI(U),
h(u),
with period
T,
. . . , .f,(~))~. i.e.,
u(t)
We wwe
= U(t + T)
for
arbitrary t. If a point on the solution curve u(t) denoted by 6(O) is given, the problem for periodic solution could become a usual initial value problem in ordinary differential equations. However, such a point and the period T are unknown. The right side of (2.1) defines a direction field in a n-dimensional region. Take a point u(‘)(O), which is close to U(O), but not necessary on the periodic solut.ion curve. Let
A(d"'(0)) n =
be the unit vector
of A(u(‘)(O)).
(24
I14~‘“‘@))l12
The pl ane with the normal
vector
n and
Periodic
Solutions of ODES
passing
the point u(O) (0) is
163
nT(u which is called the normal
- u(O)(O)) = 0,
plane.
(2.3)
We take a disc region with radius R and
center u(‘)(O) on the normal plane. When t E [0, T] the periodic solution
is a closed curve in n-dimensional
space. We suppose the curve intersects with the normal plane but does not contact it. Generally speaking, the solution curve of (2.1), which goes from point u(O)(O), is not closed. According to continuity we can suppose it will intersect with the disc again when the t increases if the point u(O) (0) is sufficiently near to the periodic solution curve. We are now in the position to solve the initial problem differential
$ = A(u) u(0) = uqo)
i=o,1,2
(2.4)
,....
The solution curve goes from UC’)(0)) which is in the the normal plane. When it comes back to the disc again increase oft, the intersection point and the correspondent t T ci) can be regarded as an ,u,‘;hz )a;z,(i), respectively. According that
of the ordinary
equations:
to the law of the least square and taking
u(T(%)), u(O) lie on the normal
disc region on along with the are denoted by approximation
into consideration
plane, let
J(4))) = ;(u(O)- u(T), ~(0) - u(T)) nT(u(T)
- r&O’(O)) = 0
nT(u(0)
- u(O) (0)) = 0.
(2.5)
When we chose u(O) on the normal plane to minimize the inner product, namely the right side of first formula in (2.5), a minimization problem in (n - l)-dimension space is formed. If the solution curve has intersection point with the normal plane, the minimization problem at least has a solution u(t) and J(u(0)) = 0. By using the u(O) as initial value of (2.4), the solution u(t) of (2.4) and T(‘) are approximations of the periodic solution and the period, respectively. We use the conjugate gradient method to solve the minimization problem (2.5), here u(O) is an independent variable on the normal plane. To this end we have to compute 6’J(u(O))/&r(O) in every step. We use the perturbation method to compute it. Suppose u(O) has a small perturbation
C. J. JIANG
164
au(O) in th e d’rsc on the normal plane. When ult=a = u(O), u(0) + &u(O), the solutions and the approximate periods of (2.4) are denoted by u, ii and T, T, respectively. Let 6u(t) = ii(t) - u(t), ST = ?; - T, SJ(u(0)) = J(u(0) + 640)) - J(u(O)), p(t) 1s . a n-dimensional vector function and We have p(T) =u(T) -u(O), a = pT(T)A(u(T))/nTA(u(T)). THEOREM 1. If u(0) h as a small perturbation plane, then the minimization problem (2.5) becomes Wu(0))
= (-p(T),WO))
Su(0)
on the normal
+ (P(T) - an,WT)).
(2.6)
PROOF. SJ(u(0))= (u(O)- u(T), Su(0) + u(T) - ii(T)) +~@u(O)+u(T)-u(T)&(O)+u(T)-s(T))
are small, by omitting the higher-order
Since both &r(O) and u(T) -U(T) small terms we obtain SJ(u(0))
= (u(O) - u(T), Su(0) + u(T) -U(F)).
Expanding A@(T)) at the point u(T) terms, we obtain also u(T)
-ii(T)=u(T)
(2.7)
and omitting the higher-order small
-ii(T)
+ii(T)
= -Su(T)
- A(a(T))6T
=-&I(T)
- A(u(T))ST.
-ii(T) (2.8)
The second formula in (2.5) nT(U( T) - u(T))
= nT(6u(T)
results in ST = -
+ A(u(T))ST)
= 0
nT6u(T) nTA(u(T))
’
place it in (2.7), (2.8) and note that
-P(T)> n;y($(u('.)))
= (-
P~(T)A(~(T)), &(T) nTA(u(T)) ’
a=~T(T)A(u(T)) nTA(u(T))
Periodic
Solutions
of ODES
We have SJ(u(0))
=
- P(T),WO) +A(u(T))
nTSu(T) nTA(G”))
= (-P(T)>WO))
+
-WT)
>
+(P(T)>wT))
-P(W(~(T))~;~($;)) (
= (-p(T),WO)) + (P(T) -an,WT)). So (2.6)
holds.
To consider
The proof is complete. (2.4),
0
the following results
d&(t) dt
are used:
aA
= dubu(t)
(2.9) = h(O)
hI(t) Here afl -
au, af2
afl -
af2
aA(U>= -au1 -au2
-a&f1
...
-af2
au,
au, ... ... ... ... afn afn ... -afn au, au2 au,
t3U
Letting
. . .
au2
(2.10)
t E [0, T], we have
THEOREMS.
Let p(t)
be a solution of the following
adjoint
differential
equations:
T(p(t) -
an)
= 0
(2.11)
p(T) = u(T) - u(O); then
aJ(U(o)) = p(O)-p(T)-
WO)
an.
(2.12)
C. J. JIANG
166
PROOF.
From (2.9) we have
(
%,p(t)-un
Integrating
(hu,P
aA(u) 3u6u,p--an
)
(
on [0, T] by parts the following
- un)IT - l’
(au(t),
=O. )
$)dt
- lT
is yielded:
(%6u(t),p
- an)dt = CL
Since (
aA(u) _)UEu,p-an)
= ([F]T(p-an),Eu)
then
(WT),
By using (2.11)
P(T)
- un) - (Su(O), P(O) - un)
we have
(WV, P(T) - 4 and putting
P(O) - un);
it into (2.6) (-p(T),
640))
= (-p(T),
640))
~J(u@)) =
= (P(O) - P(T) we obtain
= (NO),
aJ(u(o)) au(o)
+ (P(T)
- un, au(T))
+ (P(O) - un, 640)) - un, 640))
=p(O)-p(T)--an. 0
This proves the theorem.
We will look for the minimum point u(0) only on the normal plane. Both u(0) and u(T) are on the normal plane. Whether p(0) - an is on it or not, taking into consideration errors in calculation we have to compute the project of p(0) - an on the normal plane and choose the gradient vector
g = P(O) - (n, p(O))n
- p(T).
(2.13)
Periodic
Solutions
of ODES
167
Particularly, if the problem is two-dimensional, the normal plane becomes a normal line (the neighboring region is line segment) and the gradient direction coincides with the normal line direction, the procedure for solving p(0) can be omitted and can let g = u(0) -u(T). The procedure for finding out optimal u(0) is simplified as solving the minimization problem on the normal line, which is a one-dimensional problem.
3.
THE COMPUTATION
PROCEDURE
According to the preceding analysis of the algorithm, the procedure for finding the period solution can be described as below. Let u(‘)(O) approximate to u(O), the normal vector n and the normal plane are achieved from the right side of (2.1). By using u(‘)(O) as an initial value we solve If necessary we can solve (2.4); then u(t), T(O), and u(T(‘)) are obtained. the adjoint linear differential equations (2.11)by backward pattern about t, i.e., to solve p(0) from p(T(‘)). g(O) is yielded from (2.13). Let w(O) = &O). When u(‘) (0)) gci), wci) have b een obtained. The (i + l)th step solving the one-dimensional minimization problem may begin. Let
l$+l)(o) = u(qq
J(ui’+”(0))
In order to minimize
- aLw(i).
then the procedure can be completed; solution and T(ifl) is an approximate can solve the adjoint linear ordinary out P(~+')(O), let [6]
yG+l) = w(i+l)
E
(3.2)
ujitl) (t) is an approximate periodic period. If (3.2) does not hold, we differential equations (2.11) to find
(p(i+l)(o),
n)n
_
p(“+l)(T(i+l))
(gG+l), g(“+l) _ g’“‘)
z-g
(g’“‘, g(i)) (i+l) +
yG+l)w(i)
i=o,1,2,...; then turn
(3.1)
we have to find out al and solve (2.4). If
J(Uj”f”(0))<
$$“+I)= p@+l)(o)-
for
to (i + 2)th step.
(3.3)
168
C. J. JIANG
We use second-order minimization
interpolation
problem.
Denote
polynomial
for the one-dimensional
briefly
J(tJi) (0)
-
cyw(q =
J(a).
(3.4)
Given CXO,al, CYZto compute J(ao), J(cti), J(IYz), the three points can define a quadratic order interpolation polynomial of J(o),
UC?+ ba + c;
$(a) =
(3.5)
here (3.6)
The
J(Q~,
W-1),
The minimum
J(w,
CYL--~, ~1-2)
point of q(o)
Qlfl
=
denote first and second difference
of ~(a!).
is
al + QL-1
2
J(w,
-
QL-1)
2J(W,QL-13%2)’
(3.7)
is The LIAR” ’ d e t ermined by ol+i and (3.1). We take ~~~~“(0) as an initial value, then solve (2.4), and compute J(or+i). The largest one among J(al+i), J(ol), J(cr_i), J(al_2) can be omitted. So the new three points and their values are obtained again. We can stop when (3.2) holds or do it some times. We take the largest time equal to five in this paper. If the quadrative interpolation polynomial $(a) does not have minimum point or just J(al+i) J(QI), J(w-l), its zero point
is the largest one, then we take two small values among polynomial and take J(w-2) t o f orm a linear interpolation as QL+~ in order to accelerate convergence. We take the
a0 = O,cyl = -0.8(lu (i)(O) - u(~)(T(~))~~/~I(w(~)II, a2 = 0.5~1 in this paper. The numerical method, for example, the implicit Runge-Kutta method, is used for solving the initial value problem in ordinary differential equations (2.4). The discrete step of t is denoted by h. When uj+i is yielded from uj, we have to check that uj+] is across the disc on the normal plane. When u.j+i is across it, the intersect and the time t, which are denoted by u(Tci)), T(‘), can be defined by using linear interpolation. The coefficient matrix of adjoint linear differential equations (2.11) is obtained by using interpolation too. The steps of t for solving numerical solution of (2.11) are taken to coincide with ones of (2.4) and the computational procedure is taken backward pattern about t.
Periodic
Solutions
169
of ODES
TABLE ZJ(U(O)(O))
u(O) (0)
4.
Time
1.
2J(d”)(O))
T
u(i) (0) 2.806571,2.403286
2.0,2.0
4.013856e-01
5
1.480771e-11
7.1565
2.0,4.0
9.026728e-02
3
3.184155e-10
7.1565
1.684179,3.715761
2.0,6.0
1.733149e+OO
3
8.380585e-10
7.1565
0.646447,4.721644
4.0,2.0
2.397104e-01
5
2.948539e-10
7.1566
3.472134,1.551305
4.0,4.0
1.887851e+OO
3
1.659828e-10
7.1565
2.585821,2.667408
4.0,6.0
5.283966e+OO
3
3.437961e-10
7.1565
1.659904,3.743479
NUMERICAL EXAMPLE
1.
PRACTICE
The famous dx dt
=
Brusslator A
-
model in dissipation
is
(B + 1)~ + 2’1~
dy - = Bx - x’y, dt
(4.1)
where A, B are positive coefficients. When B > 1 + A2 there is a limit cycle or periodic solution [3]. We take A = 1, B = 3. Since n = 2, the normal plane degenerates into a straight line, which is the normal line. The procedure for solving adjoint differential equations can be omitted. The change of sign in the inner product (n, u(t) - u(0)) can judge if the curve is across the normal segment. So the Tci) can be computed. By using the presented method and second-order implicit Runge-Kutta method for solving (2.4) we obtain the results in Table 1 for some different initial values. All real numbers are taken single length arithmetic here. In Table 1 the third column, namely time, indicates the number of times to solve the initial problem (2.4). The last column U(~)(O) is the final approximation of u(0). Table accuracy the limit inside of
1 shows that the convergence rate is very fast and the demand on of initial values is very little. According to computational result cycle is drawn in Figure 1. There is the singular point (1, 3) in the limit cycle,
EXAMPLE 2. B. P. Belousov and A. M. Zhabotinskii noted that different organic acids and metal ions could be used in reaction and still preserve oscillation. Under some assumptions R. J. Field and R. M. Noyes [l] have given a simple model which mechanism retains quantitative features.
170
C. J. JIANG
FIG. 1. The limit cycle.
By
using some transformation
Zhabotinskii
reaction
of variable
their
model
of the Belousov-
is dx - =
s(y -
&
i(fz - y-
dt
xy +
- = dt s dz - = w(x - z),
5
-
4x2) (4.2)
xy)
dt
where s = 77.27, q = 8.375e - 6, w = f = 1.0. It is an ordinary differential equation without initial conditions. S. P. Hastings and J. D. Murray showed that the model equations have an oscillatory or periodic solution [4]. We use the present method to solve the periodic solution of this problem. Since n = 3, we have to solve the adjoint differential equations (2.11). The stiffness of the Field-Noyes model of the Belousov-Zhabotinskii reaction should be conquered by using self-adaptive fourth-order implicit Rung+ Kutta method with variable step. The idea of a self-adaptive method is described
as below:
with variable
step mentioned
here
Periodic
Solutions
of ODES
171
TABLE
g
w
i
2
CY
0
8.333739e-03 -2.00036le-02 1.87661&-01
2J(u@)(O)) 3.566940e-02
T(i) 231.503
-8,333739e-03 2.000361e-02 -1.876618e-01
1.294955e-05
-1.300868e-05
4.923687e-04
-4.922267e-04
5.186635e-05
-5.319774e-05
1
-0.799808
1.426972e-03
231.830
2
-0.399904
1.284108e-02
231.667
3
-0.999756
2.448174e-07
231.911
4
-0.799235
9.793161e-09
231.902
Suppose the step of t is h at moment of t,, we use h to compute u,+r, from u,, denote u,+r by uk+r. We use h/2 to compute two times, namely from u, to compute u,+r/z and then to compute u,+r from u,+~/Z, denote h’2 The ]]ukyr -u~+r((~ can estimate the truncation error. We u,+r by u,+r. restrict h by h,,,, which equals 0.04 in this paper, and Eis the requirement of accuracy. So the procedure is > E, th’ 1s st e P 1s not successful; we omit it, replace (1) If rlu;/;“1 -4+1/12 h by h/2, and do it again.
(2) If llu%
- ui+r]]2 < E, this step is successful. Let u,+r = uiyr and
(a) if e/10 I llu$?r - u~+r]]2 < E, the next step of t is unchanged. (b) if ll$$‘r - di+J2 it is just taken h,,,
< e/10,
the next step oft is taken as 2h, but if 2h > h,,,.
We can use the extrapolation method to improve the accuracy of u,+r C. W. Gear points out [7] that all results about convergency, stability, and
error estimate of the implicit Runge-Kutta method with fixed step is also effective to the method with variable step if the largest step is restrict. This gives a theoretical basis of self-adaptive method with variable step. The computational procedures are recorded in Table 2 and Table 3 for u(‘)(O) = (2.5, 1.65, 2.59) and (10.0, 50.0, 50.0), the final approximate values of u(0) are (2.4916579, 1.6696053, 2.4023413) and (1.0030284, 331.20989, 7651.67908), respectively. Here all real numbers are doubled precision. According to the computational result the periodic solution curve of the model is drawn in Figure 2, where the y-axis is magnified 50 times.
172
C. J. JIANG
TABLE
g
w
i
3 a
0 8.989289e+OO -2.439098efOZ
2.439098e+OZ 7.601184e+03
5.682839e-03 1.195153e+oo
-2.156741e-04
5,78646848e+07
149.009
-8.989289e+OO
-7,601184e+03
-3.724665efOl
T(i)
ZJ(u’i’(O))
1
-0.800188
2.315698e+06
2
-0.400094
2.083184e+07
199.811
3
-1.000222
1,388678e+03
228.425
4
-0.799972
5.557380efOl
231.230
5
-0.399986
4.999374e+OZ
229.853
6
-0.999953
3.334272e-02
231.900
221.021
-5.898636e-03 3.725250e+Ol -l.O12680e+OO
2.156904e-04
4,859032e-03
-4.961698e-03
1.825308e-01
-1.825280e-01 7
-0.800019
1.333710e-03
231.900
8
-0.400009
1.200338e-02
231.900
9
-1.000024
l.OlZOOZe-11
231.899
0
x-120000 FIG. 2. The periodic
solution
curve.
Periodic
The
Solutions
author
of ODES
thanks
Professor
173
Abdus Salam,
Energy Agency, and UNESCO for hospitality for Theoretical Physics, Trieste.
the International at the International
Atomic
Centre
REFERENCES 1
2 3
4
5
6 7
R. Field and R. M. Noyes, Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction, J. Chem. Phys. 160:18771884 (1974). G. Nicolis and I. Prigogina, Self-organization in nonequilibrium systems, Wiley, New York, 1977. Q. Yuaxun and Z. Xianwu, Qualitative investigation of the differential equations of Brusslator in biochemistry, A Monthly Journal of Science 25 (4):273276 (1980). S. P. Hastings and J. D. Murray, The existence of oscillatory solutions in the Field-Noyes model for the Belousov-Zhabotinski reaction, SIAM. J. Appl. Math. 28 (3):678-688 (1975). Z. Suo-Chun and G. Auchmuty, Variational principles and numerical methods for periodic solutions of autonomous ordinary differential equations, Numer. Calculation and App. Cornput., 3:141-148 (1989). W. De Ren, Solution of Nonlinear Equations and Optimization Method, The People’s Education Publishing House, Beijing, 1979. C. W. Gear, Numerical Initial Value Problems in Ordinary Difjerential Equations, Prentice-Hall, 1971.