Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

228 CHAPTER 6. INTEGRAL TRANSFORMS

The plots and integral transform packages are loaded and the 1-dimensional

diﬀusion equation entered in pde for the concentration C(x, t)ofink.

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

pde:=diff(C(x,t),t)=d*diff(C(x,t),x,x);

pde :=

∂

∂t

C(x, t)=d (

∂

∂x

C(x, t))

Since there can be no ﬂow of ink through the x = 0 end, the concentration

gradient ∂C/∂x must be zero there. This boundary condition implies that we

should use a Fourier cosine transform approach on the diﬀusion equation. The

Fourier cosine transform (fouriercos) command is applied to pde .

eq:=fouriercos(pde,x,k);

eq :=

∂

∂t

fouriercos(C(x, t),x,k)=

−

d (

√

(C)(0,t)+k

fouriercos(C(x, t),x,k)

√

π)

√

The result in eq is expressed in terms of the concentration gradient at x =0,

viz., D

(C)(0,t). This term is set equal to zero in eq2 . The notation is also

simpliﬁed by replacing fouriercos(C(x, t),x,k)withF(k, t).

eq2:=subs({fouriercos(C(x,t),x,k)=F(k,t),D[1](C)(0,t)=0},eq);

eq2 :=

∂

∂t

F(k, t)=−dk

F(k, t)

The initial concentration proﬁle f =Aδ(x − x0) of the ink is entered,

f:=A*Dirac(x-x0);

f := A Dirac(−x + x0 )

and its Fourier cosine transform taken and simpliﬁed assuming that x0 > 0.

F1:=simplify(fouriercos(f,x,k)) assuming x0>0;

F1 :=

√

2 cos(k x0 )

√

The ODE eq2 is analytically solved for F(k, t).

eq3:=dsolve(eq2,F(k,t));

eq3 := F(k, t)=

F1(k) e

(−dk

At t=0, F(k, 0)= F1(k), so F1(k)=F1 is substituted into eq3 .

eq4:=subs(_F1(k)=F1,eq3);

eq4 := F(k, t)=

√

2 cos(k x0 ) e

(−dk

√

The general desired expression for the concentration follows on applying the

Fourier cosine transform to the right-hand side of eq4 , taking us from k-space

back into x-space.

C:=fouriercos(rhs(eq4),k,x);

6.1. FOURIER TRANSFORMS 229

C :=



(−

4 td

)

cosh(

x0 x

2 td

)

The time T at which the ink concentration is an extremum at x= 0 is obtained

by evaluating C at this point, diﬀerentiating with respect to t, setting the result

equal to zero, and solving for the time t. The plot of the concentration at x=0

will conﬁrm that the time T corresponds to the maximum ink concentration at

that end of the tube.

T:=solve(diff(eval(C,x=0),t)=0,t);

T :=

2 d

The parameter values A =1, d=1, andx0= 10 are entered and the concentration

proﬁle C(x, t)andtimeT then determined.

A:=1: d:=1: x0:=10: C:=C; T:=T;

C :=



(−

100+x

4 t

)

cosh(

5 x

)

T := 50

The concentration reaches its maximum value at the x = 0 end 50 time units

after the ink drop is injected.

The concentration proﬁle is animated over the time interval t =1 to 2T , 100

frames being taken. The animation does not start at t = 0 because the initial

proﬁle is a Dirac delta function. The initial frame is shown in Figure 6.3. Run

the animation and see how the proﬁle evolves with time.

animate(C,x=0..35,t=1..2*T,frames=100,thickness=2,

numpoints=200,labels=["x","C"]);

0.05

0.1

0.2

0.25

Figure 6.3: Initial frame in the animation of the ink concentration.

230 CHAPTER 6. INTEGRAL TRANSFORMS

The concentration of ink at x = 0 is now plotted in Figure 6.4 over the time

interval t= 0 to 2 T . The maximum concentration occurs at 50 time units.

plot(eval(C,x=0),t=0..2*T,thickness=2,

labels=["t","C(0,t)"]);

0.01

0.02

C(0,t)

0.04

20 40

80 100

Figure 6.4: Time evolution of the ink concentration at x=0.

6.1.4 Diﬀusive Heat Flow

In my own time there have been inventions of this sort,...tubes for

diﬀusing warmth equally through all parts of a building...

Seneca, Roman writer, philosopher, and statesman, (5BC–65AD)

In contrast to the radiant heating referred to by Seneca, involving the con-

vective ﬂow of warm water in tubes under the ﬂoor, this recipe involves the

diﬀusive ﬂow of heat in an inﬁnitely long thin solid rod. Suppose that the

initial temperature distribution inside the rod is T (x, 0) = T

−b

.Whatis

the temperature distribution T (x, t) for times t ≥ 0? Animate the suitably

normalized temperature distribution.

Since the domain is inﬁnite, the full Fourier transform will be used. Because

the inverse will also be calculated, it is not necessary to worry about the fact

that the Maple convention for splitting the factor 1/(2π) in Fourier’s integral

theorem diﬀers from the convention adopted in this chapter. As before, the

integral transform and plots (needed for the animation) packages are loaded.

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

In order for the transforms to be explicitly calculated, the assumption that

b>0, t>0, and the heat diﬀusion coeﬃcient a>0, is supplied.

assume(a>0,b>0,t>0):

The one-dimensional heat diﬀusion equation is entered.

6.1. FOURIER TRANSFORMS 231

de:=diff(T(x,t),t)=aˆ2*diff(T(x,t),x,x);

de :=

∂

∂t

T(x, t)=a

(

∂

∂x

T(x, t))

The Fourier transform of the PDE de is taken in eq and notationally simpliﬁed

in eq2 by substituting F(k, t) for fourier(T(x, t),x,k).

eq:=fourier(de,x,k);

eq :=

∂

∂t

fourier(T(x, t),x,k)=−a

fourier(T(x, t),x,k)

eq2:=subs(fourier(T(x,t),x,k)=F(k,t),eq);

eq2 :=

∂

∂t

F(k, t)=−a

F(k, t)

Using the dsolve command, the ODE eq2 is analytically solved for F(k, t).

sol:=dsolve(eq2,F(k,t));

sol := F(k, t)=

F1(k) e

(−a

The arbitrary function F1(k) is determined by Fourier transforming the initial

temperature proﬁle, T

−b

. Its functional form is automatically substituted

into sol, the result being labeled sol2 .

_F1(k):=fourier(T0*exp(-bˆ2*xˆ2),x,k); sol2:=sol;

F1(k):=T0 e

(−

4 b

)



sol2 := F(k, t)=T0 e

(−

4 b

)



(−a

The temperature proﬁle inside the rod for t>0 follows on taking the inverse

Fourier transform (invfourier) of the right-hand side of sol2 .

T:=invfourier(rhs(sol2),k,x);

T :=

T0 e

(−

4 b

t+1

)

√

4 b

t +1

The temperature is normalized (made dimensionless) by dividing T by T0 ,and

normalized spatial (ζ) and time (τ) coordinates are introduced by substituting

x=ζ/b and t= τ/(b

Tnorm:=subs({x=zeta/b,t=tau/(bˆ2*aˆ2)},T/T0);

Tnorm :=

(−

4 τ +1

)

√

4 τ +1

The normalized temperature, Tnorm, is animated over the time interval τ =0

to 10, the spatial range being ζ =−10 to 10.

animate(Tnorm,zeta=-10..10,tau=0..10,frames=100,

thickness=2,numpoints=200,labels=["zeta","Tnorm"]);

Run the animation to see the initially Gaussian temperature proﬁle diﬀuse away.

232 CHAPTER 6. INTEGRAL TRANSFORMS

6.1.5 Deja Vu

The postmodern reply to the modern consists of recognizing that the

past, since it cannot really be destroyed, because its destruction leads

to silence, must be revisited...

Umberto Eco, Italian novelist, Reﬂections on the Name of the Rose, (1983)

If this example seems familiar, it is. We revisit Recipe 01-1-6: Mr. Dirac’s

Famous Function and, instead of using the Green’s function method, Fourier

transform the relevant diﬀerential equation and illustrate how contour integra-

tion can be used to evaluate the inverse transform.

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

which exerts a Hooke’s law restoring force (spring constant K) on the beam.

Subject to a steady point force F exerted at x = 0, the static displacement

ψ(x)=(F/K) G(x) of the beam satisﬁes the following ODE,

+ G = δ(x), (6.7)

with α =(YI/K)

1/4

, Y being Young’s modulus and I the beam’s moment of

inertia. Determine G(x) for all x.

After loading the integral transform package, the following interface com-

mand is used to set j =

√

−1, to avoid confusion with the moment of inertia.

restart: with(inttrans): interface(imaginaryunit=j):

Equation (6.7) is entered in ode and then Fourier transformed in eq .

ode:=alphaˆ4*diff(G(x),x$4)+G(x)=Dirac(x);

ode := α

(

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

eq:=fourier(ode,x,k);

eq := fourier(G(x),x,k)(α

+1)=1

eq is solved in eq2 for the Fourier transform of G(x).

eq2:=solve(eq,fourier(G(x),x,k));

eq2 :=

With Maple’s convention for the Fourier transform, the inverse Fourier trans-

form is given by G(x)=



∞

−∞

jkx

(eq2 /(2π)) dk. The integrand is now entered.

integrand:=exp(j*k*x)*eq2/(2*Pi);

integrand :=

(kxj)

(α

+1)π

The four complex k zeros of the integrand’s denominator are determined.

sol:={solve(denom(integrand)=0,k)};

6.1. FOURIER TRANSFORMS 233

sol :=

⎧

⎪

⎨

⎪

⎩

√

2 −

√

−

√

−

√

−

√

⎫

⎪

⎬

⎪

⎭

Theintegrandhasfoursimplepoles,theﬁrsttwoentriesinsol corresponding

to poles in the upper-half of the complex k-plane, the last two entries to poles in

the lower-half k-plane. The integral in G(x) is evaluated by performing a closed

contour integration along the real k-axis from −R to +R and then around a

semi-circle of radius R (where k =Re

jθ

), the limit R →∞being taken. Noting

that on the semi-circle the exponential term in the integrand takes the form

jRe

jθ

= e

jR cos(θ)x

−R sin(θ)x

, the semi-circle must be taken in the upper-half

k-plane (where sin(θ) > 0) for x>0 for the semi-circular contribution to vanish

as R →∞. In this case, the two poles in the upper-half plane will contribute

to the integral.

For x<0, on the other hand, the semi-circle must be taken in the lower-half

plane where sin(θ) < 0. The two poles in the lower-half k-plane then contribute.

A functional operator is formed for evaluating the residue of the integrand

for the ith simple pole.

res:=i->simplify(evalc(residue(integrand,k=sol[i]))):

Using the functional operator, the sequence of four residues is explicitly calcu-

lated, the ﬁrst two residues corresponding to the two poles in the upper-half

plane, the last two residues to the two poles in the lower-half plane.

eq3:={seq(res(i),i=1..4)};

eq3 :=



−

(−

√

2 α

)

√

2 (cos(

√

2 α

) − sin(

√

2 α

)+sin(

√

2 α

) j + cos(

√

2 α

) j)

απ

−

(−

√

2 α

)

√

2(−cos(

√

2 α

)+sin(

√

2 α

)+sin(

√

2 α

) j + cos(

√

2 α

) j)

απ

(

√

2 α

)

√

2 (cos(

√

2 α

)+sin(

√

2 α

) − sin(

√

2 α

) j + cos(

√

2 α

) j)

απ

(

√

2 α

)

√

2(−cos(

√

2 α

) − sin(

√

2 α

) − sin(

√

2 α

) j + cos(

√

2 α

) j)

απ



The last two residues, containing e

√

2/(2α)

, are selected in eq4 , while the ﬁrst

two residues, containing e

−x

√

2/(2α)

, are selected in eq5 .

eq4:=select(has,eq3,exp(x*sqrt(2)/(2*alpha))):

eq5:=select(has,eq3,exp(-x*sqrt(2)/(2*alpha))):

234 CHAPTER 6. INTEGRAL TRANSFORMS

For x<0 (to the left of the applied force at the origin), then G(x) is obtained

by summing the two residues in eq4 , multiplying by 2 jπ, applying the complex

evaluation command, and simplifying. The minus sign is inserted because for

the semi-circle in the lower-half k-plane, the contour is taken clockwise, rather

than counter-clockwise. The result, assigned the name GL,isthesameas

obtained previously in Recipe 01-1-6.

GL:=-simplify(evalc(2*j*Pi*(eq4[1]+eq4[2])));

GL :=

(

√

2 α

)

√

2(−sin(

√

2 α

) + cos(

√

2 α

))

The solution for x>0 is similarly obtained by using eq5 . The contour is

counter-clockwise in this case, so a minus sign is not inserted. The result GR

is identical with that previously obtained using the Green’s function method.

GR:=simplify(evalc(2*j*Pi*(eq5[1]+eq5[2])));

GR :=

(−

√

2 α

)

√

2 (cos(

√

2 α

)+sin(

√

2 α

))

6.2 Laplace Transforms

The Laplace transform of a function f(t) is deﬁned to be

L(f(t)) = F (s)=



∞

−st

f(t) dt, (6.8)

where s = x + iy is a complex variable. The inverse transform is given by the

Bromwich integral formula,

−1

(F )=f(t)=

2πi



c+i ∞

c−i ∞

F (s) ds, t > 0 (6.9)

and f(t)=0 for t<0. The integration in (6.9) is to be performed along a line

s= c + iy in the complex plane. The real number c is chosen so that s =c lies to

the right of all poles, branch points, or essential singularities. With Maple, the

laplace command is used to calculate the Laplace transform, and invlaplace

for the inverse transform.

A couple of the more important properties of Laplace transforms are:

• The Laplace transform of the nth derivative of f(t)isgivenby





= s

F (s) −

n−1



k=0





t=0

n−k−1

• Deﬁning the convolution of the functions f

(t)andf

(t)tobe

C(t)=



(t − y) f

(y) dy, then L(C(t)) = F

(s) F

(s),

where F

(s)andF

(s) are the Laplace transforms of f

(t)andf

(t).

6.2. LAPLACE TRANSFORMS 235

6.2.1 Jennifer Consults Mr. Spiegel

Histories are more full of examples of the ﬁdelity of dogs

than of friends.

Alexander Pope, English satirical poet, (1688–1744).

To illustrate the use of Laplace transforms in solving ODEs, Jennifer has cre-

ated a recipe for solving the following problem taken from Murray Spiegel’s

Advanced Mathematics. Solve the fourth order inhomogeneous ODE



(t)+2y



(t)+y(t)=sint, y(0) = 1,y



(0) = −2,y



(0) = 3,y



(0) = 0.

After loading the integral transform package, Jennifer lets the Laplace transform

of y(t) be represented by F (s) for notational convenience.

restart: with(inttrans): LT:=laplace(y(t),t,s)=F(s):

The ODE is now entered.

ode:=diff(y(t),t$4)+2*diff(y(t),t$2)+y(t)=sin(t);

ode := (

y(t)) + 2 (

y(t)) + y(t)=sin(t)

She then takes the Laplace transform of ode and substitutes LT ,

eq:=subs(LT,laplace(ode,t,s));

eq := s

F (s) − (D

(3)

)(y)(0) − s (D

(2)

)(y)(0) − s

D(y)(0) − s

y(0)

+2s

F (s) − 2D(y)(0) − 2 s y(0) + F (s)=

the result being expressed in terms of y(0) and the ﬁrst three derivatives at

t = 0. Maple has simply applied the rule for taking the Laplace transform of a

derivative to the fourth and second order derivatives. The last term, 1/(s

+1),

in eq is the Laplace transform of sin(t). If proceeding by hand, one would look

up this result in a table of Laplace transforms.

The initial conditions are now entered in ic, using the diﬀerential operator,

ic:=(y(0)=1,D(y)(0)=-2,D(D(y))(0)=3,D(D(D(y)))(0)=0):

and eq is evaluated in eq2 with the initial conditions.

eq2:=eval(eq,{ic});

eq2 := s

F (s)+4− 5 s +2s

− s

+2s

F (s)+F (s)=

Jennifer then solves eq2 for F (s).

eq3:=solve(eq2,F(s));

eq3 :=

−6 s

− 3+6s

+5s − 2 s

+ s

+3s

The solution, Y , of the original ODE then follows on applying the inverse

Laplace transform to eq3 . If proceeding by hand, as in Spiegel, one would

laboriously express eq3 as a sum of partial fractions and then use a table to

look up the inverse transforms of the various terms.

236 CHAPTER 6. INTEGRAL TRANSFORMS

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

Y := −

sin(t) (21 + t

− 16 t)+

cos(t)(8+5t)

Exactly the same answer can be obtained more easily by using the dsolve

command with the option method=laplace, as illustrated in Y2 .

Y2:=dsolve({ode,ic},y(t),method=laplace);

Y2 := y(t)=−

sin(t) (21 + t

− 16 t)+

cos(t)(8+5t)

Finally, Jennifer completes the recipe by plotting Y overthetimeintervalt=0

to 100, the solution being shown in Figure 6.5.

plot(Y,t=0..100,thickness=2,labels=["t","y"],numpoints=200);

–500

500

1000

20 40

80 100

Figure 6.5: Time-dependent solution of the ODE.

As time increases, the oscillations grow indeﬁnitely in amplitude, never settling

down to a steady-state solution.

6.2.2 Jennifer’s Heat Diﬀusion Problem

Gossip is a sort of smoke that comes from the dirty ...pipes of those

who diﬀuse it: it proves nothing but the bad taste of the smoker.

George Eliot, English novelist, (1819–80)

The Laplace transform approach can be used to solve the diﬀusion equation

with speciﬁed boundary conditions and an initial proﬁle. Jennifer will illustrate

this by solving the following temperature diﬀusion problem and animating it

overthetimeintervalt=0 to 0.05 seconds,

=4T

,T(0,t)=T (3,t)=0,T(x, 0)= 10 sin(2πx) −6sin(4πx)+2sin(6πx).

After loading the integral transform and plots packages,

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

6.2. LAPLACE TRANSFORMS 237

Jennifer enters the relevant PDE and the initial temperature proﬁle T (x, 0).

pde:=diff(T(x,t),t)=4*diff(T(x,t),x,x);

pde :=

∂

∂t

T (x, t)=4(

∂

∂x

T (x, t))

T(x,0):=10*sin(2*Pi*x)-6*sin(4*Pi*x)+2*sin(6*Pi*x);

Then pde is Laplace transformed in eq with respect to the time variable t,the

initial condition being automatically substituted. For notational simplicity the

Laplace transform of T (x, t) is replaced with F (x), the resulting ODE for F(x)

being displayed in eq2 .

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

eq2:=subs(laplace(T(x,t),t,s)=F(x),eq);

eq2 := s F (x) − 10 sin(2 πx)+6sin(4πx) − 2sin(6πx)=4(

F (x))

Since the boundary conditions are T (0,t)=T (3,t) = 0, their Laplace trans-

forms are equal to zero. So eq2 is analytically solved for F (x) subject to

F (0) = 0 and F (3) = 0.

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

sol := F (x) = ((160 sπ

+2s

+2048 π

)sin(6πx) − 13824((

+π

)sin(4πx)

−

sin(2 πx)(

+π

)) (

144

+π

))



+224 s

+12544 sπ

+147456 π

)

Finally, applying the inverse Laplace transform to the right-hand side of sol

yields the temperature proﬁle T (x, t)fort ≥ 0, which is animated over the time

interval t= 0 to 0.05.

T(x,t):=invlaplace(rhs(sol),s,t);

T (x, t):= 10sin(2πx) e

(−16 π

− 6sin(4πx) e

(−64 π

+2sin(6πx) e

(−144 π

animate(T(x,t),x=0..3,t=0..0.05,frames=100,numpoints=500,

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

–15

123

Figure 6.6: Initial frame of the animation.