Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

218 CHAPTER 5. COMPLEX VARIABLES

w(z) may be obtained by setting Z = z, solving for ω, and simplifying.

w:=simplify(solve(Z=z,w),symbolic);

w :=

(1 + e

(

2 zπ

)

) e

(−

zπ

)

Converting w to trig form and applying the combine command with the trig

option, yields the ﬁnal simpliﬁed form of w.

w:=combine(convert(w,trig),trig);

w := cosh(

zπ

)

To determine the equipotentials and electric ﬁeld lines in the strip, let’s enter

z =x + iy and set h=1 for plotting purposes. Then, w= u + iv is expressed in

terms of real and imaginary parts

z:=x+I*y: h:=1: w:=evalc(w);

w := cosh(πx) cos(πy) + sinh(πx)sin(πy) I

and the forms of u and v extracted,

u:=remove(has,w,I); v:=select(has,w,I)/I;

u := cosh(πx) cos(πy) v := sinh(πx)sin(πy)

In the w plane, the equipotentials v =const. are clearly parallel to the real

axis and the ﬁeld lines u =const. perpendicular to this axis. In the z plane,

the equipotential and ﬁeld lines can be obtained by using the contourplot

command, and selecting some ﬁxed values for the contours. A function operator

cp is created to perform this task, where V will be taken to be u and v and the

color C chosen to be blue and red, respectively.

cp:=(V,C)->contourplot(V,x=0..1,y=0..1,

[contours=seq(n,n=-4..4),-1/4,-1/2,1/4,-1/2],

grid=[100,100],color=C,thickness=2):

Using the operator cp in the following display command produces the picture

shown on the right of Figure 5.13. The conducting sheet is an equipotential so

the electric ﬁeld lines in the strip intersect the edges of the sheet perpendicularly

as expected, as well as all other equipotential curves.

display({cp(u,blue),cp(v,red)},axes=frame);

5.6 Supplementary Recipes

05-S01: Roots

Determine all the roots of z

1/9

,wherez =(1+

√

3 i). Plot the roots and z in the

same ﬁgure, superimposing them on circles centered on the origin of radii |z|

and |z|

1/9

, respectively. What is the value in radians of the principal argument?

05-S02: Fluid Flow Around a Cylinder

Consider the complex function w = u + iv=V

(z + a

/z), with z = x + iy.

(a) Determine u and v and demonstrate that they satisfy both Cauchy–

Riemann conditions and Laplace’s equation.

5.6 SUPPLEMENTARY RECIPES 219

(b) Graphically show that v may be used to represent the streamlines outside

and perpendicular to an inﬁnitely long rigid cylinder of radius a. Take

=a =1.



V = −∇u, show that the ﬂuid velocity components are

given by V

= V

(−r

+ a

cos(2 θ))/r

and V

= V

sin(2 θ)/r

,where

r is the radial distance from the center of the cylinder and θ is the angle

with respect to the x-axis.

(d) Calculate the ﬂuid speed V and show that V =0 at r = a and θ =0,π.

These two points are called stagnation points.

05-S03: Constructing f(z)

If 5 x

y −10 x

+ y

+2a

xy/((x

−y

)

+4x

) is the imaginary part of

an analytic function f(z =x + iy), determine f(z).

05-S04: Analytic or Non-analytic?

Consider the function f(z)=z

cos(z

∗

). Determine whether f is analytic or

non-analytic by (a) applying the Cauchy–Riemann conditions, (b) performing

the contour integral of f(z) along the two paths y

= x and y

= x

between

z =0 and z =1 +i.

05-S05: A Contour Integral

Evaluate the integral J =



∞

−∞

((x

− x +2)/(x

+10x

+ 9)) dx using contour

integration. Conﬁrm your answer by evaluating J with the int command.

05-S06: A Higher-order Pole

Evaluate the integral J =



∞

−∞

/(x

+ a

)

) dx,witha>0, using contour

integration. Conﬁrm your answer by evaluating J with the int command.

05-S07: Another Angular Integral

Evaluate the integral J =



(cos

(3 θ)/(5 −4 cos(θ)) dθ using contour integra-

tion. Conﬁrm your answer by evaluating J with the int command.

05-S08: A Removable Singularity

Evaluate the integral J =



∞

(sin x/x) dx using contour integration and con-

ﬁrm your answer by evaluating J with the int command. The corresponding

complex integrand sin z/z has a singularity at z = 0. Because the integrand is

ﬁnite as z → 0, the singularity is called a removable singularity.

05-S09: Another Contour Integral

Evaluate the integral J =



∞

((x

+3x

− 2)/(x

+ 1)) dx using contour inte-

gration and check the value by using the int command.

05-S10: Fluid Flow & Electric Field Around a Plate

Consider the function w = u + iv=

√

− 1, with z =x + iy.

220 CHAPTER 5. COMPLEX VARIABLES

(a) By creating a suitable ﬁgure in the x-y plane, show that the constant u

curves represent the streamlines for ﬂuid ﬂow around an inﬁnitely long

plate of ﬁnite width, lying between x=−1 and +1, inserted perpendicular

to a previously uniform ﬂuid ﬂow in the y direction. Include the plate

and the equipotentials in your ﬁgure.

(b) Taking u to represent the equipotentials, create a plot showing the electric

ﬁeld vectors and equipotentials around a thin inﬁnite conducting plate of

ﬁnite width, lying between x=−1and+1.

05-S11: Another Branch Cut

By taking a branch cut along the negative real axis from the branch point

x=0 to x=−∞, evaluate the integral J =



∞

1/3

/(1 + x

)) dx using contour

integration. Check your answer with the int command.

05-S12: Laurent Expansion

Determine the singularities of the function f =cos(πz) e

−z

/(z

(z − 1)

). De-

termine the Laurent series about each singularity, keeping eight terms in each

expansion. Calculate the residue of f at each singular point.

05-S13: Capacitor Edge Eﬀects

In elementary treatments of the parallel-plate capacitor, it is assumed that the

two plates are inﬁnitely large so that “edge eﬀects” can be ignored. Using the

conformal command, show that the conformal transformation w =z + e

when

applied to the region −7 <x<2, −π<y<πproduces the equipotentials

and electric ﬁeld lines near one edge of a ﬁnite parallel-plate capacitor. Suggest

another physical problem to which this transformation applies.

Chapter 6

Integral Transforms

Applications of the Fourier, Laplace, and Hankel transforms and their inverses

are illustrated in this chapter. The integral transform (inttrans) library pack-

age plays the key role in enabling us to perform these transforms. In some

cases, assumptions must be provided about the nature of the parameters (e.g.,

whether they are positive) in order for the transforms to be explicitly done.

6.1 Fourier Transforms

Using the Euler identity e

iθ

=cosθ + i sin θ,wherei=

√

−1, the Fourier series

representation (2.8) of the function f(x), deﬁned over the interval −L ≤ x ≤ L,

may be expressed in the complex Fourier series form

f(x)=

∞



n=−∞

inπx/L

, where c



−L

f(x) e

−inπx/L

dx. (6.1)

Setting k

=nπ/L and ∆k

≡ k

n+1

− k

=π/L, Eq. (6.1) may be written as

f(x)=

2π

∞



=−∞

∆k



−L

f(y) e

−ik

. (6.2)

In the limit that L →∞, i.e., the interval becomes inﬁnite, Eq. (6.2) yields

f(x)=

2π



∞

−∞



∞

−∞

f(y) e

ik (x−y)

dy dk, (6.3)

which is known as Fourier’s integral theorem. From (6.3), it follows that if

F (k)=

√

2π



∞

−∞

f(x) e

−ikx

dx, then f(x)=

√

2π



∞

−∞

F (k) e

ikx

dk. (6.4)

The function F (k) is called the Fourier transform of f(x), while f(x)isthe

inverse Fourier transform of F(k). It should be noted that the factor 1/(2π)

in (6.3) has been split symmetrically here. This fairly standard convention

diﬀers from that used in Maple’s fourier and invfourier commands, where

the factor is split asymmetrically, with 1 in the Fourier transform and 1/(2π)

in the inverse. Since, in solving most problems of physical interest, usually

222 CHAPTER 6. INTEGRAL TRANSFORMS

both the Fourier transform and its inverse are performed, this numerical factor

splitting issue normally need not concern us.

If f(x) is an even function, Equation (6.4) yields

(k)=





∞

f(x) cos(kx) dx, f(x)=





∞

(k) cos(kx) dk. (6.5)

(k)andf(x) are the Fourier cosine transforms of each other. On the other

hand, if f(x) is an odd function one obtains the Fourier sine transforms,

(k)=





∞

f(x)sin(kx) dx, f(x)=





∞

(k)sin(kx) dk. (6.6)

The Fourier sine and cosine transforms are useful for solving problems involving

semi-inﬁnite domains, while the “full” Fourier transforms (6.4) can be applied

to situations involving inﬁnite domains.

Some of the more important properties

of Fourier transforms are as follows:

• If the Fourier transform of f(x) (denoted by F[f (x)]) is F (k), then



f(x)



=(ik)

F (k),

provided that f(x)anditsﬁrst(n − 1) derivatives vanish at x =±∞.

• The function

C(x)=

√

2π



∞

−∞

f(x − y) g(y) dy

is called the convolution of the functions f and g (often denoted fg)

over the interval −∞ to ∞.IfF (k)andG(k) are the Fourier transforms

of f(x)andg(x), respectively, the convolution theorem is

C(x)=

√

2π



∞

−∞

F (k) G(k) e

ikx

dk.

Alternately, the theorem may be stated as follows: F(fg)=F(f) F(g).

• If one sets x =0 in the convolution theorem, the Parseval relation



∞

−∞

f(−y) g(y) dy =



∞

−∞

F (k) G(k) dk

results. Another Parseval relation is



∞

−∞

f(y) g



(y) dy =



∞

−∞

F (k) G



(k) dk.

If g = f, this last result yields Parseval’s theorem,



∞

−∞

|f(y)|

dy =



∞

−∞

|F (k)|

dk.

Similar properties hold for the sine and cosine transforms.

6.1. FOURIER TRANSFORMS 223

6.1.1 Some Fourier Transform Shapes

Spoon feeding in the long run teaches us nothing

but the shape of the spoon.

E. M. Forster, British novelist, (1879–1970)

In this recipe, the Fourier transforms of some common f(x) are examined.

Calculate the Fourier transform F (k) of the following functions, using the

symmetric factor convention:

(a) f1 = Ae

−a

; (b) f2 = Ae

−a |x|

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

Plot the f(x) and their transforms for a= A= 1. Conﬁrm Parseval’s theorem.

The inttrans library package contains the Fourier transform command,

fourier, so is loaded. To accomplish the transforms, the assumption that both

a and A arepositiveisentered.

restart: with(inttrans): assume(a>0,A>0):

A functional operator F is formed to take the Fourier transform of an arbitrary

function f and simplify the result. The result is divided by

√

2π to agree with

the symmetric factor convention.

F:=f->simplify(fourier(f,x,k)/sqrt(2*Pi)):

The Gaussian function f1 is entered, and its Fourier transform F1 calculated

using the operator F. The transform is a Gaussian proﬁle centered on k =0.

f1:=A*exp(-aˆ2*xˆ2); F1:=F(f1);

f1 := Ae

(−a

)

F1 :=

(−

4 a

)

√

The symmetric exponentially decaying function f2 is entered and its transform

calculated. The resulting proﬁle F2 is a Lorentzian line shape centered on k=0.

f2:=A*exp(-a*abs(x)); F2:=F(f2);

f2 := Ae

(−a |x|)

F2 :=

√

+ k

)

√

Making use of the Heaviside command, the rectangular proﬁle f3 is entered,

and its transform F3 determined.

f3:=A*Heaviside(x+a)-A*Heaviside(x-a); F3:=F(f3);

f3 := A Heaviside(x + a) −A Heaviside(x −a)

F3 :=

A sin(ak)

√

To plot each function and its transform, an arrow operator G is introduced to

evaluate an arbitrary function FF at a=1, A =1.

G:=FF->eval(FF,{a=1,A=1}):

Making use of G, the three given functions f1 , f2,andf3 are plotted together,

the curves being given diﬀerent colors and line styles. The resulting picture is

shown on the left of Figure 6.1.

224 CHAPTER 6. INTEGRAL TRANSFORMS

plot([G(f1),G(f2),G(f3)],x=-4..4,color=[red,blue,green],

linestyle=[1,2,3],thickness=2,labels=["x","f"]);

–4

0.8

–8

Figure 6.1: Left: the three f(x). Right: Corresponding Fourier transforms.

The three Fourier transforms F1 , F2,andF3 are plotted with matching colors

and line styles and shown on the right of the ﬁgure.

plot([G(F1),G(F2),G(F3)],k=-10..10,color=[red,blue,green],

linestyle=[1,2,3],thickness=2,labels=["x","F"]);

The oscillatory curve is the transform of the rectangle f3 , the shorter smooth

curve the transform of the Gaussian f1 , and the taller smooth curve is the

transform F2 . The correspondence is more evident on the computer screen.

Using the complex conjugate (conjugate) command, an operator L is formed

to calculate the left-hand side of Parseval’s theorem for a given function f.

L:=f->int(f*conjugate(f),x=-infinity..infinity);

L := f →



∞

−∞

(f) dx

Similarly, an operator R is created to calculate the right-hand side of Parseval’s

theorem for a speciﬁed Fourier transform FT .

R:=FT->int(FT*conjugate(FT),k=-infinity..infinity);

R := FT →



∞

−∞

(FT ) dk

Using L and R, Parseval’s theorem is now conﬁrmed for the function f1 and its

transform F1 , the results L1 and R1 being identical.

L1:=L(f1); R1:=R(F1);

L1 :=

√

2 a

R1 :=

√

2 a

It is left as an exercise for you to conﬁrm that Parseval’s theorem holds for the

other two cases.

6.1. FOURIER TRANSFORMS 225

6.1.2 A Northern Weenie Roast

Wedding: the point at which a man stops toasting a woman

and begins roasting her.

Helen Rowland, American journalist, A Guide to Men,“Syncopations” (1922)

While visiting his inlaws in Northern Ontario over the Christmas holidays,

Russell went snowmobiling with his family. After a few hours, feeling slightly

cold and hungry, they stopped at one of the warming huts along the track and

soon had a roaring ﬁre going in the ﬁreplace. Using some long thin steel rods

which had been conveniently leaning against the outer wall of the hut, they

satisﬁed their hunger by roasting weenies in the ﬁre.

Idealizing the cooking situation, let’s calculate and animate the temperature

distribution inside a thin insulated semi-inﬁnite (0 ≤ x ≤∞) steel rod which

has the x=0 end held at a ﬁxed temperature T 0 and whose interior is initially

at 0 degrees Celsius. For steel, the thermal conductivity K =46 W/(m K

◦

), the

density ρ=8 × 10

kg/m

, and the speciﬁc heat C =500 J/(kg K

◦

The plots and integral transform packages are loaded, the former needed for

the animation, the latter for the Fourier sine (fouriersin) command.

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

The one-dimensional heat diﬀusion equation is entered in de, T(x, t)being

the temperature at a distance x from the hot end of the rod at time t.The

parameter a =



K/(ρC).

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

de :=

∂

∂t

T(x, t)=a

(

∂

∂x

T(x, t))

Since the temperature is speciﬁed at x = 0, i.e., T (0,t)=T 0, we will solve the

diﬀusion equation by ﬁrst taking the Fourier sine transform of de, yielding the

time-dependent ODE shown in the output of eq .

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

eq :=

∂

∂t

fouriersin(T(x, t),x,k)=

k (

√

2T(0,t) − k fouriersin(T(x, t),x,k)

√

π)

√

The notation is simpliﬁed by replacing fouriersin(T(x, t),x,k)withF(k, t).

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

eq2 :=

∂

∂t

F(k, t)=

k (

√

2T(0,t) − k F(k, t)

√

π)

√

The boundary condition T (0,t)=T0 is substituted into eq2 .

eq3:=subs(T(0,t)=T0,eq2);

eq3 :=

∂

∂t

F(k, t)=

k (

√

2 T0 − k F(k, t)

√

π)

√

226 CHAPTER 6. INTEGRAL TRANSFORMS

Then, the ODE eq3 is analytically solved for F(k, t).

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

eq4 := F(k, t)=

√

2 T0

√

+ e

(−a

F1(k)

Since T (x, 0) = 0 for x>0, then F (k, 0)= 0. So, the arbitrary function

F1(k)

is determined by evaluating the right-hand side of eq4 at time t=0, setting the

result equal to 0, and solving for

F1(k).

_F1(k):=solve(eval(rhs(eq4),t=0),_F1(k));

F1(k):=−

√

2 T0

√

F1(k) is automatically substituted into eq4 , the result being labeled eq5 .

eq5:=eq4;

eq5 := F(k, t)=

√

2 T0

√

−

(−a

√

2 T0

√

The temperature distribution inside the rod at arbitrary time t ≥ 0 follows on

applying the following Fourier sine command to the right-hand side of eq5 ,and

simplifying with the symbolic option.

T:=simplify(fouriersin(rhs(eq5),k,x),symbolic);

T := −T0 (−1+erf(

2 a

√

))

The temperature T is expressed in terms of the error function. As a partial check

on our calculation, note that for x>0andt =0, T(x, 0) = T0 (1 − erf(∞)) =

T0 (1 − 1) = 0 as it should. The thermal conductivity (K), density (ρ), and

speciﬁc heat (C) values for steel are now entered and the diﬀusion coeﬃcient

a =



K/(ρC) calculated. In the MKS system, a has the units m/s

1/2

.For

convenience, a is converted to CGS units by multiplying by 100, since 1 m =

100 cm. The units of a are then cm/s

1/2

, the distance x being expressed in cm.

K:=46: rho:=8*10ˆ3: C:=500: a:=100*sqrt(K/(rho*C));

a :=

√

The normalized temperature distribution T/T0 is given in TT,

TT:=T/T0;

TT := 1 −erf(

5 x

√

)

which is plotted at 60 second (1 minute) time intervals up to 5 minutes.

plot([seq(eval(TT,t=60*i),i=1..5)],x=0..25,thickness=2,

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

The resulting time evolution of the temperature inside the rod is shown in

Figure 6.2, the curve furthest to the left corresponding to 1 minute, the curve

furthest to the right to 5 minutes. Even at 5 minutes the temperature at x =25

cm is a tiny fraction of the temperature at the hot end.

6.1. FOURIER TRANSFORMS 227

0.2

0.4

0.6

0.8

T/T0

Figure 6.2: Normalized temperature distribution inside rod at diﬀerent times.

Clearly, the semi-inﬁnite approximation to the long rod is not too bad. The nor-

malized temperature distribution is now animated over the same time interval,

100 frames being used.

animate(TT,x=0..25,t=0..300,frames=100,thickness=2,

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

You will have to execute the recipe on your computer to see the animation.

6.1.3 Turn Oﬀ the Boob Tube and Concentrate

Television’s perfect. You turn a few knobs, a few of those mechanical

adjustments at which the higher apes are so proﬁcient, and lean back

and drain your mind of all thought. And there you are watching the

bubbles in the primeval ooze. You don’t have to concentrate. You

don’t have to react. You don’t have to remember. You don’t miss

your brain because you don’t need it.

Raymond Chandler, American novelist, (1888–1959)

Since you are reading this book, you have turned your back (at least tem-

porarily) on watching TV. Having your undivided attention, I would like you

to concentrate on the following example involving the Fourier cosine transform.

Consider a thin semi-inﬁnite (0 ≤ x ≤∞) hollow tube, closed at the x =0

end, containing water. A drop of blue ink is injected at the point x = x0 > 0

at time t = 0. Suppose that the drop has an initial concentration proﬁle given

approximately by C(x, 0) =f(x)=Aδ(x−x0). If the diﬀusion constant of ink in

water is d, determine the concentration of ink for arbitrary time t ≥ 0. At what

time T is the concentration of ink at x = 0 a maximum? Taking the nominal

values A =1, d =1, and x0 = 10, evaluate T and animate the concentration

proﬁle C(x, t) and plot C(0,t) over a time interval up to 2T .