Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

1.1. LINEAR ODES WITH CONSTANT COEFFICIENTS 27

X := 1 +

√

2sin(2

√

2 t)+

sin(2 t)+

Y := 2 −

√

2sin(2

√

2 t)+

Z := 3 +

√

2sin(2

√

2 t) −

sin(2 t)+

The positions of the three coaches are now plotted in Figure 1.4 over the interval

t= 0 to 20, the linestyles 1, 2, 3 generating solid, dotted, and dashed curves.

plot([X,Y,Z],t=0..20,color=[red,blue,green],linestyle=

[1,2,3],thickness=2,tickmarks=[3,3],labels=["t","X,Y,Z"]);

X,Y,Z

01020

Figure 1.4: Positions of the three train coaches as a function of time.

The motion of the three coaches (represented by size 20 red boxes) is now

animated and can be observed by executing the following command line.

animate({[X,0],[Y,0],[Z,0]},t=0..100,frames=200,style=

point,symbol=box,symbolsize=20,view=[0..8,-0.1..0.1],

tickmarks=[3,2],labels=["x"," "]);

You can explore other possible motions by changing the parameter values and

initial conditions. The time range and view may have to be altered.

1.1.5 “Resonances”, A Recipe by I. M. Curious

Training is everything. The peach was once a bitter almond;

cauliﬂower is nothing but cabbage with a college education.

Mark Twain, Pudd’nhead Wilson, ch. 5, (1894)

Over the years I have been fortunate to have one or more students in each

mathematical physics class who stood out due to the quality of their work and

their innate curiosity about the subject matter. As a representative of this elite

28 CHAPTER 1. LINEAR ODES OF PHYSICS

group, let me introduce Ms. I. M. Curious who will discuss in her own words the

elegant and thorough computer algebra solution that she generated in solving

the following problem involving a familiar nonhomogeneous LODE.

A mechanical system (of unit mass) experiences a Stokes’s law drag force

Stokes

= −2 α ˙x, a Hooke’s law restoring force F

Hooke

= −ω

x, and a driving

force F

Drive

= f cos(Ω t + δ). Here α is a positive constant, ω the natural

frequency, f the force amplitude, Ω the driving frequency, and δ a phase angle.

(a) Obtain the solution of the resulting LODE for α<ω (underdamping),

subject to the initial condition x(0) = X ,˙x(0) = V . Discuss the method

used by Maple to solve the LODE and the structure of the solution.

(b) Derive an expression for the absolute value A of the steady-state ampli-

tude as a function of Ω. Determine the resonant frequency Ω

at which

A has its maximum value A

max

and evaluate A

max

x(t)andv(t)= ˙x(t) and plot the results. Approximately how long does

it take for steady-state to prevail? Evaluate Ω

and plot the resonance

curves A vs. Ω for some diﬀerent f values.

And here is I. M. to discuss her solution.

“After loading the plots package, infolevel[dsolve]:=5 is entered so as

to provide maximum information on Maple’s ODE solving method.

restart: with(plots): infolevel[dsolve]:=5:

To save on later typing, I introduce the following alias command. When a,

w, W, d, Wr,andAm are subsequently entered, the symbols α, ω,Ω,δ,Ω

,and

max

will be respectively produced in the output.

alias(alpha=a,omega=w,Omega=W,delta=d,Omega[r]=Wr,A[max]=Am);

α, ω, Ω,δ,Ω

max

The Stokes, Hooke’s law, and driving forces are entered,

F[Stokes]:=-2*a*diff(x(t),t); F[Hooke]:=-wˆ2*x(t);

F[Drive]:=f*cos(W*t+d);

Stokes

:= −2 α (

x(t)) F

Hooke

:= −ω

x(t) F

Drive

:= f cos(Ω t + δ)

which are automatically substituted into Newton’s second law in ode .

ode:=diff(x(t),t,t)-F[Stokes]-F[Hooke]=F[Drive];

ode := (

x(t)) + 2 α (

x(t)) + ω

x(t)=f cos(Ω t + δ)

Assuming that α>0, ω>0, and α<ω, ode is solved for x(t) subject to the

initial conditions and the rhs of the result taken. The list of mathematical

methods tried are given in the output, along with the solution labeled x1 .

x1:=rhs(dsolve({ode,x(0)=X,D(x)(0)=V},x(t)))

assuming(a>0,w>0,a<w);

Methods for second order ODEs:

1.1. LINEAR ODES WITH CONSTANT COEFFICIENTS 29

— Trying classiﬁcation methods —

trying a quadrature

trying high order exact linear fully integrable

trying diﬀerential order: 2; linear nonhomogeneous with symmetry [0,1]

trying a double symmetry of the form [xi=0, eta=F(x)]

Try solving ﬁrst the homogeneous part of the ODE

−> Tackling the linear ODE “as given”:

checking if the LODE has constant coeﬃcients

<− constant coeﬃcients successful

<− successful solving of the linear ODE “as given”

−> Determining now a particular solution to the nonhomogeneous ODE

building a particular solution using variation of parameters

<− solving ﬁrst the homogeneous part of the ODE successful

x1 := e

(−αt)

sin(

√

−α

t)(−2 Vω

Ω

+Ω

αX + V Ω

− Ω

f sin(δ)

+4V Ω

− Ω

fαcos(δ) − αω

f cos(δ) − 2Ω

αω

+Ωf sin(δ) ω

+4Ω

X − 2Ωf sin(δ) α

+ αω

X + Vω

)

√

− α

(ω

− 2 ω

Ω

+4Ω

+Ω

))

+ e

(−αt)

cos(

√

−α

t)(f cos(δ)Ω

− 2 fαΩsin(δ) − f cos(δ) ω

+Xω

− 2 Xω

Ω

+4X Ω

+ X Ω

)/(ω

− 2 ω

Ω

+4Ω

+Ω

)

((ω

− Ω

) cos(Ω t + δ)+2α Ωsin(Ωt + δ)) f

Ω

+(−2 ω

+4α

)Ω

+ ω

Maple has successfully solved the ODE, a task I would not have wished to do

by hand. It has recognized that ode is a second order nonhomogeneous LODE

with the corresponding homogeneous equation having constant coeﬃcients. It

has solved ode using the method of variation of parameters. The solution x1

is made up of two parts, the transient terms involving the exponentials which

vanish (since α>0) as t →∞, and the remaining steady-state contribution

which survives in this limit. Note that the steady-state part doesn’t depend on

the initial conditions. I will now extract the steady-state solution xss from x1

by removing the terms in x1 which have e

−αt

xss:=remove(has,x1,exp(-a*t));

xss :=

((ω

− Ω

) cos(Ω t + δ)+2α Ωsin(Ωt + δ)) f

Ω

+(−2 ω

+4α

)Ω

+ ω

The amplitude A of the steady-state solution can be obtained as follows. The

amplitude corresponds to maximum displacement of the oscillator at which time

the velocity is zero. I now solve for a time T at which this occurs.

T:=solve(diff(xss,t)=0,t);

T :=

−δ + arctan(

2 α Ω

− Ω

)

Ω

A is then obtained by evaluating xss at t =T and simplifying with the symbolic

30 CHAPTER 1. LINEAR ODES OF PHYSICS

option which picks a particular solution branch of the square root, the sign of

which can be either plus or minus. The magnitude of A then follows on taking

the absolute value of the result. The answer agrees with that found in Fowles

and Cassiday[FC99].

A:=simplify(eval(xss,t=T),symbolic): Amag:=abs(A);

Amag :=



√

− 2 ω

Ω

+4Ω

+Ω



To ﬁnd the resonance frequency, Amag is diﬀerentiated with respect to Ω, set

equal to zero, and the result solved for the frequency.

sol:=solve(diff(Amag,W)=0,W);

sol := 0,

√

− 2 α

, −

√

− 2 α

The resonance frequency Ω

must correspond to the non-zero positive square

root solution. Therefore, the second answer is selected in sol. The maximum

amplitude A

max

is then obtained by evaluating Amag at Ω = Ω

and again

simplifying with the symbolic option.

Wr:=sol[2]; Am:=simplify(eval(Amag,W=Wr),symbolic);

Ω

√

− 2 α

max



√

− α



To answer part (c) of the problem, the given parameter values are entered.

X:=-1: V:=2: w:=1: a:=1/5: W:=3: f:=5: d:=0:

The displacement X1 and velocity V1 are calculated, the radical expressions

being simpliﬁed with the radsimp command.

X1:=radsimp(x1); V1:=radsimp(diff(x1,t));

X1 :=

6737

9816

(−

)

sin(

√

6 t

)

√

6 −

159

409

(−

)

cos(

√

6 t

)

−

250

409

cos(3 t)+

818

sin(3 t)

V1 :=

179

9816

(−

)

sin(

√

6 t

)

√

1411

818

(−

)

cos(

√

6 t

)

750

409

sin(3 t)+

225

818

cos(3 t)

It was not speciﬁed how the results are to be plotted, so I will create a 3-

dimensional picture in gr1 of X1 vs. V1 vs. t by using the spacecurve com-

mand. The time range is taken to be t = 0 to 40, 1000 plotting points are

selected, and the plot is enclosed in a viewing box.

gr1:=spacecurve({[t,X1,V1]},t=0..40,numpoints=1000,

thickness=2,axes=box,labels=["t","x","v"]):

So that the starting point in the 3-dimensional trajectory can be easily identi-

1.1. LINEAR ODES WITH CONSTANT COEFFICIENTS 31

ﬁed, the textplot3d command is used in gr2 to place the word start (colored

red) in the vicinity of the starting point.

gr2:=textplot3d([[0,X-0.15,V,"start"]],color=red):

The graphs are superimposed with the display command. The resulting pic-

ture (Figure 1.5) is enclosed in a viewing box which may be rotated by clicking

on the computer plot and dragging with the mouse. In the upper left-hand cor-

ner of the computer screen, a small window indicates the angular coordinates

of the viewing box. Here the orientation of the box is set to be (45

◦

,45

◦

display({gr1,gr2},orientation=[45,45],tickmarks=[3,3,3]);

start

–1

–2

Figure 1.5: Space curve showing X1 vs. V1 vs. t.

The wild gyrations seen in Figure 1.5 for small t are associated with the tran-

sient, while steady-state corresponds to the uniform, cyclic behavior at larger

t. To obtain an estimate of the time beyond which steady-state prevails, the

orientation of the viewing box can be changed to (-90

◦

)soastoviewX1

vs. t. Steady-state prevails for t greater than about 25 time units.

The resonance frequency is then evaluated in W1 and the value A1 of the

absolute amplitude calculated at this frequency.

W1:=evalf(Wr); A1:=evalf(Am);

W1 := 0.9591663048 A1 := 12.75775908

The resonance frequency of about 0.96 is slightly less than the natural frequency

ω = 1. To locate the resonance frequency on the graph of the resonance curves,

I will plot a thick, vertical, dashed, blue line of height A1 at W1 in L.

L:=plot([[W1,0],[W1,A1]],color=blue,thickness=2,linestyle=3):

To free the driving frequency and the force amplitude from their previously

assigned values, they are now unassigned.

unassign(’W’,’f’):

32 CHAPTER 1. LINEAR ODES OF PHYSICS

A functional operator pl is created for plotting Amag vs. Ω for diﬀerent f

values over the frequency range Ω = 0 to 2 ω.

pl:=i->plot(eval(Amag,f=i),W=0..2*w,thickness=2):

Making use of the sequence command, resonance curves are generated for f =

1, 2, 3, 4, 5 and displayed along with L in Figure 1.6. As I expected, the resonance

curves agree with those drawn in standard texts.”

display({L,seq(pl(i),i=1..5)},tickmarks=[2,3],labels=[W,"A"]);

Omega

Figure 1.6: Resonance curves for f =1 (bottom curve) to f =5 (top).

1.1.6 Mr. Dirac’s Famous Function

The man with a new idea is a crank until the idea succeeds.

Mark Twain, Following the Equator, ch. 32, (1897).

Consider a very long uniform horizontal beam glued to an elastic foundation,

the foundation exerting a Hooke’s law restoring force. Under the inﬂuence of

a steady point force of magnitude F exerted at x = 0, the static displacement

ψ(x) of the beam is governed by the nonhomogeneous LODE [Hil57],

+ Kψ= Fδ(x) (1.3)

where Y is Young’s modulus, I the moment of inertia of the beam’s cross-

section about a horizontal axis, K the spring constant, and δ(x)theDirac delta

function. The delta function (δ function) is deﬁned by

δ(x)=0,x=0; δ(x)=∞,x=0;



−α

δ(x) dx =1, (α, β > 0). (1.4)

The δ function should be regarded as the limit of, e.g, a Gaussian function

which is made progressively taller and narrower in such a way that the area

For the Gaussian, δ(x)= lim

A→∞

(A/

√

π) e

−A

. The choice of functions is not unique.

1.1. LINEAR ODES WITH CONSTANT COEFFICIENTS 33

under the curve remains equal to 1. Changing the argument of the δ-function

from x to x − ζ shifts the location of the δ-function from x=0 to x=ζ.

It should be noted that the δ-function has the dimensions x

−1

. It has a

number of properties, the most important of which is the “selection” property.

If F (x) is a smooth function (not another δ-function), then



ζ+β

ζ−α

F (x) δ(x − ζ) dx = F (ζ), (α, β > 0). (1.5)

The Dirac delta function concept may be generalized to higher dimensions, e.g.,

δ( r ) ≡ δ(x) δ(y) δ(z) is a 3-dimensional δ-function located at x =y = z =0.

A solution of a linear ODE which has a delta function appearing in the

nonhomogeneous term is called a Green function. The goal of the following

recipe is to determine the static displacement of the beam due to the point

force, i.e., the Green function, for Equation (1.3) and plot it.

In the ODE (1.3), the moment of inertia I appears. This is a protected

symbol in Maple, representing

√

−1. Since I is a commonly used symbol in

physics, it is often desirable to unprotect it. The following interface command

is used to let j =

√

−1, thus freeing up I.

restart: interface(imaginaryunit=j):

The LODE (1.3) is inputted, the command Dirac(x) being used to enter δ(x).

ode:=Y*I*diff(psi(x),x$4)+K*psi(x)=F*Dirac(x);

ode := YI(

ψ(x)) + Kψ(x)=F Dirac(x)

The ODE ode can be simpliﬁed by introducing the parameter α =(YI/K)

1/4

and letting ψ(x)=(F/K) G(x), where G(x) is deﬁned to be the Green function.

The transformations are substituted into ode , producing ode2 ,

ode2:=subs({Y=alphaˆ4*K/I,psi(x)=F*G(x)/K},ode);

ode2 := α

K (

∂

∂x

(

F G(x)

)) + F G(x)=F Dirac(x)

which is then divided by F in ode3 , and the result expanded.

ode3:=expand(ode2/F);

ode3 := α

(

G(x)) + G(x)=Dirac(x)

For x<0andx>0, the rhs of ode3 is equal to zero. The general solution of

ode3 is obtained, assuming that x<0.

sol:=rhs(dsolve(ode3,G(x))) assuming x<0;

sol :=

C1 e

(

(1/2−1/2 j)

√

2 x

)

C2 e

(

(−1/2−1/2 j)

√

2 x

)

C3 e

(

(−1/2+1/2 j)

√

2 x

)

C4 e

(

(1/2+1/2 j)

√

2 x

)

Clearly, the general solution for x> 0 will be of a similar structure. Let’s look

at the region to the left of the origin, assigning this portion of G the name GL.

Since the beam is very long, little error is committed by taking the beam to be

34 CHAPTER 1. LINEAR ODES OF PHYSICS

of inﬁnite length. Physically, the displacement of the beam from equilibrium,

and therefore GL,shouldgotozeroasx →−∞. To satisfy this asymptotic

boundary condition, only the two exponentials which have (1/2 − 1/2 j)and

(1/2+1/2 j) in the exponents should be kept. The other two exponentials

involving (−1/2 −1/2 j)and(−1/2+1/2 j) diverge as x →−∞. The following

command selects the two desired exponential terms to keep in GL.

GL:=select(has,sol,{1/2-j/2,1/2+j/2});

GL :=

C1 e

(

(1/2−1/2 j)

√

2 x

)

C4 e

(

(1/2+1/2 j)

√

2 x

)

Similarly, applying the condition G→0asx→∞,theselect command is used

to obtain the form GR of the Green function to the right of the origin.

GR:=select(has,sol,{-1/2-j/2,-1/2+j/2});

GR :=

C2 e

(

(−1/2−1/2 j)

√

2 x

)

C3 e

(

(−1/2+1/2 j)

√

2 x

)

The four constants must now be evaluated. Since the displacement is physically

continuous at the origin, the condition GR−GL=0 at x=0 is entered in eq

eq[0]:=simplify(eval(GR-GL,x=0))=0;

:= C2 + C3 − C1 − C4 =0

An operator f is formed to simplify the nth derivative of the diﬀerence GR − GL

evaluated at x= 0. Both n and the diﬀerence s must be given.

f:=(n,s)->simplify(alphaˆn*(eval(diff(GR-GL,x$n),x=0)=s));

f := (n, s) → simplify(α

((

∂

∂x

(GR −GL))

x =0

= s))

For n =1 and 2, s =0, i.e., the ﬁrst and second derivatives of G are continuous at

x= 0. The third derivative (n=3) is discontinuous at x=0, however. Mentally

integrating ode3 from the left to the right of the δ function yields s =1/α

Making use of these results and f, the remaining three conditions are entered.

eq[1]:=sqrt(2)*f(1,0); eq[2]:=f(2,0);

eq[3]:=sqrt(2)*f(3,1/alphaˆ 4);

:= −

C2 − C2 j + C3 j − C3 +

C1 j − C1 − C4 − C4 j =0

:= ( C2 − C3 + C1 − C4 ) j =0

:= − C2 j + C2 + C3 + C3 j + C1 + C1 j − C4 j + C4 =

√

The set of four equations are solved for the four constants,

sol2:=solve({seq(eq[n],n=0..3)},{_C1,_C2,_C3,_C4});

sol2 :=



C1 =

(

−

√

C2 =

(

√

C4 =

(

√

C3 =

(

−

√



1.1. LINEAR ODES WITH CONSTANT COEFFICIENTS 35

and the solution sol2 assigned. The constants are automatically substituted into

GL and GR, which are put into real forms by applying the evalc command.

assign(sol2): GL:=evalc(GL); GR:=evalc(GR);

GL :=

√

2 e

(

√

2 x

2 α

)

cos(

√

2 x

2 α

)

−

√

2 e

(

√

2 x

2 α

)

sin(

√

2 x

2 α

)

GR :=

√

2 e

(−

√

2 x

2 α

)

cos(

√

2 x

2 α

)

√

2 e

(−

√

2 x

2 α

)

sin(

√

2 x

2 α

)

Noting that there is a symmetry in the two forms (replacing x with −x in

GL yields GR), the Green function can be put into a more compact form by

replacing x by its absolute value, |x|,inGR, then factoring, and relabeling the

result as G. The static displacement of the beam is then ψ(x)=(f/K) G(x).

G:=factor(subs(x=abs(x),GR));

G :=

√

2 e

(−1/2

√

2 |x|

)

(cos(

√

2 |x|

)+sin(

√

2 |x|

))

The Green function G is now plotted for the representative value α =1.

plot(eval(G,alpha=1),x=-10..10,thickness=2,

labels=["x","G"]);

0.35

–10

Figure 1.7: Green function G for the inﬁnite beam.

The result is shown in Figure 1.7. As expected, G and therefore the beam

displacement is localized near the origin where the point force is applied.

36 CHAPTER 1. LINEAR ODES OF PHYSICS

1.2 Linear ODEs with Variable Coeﬃcients

In this section, linear ODEs with variable coeﬃcients are considered. Our at-

tention will be focused on the second-order Sturm–Liouville (S-L) equation,



p(x)



− q(x) y + λw(x) y =0, (1.6)

where λ is a real parameter while p(x), q(x), w(x) are real functions and w(x)

is taken to be non-negative over the range x = a to b of interest. The S-L

equation, and its nonhomogeneous counterpart, plays a very important role in

the mathematical analysis of many physical problems, particularly for boundary

valueproblemswherey or its ﬁrst derivative vanish at a and b. Inthiscasethe

y ’sarereferredtoastheeigenfunctions and the λ ’s as the eigenvalues.For

certain choices of p, q, w,andλ, the general S-L equation reduces to speciﬁc

“well-known” LODEs such as Bessel’s equation, Legendre’s equation, Hermite’s

equation, etc. The standard approach to solving these equations is to seek an

inﬁnite power series solution. In some cases (e.g., for Bessel’s equation) the

solution remains as an inﬁnite series, while in other cases the inﬁnite series

have to be truncated (e.g., for Legendre’s and Hermite’s equations) to ﬁnite

polynomials to avoid divergence problems for large x. In either case, these

solutions to special forms of the S-L equation are referred to as special functions

(e.g., Bessel functions, Legendre functions (polynomials), etc.). The detailed

series derivation of these special functions will be postponed until Chapter 2. In

this section, we shall be content to see what Maple reveals about some of these

special functions and how physical problems involving these functions can be

easily dealt with. You will encounter many more members of the S-L “family”

of special functions as you progress through the book. For example, in the

Entrees you will see that special functions play an extremely important role in

the solutions of linear PDEs such as the wave and diﬀusion equations because

these equations can be “separated” into systems of S-L ODEs.

1.2.1 Introducing the Sturm–Liouville Family

My books are water; those of the great geniuses is wine.

Everybody drinks water.

Mark Twain, Mark Twains Notebooks and Journals, Notebook 26, (1887)

Jennifer, a young mathematician at the Metropolis Institute of Technology

(MIT) and a ﬁrm believer in using computer algebra to teach applied math-

ematics to engineers and physicists, has kindly consented to provide us with

some of the computer algebra recipes that she has developed. In the following

recipe, she introduces two of the more “famous” members of the S-L family,

namely the Bessel and Legendre ODEs.

Jennifer creates two functional operators, SL for generating a speciﬁc ODE

from the S-L equation (1.6) on subsequently specifying the forms of p, q, w