Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

1.2. LINEAR ODES WITH VARIABLE COEFFICIENTS 47

The probability distribution will now be explored for a speciﬁc n value, say

n= 25. Let’s check that the right form of

C2 was entered by integrating P (25)

over the entire inﬁnite range of ζ. The answer should be 1, which it is.

Prob:=int(P(25),zeta=-infinity..infinity);

Prob := 1

Classically, the energy conservation statement for the harmonic oscillator is

mω

¯hωζ

= E =(n +

)¯hω,

where v is the speed. At the turning points v =0, so that ζ = ±

√

2 n +1.

A functional operator L is now formed for calculating the magnitude of the

classical turning point for a speciﬁed n value.

L:=n->sqrt(2*n+1):

Vertical lines are plotted at the two turning points along with P (25). A list

of lists is used for each of the vertical line entries. These lines are dashed and

colored blue, while the probability curve is a solid red curve.

plot([[[L(25),0],[L(25),0.24]],[[-L(25),0],[-L(25),0.24]],

P(25)],zeta=-9..9,numpoints=2000,linestyle=[3,3,1],

color=[blue,blue,red],thickness=2,labels=["zeta","P"]);

zeta

0.1

0.2

–8

–6

–4 –2 2 4

Figure 1.10: Probability distribution for n=25.

The resulting picture is shown in Figure 1.10. From the ﬁgure, we can see that

the quantum oscillator has a non-zero probability of being outside the classical

turning points. This probability is now calculated for n = 25. The integral is

multiplied by 2, because there are two “tails” to the quantum distribution.

Prob2:=evalf(2*int(P(25),zeta=L(25)..infinity));

Prob2 := 0.04506702177

48 CHAPTER 1. LINEAR ODES OF PHYSICS

The probability of being outside the classical range is about 0.045. A bizarre

feature that occurs for odd values of n is that the probability of ﬁnding the

particle at x = 0 (the center of the potential well) is 0. You are referred to

Griﬃths for a more thorough discussion of the quantum oscillator. In that

reference, the probability distribution is drawn for n = 100. It is a trivial task

to change n in the above recipe and replot the probability distribution. Try it!

1.2.4 Going Green, the Mathematician’s Way

Colorless green ideas sleep furiously.

Noam Chomsky, American linguistic scholar illustrating that grammatical struc-

ture is independent of meaning, Syntactic Structures, (1957)

Suppose that one desires to solve a nonhomogeneous Sturm–Liouville equation,

L[y(x)] ≡



p(x)



− q(x) y = f(x). (1.11)

subject to S-L boundary conditions at x = a and b>a. The Green function

method [Zwi89] is to ﬁrst solve the corresponding Green function ODE,

L[G(x|z)] = δ(x − z), with a<(x, z) <b, (1.12)

subject to the same boundary conditions, for the Green function G(x|z). Then

the solution to (1.11) is given by

y(x)=



G(x|z) f(z) dz. (1.13)

The Green function approach, which is related ([SR66]) to the variation of pa-

rameters method, can be applied to other linear nonhomogeneous ODEs besides

the S-L equation. More importantly, it can be generalized to handle nonhomo-

geneous PDEs such as the nonhomogeneous wave and diﬀusion equations.

Here’s an introductory Green function method problem that I have often

assigned over the years in my mathematical physics course. The idea for it,

as well as for others, sprang from a table of Green functions that I stumbled

across while scanning through an old mathematical physics text by Margenau

and Murphy [MM57]. Some of these older books on mathematical techniques

for scientists are gold mines of useful information. Now to the problem.

Explicitly determine the Green function corresponding to the ODE,



+ y



− n

/x = f(x), (1.14)

where n is a positive (non-zero) integer and the boundary conditions are that

y(0) must remain ﬁnite and y(1) = 0. If n = 1 and f(x)=−x

for 0 ≤ x ≤ 1/2

and f (x)=−(1 − x)

for 1/2 ≤ x ≤ 1, use this Green function to explicitly

determine y(x) and plot the result.

The solution is as follows. First let’s note that xy



+ y



=(xy



)



,sothat

the ODE is a nonhomogeneous S-L equation with p(x)=x and q(x)=n

/x.

Now the ODE for the Green function is entered.

restart:

1.2. LINEAR ODES WITH VARIABLE COEFFICIENTS 49

ode:=x*diff(G(x),x,x)+diff(G(x),x)-nˆ2*G(x)/x=Dirac(x-z);

ode := x (

G(x)) + (

G(x)) −

G(x)

= Dirac(x − z)

The general solution of the ODE is obtained, assuming that x<z. The general

solution for x>zis of a similar mathematical structure.

sol:=rhs(dsolve(ode,G(x))) assuming x<z;

sol :=

C1 x

(−n)

C2 x

The Green function must satisfy the same boundary conditions as the solution

of the original equation. For x<z,wemusthaveG

remain ﬁnite at x =0.

Since n is a non-zero, positive integer, the term x

−n

must be removed from sol

to form GL.

GL:=remove(has,sol,xˆ(-n));

GL :=

C2 x

In some executions of the recipe, the coeﬃcient C1 appears in the above out-

put instead of

C2. To avoid diﬃculty with the solve command later, let’s

introduce our own coeﬃcient C by using the operand command to extract the

second operand (x

)inGL and multiply it by C.

GL:=C*op(2,GL);

GL := Cx

To form GR for x>z, let’s substitute C1=A and C2=B in sol.

GR:=subs({_C1=A,_C2=B}sol);

GR := Ax

(−n)

+ Bx

To satisfy the bc at x=1, GR is evaluated at that point and set equal to zero.

eq1:=eval(GR,x=1)=0;

eq1 := A + B =0

At x=z, we have the continuity condition, GL= GR .

eq2:=eval(GL=GR,x=z);

eq2 := Cz

= Az

(−n)

+ Bz

There is a discontinuity in slope at x=z,givenbyG



(z)−G



(z)=1/p(z)=1/z.

eq3:=eval(diff((GR-GL),x),x=z)=1/z;

eq3 := −

(−n)

−

The set of three equations (eq1 , eq2 , eq3 ) are solved for the set of three coeﬃ-

cients (A, B, C) and the solution sol2 is assigned.

sol2:=solve({eq1,eq2,eq3},{A,B,C}); assign(sol2):

sol2 :=



B =

2 z

(−n)

,C=

−z

(−n)

+ z

2 z

(−n)

,A= −

2 z

(−n)



On simplifying GL and GR, the Green function is completely determined.

50 CHAPTER 1. LINEAR ODES OF PHYSICS

GL:=simplify(GL); GR:=simplify(GR);

GL :=

− z

(−n)

)

2 n

GR :=

− x

(−n)

)

2 n

Using the piecewise command, the Green function can be written as a single

function G, which displays the general symmetry property G(x|z)=G(z|x).

G:=piecewise(x<=z,GL,x>=z,GR);

G :=

⎧

⎪

⎨

⎪

⎩

− z

(−n)

)

2 n

x ≤ z

− x

(−n)

)

2 n

z ≤ x

To continue with the solution, GL and GR are evaluated at n=1,

GL1:=eval(GL,n=1); GR1:=eval(GR,n=1);

GL1 :=

x (z −

)

GR1 :=

z (x −

)

and the function f(z) entered as the two pieces f1 and f2 .

f1:=-(zˆ2): f2:=-(1-z)ˆ2:

Now the integration in (1.13) has to be carried out. The integration is slightly

tricky because both G and f are piecewise functions. For x<1/2, we have

y1=



f1 dz +



1/2

f1 dz +



1/2

f2 dz.

This integration is now carried out and the result simpliﬁed.

y1:=simplify(int(f1*GR1,z=0..x)+int(f1*GL1,z=x..1/2)

+int(f2*GL1,z=1/2..1));

y1 :=

x (−6 x

− 13 + 24 ln(2))

Similarly, the integration for x ≥ 1/2 is performed and simpliﬁed.

y2:=simplify(int(f1*GR1,z=0..1/2)+int(f2*GR1,z=1/2..x)

+int(f2*GL1,z=x..1)) assuming x>0;

y2 := −

1+25x

+6x

− 32 x

+24x

ln(x)

The complete solution y to the given nonhomogeneous ODE is pieced together,

y:=piecewise(x<=1/2,y1,x>=1/2,y2);

y :=

⎧

⎪

⎨

⎪

⎩

x (−6 x

− 13 + 24 ln(2)) x ≤

−

25 x

+1+6x

− 32 x

+24x

ln(x)

≤ x

and plotted over the range x= 0 to 1.

plot(y,x=0..1,thickness=2,labels=["x","y"]);

The resulting proﬁle is shown in Figure 1.11, thus ﬁnishing the problem.

1.2. LINEAR ODES WITH VARIABLE COEFFICIENTS 51

0.005

0.01

0.015

0.02

0.2 0.4

0.6

0.8 1

Figure 1.11: Solution of the given nonhomogeneous S-L equation.

1.2.5 In Search of a More Stable Existence

If you pick up a starving dog and make him prosperous, he will not

bite you. This is the principal diﬀerence between a dog and a man.

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

Not all LODEs, even if they are of the S-L type, are readily solved in terms

of either elementary or special functions. In this case, one can easily generate

a numerical solution using the numeric option of the dsolve command. To

illustrate this approach, consider the following mechanics example.

A particle of mass m is fastened to the origin by a spring obeying Hooke’s

law whose spring constant is k(t)=k

sin(ωt). The mass is constrained to

move on a straight line through the origin. A critical frequency ω

exists such

that for ω<ω

the motion is unstable (displays unbounded growth as time

evolves), but is stable (bounded motion) for all frequencies above ω

. Setting

r = ω



m/k

, determine the critical value r

corresponding to ω =ω

and plot

representative solutions for r just below and just above r

. Take the initial

displacement to be x(0)= 1 and the initial velocity to be ˙x(0)= 0.1.

Applying Newton’s second law, the relevant LODE (dots indicating time

derivatives with respect to the time variable τ)is

m ¨x(τ)+k

sin(ωτ) x(τ )=0. (1.15)

Introducing the dimensionless time variable t =(



/m) τ , Equation (1.15)

can be rewritten (dots now indicating time derivatives with respect to t)as

¨x(t)+sin(rt) x(t)=0. (1.16)

After loading the plots package, a functional operator ode is formed for

generating Equation (1.16) for diﬀerent input values of r.

restart: with(plots):

52 CHAPTER 1. LINEAR ODES OF PHYSICS

ode:=r->diff(x(t),t,t)+sin(r*t)*x(t)=0;

ode := r → (

x(t)) + sin(rt) x (t)=0

The initial condition is entered and an unsuccessful attempt (no output is gen-

erated) is made to analytically solve the ODE for arbitrary r.

ic:=x(0)=1,D(x)(0)=0.1:

dsolve({ode(r),ic},x(t));

An arrow operator sol is now created to ﬁnd a numerical solution of the ODE

for a speciﬁed value of r. To achieve this, the option numeric

is included

in the dsolve command line. The output=listprocedure option gives the

output as a list of equations of the form variable=procedure, where the left-

hand sides are the names of the independent variable, the dependent variable(s)

and derivatives, and the right-hand sides are procedures that can be used to

compute the corresponding solution components. This output form is most

useful when the returned procedure is to be used later as will be the case here.

sol:=r->dsolve({ode(r),ic},x(t),numeric,output=listprocedure):

Using sol, a “do loop” is introduced for carrying out the repetitive numerical

calculation of x(t) at a speciﬁc time for systematically increasing values of r.

The general syntax for this common programming structure is:

for name from expression by expression to expression while expression do

statement sequence;

end do:

The statement sequence is the main body of the do loop giving the steps to be

followed in the calculation. In the following opening line of the do loop,

for i from 1 to 200 do

name is the index i, the ﬁrst expression is 1,theby expression is absent, the

third expression is 200, and there is no conditional while expression present.

Since the by expression is missing, the default is to increment i by 1 each

time the statement sequence is executed.

The statement sequence is now explained. Using sol, the following com-

mand line evaluates x(t) at arbitrary time t for r =i/100. Since i runs from

1 to 200, r will be incremented in steps of 0.01 from 0.01 to 2.00. The

concatenation operator || is used in assigning the ith result the name X||i.

This operator attaches the numbers 1 to 200 to X, so that the outputs are

assigned the names X1 ,X2 ,..., X200 . Concatenated names are useful for

employing in the sequence (seq) command.

X||i:=eval(x(t),sol(i/100));

Then each X||i is evaluated at t= 150, the result being labeled Y||i.

Y||i:=X||i(150);

A large time has been chosen so that it will be easy to distinguish between the

unstable solutions which should have grown to large amplitudes and the stable

The default numerical method is the Runge–Kutta–Fehlberg 45 ([BF89]) method (dis-

cussedintheDesserts). Other schemes are listed in the Maple Help under “dsolve,numeric”.

1.2. LINEAR ODES WITH VARIABLE COEFFICIENTS 53

solutions which remain bounded. The do loop is now ended, the output being

suppressed with a line-ending colon.

end do:

To plot the numerical output, 200 plotting points (pts) are formed using the

sequence command. The abscissa of the ith point is given by i/100, with i

running from 1 to 200. This generates the r values. Since the Y ||i for the

unstable solutions are extremely large, the log is taken of the absolute value of

the Y||i to form the ordinate values. For each plotting point, the abscissa and

ordinate values are put into a Maple list format.

pts:=seq([i/100,log(abs(Y||i))],i=1..200):

Choosing a suitable view, the points are now plotted, the default being to join

consecutive points with a straight line. The result is shown in Figure 1.12.

plot([pts],view=[0.01..2,-10..70],tickmarks=[4,3],

labels=["r",""]);

0.5

1.5

Figure 1.12: Log of the absolute value of x(t=150) versus r.

Remembering that the vertical scale is logarithmic, we can see that the region

of instability persists until r reaches slightly less than 1.5, then stability sets in

for larger r.Thatistosayr

 1.5. More precisely, looking at the numerical

output, r

lies between 1.48 (unstable) and 1.49 (stable).

The odeplot command may be used to plot the numerical solution of an

ODE. An arrow operator gr is formed to apply this command to sol(r) for a

given r value. x(t) is plotted against the time t, the time range being from 0 to

150. To obtain an accurate curve, 2000 plotting points are chosen.

gr:=r->odeplot(sol(r),[t,x(t)],0..150,numpoints=2000,

thickness=2,tickmarks=[3,3]):

Then entering gr(1.48) produces the ﬁgure shown on the left of Figure 1.13,

while gr(1.49) generates the ﬁgure on the right.

gr(1.48); gr(1.49);

54 CHAPTER 1. LINEAR ODES OF PHYSICS

–10000

10000

100 150

–10

100

150

Figure 1.13: Left: Unstable solution for r =1.48. Right: Stable for r =1.49.

For r =1.48 ever-growing oscillations occur, indicative of instability, while for

r =1.49 the oscillations remain bounded in amplitude, a sign of stability.

The motion of the particle along the x-axis for, e.g., r =1.49 can now be

animated. The coordinates C (expressed as a list with the vertical coordinate

equal to 0) of the particle at arbitrary time t are evaluated using sol(1.49).

Then, e.g., the coordinates at time t=100, are obtained

by entering C(100).

C:=eval([x(t),0],sol(1.49)): C(100);

[1.52446967145678580, 0]

An operator gr2 is formed to plot a size 20 blue circle at coordinates C(i).

gr2:=i->plot([C(i)],style=point,symbol=circle,

symbolsize=20,color=blue,labels=["x"," "]):

The command seq(gr2(i),i=0..250) plots graphs for t = i =0, 1, 2, ...250.

Displaying this sequence with the option insequence=true generates a picture

which is animated by clicking on the computer plot and then on the start arrow.

display(seq(gr2(i),i=0..250),insequence=true,

tickmarks=[3,3]);

1.3 Supplementary Recipes

01-S01: Newton’s Law of Cooling

Newton’s law of cooling is governed by the ODE,

T =−k (T − Ts), where T is

the temperature of an object at time t, k is the positive cooling constant, and

Ts is the temperature of the object’s surroundings. If T (0) = T 0, determine

T (t>0). What information does including infolevel[dsolve]:=2: provide?

Note that the numeric option of the dsolve command yields more digits in the output

than the “standard” default of 10 digits.

1.3 SUPPLEMENTARY RECIPES 55

A horseshoe is originally 100

◦

C hotter than the surrounding air. After 15

minutes, the temperature diﬀerence has fallen to 60

◦

C. How long will it take

for the horseshoe to reach a temperature 10

◦

C above the surrounding air?

01-S02: Charging a Capacitor

A capacitor C is connected in series with a resistor R and a voltage source

V = V

(t/τ)

−t/τ

,whereτ is a characteristic time.

(a) Making use of Ohm’s law for the voltage drop across R and Kirchoﬀ’s

voltage sum rule, derive the ODE for the charge q(t)onC at time t.

Then, analytically solve the ODE for q(t), given that q(0)=0.

(b) Taking R = 5 ohms, C = 2 farads, V

= 3 volts, and τ = 1 second, plot

q(t)fort = 0 to 20 seconds. At what time is q(t) a maximum and what

is the maximum charge? Answer this ﬁrst qualitatively by clicking the

mouse cursor on the maximum, and then quantitatively.

01-S03: Radioactive Chain

An important radioactive chain involves the disintegration of the unstable

238

nucleus. It decays via α-emission into

234

Th, which in turn β-decays into

234

Pa,

and so on until the stable isotope

206

Pb is created. The decay rates of the

ﬁrst three species (N is the number of atoms and λ the decay constant) in a

radioactive chain at time t can by described by

= −λ

= λ

− λ

= λ

− λ

(a) If N

(0)= N, N

(0)= 0, and N

(0)= 0, determine N

(t) N

(t), and N

(t).

(b) The λ values for the uranium chain are vastly diﬀerent, so for plotting

purposes consider a hypothetical radioactive chain for which N = 1000,

=1, λ

=2, and λ

=3. Plot N

(t) N

(t), and N

(t) in the same ﬁgure

for t =0..5 using diﬀerent colors and line styles for each curve. At what

time is N

a maxima? What is N

at this time?

01-S04: Stokes and Newton Join Forces

The drag force ([FC99]) on a sphere of diameter d moving with speed v is, in

general, given in SI units by F

drag

= −av − bv

,witha =1.55 × 10

−4

d and

b=0.22 d

. I.e., it is a combination of Stokes’s and Newton’s resistance laws

(a) A spherical mass m is dropped from rest. Derive and solve the relevant

nonlinear ODE for v(t). What method does Maple use to solve the ODE?

(b) A basketball (d =0.25, m =0.60), a raindrop (d =10

−4

, m =0.52 × 10

−9

and a soap bubble (d=10

−2

, m =10

−7

), are all dropped from rest. Taking

the gravitational acceleration g =9.8m/s

, determine v(t) for each body.

mine V for each body and the time it takes to come within 1% of V ?

01-S05: Exploring the RLC Series Circuit

(a) Derive the ODE for a series circuit consisting of a resistor R, an inductor

L, and a capacitor C and cast it into the form ¨q +2α ˙q + ω

q =0, where

q(t) is the charge on C at time t, α=R/(2 L)andω =1/

√

LC.

56 CHAPTER 1. LINEAR ODES OF PHYSICS

(b) Given q(0) = A,˙q(0) = B, solve the ODE for overdamping (α>ω),

underdamping (α<ω), and critical damping (α = ω). What fact does

Maple use in successfully arriving at each solution? Calculate the currents.

and the currents for the three cases. Plot the q(t) together using diﬀerent

colors and linestyles for each case. Do the same for the currents.

01-S06: The Whirling Bar Revisited

Suppose that the rotating bar of recipe 01-1-3 has hinged-end (y and y



(the

bending moment) are both zero) boundary conditions at each end. Obtain the

ﬁrst four critical frequencies and corresponding shapes and plot the latter in a

single ﬁgure. Compare the lowest critical frequency to that in 01-1-3.

01-S07: Driven Couple Oscillators

Masses m

and m

, with equilibrium positions at x = 1 and x =5, are free to

move horizontally on a smooth surface (the x-axis). m

is connected to a ﬁxed

wall on its left by a linear spring (spring constant k) and on its right to m

with an identical spring. A driving force F =f sin(ωt)actstotherightonm

Air resistance is present, given by Stokes’s drag law, F

drag

=−av.

(a) Derive the governing ODEs for the displacements x

(t)andx

(t)ofm

and m

from equilibrium. Taking m

=2, m

=1, k =1, a =1/100, ω =2,

f =2, x

(0) = 1/10, x

(0) = 0, ˙x

(0) = 0, and ˙x

(0) = 0, solve the ODEs

using the Laplace transform option. Express x

(t)andx

(t) in real forms.

(b) Extract the steady-state parts of x

and x

and plot them. Discuss the

results. Plot the transient expressions over a time interval for which the

transients become small. This will reveal the envelope of the transients.

and m

about their equilibrium positions over

a time interval for which the transients become small.

01-S08: Some Properties of the Delta Function

(a) Two common representations of the δ function are G(x)=(α/

√

π) e

−α

and F (x)=sin(αx)/(πx)withα →∞. Conﬁrm that these representa-

tions are reasonable by (i) plotting G(x)andF (x) for increasing α,

(ii) showing that



∞

−∞

G(x) dx =1 and



∞

−∞

F (x) dx=1 for α>0.

(b) Using Dirac(x-a), conﬁrm the following: (i)



∞

−∞

δ(x − a) dx =1,

(ii)



∞

−∞

f(x) δ(x − a) dx =f(a), (iii)



∞

−∞

f(x) δ



(x − a) dx=−f



(a),

(iv)



∞

−∞

δ(b(x − a)) dx =1/|b|.

andwithedgesatx = ±1. Plot h(x)forx = −3 to 3, using constrained

scaling. Diﬀerentiate h(x) with respect to x and interpret the result.

01-S09: A Green Function

Derive the Green function G which is the solution to G



+ k

G = −δ(x − ζ),

subject to the boundary conditions G(−1)= G(1) and G



(−1)= G



(1). Simplify