Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

248 CHAPTER 6. INTEGRAL TRANSFORMS

BI := −



∞

sin(a

√

r) e

(−rt)

rπ

dr +1

Applying the value command to BI yields the answer to the ﬁrst part of the

question, the result being expressed in terms of the error (erf) function or,

on applying convert(BI,erfc), in terms of the complimentary error (erfc)

function.

BI:=value(BI); BI:=convert(BI,erfc);

BI := −erf(

√

)+1 BI := erfc(

√

)

So, L

−1

−a

√

/s)=erfc(a/(2

√

t)).

As a check, I. M. will derive exactly the same result by directly applying

the inverse Laplace transform command. She mentally replaces a

√

s with

√

and t with t/a

in the Bromwich integral, so that the function to be inverse

transformed is F2 = e

−

√

/y. Loading the integral transform package,

F2:=exp(-sqrt(y))/y; with(inttrans):

F2 :=

(−

√

she inverse transforms F2 and applies the radical simpliﬁcation command.

check:=radsimp(invlaplace(F2,y,t/aˆ2));

check := erfc(

√

)

The result is exactly the same as in BI .

To answer the second part of the question, I. M. notes that e

−a

√

s =

−

−a

√

/s)=−

dF (s)

. Thus, the second inverse Laplace transform, ILT2 ,

follows on diﬀerentiating BI with respect to a and multiplying by −1.

ILT2:=-diff(BI,a);

ILT2 :=

(−

4 t

)

√

Once again I. M. conﬁrms her answer by making the same variable transforma-

tion as above, forming the new function F3 = e

−

√

/(a

√

y), and applying the

inverse Laplace transform to F3 and simplifying.

F3:=exp(-sqrt(y))/(a*sqrt(y));

F3 :=

(−

√

check2:=radsimp(invlaplace(F3,y,t/aˆ2));

check2 :=

(−

4 t

)

√

πt

Theanswer,ofcourse,isthesameasinILT2 .

6.3. BROMWICH INTEGRAL AND CONTOUR INTEGRATION 249

6.3.3 Cooling That Weenie Rod

Write while the heat is in you. The writer who postpones the record-

ing of his thoughts uses an iron which has cooled to burn a hole with.

He cannot inﬂame the minds of his audience.

Henry David Thoreau, American author, philosopher, naturalist (1817–62)

Returning to our “tale” of the Northern weenie roast, having ﬁnished with

the long iron rods that they used for cooking their weenies, Russell places them

outdoors again to cool in the 0

◦

temperature. This recipe is an idealization of

that cooling process.

AverylongcircularrodofradiusR, initially having a constant temperature

throughout, has a constant temperature of 0

◦

applied to its surface for times

t ≥ 0. Determine the temperature at an arbitrary point inside the rod after the

cooling process has begun. Animate the temperature proﬁle inside the rod.

If the rod were of ﬁnite length, one would use cylindrical coordinates (r, θ, z)

where r is the radial distance from the cylinder axis, θ the angle around the cir-

cumference of the cylinder, and z the distance along the cylinder axis. Treating

the rod as being inﬁnitely long (which might have presented a “slight” problem

for the weenie roast!) removes the z-dependence from the problem. Since the

rod initially has a constant temperature throughout and the boundary condi-

tion at the surface has no angular dependence, the temperature T inside the

rod depends only on the radial distance r, i.e., T =T (r, t). T (r, t) then satisﬁes

the following heat diﬀusion equation (a

being the heat diﬀusion coeﬃcient) for

0 <r<Rand t>0,

∂T

∂t

= a



∂

∂r

∂T

∂r



,T(R, t)=0,T(r, 0) = T

. (6.14)

The problem can be made dimensionless by setting u ≡ T/T

, x ≡ r/R,and

τ ≡ a

t/R

,sothat

∂u

∂τ

∂

∂x

∂u

∂x

,u(1,τ)=0,u(x, 0) = 1. (6.15)

After loading the plots and integral transform packages, the normalized heat

ﬂow equation given in (6.15) is entered in pde .

restart: with(plots): with(inttrans):

pde:=diff(u(x,tau),tau)=diff(u(x,tau),x,x)+diff(u(x,tau),x)/x;

pde :=

∂

∂τ

u(x, τ )=(

∂

∂x

u(x, τ )) +

∂

∂x

u(x, τ )

For later notational convenience, the addtable command is used to replace the

Laplace transform of u(x, τ ), with respect to τ,withthesymbolF (x, s).

addtable(laplace,u(x,tau),F(x,s),tau,s):

250 CHAPTER 6. INTEGRAL TRANSFORMS

The initial condition u(x, 0) = 1 is entered, and pde Laplace transformed with

respect to τ in de. The initial condition is automatically substituted into de.

u(x,0):=1; de:=laplace(pde,tau,s);

u(x, 0) := 1

de := s F (x, s) − 1=(

∂

∂x

F (x, s)) +

∂

∂x

F (x, s)

The spatially dependent ODE de is analytically solved for F (x, s)insol,the

answer involving zeroth-order Bessel functions of the ﬁrst and second kinds.

sol:=dsolve(de,F(x,s));

sol := F (x, s) = BesselJ(0,

√

−sx) F2 (s) + BesselY(0,

√

−sx) F1 (s)+

Since Y

(

√

−sx) diverges at x = 0 (the central axis of the rod), it is removed

from the right-hand side of sol, F (x, s) then taking the form shown in sol2 .

sol2:=remove(has,rhs(sol),BesselY);

sol2 := BesselJ(0,

√

−sx) F2 (s)+

From the boundary condition u(1,τ) = 0, one must have F (x =1,s)=0. This

transformed boundary condition (bc) is now applied to sol2 .

bc:=eval(sol2,x=1)=0;

bc := BesselJ(0,

√

−s) F2 (s)+

The boundary condition bc is then solved for the arbitrary function

F2 (s)

which, because it is assigned, is automatically substituted into sol2 ,orF (x, s).

_F2(s):=solve(bc,_F2(s)); F(x,s):=sol2;

F2 (s):=−

BesselJ(0,

√

−s) s

F (x, s):=−

BesselJ(0,

√

−sx)

BesselJ(0,

√

−s) s

To determine u(x, τ ), we must calculate the inverse Laplace transform of F (x, s).

As you may verify, using the inverse Laplace transform command in the inte-

gral transform package will not work here because, in addition to the simple

pole at s=0, F (x, s) has an inﬁnite number of simple poles associated with the

zeros of J

(

√

−s). This implies that u(x, τ ) will consist of an inﬁnite series. To

determine the series form, the Bromwich integral must be evaluated.

Choosing the “basic” contour Γ on the left of Figure 6.8, all the poles lie

inside Γ in the limit that the radius R of the circular arc goes to inﬁnity. Since

the circular arc contribution to the line integral vanishes in this limit, the inverse

Laplace transform will be the sum of the residues of the poles inside Γ. We will

now calculate the residues.

For the p ole at s =0, applying the residue command to e

sτ

F (x, s) yields 0.

6.3. BROMWICH INTEGRAL AND CONTOUR INTEGRATION 251

residue(exp(s*tau)*F(x,s),s=0);

Now the term J

(

√

−s) in the denominator of the ﬁrst term of F (x, s)has

simple zeros at, say, λ

, λ

,....Thus the integrand corresponding to the ﬁrst

term has simple poles at s=s

= −λ

, n=1, 2, ... Theﬁrst15valuesofs

are

calculated in S using the BesselJZeros command, and S is assigned.

S:=[seq(s[n]=-(BesselJZeros(0,n))ˆ2,n=1..15)]; assign(S):

S := [s

= −BesselJZeros(0, 1)

= −BesselJZeros(0, 2)

, ···

Care must be taken in evaluating the residues at the zeros of J

(

√

−s). The

residue of the integrand at s=s

is given by

lim

s→s

(s − s

)



−

sτ

(

√

−sx)

(

√

−s)



= −



sτ

(

√

−sx)

(

√

−s)



s=s

where use has been made of L’Hospital’s rule. A functional operator U is

created to implement this latter form for the residue of the integrand at an

arbitrary zero of J

U:=sn->simplify(-eval(exp(s*tau)*BesselJ(0,sqrt(-s)*x)

/(s*diff(BesselJ(0,sqrt(-s)),s)),s=sn));

U := sn → simplify



−

(sτ)

BesselJ(0,

√

−sx)

s (

BesselJ(0,

√

−s))

s = sn



Making use of U, the formal series solution for u is given out to 15 terms.

u:=Sum(U(s[n]),n=1..15);

u :=



n=1



2 e

τ)

BesselJ(0,

√

−s

BesselJ(1,

√

−s

)

√

−s



The explicit series form of the normalized temperature distribution is now de-

termined in u2 , only the ﬁrst few terms being displayed here in the text.

u2:=evalf(value(u));

u2 := 1.601974697 e

(−5.783185964 τ )

BesselJ(0., 2.404825558 x)

− 1.064799259 e

(−30.47126234 τ )

BesselJ(0., 5.520078110 x)

+0.8513991928 e

(−74.88700679 τ )

BesselJ(0., 8.653727913 x)

− 0.7296452400 e

(−139.0402844 τ )

BesselJ(0., 11.79153444 x)+···

For τ>0, the terms quickly drop oﬀ in magnitude, but for τ = 0 many terms

must be kept to reduce the “ringing” associated with the initial step-function

temperature proﬁle.

The series solution u2 is animated over the time interval τ =0.001 to 2. The

animation is started slightly past τ = 0 to eliminate the worst of the ringing.

The initial frame of the animation is shown in Figure 6.9.

252 CHAPTER 6. INTEGRAL TRANSFORMS

animate(u2,x=0..1,tau=0.001..2,numpoints=100,frames=50,

labels=["x=r/R","u=T/T0"]);

0.2

0.4

u=T/T0

0.8

0.2 0.4

x=r/R

0.8 1

Figure 6.9: Initial frame of animated temperature distribution inside the rod.

By the time the animation reaches τ = 2, the normalized temperature proﬁle

inside the rod has eﬀectively decayed to 0.

6.4 Other Transforms

While the Fourier and Laplace transforms are the most recognized and useful

of integral transforms, there exist other transforms which can be applied in

particular geometries. We will conclude this chapter with an example involving

the Hankel transform applied to a problem having cylindrical symmetry.

The Hankel transform (F

(k)) of order m>−1/2 of a function f(r), and

the inverse transform, are given by

(k)=



∞

f(r) J

(kr) rdr, f(r)=



∞

(k) J

(kr) kdk, (6.16)

where J

is the mth order Bessel function of the ﬁrst kind. The reader is re-

ferred to standard mathematical physics texts for a discussion of the properties

of Hankel transforms. However, it should be mentioned that there is no simple

convolution theorem for the Hankel transform.

The Maple commands for performing the transform operations in equa-

tion (6.16) are hankel(f(r),r,k,m) and hankel(F

(k),k,r,m), respectively.

The integral transform package must be loaded to use these commands. In

some cases, as the following recipe illustrates, it may be easier to make use of

the basic deﬁning relations (6.16) than Maple’s Hankel transform commands.

6.4. OTHER TRANSFORMS 253

6.4.1 Meet the Hankel Transform

Reason transformed into prejudice is the worst form of prejudice,

because reason is the only instrument for liberation from prejudice.

Allan Bloom, American educator, author, (1930–1992)

We wish to determine the subsequent transverse vibrations of a thin, inﬁ-

nite, elastic membrane which is initially at rest and has an initial displacement

ψ(r, θ, t =0)=A/



1+r

. The solution ψ(r, θ, t > 0) will then be animated

in three dimensions, taking A=1, a=10, and wave speed c =1.

The plots, integral transform, and VectorCalculus packages are loaded, being

required for the animate, hankel,andLaplacian commands, respectively.

restart: with(plots): with(inttrans): with(VectorCalculus):

It is assumed that c>0, t>0, A>0, a>0, and r>0. The initial condition

for the shape is also entered.

assume(c>0,t>0,A>0,a>0,r>0): ic:= A/sqrt(1+rˆ2/aˆ2);

ic :=



Since the initial shape has no angular (θ) dependence, the subsequent membrane

displacement must also be independent of θ, i.e., ψ = ψ(r, t). This also suggests

that we set m = 0 when applying the Hankel transform and its inverse. The

relevant wave equation for ψ(r, t) is now entered in polar coordinates in pde .

pde:=expand(Laplacian(psi(r,t),’polar’[r,theta]))

-(1/cˆ2)*diff(psi(r,t),t,t)=0;

pde :=

∂

∂r

ψ(r, t)

∂

∂r

ψ(r, t)) −

∂

∂t

ψ(r, t)

Then pde is Hankel transformed with respect to the radial coordinate r,with

m= 0, and multiplied by −c

. A time-dependent ODE results in eq .

eq:=-cˆ2*hankel(pde,r,k,0);

eq := (

∂

∂t

hankel(ψ(r, t),r,k,0)) + k

hankel(ψ(r, t),r,k,0) c

To simplify the notation, the Hankel transform of ψ(r, t) with respect to r with

m= 0 is replaced with F (t)ineq .

eq2:=subs(hankel(psi(r,t),r,k,0)=F(t),eq);

eq2 := (

F (t)) + k

F (t) c

Since ∂ψ(r, t)/∂t =0 at t =0, eq2 is analytically solved using dsolve for F (t),

subject to the initial condition dF (t)/dt =0 at t =0.

sol:=dsolve({eq2,D(F)(0)=0},F(t));

sol := F(t)=

C2 cos(ckt)

254 CHAPTER 6. INTEGRAL TRANSFORMS

At t = 0, one has from eq2 ,

C2=F (0). So we need to calculate the Hankel

transform (with m=0) of the initial membrane shape. If the hankel command

was used,

C2 would be expressed as a combination of generalized hypergeomet-

ric functions, which are diﬃcult to reduce to a simpler form. This then would

lead to further diﬃculty in performing the inverse Hankel transform. Instead,

let’s determine

C2 by using the basic deﬁning relation (6.16) and calculating



∞

ic r J

(kr) dr. The result is then simpliﬁed.

_C2:=simplify(int(ic*r*BesselJ(0,k*r),r=0..infinity));

C2 :=

(−ak)

With the form of

C2 automatically substituted, F(t)isgivenbytherhsofsol.

F(t):=rhs(sol);

F (t):=

(−ak)

a cos(ckt)

Applying the inverse Hankel transform command to F (t) would again generate

an answer in terms of generalized hypergeometric functions which is diﬃcult to

simplify and express in real form. So, the basic deﬁning relations will be used.

First, however, it’s necessary to convert F (t) to an exponential form.

F2:=simplify(convert(F(t),exp));

F2 :=

Aa(e

(k (−a+ctI))

+ e

(−k (a+ctI))

)

Then ψ is obtained by calculating the integral



∞

F2 k J

(kr) dk.

psi:=int(F2*k*BesselJ(0,k*r),k=0..infinity);

ψ :=

2 r



(a − ctI)

2 r



(a + ctI)

The displacement ψ is expressed in terms of I ≡

√

−1. As the initial proﬁle was

real, ψ should be real, however, not complex. That it is real may be conﬁrmed

by applying the complex evaluation command to ψ and simplifying.

psi2:=simplify(evalc(psi));

ψ2:= Aa

√



√

+2r

−2 r

+2a

−c



√

+2r

−2 r

+2a

)

To animate the vibrations, the given parameter values are substituted into ψ2,

which is also converted to Cartesian coordinates by setting r=



+ y

psi3:=subs({A=1,a=10,c=1,r=sqrt(xˆ2+yˆ2)},psi2);

Finally, ψ3isanimatedwiththeanimate command. A triangular grid style is

used to smooth out the wave form.

animate(plot3d,[psi3,x=-100..100,y=-100..100],t=0..90,

frames=100,axes=boxed,gridstyle=triangular,orientation

=[25,65],tickmarks=[3,3,2],labels=["x","y","psi"]);

Execute the recipe to see what happens.

6.5 SUPPLEMENTARY RECIPES 255

6.5 Supplementary Recipes

06-S01: Verifying the Convolution Theorem

Using the symmetric factor convention, verify that the convolution theorem is

satisﬁed for the functions f1 = Ae

−a

and f2 = Ae

−a |x|

with a>0, A>0.

06-S02: Bandwidth Theorem

An approximately monochromatic plane wave packet in one dimension has the

instantaneous form u(x, 0) = f(x) e

, with f(x) the envelope function and

the central wave number. Consider:

(a) f

=Ae

−a

; (b) f

=Ae

−a|x|

; (c) f3 =A for |x| <a, 0for|x| >a.

Foreachofthef(x), explicitly evaluate the root mean square deviations from

the means, i.e., ∆x =

√

> − <x>

and ∆k =

√

> − <k>

.The

means <>are deﬁned in terms of the respective intensities |u(x, 0)|

= |f(x)|

and |F (k)|

, respectively, where F (k) is the Fourier transform of f(x). Show

that in each case the bandwidth theorem (the optical analogue of the uncertainty

principle)∆x ∆k ≥ 1/2 is satisﬁed.

06-S03: Solving an Integral Equation

Solve the following integral equation for f(x), and plot the solution:



∞

f(x) cos(αx) dx = {1 − α for 0 ≤ α ≤ 1, 0forα>1}.

06-S04: Verifying Parseval’s Theorem

For f (x)=x

sin(x) e

−2x

, calculate the Fourier transform F (k)=F(f(x)) using

the symmetric convention. Plot the intensities |f(x)|

and |F (k)|

.Verifythat

Parseval’s theorem is satisﬁed.

06-S05: Heat Diﬀusion in a Copper Rod

Suppose that the initial temperature distribution inside an inﬁnitely long, thin,

insulated, copper rod is T = T 0for−x0 ≤ x ≤ x0 and zero otherwise. By

Fourier transforming the heat diﬀusion equation with respect to x, solving the

resulting time-dependent ODE with the initial condition, and performing the

inverse transform, determine the temperature distribution inside the rod for

arbitrary time t>0. For copper, the thermal conductivity K = 386 W/m·K

◦

the density ρ = 8900 kg/m

, and the speciﬁc heat capacity C = 390 J/kg· K

◦

If x0 = 1 m and T 0 = 100

◦

, animate the temperature proﬁle over the time

interval t = 0 to the time it takes for the temperature at the origin to drop to

◦

. Determine the temperature at x =2 m at the time the temperature at the

origin has dropped to 50

◦

06-S06: Solving Another Integral Equation

Using the convolution theorem, solve the following integral equation for y(x):



∞

−∞

(y(u)/((x − u)

+ a

)) du =1/(x

+ b

)

, 0 <a<b.

Verify the solution by direct substitution of y(u) back into the integral equation.

Taking a= 1 and b=2, plot y(x) over the range x=−5to+5.

256 CHAPTER 6. INTEGRAL TRANSFORMS

06-S07: Free Vibrations of an Inﬁnite Beam

Consider a horizontal beam which is initially at rest (∂ψ(x,0)/∂t =0) with the

shape ψ(x,0) = Ae

−b

, A and b being positive constants. When released, the

transverse vibrations of the beam are governed by ψ

+(1/a

) ψ

=0, with

a positive and the subscripts denoting partial derivatives. By Fourier trans-

forming the beam PDE with respect to x, solving the resulting time-dependent

ODE with the initial conditions, and performing the inverse transform, deter-

mine ψ(x, t)fort>0. Taking A= b= a=1, animate the motion of the beam.

06-S08: A Potential Problem

The potential V (x, y) (in volts) in the region x>0, y>0 has the boundary

conditions V (0,y)=0 and V (x, 0) = 1. Using an appropriate Fourier transform

approach, determine V (x, y). Plot the equipotentials in 0.1 volt increments.

06-S09: Solving an ODE

Use the Laplace transform method to solve the ODE



(t)+8y(t)=32t

− 16 t, y(0) = y



(0) = y



(0) = 0.

Verify the solution by directly solving the ODE using the dsolve command

with the Laplace transform option. Plot y(t) over the range t = 0 to 1.6, and

determine qualitatively and quantitatively the time at which the ﬁrst minimum

in the solution curve occurs.

06-S10: Impulsive Force

A unit mass attached to a spring and lying on a ﬂat rough table experiences

an impulsive force F (t)=5t

sin(t)fort ≤ π and = 0 for t ≥ π, the ODE being



(t)+2y



(t)+y(t)=F (t),y(0) = 0,y



=1.

Use the Laplace transform method to determine the displacement y(t)ofthe

mass. Conﬁrm your result with the dsolve command, using the Laplace trans-

form option. Plot the solution over the range t=0 to 10. Determine the maxi-

mum displacement and the time at which it occurs.

06-S11: Bromwich Integral

Consider F (s)=1/(s

− 6s

+9s

+8s

− 24s

+ 16). Identify the singular

points of F (s). Calculate L

−1

(F (s)) by replacing the Bromwich integral with an

appropriate closed contour integral and using Cauchy’s residue theorem. Verify

the inverse Laplace transform by directly using the invlaplace command.

06-S12: Branch Point

Consider F (s)=

√

s/(s − 1). Identify the singular points. Plot a modiﬁed

Bromwich contour which can be used to evaluate L

−1

(F (s)). Use this contour

to explicitly carry out the evaluation. Verify your answer by directly using the

invlaplace command.

Chapter 7

Calculus of Variations

Although one could solve variational calculus problems with computer algebra

by completely mimicking a hand calculation, in all the recipes of this chapter

we will use the VariationalCalculus package to ease the work.

7.1 Euler–Lagrange Equation

Among all functions y(x) with ﬁxed points y(x

)=y

and y(x

)=y

at two

distinct points A and B respectively, ﬁnd the y(x) which gives an extremum

(minimum or maximum) or stationary value to the integral

I[y]=



F (x, y, y



) dx. (7.1)

The value of the integral I depends on what functional form is chosen for y(x).

In a given problem, F is a known function of x, y,andy



≡

.They which

gives an extremum value to I is found by solving the Euler–Lagrange equation,

∂F

∂y

−



∂F

∂y





=0. (7.2)

In physical problems, it is usually evident whether y minimizes or maximizes I.

If F does not explicitly depend on x, Equation (7.2) may be written as



F −y



∂F

∂y





=0, (7.3)

so that a ﬁrst integral of Equation (7.2) is F − y



(∂F/∂y



) = constant.

7.1.1 Betsy’s In A Hurry

Work expands so as to ﬁll the time available for its completion.

C. Northcote Parkinson, English writer, Parkinson’s Law (1958)

One of the oldest variational examples is the brachistochrone

problem pro-

posed by John Bernoulli in 1696 and independently solved by him and his

brother James, Gottfried Leibniz, and Isaac Newton.

brachistos ≡shortest, chronos ≡time