Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

8.3. NONLINEAR ODES: APPROXIMATE METHODS 309

saddle point at the origin. The fourth singular point at H =0, S =−9000 plays

no role in the “ﬂow” of the trajectories. The phase-plane portrait doesn’t show

the times involved, but of course the other plot on the rhs of Figure 8.7 does.

If every squid was removed from the area occupied by the herring and

from all surrounding areas, H would asymptotically approach an upper limit of

110,000 herring (the H value of the saddle point to which the trajectory would

be attracted). To answer the last part of part (b), we numerically solve the

system of NLODEs with the initial condition H(0) = 110, 000, S(0) = 2, giving

the output as a listprocedure.

sol2:=dsolve({sys,H(0)=110000,S(0)=2},{H(t),S(t)},

type=numeric,output=listprocedure):

The number of herring and squid at an arbitrary time t years later is now

determined using sol2.

H:=eval(H(t),sol2): S:=eval(S(t),sol2):

The number (per acre of seabed) of herring and squid two years later is now

calculated, and (using the round command) found to be about 108,794 and 26.

H2:=round(H(2)); S2:=round(S(2));

H2 := 108794 S2 := 26

8.3 Nonlinear ODEs: Approximate Methods

The perturbation and Krylov–Bogoliubov approximation methods, which may

be applied when the nonlinear terms in the NLODEs are small, are illustrated,

along with the Ritz trial function method which can be used when they are not.

8.3.1 Poisson’s Method Isn’t Fishy

Let a man get up and say, Behold, this is the truth, and instantly

I perceive a sandy cat ﬁlching a piece of ﬁsh in the background.

Look, you have forgotten the cat, I say.

Virginia Woolf, British novelist, (1882–1941)

Perturbation methods are applicable to NLODEs when the nonlinear terms

are small. As a simple example, consider the motion of a unit mass moving in

a viscous medium characterized by a Hooke’s law restoring force F

Hooke

= −y

and a drag force F

drag

= −av− v

= −a (dy/dt) −  (dy/dt)

where y is the

displacement from equilibrium, v is the velocity, a is a positive coeﬃcient, and

 is a small positive parameter. The governing NLODE then is

¨y + a ˙y +  ˙y

+ y =0. (8.13)

For  = 0, (8.13) reduces to the linearly damped simple harmonic oscillator

equation. Our goal is to derive an approximate analytic solution when  is

small, but not zero, for the initial condition y(0)= 1, ˙y(0)= 0. Two values of a

will be considered, ﬁrst a =1/2 which is entered below, and later, a =0. The

Poisson perturbation method is to assume a series solution in powers of , viz.,

310 CHAPTER 8. NLODES & PDES OF PHYSICS

y(t)=y

(t)+y

(t)+

(t)+···+ 

(t). In our recipe, let’s take N =2,

i.e., we will stop at second order in  in our expansion.

restart: with(plots): N:=2: a:=1/2:

A functional operator is introduced for generating the governing NLODE.

ode:=y->diff(y(t),t,t)+a*diff(y(t),t)

+epsilon*(diff(y(t),t))ˆ3+y(t)=0;

ode := y → (

y(t)) + a (

y(t)) +  (

y(t))

+ y(t)=0

The Poisson perturbation expansion is now entered out to order N.

Y(t):=sum(y[n](t)*epsilonˆn,n=0..N);

Y (t):=y

(t)+y

(t)  + y

(t) 

The perturbation expansion is automatically substituted into the NLODE in

the following command line and the result expanded. The very lengthy output

is suppressed with a line-ending colon.

ode2:=expand(ode(Y)):

Powers of  are now collected in ode2 , all powers of order 

(N+1)

and higher

being set equal to zero. The result is then simpliﬁed.

eq1:=simplify(collect(ode2,epsilon),{epsilonˆ(N+1)=0}):

Since  is arbitrary (but small), the coeﬃcient of each power of  must be equal

to zero. This will generate a linear ODE for each order of the perturbation

expansion. These equations are produced in the following do loop.

fornfrom0toNdo

Eq[n]:=coeff(lhs(eq1),epsilon,n)=0;

end do;

:= (

(t)) +

(

(t)) + y

(t)=0

:= (

(t)) + (

(t))

+ y

(t)+

(

(t)) = 0

:= (

(t)) + 3 (

(t))

(

(t)) +

(

(t)) + y

(t)=0

The zeroth order equation, Eq

, is just the linearly damped simple harmonic

oscillator equation. It is now solved subject to the initial condition y

(0) = 1,

˙y

(0) = 0, and the solution is assigned.

sol[0]:=dsolve({Eq[0],y[0](0)=1,D(y[0])(0)=0},y[0](t));

assign(sol[0]):

sol

:= y

(t)=

√

15 e

(−

)

sin(

√

15 t

)+e

(−

)

cos(

√

15 t

)

Since the zeroth-order solution satisﬁes the initial condition, all higher-order

equations must be solved subject to y

(0) = ˙y

(0) = 0, where n =1, ···,N.

8.3. NONLINEAR ODES: APPROXIMATE METHODS 311

This is done in the following do loop, using the Laplace transform option of the

dsolve command. Each solution is assigned and the do loop ended.

fornfrom1toNdo

sol[n]:=dsolve({Eq[n],y[n](0)=0,D(y[n])(0)=0},y[n](t),

method=laplace); assign(sol[n]):

end do:

The lengthy perturbation expansion is now written out

in Y ,

Y:=sum(’epsilonˆn*rhs(sol[n])’,’n’=0..N):

and  and exponential terms sequentially collected.

Y:=collect(Y,[epsilon,exp]);

Y := (

sin(

√

15 t

)

√

15+cos(

√

15 t

)) e

(−

)

+ ((−

305

cos(

√

15 t

116

13725

sin(

√

15 t

)

√

15+

cos(

√

15 t

)

sin(

√

15 t

)

√

15)e

(−

3 t

)

+(−

cos(

√

15 t

305

sin(

√

15 t

)

√

15) e

(−

)

+

((

127494

161101

cos(

√

15 t

) −

10818

805505

sin(

√

15 t

)

√

15) e

(−

)

−(

108

cos(

√

15 t

1525

sin(

√

15 t

)

√

15 +

1525

sin(

√

15 t

)

√

15) e

(−

3 t

)

+(−

648

211975

cos(

√

15 t

) −

1544

3179625

√

15sin(

√

15 t

) −

486

28975

cos(

√

15 t

)

2306

86925

sin(

√

15 t

)

√

15+

28944

28975

cos(

√

15 t

6896

86925

sin(

√

15 t

)

√

15)e

(−

5 t

)

Clearly, deriving the above solution by hand would be a tedious task. Let’s now

take a speciﬁc value for the parameter, say =0.4.

epsilon:=0.4:

The perturbation solution Y is plotted over the time interval t =0 to 20 in gr1.

gr1:=plot(Y,t=0..20,labels=["t","y"],thickness=2):

For comparison, the NLODE is now solved numerically, subject to the given

initial condition,

X:=dsolve({ode(x),x(0)=1,D(x)(0)=0},x(t),numeric,

output=listprocedure):

and plotted in gr2. A point style is chosen, with ﬁfty size 12 black circles.

gr2:=odeplot(X,[t,x(t)],0..20,style=point,symbol=circle,

symbolsize=12,numpoints=50,color=black):

The two graphs are superimposed with the display command,

display([gr1,gr2]);

Note the enclosure of the summand and summation index in “right quotes” here. This is

required for the sum to be carried out. See Maple’s Help on sum for a discussion of this issue.

312 CHAPTER 8. NLODES & PDES OF PHYSICS

–6

Figure 8.8: Numerical (circles) and perturbation (line) solutions for a =1/2

(left) and a=0 (right).

the resulting picture being shown on the left of Figure 8.8. In this case, the

perturbation expansion ﬁts the numerical solution very well.

However, if the recipe is re-executed with a = 0, trouble occurs. Y then

takes the following form, containing secular terms involving powers of t.

Y = cos(t)+ (

sin(3 t)+

sin(t) −

t cos(t))

+ 

(−

t sin(t)+

128

cos(t) −

256

t sin(3 t) −

1024

cos(5 t)+

1024

cos(t)).

Referring to the right hand side picture in Figure 8.8, the perturbation solution

agrees with the exact numerical result for short times, but eventually the secular

terms cause the perturbation expansion to diverge from the exact solution, the

oscillations in fact increasing with time, rather than decreasing. Going to higher

order in the perturbation expansion will not “cure” the situation, as even higher

powers of t then occur. Their are a variety of methods (see Zwillinger [Zwi89])

for dealing with the presence of secular terms in the perturbation solution of a

NLODE. In the following recipe, we illustrate one method which can be applied

to a NLODE having a periodic solution.

8.3.2 Lindstedt Saves the Day

I like to do all the talking myself.

It saves time, and prevents arguments.

Oscar Wilde, Anglo-Irish playwright, author, (1854–1900)

Consider a spring whose displacement from equilibrium satisﬁes the NLODE,

¨y(t)+y(t)+y(t)

=0, (8.14)

8.3. NONLINEAR ODES: APPROXIMATE METHODS 313

with >0 and small. Our goal is to derive a second-order perturbation solution

to equation (8.14), subject to the initial condition y(0) = 1, ˙y(0) = 0. The

solution then will be compared to the exact numerical solution for  =1/2. If

the Poisson method is applied, secular terms will occur, which would destroy

the expected periodic solution. So what do we do?

For  =0, (8.14) is just the simple harmonic oscillator equation which has the

periodic solution cos(t), with frequency 1. The Lindstedt perturbation method

is to assume not only a perturbation expansion for y, but that the frequency

changes slightly to a new frequency Ω, which can also be represented by a per-

turbation series. Introducing a new time variable T =Ωt, and setting X ≡ y,

equation (8.14) can be rewritten as

Ω

X(T )+X(T )+X(T )

=0. (8.15)

Lindstedt’s procedure is to assume that

X(T )=x

(T )+x

(T )+

(T )+···, and Ω = 1+ω

+

+···. (8.16)

The secular terms are removed in each order by imposing the periodicity condi-

tion X(T +2π)=X(T ), i.e., x

(T +2π)=x

(T ), for n=0, 1, 2, .... The initial

conditions are x

(0)= 1, ˙x

=0, and x

(0)= 0, ˙x

(0)= 0 for n=1, 2, ....

The following recipe carries out Lindstedt’s method for the given NLODE.

The maximum power, N =2, of  to be kept in the expansion is entered. The

following alias command produces the output symbol  when e is entered.

restart: with(plots): N:=2: alias(epsilon=e):

The NLODES (8.14) and (8.15) are entered in ode1a and ode1b, respectively.

ode1a:=diff(y(t),t,t)+y(t)+e*y(t)ˆ5=0;

ode1a := (

y(t)) + y(t)+ y(t)

ode1b:=Omegaˆ2*diff(X(T),T,T)+X(T)+e*X(T)ˆ5=0;

ode1b := Ω

(

X (T )) + X (T )+ X (T )

The perturbation expansions for X(T ) and Ω, given in (8.16), are inputted.

X(T):=sum(’x||n(T)*eˆn’,’n’=0..N);

X (T ):=x0 (T )+x1 (T )  + x2 (T ) 

Omega:=1+sum(’omega||n*eˆn’,’n’=1..N);

Ω:=1+ω1  + ω2 

The above series are automatically substituted into ode1b, which is expanded.

ode2:=expand(ode1b):

Powers of  are collected in ode2 , and terms higher than 

set equal to 0.

ode3:=simplify(collect(ode2,e),{eˆ(N+1)=0}):

A functional operator eq is formed to set the coeﬃcient of the nth power of 

on the lhs of ode3 equal to zero.

eq:=n->coeff(lhs(ode3),e,n)=0:

314 CHAPTER 8. NLODES & PDES OF PHYSICS

Using the functional operator in the following sequence command produces the

linear ODEs eq0 , eq1 ,andeq2 , corresponding to 

, 

,and

, respectively.

eqs:=seq(eq||n=eq(n),n=0..N); assign(eqs);

eqs := eq0 =((

x0 (T )) + x0 (T )=0),

eq1 =(x1 (T )+(

x1 (T )) + x0 (T )

+2ω1(

x0 (T )) = 0),

eq2 =((

x2 (T )) + 5 x0 (T )

x1 (T )+2ω1(

x1 (T )) + x2 (T )

+ ω1

(

x0 (T )) + 2 ω2(

x0 (T )) = 0)

An operator sol is created to analytically solve the nth order ODE, using the

Laplace transform method. The order n and initial value A of x

must be given.

sol:=(n,A)->dsolve({eq||n,x||n(0)=A,D(x||n)(0)=0},x||n(T),

method=laplace):

Setting A=1, the zeroth-order ODE is solved and the solution assigned.

sol0:=sol(0,1); assign(sol0):

sol0 := x0 (T ) = cos(T )

The1stODEissolvedwithA =0, trig terms are combined, and sol1 assigned.

sol1:=combine(sol(1,0),trig); assign(sol1);

sol1 := x1 (T )=

128

cos(3 T ) −

cos(T )+

384

cos(5 T )

−

T sin(T )+T sin(T ) ω1

The secular terms, −

T sin(T )+T sin(T ) ω1, occur and must be removed

from x1 (T ) to maintain periodicity, i.e., to have x1 (T +2π)=x1 (T ). This

implies that we must have ω1=5/16. This frequency is extracted and the

secular term removed from x1 (T ) in the following command lines.

omega||1:=solve(coeff(x||1(T),sin(T))=0,omega||1);

x||1(T):=subs(sin(T)=0,x||1(T));

ω1:=

x1 (T ):=

384

cos(5 T )+

128

cos(3 T ) −

cos(T )

Then, the second-order ODE is solved, and the solution assigned.

sol2:=combine(sol(2,0),trig); assign(sol2):

sol2 := x2 (T )=

98304

cos(9 T )+

294912

cos(7 T ) −

192

cos(3 T )

3791

147456

cos(T )+T sin(T ) ω2+

215

3072

T sin(T )

Secular terms occur in sol2 , which are removed by choosing ω2=−215/3072.

omega||2:=solve(coeff(x||2(T),sin(T))=0,omega||2);

8.3. NONLINEAR ODES: APPROXIMATE METHODS 315

x||2(T):=subs(sin(T)=0,x||2(T));

ω2:=

−215

3072

x2 (T ):=

98304

cos(9 T )+

294912

cos(7 T ) −

192

cos(3 T )+

3791

147456

cos(T )

The perturbation solution is now completely determined to second order in .

Setting =1/2andT =Ωt, the Lindstedt perturbation solution X is displayed.

This would be a tedious result to derive by hand.

e:=1/2: X:=subs(T=Omega*t,x||0(T)+e*x||1(T)+eˆ2*x||2(T));

X :=

581327

589824

cos(

13993 t

12288

768

cos(

69965 t

12288

384

cos(

13993 t

4096

)

393216

cos(

41979 t

4096

1179648

cos(

97951 t

12288

)

AgraphofX,aswellasthe = 0 solution, is created over the time interval

t= 0 to 20.

gr1:=plot({X,cos(t)},t=0..20):

The original NLODE, ode1a, is numerically solved and plotted over the same

time interval, size 12 black circles being used for the numerical curve.

Y:=dsolve({ode1a,y(0)=1,D(y)(0)=0},y(t),numeric,

output=listprocedure):

gr2:=odeplot(Y,[t,y(t)],0..20,style=point,symbol=circle,

symbolsize=12,numpoints=120,color=black):

The two graphs are superimposed, the resulting picture being shown in Fig-

ure 8.9. The Lindstedt curve ﬁts the numerical data quite well, and diﬀers

appreciably from the =0 solution.

display({gr1,gr2});

–1

Figure 8.9: Lindstedt curve lies on numerical circles. Other curve: =0 solution.

316 CHAPTER 8. NLODES & PDES OF PHYSICS

8.3.3 Krylov–Bogoliubov Have A Say

He who says he hates every kind of ﬂattery, and says it in earnest,

certainly does not yet know every kind of ﬂattery.

G. C. Lichtenberg, German physicist and philosopher, (1742–99)

For a =0, equation (8.13) reduces to

¨y + y +  ˙y

=0, (8.17)

the initial condition being y(0)= 1, ˙y(0)= 0. For  small, we saw that the Poisson

perturbation expansion broke down for this NLODE, generating secular terms

which grew with time. Since the solution of this equation is not periodic, the

Lindstedt procedure cannot be used to remove the secular terms.

A procedure developed by Krylov and Bogoliubov [KB43] can be employed

to solve NLODEs of the general structure

¨x + ω

x + f(x, ˙x)=0, (8.18)

with  a small parameter. Equation (8.17) is of this structure with x = y,

= 1 and f =˙x

.For = 0, a general solution of (8.18) would be of the form

x=a sin(ω

t + φ), with a and φ constants determined by the initial conditions.

If x is the displacement, the velocity would be ˙x = aω

cos(ω

t + φ). In the

Krylov–Bogoliubov (KB) method

, exactly the same forms are assumed for x

and ˙x, but the amplitude a and phase factor φ are allowed to become time-

dependent. Provided that a(t)andφ(t) are “slowly varying”, their structures

are determined by solving [EM00] [Zwi89],

˙a = −



2πω



2π

f(a sin ψ, aω

cos ψ)cosψdψ,

φ =



2πaω



2π

f(a sin ψ, aω

cos ψ)sinψdψ,

(8.19)

where ψ ≡ ω

t + φ. The criterion for being slowly varying is that



˙a



2 π

 1, and



 1. (8.20)

Recalling that ω

= 1, for equation (8.17) we have f =(a cos ψ)

. Since the

relevant integrand cos

ψ sin ψ is an odd function of ψ,

φ =0,soφ is a constant.

We need only determine the form of a(t) and make sure that the solution

x = a(t)sin(t + φ) satisﬁes the initial condition. We shall now do this for

the same initial condition as in Recipe 08-3-1, taking x(0) = 1 and ˙x(0) = 0.

The KB solution will be compared to the numerical solution for the same value

of  as in that recipe, viz., =0.4.

Use of the following alias command allows us to enter e in the input and

generate  in the output.

Also referred to as the method of averaging [Zwi89].

8.3. NONLINEAR ODES: APPROXIMATE METHODS 317

restart: with(plots): alias(epsilon=e):

The form of f isenteredintermsofψ,

f:=(a(t)*cos(psi))ˆ3;

f := a(t)

cos(ψ)

and the integration in the ˙a equation in (8.19) carried out.

eq:=diff(a(t),t)=-(e/(2*Pi))*int(f*cos(psi),psi=0..2*Pi);

eq :=

a(t)=−

 a(t)

The resulting ﬁrst-order ODE eq is analytically solved for a(t),

sol:={dsolve(eq,a(t))};

sol := {a(t)=

√

3 t+4 C1

, a(t)=−

√

3 t+4 C1

}

generating two answers with

C1 an arbitrary integration constant. The rhs

of the positive square root solution (the ﬁrst one here) is selected to form the

solution x. The two constants

C1 and φ remain to be determined.

x:=rhs(sol[1])*sin(t+phi);

x :=

2sin(t + φ)

√

3 t+4 C1

Taking =0.4, we set x(0)=1 in initial condition ic1 ,

e:=0.4: ic1:=eval(x,t=0)=1;

ic1 :=

√

4sin(φ)

√

and ˙x(0) = 0 in ic2 .

ic2:=eval(diff(x,t),t=0)=0;

ic2 := −

0.07500000000

√

4sin(φ)

(3/2)

√

4 cos(φ)

√

Then ic1 and ic2 are numerically solved for φ and

C1 . The solution sol2 is

assigned, and the ﬁnal KB solution, x,displayed.

sol2:=fsolve({ic1,ic2},{phi,_C1}); assign(sol2): x:=x;

sol2 := {

C1 =0.9769696007,φ= −11.14792061}

x :=

2sin(t − 11.14792061)

√

1.2 t +3.907878403

The KB solution is plotted over the time interval t= 0 to 20.

gr1:=plot(x,t=0..20,labels=["t","x"]):

The NLODE is entered,

ode:=diff(y(t),t,t)+y(t)+e*diff(y(t),t)ˆ3=0;

ode := (

y(t)) + y(t)+0.4(

y(t))

solved numerically subject to the same initial condition,

318 CHAPTER 8. NLODES & PDES OF PHYSICS

Y:=dsolve({ode,y(0)=1,D(y)(0)=0},y(t),numeric,

output=listprocedure):

and plotted using 50 size 12 black circles to represent the numerical curve.

gr2:=odeplot(Y,[t,y(t)],0..20,style=point,symbol=circle,

symbolsize=12,numpoints=50,color=black):

The KB and numerical solutions are superimposed with the display command,

display([gr1,gr2]);

the resulting picture being shown in Figure 8.10.

–0.6

Figure 8.10: Circles: numerical solution. Curve: KB solution.

The KB solution is in quite good agreement with the exact (numerical) solution.

8.3.4 A Ritzy Approach

We must never confuse elegance with snobbery.

Yves Saint Laurent, French couturier, Ritz, no. 85 (London, 1984)

The perturbation and Krylov–Bogoliubov methods can be applied to a NLODE

when the nonlinear terms are small compared to the linear parts of the equation.

The Ritz method is more general, not being restricted to small nonlinearities.

Consider a general ODE,

f(x, ˙x, ¨x, ..., t)=0, (8.21)

where f is some nonlinear function of its arguments. In the spirit of the

Rayleigh–Ritz method for estimating eigenvalues, we try to guess at a trial

wave function Φ which comes “close” to satisfying (8.21). Since, in general, Φ

will not be an exact solution, substitution of Φ into (8.21) will not generate

0 on the right-hand side, but some time-dependent contribution e(t), referred

to as the residual, i.e., f(Φ,

Φ,

Φ, ... t)=e(t). Paralleling the method of least

squares for ﬁtting experimental data, the Ritz method calculates the total error