Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

1.3 SUPPLEMENTARY RECIPES 57

G as much as possible and write it as a piecewise function. Check that G obeys

the general symmetry property G(x|ζ)=G(ζ|x). At what k values does G have

singularities where G diverges?

01-S10: A Potpourri of General Solutions

Classify the following LODEs as to order, homogeneneous or nonhomogeneous,

and constant or variable coeﬃcients. Obtain the general solution for each and

identify which part of the solution is the particular solution and which part is the

solution to the corresponding homogeneous ODE. Identify any non-elementary

functions in the solution. What method of attack does Maple use?

(a) xy



+(1+x) y = e

;

(b) y



+ xy =0;



+3y



+2y = e

;

(d) x



− 6 y = x

ln x;

(e) y



− (3/x) y



+ (15/(4 x

√

x) y =0;

(f) x



+(1− 2 α) xy



+(k

2 β

+(α

− p

)) y =0.

01-S11: Chebyshev Polynomials

Taking p =

√

1 − x

, q =0, w =1/

√

1 − x

,andλ = n

, show that the S-L

ODE yields the Chebyshev equation,(1− x

) y



− xy



+ n

y =0. For n a

non-negative integer, Chebyshev’s ODE is known to have polynomial solutions

(the Chebyshev polynomials T

(x)) deﬁned over the range x =−1 to +1. Using

the command expand(ChebyshevT(n,x)), generate the T

(x)forn = 0 to 10.

Then use the dsolve command to obtain the general solution of the Chebyshev

ODE and show that it can be expressed in terms of the T

(x). Demonstrate

that, e.g.,



−1

(x) T

(x) w(x) dx=0, i.e., the Chebyshev polynomials form a

set of orthogonal functions. Plot the ﬁrst ﬁve polynomials for x =−1to1.

01-S12: The Growing Pendulum

Consider a pendulum which consists of a point mass m at the bottom end of

a light supporting rod of length L which is allowed to move in a vertical plane

about a pivot point at its top end. Suppose that L increases at a steady rate,

i.e., L = L

+ vt,wherev>0 is a constant speed and t the time.

(a) Letting the rod make an angle θ(t) with the vertical and neglecting drag,

use Newton’s second law for angular motion (d(Iω)/dt = T ,whereI is

the moment of inertia, ω the angular velocity, and T the torque) to derive

the LODE for small θ. Solve the ODE for θ(t), given θ(0) = Θ,

θ(0) = 0.

Simplify θ(t) and identify the functions which occur.

(b) Taking L

=1 m, g =10 m/s

,Θ=π/6 rads, and v =0.5 m/s, plot θ(t)for

t = 0 to 100 seconds. Describe the behavior of the period. Then animate

the motion of the pendulum arm (representing it as a thick line) over this

time interval taking 100 frames and using constrained scaling.

58 CHAPTER 1. LINEAR ODES OF PHYSICS

01-S13: Another Green Function

Assuming that m = 0, derive the Green function G which is the solution to

(1 − x

) G



− 2 xG



− m

G/(1 − x

)=δ(x − z), subject to the boundary

conditions that G remains ﬁnite at the end points x= ±1 of the range. Simplify

G as much as possible and write it as a piecewise function. Plot G over the

range x=−1to1,takingm =3 and z =1/2.

01-S14: Going Green, Once Again

A light elastic string, under tension T andﬁxedatx=±L, is horizontal in the

absence of any applied forces. The string is embedded in an elastic membrane

which exerts a Hooke’s law (spring constant k) restoring force if the string is

displaced from equilibrium. If a force fe

−βx

is applied to the string, derive

the ODE for its shape y(x), assuming that y(x) is small. Using the Green

function method, determine y(x). The answer will involve the error function,

erf(z) ≡(2/

√



−t

dt. Check your result by solving the original ODE with

Maple’s dsolve command, subject to the given boundary condition. Taking

T =100 N, k =10 N/m, L=1 m, β =1/100 m

−2

,andf =50 N, plot y(x).

Chapter 2

Applications of Series

In this chapter, it is assumed that you are already familiar with the basic

concepts of inﬁnite series covered in standard calculus texts such as Calculus

by James Stewart [Ste87]. The emphasis here is on applications of series, the

topics being Taylor series, series solutions of LODEs, Fourier series, Legendre

and Bessel series, and summing series. Laurent series are covered in Chapter 5.

2.1 Taylor Series

The Taylor expansion of a function f(x)=f(x

+ h)aboutx

is given by

f(x)=f(x

)+hf







)+···. (2.1)

Thus, e.g., taking x

=0, so h=x − x

=x,

sin(x)=x −

−

+ ···=

∞



n=0

(−1)

2n+1

(2n + 1)!

≡

∞



n=0

Applying the ratio test, r ≡|u

n+1

|= x

/((2n + 3)(2n + 2)) → 0asn →∞.

Since r<1, the series converges for all x, i.e., the radius of convergence is ∞.

The extension to functions of more than one variable is straightforward, e.g.,

for a function f(x, y)=f(x

+ h, y

+ k),theTaylorseriesisoftheform,

f(x, y)=f(x



∂

∂x

∂

∂y



f(x



∂

∂x

∂

∂y



f(x

)+···.

2.1.1 Polynomial Approximations

Every man sees in his relatives, and especially in his cousins,

a series of grotesque caricatures of himself.

H. L. Mencken, American journalist, (1880–1956)

In this recipe, I will demonstrate how the Taylor series expansion may be used

to obtain the smallest polynomial approximation to the integral

f =



du, valid within ±0.000001 over the range 0 ≤ x ≤ 1.

60 CHAPTER 2. APPLICATIONS OF SERIES

The desired accuracy A is entered, and then the integral f is displayed and

evaluated in the following command line.

restart: A:=0.000001:

f:=Int(uˆ4*exp(uˆ2),u=0..x)=int(uˆ4*exp(uˆ2),u=0..x);

f :=



)

du =

)

−

)

−

√

π erf(xI)

Theanswerisgivenintermsoftheerror function,erf(z)=(2/

√

π)



−u

du,

with z = Ix=

√

−1 x. The result can be expressed in terms of real quantities by

converting the rhs to an imaginary error function,erﬁ(x)=(2/

√

π)



du.

f:=convert(rhs(f),erfi);

f :=

)

−

)

√

π erﬁ(x)

Then f is simpliﬁed by collecting exponential terms.

f:=collect(f,exp);

f := (

−

x) e

)

√

π erﬁ(x)

Although f can be evaluated numerically, a polynomial approximation to f is

now obtained. An operator F is formed which Taylor expands f in powers of

x about x=0, dropping terms of order x

N+1

F:=N->taylor(f,x=0,N+1):

A second arrow operator F2 is created to remove the order of term in the series.

F2:=N->convert(F(N),polynom):

Using F2(2*n+1), a sequence of polynomial approximations to the integral is

generated in T for n = 0 to 10 (not all polynomials are shown here). Note that

the Taylor series expansion of the integral generates only odd order polynomials,

hence the argument 2 n + 1. For plotting purposes, T is assigned.

T:=seq(t||(2*n+1)=F2(2*n+1),n=0..10); assign(T):

T := t1 =0, t3 =0, t5 =

, t7 =

, t9 =

t11 =

, ............................................................,

t21 =

312

1800

12240

95760

846720

The diﬀerence between the exact integral f and each polynomial for n =2

(yielding t5 ) and 10 (yielding t21 ) is plotted over the range x =0 to 1.

plot([seq(f-t||(2*n+1),n=2..10)],x=0..1,view=[0..1,0..A],

axes=box,thickness=2);

Since the diﬀerences turn out to be positive, the vertical view has been taken to

be between 0 and A, the upper limit setting the desired accuracy. The sequence

2.1. TAYLOR SERIES 61

of diﬀerences is shown in Figure 2.1, the diﬀerence corresponding to n =2 on

the far left and that corresponding to n = 10 on the far right. Since the ﬁrst

diﬀerence curve that remains completely within the vertical viewing range is

the curve corresponding to n =10, the smallest polynomial which is within the

desired accuracy over the range x=0to1isgivenbyt21 .

2e–07

4e–07

6e–07

8e–07

1e–06

0.2 0.4

0.6

0.8 1

Figure 2.1: Diﬀerences between exact integral and polynomial approximations.

2.1.2 Finite Diﬀerence Approximations

I’m not a teacher: only a fellow-traveler of whom you asked the way.

I pointed ahead–ahead of myself as well as you.

George Bernard Shaw, Anglo-Irish playwright (1856–1950)

Replacing ordinary and partial derivatives by their ﬁnite diﬀerence approxi-

mations (FDAs) is useful in numerically solving ODEs and PDEs as will be

illustrated in Chapter 9. These approximations are based on a Taylor series

approach. For example, an FDA to y



(x) can be obtained by Taylor expanding

y(x ± h)aboutx for small h, viz.,

y(x ± h)=y(x) ± hy



(x)+



(x) ±



(x)+O(h

), (2.2)

where O(h

) indicates that terms of order h

(and higher) are being neglected.

Forming the sum y(x + h)+y(x − h), and solving for y



, yields the FDA,



(x)=[y(x + h)+y(x − h) − 2y(x)]/h

+ O(h

). (2.3)

Eq. (2.3), known as the central diﬀerence approximation to the second deriva-

tive, was easily derived by hand. However, a computer algebra approach can

sometimes prove useful in deriving other FDAs as our mathematical physics

student, I. M. Curious, will now illustrate in answering the following problem.

(a) Show that to O(h

) an alternate FDA to the second derivative is given by



(x)=[−y(x+2h)+16y(x+h)−30y(x)+16y(x−h)−y(x−2h)]/(12h

62 CHAPTER 2. APPLICATIONS OF SERIES

Taking 12 digits accuracy, plot the diﬀerence between the FDA and the

exact y



for y = x

5.1

over the range x=0.1to10forh=0.1andh=0.05.

(b) Show that an FDA for f

xxyy

≡ ∂

f(x, y)/∂x

∂y

is given by

xxyy

=[f (x+h, y + k)+f(x − h, y+k)+f (x+h, y − k)+f(x − h, y − k)

−2(f(x+h, y)+f(x −h, y)+f (x, y+k)+f(x, y −k))+4 f (x, y)]/(h

and produce a 3-dimensional plot of the diﬀerence between the FDA and

the exact fourth derivative for f =sin(x − y) e

−x

. Take h = k =0.05,

12 digits accuracy, and a plotting range x= −3...3, y =−3...3.

I. M. begins her recipe for part (a) by setting the accuracy to 12 digits.

restart: Digits:=12:

Noting that the wording of the problem implies that terms of order h

have

been dropped in the Taylor series, she creates an arrow operator t to Taylor

expand y(x + ah)inpowersofh about an arbitrary point x, dropping terms of

O(h

). The order of term is removed by enclosing the taylor command with

convert( ,polynom). The quantity a is an integer which must be supplied.

t:=a->y(x+a*h)=convert(taylor(y(x+a*h),h,6),polynom):

Asatest,I.M.choosesa =2 and generates the Taylor series for y(x +2h).

t(2);

y(x +2h)=y(x)+2D(y)(x) h +2(D

(2)

)(y)(x) h

(3)

)(y)(x) h

(4)

)(y)(x) h

(5)

)(y)(x) h

In the output, (D

(4)

)(y)(x), for example, stands for the 4th derivative, y



(x).

Making use of the operator t, the proposed FDA is entered in eq1 .

eq1:=-t(2)+16*t(1)-30*t(0)+16*t(-1)-t(-2);

eq1 := −y(x +2h)+16y(x + h) −30 y(x)+16y(x − h) − y(x − 2 h)

=12(D

(2)

)(y)(x) h

On the rhs of the output, only the second derivative term survives, all other

terms up to O(h

) having canceled. Using the isolate command, the second

derivative (entered as D[1,1](y)(x))isisolatedtothelhsofeq2 ,

eq2:=isolate(eq1,D[1,1](y)(x));

eq2 := (D

(2)

)(y)(x)

= −

y(x +2h) − 16 y(x + h)+30y(x) − 16 y(x − h)+y(x − 2 h)

thus conﬁrming the FDA.

An arrow operator y is formed to evaluate y = x

5.1

at an arbitrary x value.

So that the FDA can be evaluated for diﬀerent h values, the rhs of eq2 is turned

into a functional operator A, depending on h, using the unapply command.

y:=x->xˆ5.1: A:=unapply(rhs(eq2),h):

The exact y



(x) is calculated and the FDAs generated for h=0.1andh=0.05.

2.1. TAYLOR SERIES 63

exact:=diff(y(x),x,x); approx1:=A(0.1): approx2:=A(0.05):

exact := 20.91 x

3.1

The diﬀerences between the FDAs approx1 and approx2 and exact are plotted

for x=0.1 to 10 as solid red and dashed green curves, respectively, the resulting

plot being shown in Figure 2.2.

plot([approx1-exact,approx2-exact],x=0.1..10,color=

[red,green],thickness=2,linestyle=[1,3],tickmarks=[2,2]);

–5e–05

510

Figure 2.2: FDA minus the exact result for h=0.1 (solid) and h=0.05 (dashed).

As expected, the accuracy of the FDA is improved with the smaller value of h.

Both curves have greatest accuracy over a limited range of x.

Now, I. M. tackles part (b) of the problem. The plots library package is

loaded and 12 digits accuracy entered.

restart: with(plots): Digits:=12:

The FDA involves a Taylor expansion in two variables. I. M. uses the multivari-

ate Taylor series command, mtaylor, to create an operator t for performing the

expansion of f(x + ah, y + bk)inpowersofh and k, neglecting terms of sixth

order. mtaylor does not generate the order of term, so convert( ,polynom)

needn’t be applied. Two integers a and b must be entered as arguments in t.

t:=(a,b)->f(x+a*h,y+b*k)=mtaylor(f(x+a*h,y+b*k),[h,k],6):

The given FDA is entered in eq1 .

eq1:=t(1,1)+t(-1,1)+t(1,-1)+t(-1,-1)

-2*(t(1,0)+t(-1,0)+t(0,1)+t(0,-1))+4*t(0,0);

eq1 := f (x + h, y + k)+f (x − h, y + k)+f (x + h, y −k)+f (x − h, y −k)

− 2 f (x + h, y) − 2 f (x − h, y) − 2 f (x, y + k) − 2 f (x, y − k)

+4f (x, y)=h

1, 1, 2, 2

(f)(x, y)

64 CHAPTER 2. APPLICATIONS OF SERIES

Only a single term involving D

1, 1, 2, 2

results on the rhs of the output, the

subscripts 1 and 2 denoting derivatives with respect to x and y, respectively.

Dividing eq1 by h

, and simplifying, yields the desired FDA.

eq2:=simplify(eq1/(hˆ2*kˆ2));

eq2 := (f(x + h, y + k)+f(x −h, y + k)+f(x + h, y − k)+f(x − h, y −k)

− 2f(x + h, y) − 2f(x −h, y) − 2f(x, y + k) − 2f(x, y −k)

+4f(x, y))/(h

)=D

1, 1, 2, 2

(f)(x, y)

The function f(x, y)=sin(x − y) e

−x

is entered, as well as h= k =0.05.

f:=(x,y)->sin(x-y)*exp(-xˆ2*yˆ2): h:=0.05: k:=0.05:

The exact 4th derivative is calculated and the FDA given by the lhs of eq2 .

exact:=diff(f(x,y),x,x,y,y); approx:=lhs(eq2):

The plot3d command is used to create a 3-dimensional plot of approx minus

exact for x = −3...3, y = −3...3. The grid is taken to be 25 ×25, the shad-

ing to be XYZ (i.e., the color varies in the 3 directions), the axes boxed, the

plot “illuminated” with a light source at a certain orientation given by enter-

ing lightmodel=light2, the number of plotting points to be 2000, and the

orientation speciﬁed. The resulting picture is shown in Figure 2.3.

plot3d(approx-exact,x=-3..3,y=-3..3,grid=[25,25],

shading=XYZ,axes=box,lightmodel=light2,numpoints=2000,

tickmarks=[3,3,3],orientation=[20,40]);

–2

–0.4

–0.2

0.2

0.4

Figure 2.3: Diﬀerence between FDA and exact 4th derivative.

The regions where the FDA and exact result diﬀer appreciably can be seen as

a series of ridges and valleys. The 3-dimensional plot can be rotated on the

computer screen by dragging with the mouse.

2.2. SERIES SOLUTIONS OF LODES 65

2.2 Series Solutions of LODEs

Consider an nth-order homogeneous linear ODE with variable coeﬃcients,

+ a

n−1

(x)

n−1

+ ···+ a

(x)

+ a

(x) y =0. (2.4)

If a

(x), ..., a

n−1

(x)areregular (single-valued and analytic) at a point x = x

then x

is referred to as an ordinary point of the ODE. Near x

, the general

solution of the LODE can be written as a Taylor series, y =



∞

m=0

(x −x

)

whose radius of convergence is the distance to the nearest singular point (a

non-ordinary point) of the ODE. The coeﬃcients c

are found by substituting

the series into the ODE and setting the coeﬃcients of each power of x to zero.

If x

is not an ordinary point, but (x − x

) a

n−1

(x), (x − x

)

n−2

(x),...,

(x −x

)

(x) are regular at x

,thenx

is a regular singular point.Nearsuch

an x

,aFrobenius series solution of the form y =(x−x

)



∞

m=0

(x−x

)

with c

= 0 and s not necessarily an integer, can always be found. The radius

of convergence is again the distance to the nearest singular point outside x

In this section, series solutions to two familiar Sturm–Liouville ODEs are

obtained, illustrating the two types of expansions mentioned above.

2.2.1 Jennifer Renews an Old Acquaintance

We need two kinds of acquaintances, one to complain to,

while to the others we boast.

Logan Pearsall Smith, American essayist, Afterthoughts,“Other People”, 1931

In recipe 01-2-1, Jennifer introduced the Legendre functions. In this recipe,

she renews her acquaintance with these important special functions, showing

us how they arise as a series solution to the Legendre ODE,

(1 − x

) y



− 2 xy



+ n (n +1)y =0. (2.5)

Casting (2.5) into the “standard” form (2.4), Equation (2.5) has regular singular

points at x = ±1. Jennifer seeks a Taylor series solution about the ordinary

point x

= 0. She decides that it will suﬃce for calculation purposes to work

with a ﬁnite number of terms, say N = 7, in the series.

restart: N:=7:

The left-hand side of Legendre’s equation is now entered.

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

ode := (1 − x

)(

y(x)) −2 x (

y(x)) + n (n +1)y(x)

The power series y(x)=



m=0

is inputted using the add command.

y(x):=add(c||m*xˆm,m=0..N);

y(x):=c0 + c1 x + c2 x

+ c3 x

+ c4 x

+ c5 x

+ c6 x

+ c7 x

y(x) is automatically substituted into ode and powers of x are collected.

66 CHAPTER 2. APPLICATIONS OF SERIES

ode2:=collect(ode,x);

ode2 := (−56 c7 + n (n +1)c7 ) x

+(−42 c6 + n (n +1)c6 ) x

+ (42 c7 − 30 c5 + n (n +1)c5 ) x

+(n (n +1)c4 − 20 c4 +30c6 ) x

+(−12 c3 + n (n +1)c3 +20c5 ) x

+(−6 c2 +12c4 + n (n +1)c2 ) x

+(6c3 + n (n +1)c1 − 2 c1 ) x +2c2 + n (n +1)c0

Since ode2 represents the lhs of the Legendre equation, it must be set equal

to zero to form the complete equation. But x is arbitrary, so the coeﬃcient of

each power of x must separately be equal to zero. A functional operator eq is

created to set the coeﬃcient of x

in ode2 equal to zero.

eq:=m->coeff(ode2,x,m)=0:

The set of equations eq(m) is solved for the c||(m+2) for m=0 to N −2.

sol:=solve({seq(eq(m),m=0..N-2)},{seq(c||(m+2),m=0..N-2)}):

The solution sol is assigned, the coeﬃcients factored and the series formed,

assign(sol): y:=add(factor(c||m)*xˆm,m=0..N):

and the coeﬃcients of c0 and c1 collected in y.

y:=collect(y,{c||0,c||1});

y := (1 −

n (n +1)x

n (n − 2) (n +3)(n +1)x

−

n (n − 2) (n − 4) (n +5)(n +3)(n +1)x

720

)c0

+(x −

(n +2)(n − 1) x

(n − 1) (n − 3) (n +4)(n +2)x

120

−

(n − 1) (n − 3) (n − 5) (n +6)(n +4)(n +2)x

5040

)c1

y gives the ﬁrst few term in the general series solution of Legendre’s ODE.

It consists of an even (coeﬃcient of c0 ) and odd (coeﬃcient of c1 ) series in

x. Noting that the denominator of the x

term is just m! (e.g., for m =6,

6! = (6)(5)(4)(3)(2)(1) = 720), the structure of higher order terms in the inﬁnite

series is easy to deduce. In deriving y, Jennifer has mimicked a hand calculation.

The same result could be more easily obtained by using the series option in

the dsolve command, as Jennifer will now demonstrate. She unassigns y and

sets the order of the ﬁrst term to be neglected in the series solution. Since she

has taken N =7,theorderhereis8,i.e.,termsofO(x

)aredroppedinthe

series. If the order is not speciﬁed, the default is to neglect terms of O(x

unassign(’y’): Order:=N+1;

Order := 8

Then ode is solved for y(x), subject to the initial condition y(0) = c0 and



(0)= c1 , using the series option in the dsolve command.

Y:=dsolve({ode,y(0)=c||0,D(y)(0)=c||1},y(x),series);