Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

288 CHAPTER 7. CALCULUS OF VARIATIONS

wave equation for small vibrations of a ﬂexible string, ﬁxed at its end points

x= 0 and a, and under tension τ.

07-S12: Betsy Bug’s Ride

Betsy bug is perched on a bead which slides along a smooth wire bent in the

shape of a stylized “W”. I.e., the bead’s vertical coordinate at time t is z(t)=

−c

r(t)

+ c

r(t)

,withc

> 0, c

> 0, and r(t) the radial distance from the

vertical symmetry axis. The wire is rotating about it’s symmetry axis with

angular velocity ω radians/s, so the bead’s horizontal coordinates at time t are

x(t)=r(t) cos(ωt), y(t)=r(t)sin(ωt).

(a) Using the Lagrangian approach, ﬁnd the equation of motion for r(t).

(b) Show that there is critical frequency ω



2 g (2 c

r(0)

− c

), where

g is the acceleration due to gravity, for which the bead will not move

relative to the wire.

=1 m

−1

, c

=1/2m

−3

, gravitational acceleration g =9.8m/s

r(0)= 2 m, ˙r(0)=0 m/s, determine the critical frequency.

(d) Taking ω =0.1 ω

, numerically solve the equation of motion for r(t).

(e) Animate the motion of the bead and the wire.

Part III

THE DESSERTS

The reasonable man adapts himself to the world;

the unreasonable one persists in trying to adapt the

world to himself. Therefore, all progress depends

on the unreasonable man.

George Bernard Shaw, Anglo-Irish playwright, Man and Superman, (1903)

Perfection is achieved not when you have nothing more

to add, but when you have nothing left to take away.

Antoine de Saint-Exupery, French pilot and author, (1900–1944)

I may not have gone where I intended to go,

but I think I have ended up where I intended to be.

Douglas Adams, English humorist and science ﬁction writer, (1952–2001)

Chapter 8

NLODEs & PDEs of Physics

So far, this gourmet selection of computer algebra recipes has emphasized phys-

ical examples governed by linear ODEs and PDEs. However, much of modern

scientiﬁc research involves phenomena described by nonlinear diﬀerential equa-

tions, i.e, equations which are not linear in the dependent variable(s). In the

Desserts, I will give you a small taste of this intellectually delectable area of

mathematical physics. If you enjoy the nonlinear recipes that follow, and crave

more, you should consult [EM00] for a much deeper treatment of the subject.

This chapter illustrates a few of the basic analytic (exact and approximate)

and graphical methods for solving nonlinear ODEs and PDEs (NLODEs and

NLPDEs). Basic numerical methods will be covered in Chapter 9. A compre-

hensive survey of analytic and numerical approaches may be found in Daniel

Zwillinger’s Handbook of Diﬀerential Equations [Zwi89].

8.1 Nonlinear ODEs: Exact Methods

In Chapter 1, we encountered a simple nonlinear ODE, describing the motion

of a falling badminton bird acted upon by a nonlinear drag force, which had an

exact analytic solution. Only a handful of NLODEs of physical interest can be

exactly solved using elementary techniques. Here are a few more examples of

increasing mathematical complexity.

8.1.1 Jacob Bernoulli and the Nonlinear Diode

I recognize the lion by his paw.

Jacob Bernoulli, Swiss mathematician, after reading an anonymous solution to

a problem that he realized was Newton’s solution, (1654–1705)

What does Jacob Bernoulli who lived some 300 years ago have to do with

a modern nonlinear diode? As you will shortly see, the governing circuit equa-

tion in the following recipe is the nonlinear Bernoulli equation which Jacob

discovered in another context in 1690 and solved in 1696. Bernoulli’s equation

292 CHAPTER 8. NLODES & PDES OF PHYSICS

is a ﬁrst order nonlinear ODE of the form



+ f

(x) y = f

(x) y

(8.1)

which can be cast into a linear ODE, and thus solved, by introducing a new

dependent variable p =1/y

n−1

. Bernoulli’s equation is known to Maple and

thus can be easily solved without having to carry out the details of the variable

transformation. Now let’s formulate the nonlinear diode problem.

A linear capacitor C is connected in series to a nonlinear diode, a circuit

element which has a current (i)–voltage (v) relation of the form i = av + bv

where a and b are positive constants and n=2, 3, 4, ···Given that the voltage

across the capacitor at time t=0 is v(0) =V , what is the voltage v(t)fort>0?

Express the solution in dimensionless form in terms of a single parameter β,

which remains to be identiﬁed. Plot the analytic result for n =2to5andβ =2.

To solve this problem, the PDEtools library package is loaded. This contains

the dchange command which will be used to cast the governing circuit equation

and the voltage into dimensionless form.

restart: with(PDEtools):

From the deﬁnition of capacitance, the charge q on the capacitor at time t is

q = Cv(t). So, the current in the circuit is i = −dq/dt = av(t)+bv(t)

These two basic relations are entered, the expression for q being automatically

substituted into i.

q:=C*v(t);

q := C v(t)

i:=-diff(q,t)=a*v(t)+b*v(t)ˆn;

i := −C (

v(t)) = a v(t)+b v(t)

The above ODE is clearly of Bernoulli’s form with y ≡ v, f

≡ a/C and

≡−b/C. It can be transformed into dimensionless form by introducing a

new time variable τ =a t/C and a dimensionless voltage variable y(τ )=v(t)/V .

tr:={t=C*tau/a,v(t)=V*y(tau)};

tr := {t =

Cτ

, v(t)=V y(τ )}

The above transformation is applied to i in the following dchange command.

ode:=dchange(tr,i,[tau,y(tau)]);

ode := −aV (

dτ

y(τ)) = aV y(τ)+b (V y(τ ))

Then ode is divided by −aV and the result expanded.

ode:=expand(ode/(-a*V));

ode :=

dτ

y(τ)=−y(τ) −

b (V y(τ ))

Introducing a dimensionless parameter β through the substitution b = aβ/V

(n−1)

and simplifying with the symbolic option produces the desired dimensionless

nonlinear ODE in ode2 .

8.1. NONLINEAR ODES: EXACT METHODS 293

ode2:=simplify(subs(b=a*beta/Vˆ(n-1),ode),symbolic);

ode2 :=

dτ

y(τ)=−y(τ) − β y(τ)

Note that this equation now only depends on one parameter, β, and that the

initial condition at τ =t=0 is y(τ )=1. Instead of mimicking a hand calculation,

and making the suggested dependent variable substitution, I will directly use

the dsolve command. First let’s enter the following infolevel command to

give us information on what approach Maple takes in solving ode2 .

infolevel[dsolve]:=2:

Then, ode2 is solved for y(τ), subject to the initial condition y(0)= 1.

sol:=simplify(dsolve({ode2,y(0)=1},y(tau)));

Methods for ﬁrst order ODEs:

— Trying classiﬁcation methods —

trying a quadrature

trying 1st order linear

trying Bernoulli

<− Bernoulli successful

sol := y(τ)=(−β + e

((n−1) τ )

+ e

((n−1) τ )

β)

(−

n−1

)

The NLODE is identiﬁed as a Bernoulli equation and successfully solved. To

plot the solution, let’s enter the suggested value β =2.

beta:=2:

The rhs of sol is now plotted for n=2, 3, 4, 5, a diﬀerent color being assigned

to each curve. The resulting graph is shown in Figure 8.1, the bottom curve

corresponding to n = 2, the next highest curve to n = 3, and so on. Because

y(τ) < 1asτ increases, the curves decrease less rapidly as n is increased.

plot([seq(rhs(sol),n=2..5)],tau=0..2,color=[red,green,blue,

black],thickness=2,labels=["tau","y"],tickmarks=[2,2])

0.5

01 2

tau

Figure 8.1: Normalized voltage as a function of time.

294 CHAPTER 8. NLODES & PDES OF PHYSICS

8.1.2 The Chase

I wouldn’t think of asking you to lie;

you haven’t the necessary diplomatic training.

AlinespokenbyoneofthecharactersinthemovieSea Chase (1955)

A classic problem [Dav62] in the history of nonlinear ODEs involves the curve

of pursuit. This is the trajectory generated by a point P which moves in such

a way that its direction of motion is always towards a second point P



, con-

strained to move along a prescribed path. This type of problem originated with

Leonardo da Vinci in the 15th century, but the example in the following recipe

is due to the French hydrographer Pierre Bouguer who published his solution in

1732. Bouguer’s formulation is expressed in terms of one ship pursuing another,

the latter moving along a straight line.

I will phrase the problem somewhat diﬀerently. Referring to the lhs of

Figure 8.2, rascally Roger Rabbit is being pursued by ferocious Freddy Fox.

freddy fox

roger rabbit

capture?

freddy’s

path

(x,y)

(a,Y)

a-x

Y-y

0.5

0.2 0.4

0.8 1

Figure 8.2: Left: Geometry for pursuit. Right: Freddy’s calculated path.

Roger, who is initially at the point x= a, y = 0, runs at a constant speed along

the vertical line x = a. Freddy Fox, who is initially at the origin O, pursues

Roger by constantly aiming at Roger’s current position. Freddy’s constant

speed is n times that of Roger’s, with n>1. The problem is to derive the

equation y(x) of Freddy’s path (indicated schematically in the picture) and

determine the theoretical point of capture. If a = 1 km and n =3/2, will Roger

escape capture by diving into his hole which is located 1195 meters from Roger’s

starting point. Plot Freddy’s path up to the theoretical point of capture. So

let the chase begin!

The assumption that a>0, x ≥ 0, x<a,andn>1 is entered.

restart: assume(a>0,x>=0,x<a,n>1):

From the ﬁgure, the slope of Freddy’s trajectory when he is at the point

(x, y(x)) and Roger is at (a, Y )isgivenbydy/dx =(Y − y(x))/(a − x).

8.1. NONLINEAR ODES: EXACT METHODS 295

slope:=diff(y(x),x)=(Y-y(x))/(a-x);

slope :=

y(x)=

Y − y(x)

a − x

The slope equation is solved for Y ,

Y:=solve(slope,Y);

Y := (

y(x)) a − (

y(x)) x + y(x)

and the x derivative of Y then taken, and factored.

Yder:=factor(diff(Y,x));

Yder := (

y(x)) (a − x)

The arclength along Freddy’s path is ds=



(dx)

+(dy)

and his speed, ds/dt,

is n times Roger’s speed, dY/dt.So,n (dY /dx)=



1+(dy/dx)

is entered.

ode:=n*Yder=sqrt((diff(y(x),x))ˆ2+1);

ode := n (a − x)(

y(x)) =



(

y(x))

Freddy’s path is described by a second order NLODE which, amazingly, can

be solved exactly. Although the dsolve command could be applied directly to

ode , the subsequent manipulations to put the solution in a compact form are

messy. Instead, let’s mimic the approach that we would undertake by hand.

Making the substitution dy(x)/dx = p(x) reduces ode to a ﬁrst order ODE in

p(x) displayed in ode2 .

ode2:=subs(diff(y(x),x)=p(x),ode);

ode2 := n (a − x)(

p(x)) =



p(x)

Freddy’s initial slope is zero, i.e., p(0) = 0. ode2 is now analytically solved for

p(x), subject to the initial slope condition. The right-hand side of the solution

is then taken and the result equated to dy/dx. Occasionally, on executing the

worksheet, the form of the output will diﬀer in ode3 from what is reproduced

here in the text. This will not aﬀect the shape of Freddy’s trajectory or the

theoretical point of capture.

ode3:=diff(y(x),x)=rhs(dsolve({ode2,p(0)=0},p(x)));

ode3 :=

y(x)=−sinh(

ln(a − x) −ln(a)

)

ode3 is analytically solved for y(x), subject to the initial condition y(0) = 0.

Again, the right-hand side of the solution is taken. The equation of Freddy’s

path is then displayed in the output of the following command line.

y:=rhs(dsolve({ode3,y(0)=0},y(x)));

296 CHAPTER 8. NLODES & PDES OF PHYSICS

y := −

nx(a − x)

(

)

2(n +1)a

(

)

an(a − x)

(

)

2(n +1)a

(

)

(−

2(n − 1)

) a

(

)

(a − x)

(

)

−

2(n +1)

2(n − 1)

Next, let’s substitute x= a − X

into y and expand the result.

y:=expand(subs(x=a-Xˆn,y));

y :=

n (X

)

(

)

2(n +1)a

(

)

−

(

)

2(X

)

(

)

(n − 1)

−

2(n +1)

2(n − 1)

The previous substitution is reversed by setting X =(a − x)

1/n

in y. Applying

the combine command yields a compact equation describing Freddy’s trajectory.

y:=combine(subs(X=(a-x)ˆ(1/n),y));

y :=

(a − x)

(

+1)

(−

)

2(n +1)

−

(a − x)

(1 −

)

(

)

2(n − 1)

−

2(n +1)

2(n − 1)

The theoretical point of capture, Yc, follows on evaluating y at x=a.

Yc:=simplify(eval(y,x=a));

Yc :=

− 1

Roger would be captured at Yc =(an)/(n

− 1), unless he gets to his hole

ﬁrst. So, does he avoid capture? Substituting a = 1 and n =3/2intoYc and

multiplying by 1000 to convert the result into meters,

Yc:=evalf(subs({n=3/2,a=1},Yc))*1000;

Yc := 1200.000000

yields a value for Yc of 1200 meters. Fortunately, for Roger, his hole was only

1195 meters from his starting point, so he survives for another day. To plot

Freddy’s path up to the point of theoretical capture, the values n =3/2and

a= 1 are substituted into y, yielding the result shown in y2.

y2:=subs({n=3/2,a=1},y);

y2 :=

3(1− x)

(5/3)

−

3(1− x)

(1/3)

Freddy’s trajectory y2 is then plotted up to the theoretical point of capture,

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

the result being shown on the right-hand side of Figure 8.2. Of course, if Roger’s

hole were not conveniently close, he would have had to take evasive action, lead-

ing to a much more diﬃcult pursuit problem to solve.

8.1. NONLINEAR ODES: EXACT METHODS 297

8.1.3 Not As Hard As It Seems

Often, the less there is to justify a traditional custom,

the harder it is to get rid of it.

Mark Twain, American author, Tom Sawyer, ch. 5, 1876

When a spring is stretched from its equilibrium length by an amount x which

is no longer small, Hooke’s law must be modiﬁed so as to include nonlinear or

anharmonic terms in x. For symmetric oscillations of a nonlinear spring about

equilibrium, the next terms retained in the Taylor expansion of the force must

be cubic in x. In this case, the force F required to stretch the nonlinear spring

will be of the form F =k

x + k

(x + ax

), with k

> 0anda ≡ k

If a>0, one has a hard spring because it is harder to stretch than a linear

spring (a = 0). If a<0, one has a soft spring. Inthisrecipe,wewillsolvethe

equation of motion for the oscillations of a mass m attached to a (light) hard

spring allowed to move horizontally on a smooth horizontal surface. Applying

Newton’s second law, the equation of motion of m at arbitrary time τ is

m ¨x(τ)+k

(x(τ)+ax(τ )

)= 0 (8.2)

or, on introducing a new time variable t=



/m τ ,

¨x(t)+x(t)+ax(t)

=0. (8.3)

Unlike the situation for the linear spring where the solution is in terms of ele-

mentary functions (sine or cosine), the solution for the hard spring will involve

another “special” function, the Jacobian elliptic function. As you will see, by

using computer algebra, deriving the solution is not as hard as it seems.

Taking, say, a=1, the equation of motion (8.3) is entered in ode .

restart: a:=1:

ode:=diff(x(t),t,t)+x(t)+a*x(t)ˆ3=0;

ode := (

x(t)) + x (t)+x (t)

The dsolve command is applied to ode .Twoimplicit solutions are generated,

with the time t giveninintegralform(

a is the integration variable) and two

arbitrary constants,

C1 and C2 , present. The positive square root solution

must be selected in order to produce a positive period for the oscillatory motion.

sol:=dsolve(ode,x(t)); #choose positive root

sol :=



x (t)

√

−4 a

− 2 a

+4 C1

a − t − C2 =0,



x (t)

−

√

−4 a

− 2 a

+4 C1

a − t − C2 =0

The positive square root solution (the ﬁrst one here) is selected and diﬀerenti-

ated with respect to t. This removes the constant

C2 .

eq1:=diff(sol[1],t);

298 CHAPTER 8. NLODES & PDES OF PHYSICS

eq1 :=

x(t))



−4 x (t)

− 2 x (t)

+4 C1

− 1=0

eq1 is solved in eq2 for dx(t)/dt, i.e., the speed.

eq2:=solve(eq1,diff(x(t),t));

eq2 :=



−4 x (t)

− 2 x (t)

+4 C1

In order to perform the subsequent integrations, the dependent variable x(t)is

replaced in eq2 with a time-independent symbol, say y .

eq3:=subs(x(t)=y,eq2);

eq3 :=



−4 y

− 2 y

+4 C1

Let’s call the amplitude of the oscillations A.Wheny ≡ x(t)=A, the speed of

the mass m will be zero. This condition is expressed in eq4 .

eq4:=subs(y=A,eq3)=0;

eq4 :=

√

−4 A

− 2 A

+4 C1

eq4 is solved for the constant

C1 , which is automatically substituted into eq3

which is relabeled as eq5 .

_C1:=solve(eq4,_C1); eq5:=eq3;

C1 := A

eq5 :=



−4 y

− 2 y

+4A

+2A

Since dx(t)/dt ≡ dy/dt = eq5 ,theperiodT (time for a complete oscillation)

is given by T =2



−A

(1/eq5 ) dy. This integral is entered. To accomplish the

integration, the assumption that A>0 must be included.

T:=2*int(1/eq5,y=-A..A) assuming A>0;

T :=

√

2 EllipticK(

√

2 A

√

1+A

)

√

2+2A

Highlighting EllipticK in the output and opening the Help page reveals that the

period T involves the complete elliptic integral of the ﬁrst kind. The incomplete

elliptic integral of the ﬁrst kind is deﬁned [AS72] by

F (φ\α)=



(1 − sin

α sin

θ)

−1/2

dθ, (8.4)

or, on setting m ≡ sin

α, y ≡ sin θ,andx ≡ sin φ,

F (φ |m)=



[(1 − y

)(1− my

)]

−1/2

dy. (8.5)

The complete elliptic integral K ≡ K(m) corresponds to setting φ = π/2in

F (φ\α), or x =1 in F (φ |m). The command EllipticF(sin(phi),sqrt(m))

generates F (φ |m), while EllipticK(sqrt(m)) produces K(m). The period is

now plotted over the range A=0 to 5, the result being shown in Figure 8.3.