Enns R.H. Computer Algebra Recipes for Mathematical Physics

Подождите немного. Документ загружается.

9.1. ORDINARY DIFFERENTIAL EQUATIONS 339

The time at the end of the do loop is recorded and the beginning time subtracted

to give us the elapsed cpu time for the algorithm.

cpu_time:=time()-begin;

cpu

time := 0.980

It took

about 1 s to generate all the data points needed for plotting. A 3-

dimensional plot of t vs. x vs. y is produced, using the sequence of n data

points, with the pointplot3d command. The orientation [-90,0] produces x vs.

t as shown on the left of Figure 9.1. Each cross represents a data point. In the

slowly varying regions, the data points are so thick that they blend into a solid

line. If you wish to see the phase plane portrait, y vs. x, simply rotate the

computer plot by clicking on it and dragging with your mouse.

pointplot3d([seq(pt[j],j=0..n)],axes=normal,symbol=cross,

symbolsize=8,color=red,labels=["t","x","y"],

orientation=[-90,0]);

–2

–1

–2

–1

20 40

Figure 9.1: Euler solution of VdP eq. Left:“First principles”; Right:“Dial-up”.

Maple has built-in options in its dsolve command to implement some common

speciﬁc numerical algorithms, the (forward) Euler method being one of them.

Let’s now illustrate this “dial-up” approach, entering the system of ﬁrst-order

ODEs for the VdP oscillator.

sys:=diff(X(t),t)=Y(t),diff(Y(t),t)=F(X(t),Y(t));

sys :=

X (t)=Y (t),

Y (t)=20(1− X (t)

) Y (t) −X (t)

The dependent variables are speciﬁed along with the initial condition.

vars:={X(t),Y(t)}: ic:=X(0)=0.2,Y(0)=0.1:

All cpu times quoted in this chapter are for a 1 GHz personal computer.

340 CHAPTER 9. NUMERICAL METHODS

The dsolve command is applied to sys, subject to the initial condition, a numer-

ical solution being sought with stepsize h. To obtain the Euler approximation,

rather than Maple’s default scheme, the option method=classical[foreuler]

is speciﬁed. The output is given as a listprocedure, so that the numerical data

can be extracted for comparison with the ﬁrst-principles calculation.

sol:=dsolve({sys,ic},vars,numeric,stepsize=h,

method=classical[foreuler],output=listprocedure):

The numerical solution, X(t)vs. t, is then plotted over the time interval t=0 to

hn= 60, using the odeplot command. The same number of points is taken as

in the ﬁrst-principles calculation. The result is shown on the right of Figure 9.1,

a line style having been selected.

odeplot(sol,[t,X(t)],0..h*n,axes=normal,style=line,

numpoints=n,labels=["t","x"]);

The two pictures appear to be identical, but we can conﬁrm the agreement

quantitatively by considering a particular time, say, t = 1. The solution sol is

used to evaluate X at arbitrary time t in XX.UsingXX, the dial-up value of x at

t = 1 is given in x1 . The ﬁrst-principles answer is found by entering pt[100],

the x value at 0.01 × 100 = 1.00 s being the second entry in the output.

XX:=eval(X(t),sol): x1:=XX(1); pt[100];

x1 := 1.64892899170276830

[1.00, 1.64892899170276830, −0.0478758403582165596]

The dial-up and ﬁrst-principles values of x agree to 18 digits!

Because  was increased from 5 in our earlier recipe (08-2-1)to20here,the

relaxation oscillations in Figure 9.1 are characterized by even longer periods of

relative inaction between the abrupt temporal changes. It should be mentioned

that an example of a relaxation oscillator in nature (although not governed by

the VdP equation) is the famous geyser Old Faithful in Yellowstone Park.

9.1.2 Survival of the Fittest

The price which society pays for the law of competition ... is great.

But, whether the law be benign or not, ... it is here; we cannot evade

it; no substitutes for it have been found; and while the law may be

sometimes hard for the individual, it is best for the race, because it

ensures the survival of the ﬁttest in every department.

Andrew Carnegie, American industrialist, philanthropist, (1835–1919)

In terms of a Taylor series expansion, the Euler method is of O(h) accurate.

This necessitates taking h very small, which leads to a longer computation time

than for an algorithm which is of higher order accuracy in h, thus permitting a

bigger h to be used. In this recipe, the Volterra–Lotka (V-L) equations will be

solved with the modiﬁed Euler method, which is of O(h

) accuracy, as well as

with the less accurate Euler method.

9.1. ORDINARY DIFFERENTIAL EQUATIONS 341

The V-L competition equations are used in mathematical biology to describe

the interaction of biological populations, such as, e.g., rats and ferrets. Letting

r(t)andf(t) be the rat and ferret population numbers per unit area at time t,

the V-L equations are

˙r = R(r, f )=a

r −a

rf,

f = F (r, f )=−b

f + b

rf, (9.7)

where the coeﬃcients a

, a

, b

,andb

are all positive. In (9.7), it is assumed

that the ferrets survive by eating only rats, which munch on garbage. If the

rats and ferrets did not interact (a

= b

= 0), the rat (ferret) equation would

become ˙r = a

r (

f = −b

f), which has an exponentially growing (decaying)

solution for the rat (ferret) number. Inclusion of the nonlinear interaction terms

in each ODE will dramatically change the nature of the solution.

Setting t

k+1

= t

+ h, the modiﬁed Euler algorithm for (9.7) is

k+1

= r

+ h (R

[k]+R

[k])/2,f

k+1

= f

+ h (F

[k]+F

[k])/2, (9.8)

with R

[k] ≡ R(r

), F

[k] ≡ F (r

), R

[k] ≡ R(r

+ hR

[k],f

+ hF

[k]),

and F

[k] ≡ F (r

+ hR

[k],f

+ hF

[k]).

That this scheme is of O(h

) accuracy can be conﬁrmed by expanding (9.8)

in powers of h and showing that the resulting equations agree to order h

with

the Taylor expansions of r

k+1

= r(t

k+1

= t

+ h), f

k+1

= f(t

k+1

As an illustrative example, let’s take the nominal values a

=2, a

=0.02,

=1, b

=0.01, with r

=300 and f

=100 at t

=0.

restart: with(plots):

a1:=2: a2:=0.02: b1:=1: b2:=0.01:

t[0]:=0: r[0]:=300: f[0]:=100:

10 digits accuracy is taken and for the modiﬁed Euler method the stepsize

h=0.06 and n=500 iterations are considered. For the Euler method, a smaller

stepsize H = h/3=0.02 is chosen and a correspondingly larger number N =

3 n =1500 iterations taken to make the total time for both solutions the same.

Digits:=10: h:=0.06: n:=500: H:=h/3: N:=3*n:

Functional operators are formed to calculate R and F , i.e., the rhs of eqs. (9.7).

R:=(r,f)->a1*r-a2*r*f: F:=(r,f)->-b1*f+b2*r*f:

Using spacecurve, an operator gr is formed to plot the numerical solution for

either method as a red curve in 3 dimensions. The numerical points p must be

speciﬁed as well as the upper limit n in the sequence command. The orientation

is taken to be [0,90], so that a phase-plane portrait, f vs. r, results. Of course,

the 3-dimensional plots can be rotated on the computer screen.

gr:=(p,n)->spacecurve([seq(p[k],k=0..n)],color=red,axes=

normal,labels=["t","r","f"],orientation=[0,90]):

The Euler solution will be obtained ﬁrst. The start time for the algorithm is

recorded in begin1,

begin1:=time():

342 CHAPTER 9. NUMERICAL METHODS

and the Euler method applied to the V-L equations, the stepsize being H and

N the number of iterations .

forkfrom0toNdo

t[k+1]:=t[k]+H;

r[k+1]:=r[k]+H*R(r[k],f[k]);

f[k+1]:=f[k]+H*F(r[k],f[k]);

pt1[k]:=[t[k],r[k],f[k]];

end do:

The cpu time for the above algorithm is determined.

cpu_time1:=time()-begin1;

cpu

time1 := 0.208

It took about 0.2 seconds cpu time to execute the Euler method. Using gr and

the numerical points pt1, the Euler solution of the V-L equations is plotted,

the resulting picture being shown on the left of Figure 9.2.

gr(pt1,N);

100

200

300

0 100 300

500

100

200

100 200

300

Figure 9.2: Left: Euler solution. Right: Modiﬁed Euler solution.

From this picture, one might conclude that the trajectory is an unstable spiral,

unwinding from the starting point r

= 300, f

= 100. Let’s see what happens

when the modiﬁed Euler method is now applied to the V-L equations. Again

the starting and ending times for the algorithm are recorded.

begin2:=time():

forkfrom0tondo

t[k+1]:=t[k]+h;

R1[k]:=R(r[k],f[k]); F1[k]:=F(r[k],f[k]);

R2[k]:=R(r[k]+h*R1[k],f[k]+h*F1[k]);

9.1. ORDINARY DIFFERENTIAL EQUATIONS 343

F2[k]:=F(r[k]+h*R1[k],f[k]+h*F1[k]);

r[k+1]:=r[k]+h*(R1[k]+R2[k])/2;

f[k+1]:=f[k]+h*(F1[k]+F2[k])/2;

pt2[k]:=[t[k],r[k],f[k]];

end do:

cpu_time2:=time()-begin2;

cpu

time2 := 0.164

Because a larger stepsize was used with fewer iterations, the cpu time for the

modiﬁed Euler method is slightly less here than for the Euler algorithm. Again

using gr, and considering the numerical points pt2, the resulting trajectory is

asshownontherightofFigure9.2.

gr(pt2,n);

Because the cumulative numerical error for the modiﬁed Euler method is less,

a more accurate curve is obtained. It appears that the correct solution is a

cyclic variation in the rat and ferret population numbers. This is in dramatic

contrast to the situation where there is no interaction. Although the parameter

values were artiﬁcial and there are limitations

to the V-L equations, such cyclic

behavior has been observed in a wide variety of competing biological species.

Since the modiﬁed Euler and Euler algorithms are actually equivalent to

ﬁnite Taylor expansions in h, it would not be surprising that the algorithms

break down for large values of h. What may surprise you is that the breakdown

can occur for values of h considerably less than 1, the breakdown showing up as

a numerical instability with the solution diverging to inﬁnity. Without changing

any of the other parameters you should see what happens in this recipe as you

increase the value of h. You may have to cut down on n to avoid numerical

overﬂow. Numerical instabilities for suﬃciently large h are an inherent feature

of explicit schemes.

9.1.3 A Chemical Reaction

What a man calls his “conscience” is merely the mental action that

follows a sentimental reaction after too much wine or love.

Helen Rowland, American journalist, A Guide to Men, 1922

The Euler and modiﬁed Euler methods are the simplest examples of single-step

explicit methods for numerically solving ordinary diﬀerential equations. They

represent the lowest approximations in terms of accuracy of the Rung´e–Kutta

(RK) methods ([LS71][PFTV90]). Historically, the fourth-order (accurate to or-

der h

) RK method has been the “work horse” of these methods as it combines

reasonable accuracy with reasonable speed.

For example, if there were no interaction, the rat population would not grow indeﬁnitely

if conﬁned to a given area because of limited resources, diseases, use of poisons, etc. Further,

ferrets eat other creatures such as rabbits.

344 CHAPTER 9. NUMERICAL METHODS

Given the general ODE, dx/dt= f(t, x), the 4th-order RK algorithm is

k+1

= x

+[K

+2K

+ K

]/6, (9.9)

with K

=hf(t

), K

=hf(t

+h/2,x

/2), K

=hf(t

+h/2,x

/2),

and K

=hf(t

+h, x

). The algorithm is easily generalized to higher-order

ODEs or systems of ﬁrst-order ODEs.

As a simple example, consider the irreversible chemical reaction

+2H

O+3S → 4KOH+2Cr

+3SO

with initially N

molecules of potassium dichromate (K

), N

molecules

of water (H

O), and N

atoms of sulphur (S). The number x of potassium hy-

droxide (KOH) molecules at time t seconds is given by the rate equation

dx/dt = k(2N

− x)

(2N

− x)

(4N

/3 − x)

with k =1.64 × 10

−20

−1

and x(0) = 0.

Taking N

= 1000, N

= 2000, N

= 5000, stepsize h =0.001 s, and n = 200

iterations, we will solve the rate equation using the fourth-order RK algorithm

and graphically compare the solution with that obtained using Maple’s dsolve

command with the option method=classical[rk4]. We will also determine

the number of KOH molecules at 0.1 s.

After loading the plots package and setting the accuracy at 10 digits,

restart: with(plots): Digits:=10:

an operator F is created for calculating the rhs of the rate equation.

F:=x->k*((2*N1-x)ˆ2)*((2*N2-x)ˆ2)*((4*N3/3-x)ˆ3);

F := x → k (2 N1 − x)

(2 N2 − x)

(

4 N3

− x)

The various parameter values are entered in the next two command lines.

k:=1.64*10ˆ(-20): h:=0.001: n:=200:

x[0]:=0: t[0]:=0: N1:=1000: N2:=2000: N3:=5000:

After recording the start time, the 4th-order RK algorithm is applied to the

rate equation, making use of F . A plotting point, pt,isformedoneachstep.

begin:=time():

forjfrom0tondo

K1[j]:=h*F(x[j]);

K2[j]:=h*F(x[j]+K1[j]/2);

K3[j]:=h*F(x[j]+K2[j]/2);

K4[j]:=h*F(x[j]+K3[j]);

x[j+1]:=x[j]+(K1[j]+2*K2[j]+2*K3[j]+K4[j])/6;

t[j+1]:=t[j]+h;

pt[j]:=[t[j],x[j]];

end do:

9.1. ORDINARY DIFFERENTIAL EQUATIONS 345

The cpu time for executing the algorithm is determined, and a graph of the

numerical points created in gr1. The points are represented by size 12 blue

circles. To avoid a messy overlap, only every third point is plotted.

cpu_time:=time()-begin;

cpu

time := 0.153

gr1:=pointplot([seq(pt[3*j],j=0..n/3)],symbol=circle,

symbolsize=12,color=blue):

To determine the number of KOH molecules at 0.1 s, the plotting point at

n/2 = 100 is determined and the round command applied to the second entry.

Pt:=pt[n/2]; N:=round(Pt[2]);

Pt := [0.100, 1586.746962] N := 1587

At 0.1 s, there are 1587 KOH molecules. Now the rate equation is solved using

the dsolve command. The relevant ODE is entered, again making use of F .

ode:=diff(X(t),t)=F(X(t));

Taking the same stepsize h and initial condition X(0) = 0, ode is numerically

solved for X(t), with the option method=classical[rk4].

sol:=dsolve({ode,X(0)=0},X(t),numeric,stepsize=h,

method=classical[rk4]):

This solution is then plotted in gr2 with odeplot, using a line style.

gr2:=odeplot(sol,[t,X(t)],0..h*n,numpoints=n,style=line):

The two graphs, gr1 and gr2, are then displayed in Figure 9.3, the circles lying

exactly on top of the solid line as expected.

display({gr1,gr2},labels=["t","x"]);

800

1600

0.1

0.2

Figure 9.3: Circles: “First principles” solution. Line: “Dial-up” solution.

346 CHAPTER 9. NUMERICAL METHODS

9.1.4 Parametric Excitation

It is the unknown that excites the ardor of scholars, who,

in the known alone, would shrivel up with boredom.

Wallace Stevens, American poet, (1879–1955)

In numerically solving some nonlinear ODEs, regions of solution space can

be encountered where the solution is changing relatively slowly and a larger

stepsize could be used while maintaining a reasonable accuracy, while in other

regions the solution is rapidly changing and a smaller stepsize should be used.

In such situations, an algorithm which adapts its stepsize according to the so-

lution “terrain” is desirable. One of the most popular adaptive step schemes is

the Rung´e–Kutta–Fehlberg [Feh70] fourth-ﬁfth (RKF 45) method, which uses

both the fourth- and ﬁfth-order RK algorithms and adapts the stepsize accord-

ing to the diﬀerence between the two solutions. Fortunately, one of the options

in the numerical dsolve command is precisely this scheme. To invoke it, simply

include the option method=rkf45.

As an illustrative example of using the rkf45 option, consider the follow-

ing interesting problem. The pivot point for the simple pendulum (a mass m

attached to a light rod of length r making an angle θ with the vertical) is under-

going vertical oscillations given by A sin(ωt). Using the Lagrangian approach,

show that the relevant equation of motion is

θ +[W

− Aω

sin(ωt)/r]sinθ =0,

with W =



g/r,whereg is the gravitational acceleration. This nonlinear

ODE with a time-dependent coeﬃcient is an example of a parametric excitation.

Taking A=1 m, r=9.8m,g =9.8m/s

, θ(0) =2 π/3 rads, and

θ(0)= 0 rads/s,

numerically solve the equation of motion using the RKF 45 method for ω =0.1

and 1.2 rads/s. Plot θ(t) over the time interval t= 100 to 300 s.

The plots and VariationalCalculus packages are loaded,

restart: with(plots): with(VariationalCalculus):

and the horizontal (x) and vertical (y) coordinates of m entered. The coordi-

nates are measured from the equilibrium position of the suspended mass.

x:=r*sin(theta(t)); y:=r*(1-cos(theta(t)))+A*sin(omega*t);

x := r sin(θ(t)) y := r (1 − cos(θ(t))) + A sin(ωt)

The kinetic energy T =

m(˙x

+˙y

) is calculated and simpliﬁed.

T:=simplify((m/2)*(diff(x,t)ˆ2+diff(y,t)ˆ2));

T :=

m(r

(

θ(t))

+2r sin(θ(t)) (

θ(t)) A cos(ωt) ω + A

cos(ωt)

)

The potential energy V =mgy is entered and the Lagrangian L=T −V formed.

V:=m*g*y: L:=T-V;

9.1. ORDINARY DIFFERENTIAL EQUATIONS 347

L :=

m(r

(

θ(t))

+2r sin(θ(t)) (

θ(t)) A cos(ωt) ω

+ A

cos(ωt)

) − mg(r (1 − cos(θ(t))) + A sin(ωt))

The Euler–Lagrange command is applied to L, simpliﬁed, and equated to 0.

eq:=simplify(EulerLagrange(L,t,theta(t))[1])=0;

eq := −mr(g sin(θ(t)) + r (

θ(t)) − sin(θ(t)) A sin(ωt) ω

)=0

The output in eq is the equation of motion. Because L contained an explicit

time-dependence, no ﬁrst integral was produced here. Multiplying eq by −1

and dividing by mr

produces the result in eq2 .

eq2:=expand(-eq/(m*rˆ2));

eq2 :=

g sin(θ(t))

θ(t)) −

sin(θ(t)) A sin(ωt) ω

Next, we substitute g =rW

into eq2 .

eq3:=subs(g=r*Wˆ2,eq2);

and obtain the desired NLODE by collecting the sin(θ(t)) terms in eq3 .

eq4:=collect(eq3,sin(theta(t)));

eq4 := (W

−

A sin(ωt) ω

)sin(θ(t)) + (

θ(t)) = 0

The given parameter values are entered and the frequency W =



g/r calculated.

A:=1: r:=9.8: g:=9.8: W:=sqrt(g/r);

W := 1.000000000

The initial condition is entered.

ic:=theta(0)=2*Pi/3,D(theta)(0)=0;

ic := θ(0) =

2 π

, D(θ)(0) = 0

Applying the unapply command to eq4 turns the ODE into a functional oper-

ator, where the frequency ω must be supplied.

ode:=unapply(eq4,omega):

An operator sol is formed for numerically solving the ODE for θ(t), using the

RKF 45 method, for a given ω and subject to the initial condition.

sol:=omega->dsolve({ode(omega),ic},theta(t),numeric,

method=rkf45):

An operator p is created to plot the numerical solution (θ(t)vs. t) for a speciﬁed

value of ω. The time range is taken from t = 100 to 300, the lower limit being

such as to eliminate the transient. To avoid an incorrect set of ﬁgures, 3000

plotting points are chosen.

p:=omega->odeplot(sol(omega),[[t,theta(t)]],100..300,

numpoints=3000,thickness=2):

348 CHAPTER 9. NUMERICAL METHODS

Entering p(0.1) and p(1.2) yields the pictures on the left and right of Fig. 9.4.

p(0.1); p(1.2);

–2

–1

theta

150

200

250

300

–20

theta

150 200 250 300

Figure 9.4: RKF 45 solution for ω =0.1 (left) and ω =1.2 (right).

For ω =0.1, a periodic response of the pendulum is obtained, while for ω =1.2

the solution appears to be highly irregular or chaotic. Other values of ω yield

diﬀerent results.

9.1.5 A Stiﬀ System

A stiﬀ apology is a second insult ....

The injured party does not want to be compensated because he has

been wronged; he wants to be healed because he has been hurt.

G. K. Chesterton, British author, (1874–1936)

A stiﬀ ODE, or ODE system, is one for which there are two or more very

diﬀerent time or spatial scales for the independent variable. The shortest char-

acteristic time τ (or distance) acts as an approximate boundary between nu-

merical stability and instability for any explicit numerical scheme with a ﬁxed

stepsize h.Whenh is less than τ , the scheme is stable, but as h is increased

above τ the scheme becomes increasingly unstable, with ﬂoating point overﬂow

eventually occuring. Numerical algorithms having an adaptive stepsize, such

as RKF45, circumvent this problem by adjusting the stepsize appropriately. In

this recipe, we shall ﬂesh these ideas out by numerically solving a stiﬀ system

of linear ODEs with both the fourth-order RK and the RKF 45 algorithms.

This system will also be solved analytically so as to illustrate the existence of

two widely diﬀerent characteristic times and to provide us with a check on the

accuracy of our numerical algorithms.