Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

238 CHAPTER 6. INTEGRAL TRANSFORMS

The initial frame of the animation is shown in Figure 6.6, with 100 frames

being taken and 500 plotting points used to produce a smooth ﬁgure. The

temperature distribution becomes eﬀectively zero throughout the entire range

x= 0 to 3 on completion of the animation.

6.2.3 Daniel Strikes Yet Again: Mr. Laplace Appears

Propaganda is a soft weapon; hold it in your hands too long,

and it will move about like a snake, and strike the other way.

Jean Anouilh, French playwright, (1910–1987)

In Example 4.1.3, we imagined a scenario where young Daniel struck a light,

elastic, initially horizontal, string of length L ﬁxed at both ends. Our approach

to determining the subsequent motion of the string was to build up a solution

using a Fourier series approach. In the present recipe, a similar struck string

is considered but now solved for the transverse displacement U(x, t) using the

Laplace transform method. In particular, we take L = π m, wave speed c =2

m/s, and an initial transverse velocity ∂U(x, 0)/∂t =3sin(2x) − 2sin(5x).

After solving for U (x, t), the motion of the string is animated over the time

interval t=0 to T =5 s.

After loading the integral transform and plots packages,

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

the values of c, L,andT are entered.

c:=2: L:=Pi: T:=5:

To simplify the notation, the addtable command is used to let F (x) represent

the Laplace transform of U(x, t) with respect to t in subsequent outputs.

addtable(laplace,U(x,t),F(x),t,s):

The wave equation for U(x, t) is entered in pde ,thevalueofc being automati-

cally substituted.

pde:=diff(U(x,t),x,x)=(1/cˆ2)*diff(U(x,t),t,t);

pde :=

∂

∂x

U (x, t)=

(

∂

∂t

U (x, t))

The Laplace transform of pde with respect to time t is taken in eq.

eq:=laplace(pde,t,s);

eq :=

F (x)=

F (x) −

(U)(x, 0) −

s U (x, 0)

The output in eq is given in terms of the initial proﬁle (U (x, 0)) and the initial

transverse velocity (D

(U)(x, 0), with D

indicating a time derivative). These

are supplied in the initial condition, ic, which is substituted into eq in eq2 .

ic:=U(x,0)=0,D[2](U)(x,0)=3*sin(2*x)-2*sin(5*x);

ic := U (x, 0) = 0, D

(U)(x, 0) = 3 sin(2 x) − 2sin(5x)

eq2:=subs(ic,eq);

6.2. LAPLACE TRANSFORMS 239

eq2 :=

F (x)=

F (x) −

sin(2 x)+

sin(5 x)

The ODE eq2 is analytically solved for F (x), subject to the two boundary

conditions, F (0)= 0 and F (L)=0

sol:=dsolve({eq2,F(0)=0,F(L)=0},F(x));

sol := F (x)=−2(16 cos(x)

+ 256 cos(x)

− 192 cos(x)

− 12 cos(x)

−3 s

cos(x) − 300 cos(x)+s

+ 16)sin(x)/(s

+ 116 s

+ 1600)

The solution U of the wave equation follows on taking the inverse Laplace

transform of the right-hand side of sol.

U:=invlaplace(rhs(sol),s,t);

U :=

sin(x) (15 cos(x)sin(4t) + 2 sin(10 t) (4 cos(x)

− (4 cos(x)

− 1)

))

U can be simpliﬁed by applying the combine command with the trig option.

U:=combine(U,trig);

U :=

cos(5 x+10 t) −

cos(−5 x+10 t) −

cos(2 x+4t)+

cos(−2 x+4 t)

The motion of the string is now animated over the time interval t=0 to T ,with

100 frames being taken and 500 plotting points used.

animate(U,x=0..L,t=0..T,frames=100,numpoints=500,

thickness=2,labels=["x","U"]);

You will have to run the animation to see how the string vibrates.

6.2.4 Inﬁnite-medium Green’s Function

Private information is practically the source

of every large modern fortune.

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

In this recipe we consider the 3-dimensional inﬁnite-medium diﬀusion equa-

tion (diﬀusion coeﬃcient d) with a point source of unit magnitude located at

r =(x, y, z) = 0 and turned on at the instant t = 0. The relevant PDE for the

Green’s function G(r, t) is of the form

∂G

∂t

− d ∇

G = δ(r) δ(t). (6.10)

Prior to the switching on of the source G = 0 everywhere. Our goal is to

determine G(x, y, z, t)fort>0 and animate the solution.

In addition to the usual integral transform and plots packages that have

appeared regularly in the recipes of this chapter, the VectorCalculus package is

also loaded because the Laplacian command will be used to enter ∇

in (6.10).

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

Eq. (6.10) is entered in Cartesian coordinates, noting that δ(r)=δ(x) δ(y) δ(z).

240 CHAPTER 6. INTEGRAL TRANSFORMS

pde:=diff(G(x,y,z,t),t)-d*Laplacian(G(x,y,z,t),’cartesian’

[x,y,z])=Dirac(x)*Dirac(y)*Dirac(z)*Dirac(t);

pde := (

∂

∂t

G(x, y, z, t)) − d ((

∂

∂x

G(x, y, z, t)) + (

∂

∂y

G(x, y, z, t))

∂

∂z

G(x, y, z, t))) = Dirac(x) Dirac(y) Dirac(z) Dirac(t)

The Green’s function G(x, y, z, t) is equal to 0 up to the instant that the source

is switched on. The Laplace transform of the ﬁrst time derivative of G will be

expressed in terms of G(x, y, z, 0), which is set equal to 0.

G(x,y,z,0):=0:

A functional operator F is formed for calculating the 1-dimensional Fourier

transform of an expression e with respect to an arbitrary variable v, the result

being given in terms of a second variable k. A similar functional operator K

is introduced to perform the 1-dimensional inverse Fourier transform of e from

k-space back to v-space.

F:=(e,v,k)->fourier(e,v,k): K:=(e,k,v)->invfourier(e,k,v):

Using F, a 3-dimensional Fourier transform of pde is taken in the spatial coordi-

nates x, y,andz. This is implemented by nesting 3 one-dimensional transforms,

the result being expressed in terms of kx, ky,andkz . A Laplace transform in

time is then performed, and the result factored.

eq1:=factor(laplace(F(F(F(pde,x,kx),y,ky),z,kz),t,s));

eq1 := laplace(fourier(fourier(fourier(G(x, y, z, t),x,kx),y,ky),z,kz ),t,s)

(d kz

+ s + d kx

+ d ky

)=1

eq1 is then solved for the Laplace-Fourier transform of G(x, y, z, t).

eq2:=solve(eq1,laplace(F(F(F(G(x,y,z,t),x,kx),y,ky),z,kz),t,s));

eq2 :=

s + d kz

+ d ky

+ d kx

Applying the inverse Laplace transform to eq2 putsusbackintot-space.

eq3:=invlaplace(eq2,s,t);

eq3 := e

(−(d kz

+d ky

+d kx

) t)

Finally, by performing three 1-dimensional inverse Fourier transforms on eq3 ,

the desired Green’s function solution G to Equation (6.10) is obtained. It is

necessary to assume that both d and t are greater than zero for the transforms

to be explicitly calculated.

G:=K(K(K(eq3,kx,x),ky,y),kz,z) assuming d>0,t>0;

G :=

(−

4 td

)

tdπ

(3/2)

√

G is expressed in terms of the radial distance r by making the algebraic substi-

tution x

+ y

+ z

6.2. LAPLACE TRANSFORMS 241

G:=algsubs(xˆ2+yˆ2+zˆ2=rˆ2,G);

G :=

(−

4 td

)

tdπ

(3/2)

√

As expected for a point source, G has spherical symmetry, depending only on

the radial distance. The above form of G is valid for t>0, with G=0 for t<0.

Taking the diﬀusion coeﬃcient d=1, G is animated over the radial range r =0

to 2.5 and time interval t =0.1 to 3, the opening frame of the animation being

shown in Figure 6.7.

animate(eval(G,d=1),r=0..2.5,t=0.1..3,frames=100,

numpoints=200,thickness=2,labels=["r","G"]);

0.7

Figure 6.7: Initial frame of animation of the inﬁnite-medium Green’s function.

As discussed in Wallace [Wal84], Green’s functions can be similarly derived for

semi-inﬁnite, and even ﬁnite media, with speciﬁed boundary conditions. The

Green’s functions can in turn be used to solve problems involving distributed

time-dependent sources with the same boundary conditions. The interested

reader is referred to Wallace’s text.

6.2.5 Our Field of Dreams

A man is not old until regrets take the place of dreams.

John Barrymore, American actor, (1882–1942)

One summer I worked, along with other students, at a uranium mine on the

Arctic circle to earn suﬃcient money to pay my college expenses for the next

year. The students and the permanent workers would often play baseball under

the midnight sun on a ﬁeld which consisted of processed mine tailings. On and

around this “ﬁeld of dreams”, the students would chat about what they planned

to do in life. From this ﬁeld of dreams came future engineers, chemists, and at

least one physicist. This recipe is dedicated to those friends of long ago.

242 CHAPTER 6. INTEGRAL TRANSFORMS

Radioactive gas is diﬀusing at a steady rate into the atmosphere from leveled

mine tailings. Treating the ground and atmosphere as semi-inﬁnite media, the

X-axis is taken vertically upwards with the ground-atmosphere boundary at

X = 0. Letting d be the diﬀusion constant and λ the decay rate of the radioactive

gas (only one species is considered), the relevant modiﬁed diﬀusion equation for

the gas concentration C and boundary condition (Fick’s law) are of the form

∂C(X, T )

∂T

= d

∂

C(X, T )

∂X

− λC(X, T ), −d

∂C(0,T)

∂X

= k, (6.11)

with T thetimeandk a positive constant (units of kg/(m

· s)). The various

constants may be removed from the equations by setting t = λT, x =



λ/d X,

and c(x, t)=(

√

dλ/k) C(X, T ). Then Equation (6.11) becomes

∂c(x, t)

∂t

∂

c(x, t)

∂x

− c(x, t),

∂c(0,t)

∂x

= −1. (6.12)

Assuming that the concentration of radioactive gas in the atmosphere is initially

zero, our goal is to determine the normalized concentration c(x, t) in the region

x>0 at arbitrary time t>0 and animate the result.

In addition to the integral transform package, the DEtools package is loaded.

Using the latter will enable us to seek exponential solutions to the Laplace

transformed PDE which are convenient for solving the problem.

restart: with(inttrans): with(DEtools):

The normalized diﬀusion equation in (6.12) is entered in pde .

pde:=diff(c(x,t),t)=diff(c(x,t),x,x)-c(x,t);

pde :=

∂

∂t

c(x, t)=(

∂

∂x

c(x, t)) − c(x, t)

The addtable command is used to simplify the notation. The Laplace trans-

form of c(x, t) will be given as F(x).

addtable(laplace,c(x,t),F(x),t,s):

The initial normalized concentration, c(x, 0), of radioactive gas in the atmo-

sphere is set equal to zero, and the Laplace transform is applied to pde .

c(x,0):=0: eq:=laplace(pde,t,s);

eq := s F (x)=(

F (x)) − F (x)

The expsols command is used to obtain exponential solutions for F (x)ineq .

sol:=expsols(eq,F(x));

sol := [e

(

√

s+1 x)

(−

√

s+1 x)

]

As x →∞, the concentration of radioactive gas in the atmosphere must go to

zero. Therefore the exponential solution with the minus sign (the second result

in sol) is selected and multiplied by an arbitrary constant B.

sol2:=F(x)=B*sol[2];

6.3. BROMWICH INTEGRAL AND CONTOUR INTEGRATION 243

sol2 := F (x)=Be

(−

√

s+1 x)

The boundary condition at x=0 in (6.12) is applied as follows. The right-hand

side of sol2 is diﬀerentiated with respect to x and evaluated at x =0. This

result is equated to the Laplace transform of −1.

bc:=eval(diff(rhs(sol2),x),x=0)=laplace(-1,t,s);

bc := −B

√

s +1=−

The Laplace transformed boundary condition, bc, is solved for B and the form

of the transformed concentration, F (x), displayed.

B:=solve(bc,B); sol2;

B :=

√

s +1s

F (x)=

(−

√

s+1 x)

√

s +1s

To facilitate the calculation of the inverse transform, let’s make a variable

change, substituting s = y − 1ontherhsofsol2 and multiplying the result

by e

−t

to keep the Bromwich integral (6.9) the same.

F2:=exp(-t)*subs(s=y-1,rhs(sol2));

F2 :=

(−t)

(−

√

yx)

√

y (y − 1)

Taking the inverse Laplace transform of F2 and assuming that x>0 produces

the concentration c, the answer given in terms of the complimentary error func-

tion,erfc(z)=1− erf(z).

c:=invlaplace(F2,y,t) assuming x>0;

c := −

erfc(

x +2t

√

) e

erfc(−

−x +2t

√

) e

(−x)

The concentration c is now animated over the time interval t = 0 to 10, 100

frames being taken.

plots[animate](c,x=0..10,t=0..10,frames=100,thickness=2,

numpoints=200,color=blue,labels=["x","c"]);

Execute the recipe to see how the concentration of radioactive gas in the atmo-

sphere builds up from zero to the steady-state concentration, where the diﬀusion

of radioactive gas molecules is balanced by their radioactive decay.

6.3 Bromwich Integral and Contour Integration

In performing inverse Laplace transforms in the previous section, we bypassed

evaluating the Bromwich integral by simply using the invlaplace command.

This approach is equivalent to consulting tables of Laplace transforms and their

inverses in the pre-computer age. Occasionally, Maple may stumble, unless

carefully guided, in performing the inverse Laplace transform, so it is instructive

to look at a few examples of evaluating the Bromwich integral.

244 CHAPTER 6. INTEGRAL TRANSFORMS

Recall that the inverse Laplace transform L

−1

(F (s)) = f(t)isgivenby,

f(t)=

2πi



c+i ∞

c−i ∞

F (s) ds, t > 0, and f(t)=0,t<0. (6.13)

The integration involved in determining f(t)fort>0 is to be performed along

a line s=c + iy in the complex s-plane, c being chosen so that s=c lies to the

right of all the poles of F (s).

The central idea is to replace the integral in (6.13) with the closed con-

tour integral

2πi



F (s) ds, and apply Cauchy’s residue theorem. The basic

Bromwich contour Γ is shown on the left of Figure 6.8, consisting of a vertical

“leg” A → B between s=c −iLand s= c + iL, and a circular arc B → C → A,

centeredontheoriginO, of radius R.

Re(s)

Im(s)

c+iL

c-i L

cut

Figure 6.8: Left: Basic contour. Right: Contour modiﬁed with branch cut.

For L →∞, the integral contribution A → B will be the Bromwich integral.

By Cauchy’s residue theorem, assuming that the circular arc (B → C → A)

contribution approaches 0 as R →∞, the Bromwich integral will be equal to

the sum of the residues of e

F (s)atthepolesofF (s)insideΓ.

A suﬃcient condition for



B→C→A

F (s) ds → 0asR →∞is that

|F (s)| < constant/R

,withk>0. This condition always holds if F (s)=

P (s)/Q(s), where P (s)andQ(s) are polynomials and the degree of P is less

than that of Q. The suﬃciency condition still holds if F (s) has branch points.

However, then the basic Bromwich contour must be modiﬁed. For example,

if F (s) has a branch point at the origin, the contour is modiﬁed as shown on the

right of Figure 6.8, a branch cut being introduced from −∞ to 0. The contour

is then given by A → B → C → D → E → F → A, the horizontal legs C → D

and E → F ultimately being shrunk onto the Re(s) axis and the small circular

arc D → E shrunk onto the origin.

6.3. BROMWICH INTEGRAL AND CONTOUR INTEGRATION 245

6.3.1 Spiegel’s Transform Problem Revisited

When one pays a visit it is for the purpose of wasting

other people’s time, not one’s own.

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

In Recipe 06-2-1, we successfully solved a fourth order inhomogeneous ODE

with speciﬁed initial conditions, suggested by Spiegel, using the laplace and

invlaplace commands in the integral transform package. In this recipe, we

shall solve the inverse Laplace transform of the previously obtained transform

F (s)=(s

− 2 s

+6s

− 6 s

+5s − 3)/(s

+3s

+1),

by evaluating the Bromwich integral using contour integration. In this case,

F (s) ≡ P (s)/Q(s) involves the ratio of two polynomials, the largest exponent of

s in P being one less than in Q. On the circular arc (radius R)oftheBromwich

contour Γ, s = Re

iθ

,so|F (s)|∼1/R for suﬃciently large R. Applying the

suﬃciency condition, the circular arc contribution will be 0 for R →∞.So,

taking the vertical leg to the right of all poles of F (s), the Bromwich integral

is equal to the sum of the residues of e

F (s)atthepolesofF (s)insideΓ.

F (s) is now entered and the Bromwich integrand e

F (s) formed in f.

restart:

F(s):=(sˆ5-2*sˆ4+6*sˆ3-6*sˆ2+5*s-3)/(sˆ6+3*sˆ4+3*sˆ2+1);

F (s):=

− 2 s

+6s

− 6 s

+5s − 3

+3s

f:=exp(s*t)*F(s);

f :=

(st)

− 2 s

+6s

− 6 s

+5s − 3)

+3s

The denominator of f is factored in d, revealing that F (s) has third-order poles

at s

+ 1=0, i.e., at s = ±i.

d:=factor(denom(f));

d := (s

+1)

Setting the denominator of f equal to 0, and solving for s yields the explicit

locations (±I) of the poles, needed for evaluating the residues of F (s). Each

root is reproduced three times, again indicating that each pole is of third order.

pole:=solve(denom(f)=0,s);

pole := −I, I, −I, I, −I, I

Creating a functional operator r to evaluate the residue at an arbitrary pole,

r:=v->residue(f,s=pole[v]):

the residues of the two poles are determined in r1 and r2 .

r1:=r(1); r2:=r(2);

r1 :=

I (8 t − 4 It−

+8I −

I (2 t − 16 + 5 I))

(tI)

246 CHAPTER 6. INTEGRAL TRANSFORMS

r2 :=

I(−8 te

(tI)

− 4 Ite

(tI)

+8Ie

(tI)

−

(tI)

(−2 t +16+5I))

The inverse Laplace transform, ILT , (Bromwich integral), then is equal to the

sum of the two residues. Applying the complex evaluation command to the

sum and simplifying yields the following result:

ILT:=simplify(evalc(r1+r2));

ILT := cos(t)+

t cos(t) −

sin(t) −

sin(t)+2t sin(t)

The answer is further simpliﬁed by collecting sine and cosine terms.

ILT:=collect(%,[sin,cos]);

ILT := (2 t −

−

)sin(t)+(1+

5 t

) cos(t)

As a check, the inverse Laplace transform of F (s)=f/e

is calculated with the

invlaplace command, yielding the same analytic result.

check:=inttrans[invlaplace](f/exp(s*t),s,t);

check := −

sin(t)(−16 t + t

+ 21) +

cos(t)(8+5t)

6.3.2 Ms. Curious’s Branch Point

Propaganda is that branch of the art of lying which consists in nearly

deceiving your friends without quite deceiving your enemies.

F. M. Cornford, British author, poet, (1874–1943)

Ms. I. M. Curious has reached an important branch point in her educational

career. Although, she has done very well to this point in her course work,

she has to decide on whether to continue as a physics major and eventually

seek an academic career, or would she better oﬀ to switch into the engineering

physics program and go into industry. On seeking Jennifer’s advice, Jennifer

has temporarily deﬂected I. M.’s question by suggesting that, whatever decision

I. M. makes, this week’s homework is due tomorrow and no extensions will be

allowed. Ironically, one of the questions on this week’s quiz also deals with a

branch point, albeit of the mathematical variety.

Speciﬁcally, I. M. is requested to calculate L

−1

−a

√

/s), with a>0, using

contour integration. Then, employing this answer, determine L

−1

−a

√

s).

Here is the recipe that I. M. has created. She begins by entering F(s)=

−a

√

/s, which has a branch point at s = 0, and then forming the Bromwich

integrand f =e

F (s)/(2πi).

restart:

F(s):=exp(-a*sqrt(s))/s; f:=exp(s*t)*F(s)/(2*Pi*I);

6.3. BROMWICH INTEGRAL AND CONTOUR INTEGRATION 247

F(s):=

(−a

√

f :=

−1

(st)

(−a

√

sπ

Choosing the modiﬁed contour Γ on the right of Figure 6.8, I. M. notes that the

only singular point in F(s)isats= 0, which is outside the contour. Therefore,

by Cauchy’s theorem, the closed line integral



fds=0.

I. M. takes the radius of the outer circular arc (B → C and F → A in

Figure 6.8) to be R and the radius of the inner arc (D → E)tobe. The limits

R →∞and  → 0 are to be taken. As R →∞, the outer arc contribution

(



B→C



F →A

) fds → 0, and



A→B

fds is just the Bromwich integral, BI .

Then, BI =−(



C→D



D→E



E→F

) fds, the three integrals to be evaluated.

To determine IDE =



D→E

fds,I.M.setss = s1 = e

iθ

. To evaluate

ICD =



C→D

fds,I.M.notesthatθ =π for this leg, so that s= re

iπ

=−r and

s2 ≡

√

iπ/2

√

r.ForIEF =



E→F

fds, θ = −π so that s=re

−iπ

=−r

still, but s3 ≡

√

−iπ/2

=−i

√

r.I.M.nowenterss1 , s2 ,ands3 .,

s1:=epsilon*exp(I*theta); s2:=sqrt(r)*exp(I*Pi/2);

s3:=sqrt(r)*exp(-I*Pi/2);

s1 := e

(θI)

s2 :=

√

rI s3 := −I

√

The integral IDE is given by lim

→0



−π

(fis)|

s=s1

dθ, which is now evaluated.

IDE:=int(limit(subs(s=s1,f*I*s),epsilon=0),theta=Pi..-Pi);

IDE := −1

To evaluate ICD and IEF, a functional operator G is introduced to calculate



∞

(f)|

√

s=v,s=−r

dr,wherev is equal to s2 for ICD,ands3 for IEF .Theinert

form of the integral command is used here, as Maple is not able to actually

perform the integrations.

G:=v->Int(subs({sqrt(s)=v,s=-r},f),r=0...infinity):

Using G, ICD and IEF take the following forms. For the latter, I. M. inserts a

minus sign, since the integral IEF is in the opposite direction to ICD.

ICD:= G(s2); IEF:=-G(s3);

ICD :=



∞

(−rt)

(−Ia

√

rπ

dr IEF := −



∞

(−rt)

√

rI)

rπ

The Bromwich integral BI then is equal to −(IDE + ICD + IEF ). This sum is

now converted to a trig form and the combine command applied.

BI:=-combine(convert(IDE+ICD+IEF,trig));

BI := −



∞

cosh(rt)sin(a

√

r) − sinh(rt)sin(a

√

rπ

dr +1

The integrand of BI can be simpliﬁed by making the algebraic substitution

cosh(rt) − sinh(rt)=e

−rt

BI:=algsubs(cosh(r*t)-sinh(r*t)=exp(-r*t),BI);