Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

9.2. PARTIAL DIFFERENTIAL EQUATIONS 359

the end points are included in the list of plotting points. The points will be

joined by straight lines.

for n from 1 to numplots do

T:=Aˆ(20*n) . v;

p[n]:=plot([[0,0],seq([i/(N+1),T[i]],i=1..N),[1,0]],

tickmarks=[2,2],labels=["x","y"]);

end do:

The cpu time for the do loop is about 20 s.

cpu:=time()-begin;

cpu := 20.037

Making use of the display command with the insequence=true option, the

heat diﬀusion process in the rod is animated. The opening frame of the anima-

tion is shown in Figure 9.10. Execute the recipe to see what happens.

plots[display]([seq(p[i],i=1..numplots)],insequence=true,

view=[0..L,0..vmax+1]);

Figure 9.10: Opening frame of the heat diﬀusion animation.

9.2.3 Von Neumann Stability Analysis

No civilization would ever have been possible

without a framework of stability....

Hannah Arendt, American political philosopher, (1906–75)

That the explicit scheme used in the heat diﬀusion equation recipe 09-2-2 is

conditionally stable, i.e., stable for r =k/h

≤ 1/2, could be proved by executing

the recipe with varying values of r and looking at the numerical output. Alter-

nately, one can use the Von Neumann stability analysis as will be demonstrated

in the following recipe.

360 CHAPTER 9. NUMERICAL METHODS

Recall that the relevant algorithm is of the structure

m,n+1

= rT

m−1,n

+(1− 2 r) T

m,n

+ rT

m+1,n

(9.21)

with m and n indexing the spatial and time steps, respectively. Assume that

m,n

= T

m,n

+ U

m,n

,whereT

is the exact solution of the diﬀerence scheme

and U represents a small numerical error due to roundoﬀ, etc. Substituting this

form into (9.21) yields an identical ﬁnite diﬀerence scheme for U . The question

is “Does U grow or decay with time?”, being unstable (stable) if it grows (de-

cays). To answer this question, we assume that U

m,n

can be represented by the

following exponential form,

m,n

= e

imθ

inλ

, (9.22)

with θ real and λ= α + iβ. Itisassumedthatbothα and β are real.

restart: assume(alpha::real,beta::real):

An operator U is formed to calculate U

m,n

for arbitrary subscripts m and n.

U:=(m,n)->exp(I*m*theta)*exp(I*n*lambda);

U := (m, n) → e

(mθI)

(nλI)

Eq. (9.21) is entered in terms of U, each term being automatically evaluated.

eq:=U(m,n+1)=r*U(m-1,n)+(1-2*r)*U(m,n)+r*U(m+1,n);

eq := e

(mθI)

((n+1) λI)

= re

((m−1) θI)

(nλI)

+(1 − 2 r) e

(mθI)

(nλI)

+ re

((m+1) θI)

(nλI)

Then, eq is divided by U

m,n

and simpliﬁed.

eq2:=simplify(eq/U(m,n));

eq2 := e

(λI)

=2r cos(θ)+1− 2 r

We substitute λ= α + iβ into eq2 .

eq3:=subs(lambda=alpha+I*beta,eq2);

eq3 := e

((α+βI) I)

=2r cos(θ)+1− 2 r

The absolute values of the lhs and rhs of eq3 are equated, and simpliﬁed.

eq4:=simplify(abs(lhs(eq3)))=abs(rhs(eq3));

eq4 := e

(−β)

= |2 r cos(θ)+1− 2 r|

For stability, β ≥ 0 is required, because if β<0, then U

m,n

∼ e

|β|n

diverges as

n (i.e., time) increases. For β ≥ 0, e

−β

≤ 1, so we must have the rhs of eq4 less

than or equal to 1 for stability. The rhs of eq4 is turned into an operator g in

terms of r by using the unapply command.

g:=unapply(rhs(eq4),r);

g := r →|2 r cos(θ)+1− 2 r|

Using g,therhsofeq4 is plotted over the range θ =−2.5 π to 2.5 π for r =0.49

and 0.51, being represented on the computer screen by red and blue curves,

respectively. A horizontal green line is also plotted at a vertical height 1. The

corresponding black and white picture is shown in Figure 9.11.

9.2. PARTIAL DIFFERENTIAL EQUATIONS 361

plot([g(.49),g(.51),1],theta=-2.5*Pi..2.5*Pi,

color=[red,blue,green],numpoints=1000);

–6

theta

Figure 9.11: The r =0.51 curve > 1forsomeθ,ther =0.49 curve doesn’t.

Since the r =0.51 curve exceeds 1 for some ranges of θ, the numerical scheme is

unstable for this r value. The r =0.49 curve doesn’t exceed 1, so the algorithm

is stable. By varying r between 0.49 and 0.51, you may use the recipe to check

that the critical value of r for numerical stability or instability is 1/2.

9.2.4 Sometimes It Pays to be Backwards

Diaper backward spells repaid. Think about it.

M. McLuhan, Canadian communications theorist, Vancouver Sun, June 1969

An unconditionally stable algorithm can be obtained for ∂T/∂y = ∂

T/∂x

approximating the y derivative with the BWDA. In this case, on still replacing

the x derivative with the standard CDA, the numerical algorithm becomes

(1 + 2 r) T

i,j

− rT

i+1,j

− rT

i−1,j

= T

i,j−1

. (9.23)

Let’s use equation (9.23) to resolve the example given in recipe 09-2-2,which

involved determining the evolution of the temperature distribution in a thin

rodoflengthL = 1 m, with boundary conditions T (0,y)=T (L, y) = 0, and an

initial proﬁle T (x, 0)= f(x) =200 x

(L − x).

We take M = 30 spatial steps, so h = L/M =1/30. Choosing the time (y)

step to be k =0.002 produces a ratio r = k/h

=1.8. Not only is r larger than

1/2, it’s greater than 1! Nevertheless, as you will see, the scheme remains stable,

although as r is increased the algorithm becomes increasingly less accurate. The

total number of time steps is taken to be N =250.

restart: with(plots):

L:=1: M:=30: h:=L/M; k:=0.002; r:=k/hˆ2; N:=250:

362 CHAPTER 9. NUMERICAL METHODS

h :=

k := 0.002 r := 1.800

Thevaluesofx

=ihare calculated in Xcoords for i=0 to M, and the sequence

of numbers assigned.

Xcoords:=seq(x[i]=i*h,i=0..M): assign(Xcoords):

A functional operator f is introduced for generating the initial temperature

value at the ith spatial mesh point.

f:=i->evalf(200*(x[i]ˆ3*(L-x[i]))):

Using f, the sequence of initial T values at the grid points is calculated.

ic:=seq(T[i,0]=f(i),i=0..M):

The two boundary conditions are entered,

bc1:=seq(T[0,j]=0,j=0..N): bc2:=seq(T[M,j]=0,j=0..N):

and are assigned, along with the initial condition.

assign(ic,bc1,bc2):

An operator F is formed for implementing the algorithm (9.23).

F:=(i,j)->(1+2*r)*T[i,j]-r*T[i+1,j]-r*T[i-1,j]-T[i,j-1];

F := (i, j) → (1 + 2 r) T

i, j

− rT

i+1,j

− rT

i−1,j

− T

i, j−1

The start time for executing the algorithm in the do loop is recorded.

begin:=time():

We now generate and numerically solve the algorithm equations, F (i, j)= 0, at

each of the internal grid points on each of the j time steps, for j =1 to N .

forjfrom1toNdo

eqs[j]:={seq(F(i,j)=0,i=1..M-1)};

sol[j]:=fsolve(eqs[j],{seq(T[i,j],i=1..M-1)});

assign(sol[j]):

end do:

The cpu time for the do loop is now displayed.

cpu:=time()-begin;

cpu := 14.110

A graphing operator gr is created for producing the sequence (i =0 to M )of

plotting points (x

, T

i,j

)onthejth time step and plotting them.

gr:=j->plot([seq([x[i],T[i,j]],i=0..M)]):

The temperature diﬀusion process is animated with the display command,

the insequence=true option being included. Click on the computer plot and

the start arrow in the tool bar to see the animation.

display([seq(gr(j),j=0..N)],insequence=true,labels=["x","T"]);

9.2. PARTIAL DIFFERENTIAL EQUATIONS 363

9.2.5 Daniel Still Separates, I Now Iterate

The ultimate result of shielding men from the eﬀects of folly

is to ﬁll the world with fools.

Herbert Spencer, British philosopher, sociologist, (1820–1903)

To numerically solve the 1-dimensional wave equation, we can replace each

second derivative with a CDA, viz.,

∂

∂x

∂

∂t

=⇒

i+1,j

− 2 U

i,j

+ U

i−1,j

i,j+1

− 2 U

i,j

+ U

i,j−1

=⇒ U

i,j+1

=2(1−r) U

i,j

+ r (U

i+1,j

+ U

i−1,j

) −U

i,j−1

,r=(ck/h)

. (9.24)

In this case, the unknown U

i,j+1

depends explicitly on U values on the previous

two steps, as shown schematically in Figure 9.12, where x is horizontal and

t vertical. To start the iteration scheme (9.24), we must therefore know the

UU U

i−1,j i,j i+1,j

i, j+1

i, j −1

Figure 9.12: Mesh points involved in wave equation algorithm (9.24).

values of U on both the j = 0 and j = 1 time steps. To see how we deal with

this issue, let’s consider the following possible scenario.

Young Daniel is playing with a light horizontal string of length L ﬁxed at

x= 0 and L, i.e., U(0,t)=U(L, t)= 0. Before he manages to separate it into its

individual strands, his older cousin Justine cleverly manages to give it an initial

proﬁle U(x, 0) = f(x)=Ax

(L − x) and transverse velocity

U(x, 0) = g(x)=

Bx(L −x)

. Our task is to iterate the numerical algorithm (9.24) to determine

the subsequent motion of the string and then animate the string vibrations,

taking L=1 m, A=10 m

−3

, B =5 m

−2

,andc=3 m/s.

On the zeroth time row, we have U

i,0

=f(x

). Using the forward diﬀerence

approximation, the U values on the ﬁrst time row, i.e., U

i,1

, are found thus:

364 CHAPTER 9. NUMERICAL METHODS

∂U(x

, 0)

∂t

i,1

− U

i,0

= g(x

), so U

i,1

= f(x

)+kg(x

)+O(k

). (9.25)

An even better approximation for U

i,1

, which we shall use in our recipe, is

i,1

= f(x

)+kg(x

[f(x

i+1

) − 2 f(x

)+f(x

i−1

)] + O(k

). (9.26)

To begin the recipe, the values of L, c, A,andB are entered. The total time is

taken to be T =2 s, and the number of time and spatial steps are N = 200 and

M =20, respectively.

restart: with(plots):

L:=1: T:=2: M:=20: N:=200: c:=3: A:=10: B:=5:

The spatial step h= L/M, time step k = T/N,andr =(ck/h)

are calculated.

h:=evalf(L/M); k:=evalf(T/N); r:=evalf((c*k/h)ˆ2);

h := 0.05000000000 k := 0.01000000000 r := 0.3600000000

Here, we have r =0.36. Using the Von Neumann stability analysis, it can be

shown that the condition for stability is that r ≤ 1, referred to as the Courant

stability condition.

The boundary conditions U

0,j

= 0 and U

M,j

=0 for j =0 to N are entered

in bc1 and bc2, which are both assigned.

bc1:=seq(U[0,j]=0,j=0..N): assign(bc1):

bc2:=seq(U[M,j]=0,j=0..N): assign(bc2):

Functional operators f and g are introduced for calculating f(x)andg(x)at

the internal grid points x

=ih.

f:=x->evalf(A*xˆ3*(L-x)): g:=x->evalf(B*x*(L-x)ˆ2):

The sequence of values U

i,0

= f(ih) are generated in ic1 for i =1...M − 1and

the result assigned. Similarly, the sequence of values U

i,1

are generated using

(9.26) in ic2, and the result also assigned.

ic1:=seq(U[i,0]=f(i*h),i=1..M-1): assign(ic1):

ic2:=seq(U[i,1]=f(i*h)+k*g(i*h)

+(r/2)*(f((i+1)*h)-2*f(i*h)+f((i-1)*h)),i=1..M-1):

assign(ic2):

Thetimeatthebeginningofthedoloopisrecorded.

begin:=time():

The algorithm (9.24) is iterated from j =1 to N −1 and the solution assigned.

for j from 1 to N-1 do

sol:=seq(U[i,j+1]=2*(1-r)*U[i,j]+r*(U[i+1,j]+U[i-1,j])

-U[i,j-1],i=1..M-1);

assign(sol):

end do:

The cpu time for the do loop is found to be about 0.4 s.

cpu:=time()-begin;

9.2. PARTIAL DIFFERENTIAL EQUATIONS 365

cpu := 0.404

An operator gr is formed for plotting the spatial proﬁle on the jth time step.

gr:=j->plot([seq([i*h,U[i,j]],i=0..M)],labels=["x","U"]):

The motion is animated by including insequence=true in the display com-

mand. The opening frame of the animation is shown in Figure 9.13.

display(seq(gr(j),j=0..N),insequence=true);

–1

–0.5

0.5

0.2 0.4

0.6

0.8 1

Figure 9.13: Opening frame of the vibrating string animation.

9.2.6 Interacting Laser Beams

The exercise of power is determined by thousands of interactions

between the world of the powerful and that of the powerless,

al l the more so because these worlds are never divided by

a sharp line: everyone has a small part of himself in both.

Vclav Havel, Czech playwright, president, (b. 1936)

The interaction of two intense laser pulses as they pass through each other in

opposite directions in a certain resonant absorbing ﬂuid can be described [RE76]

by the following normalized PDEs for the laser intensities u and v:

∂u

∂x

∂u

∂y

= −g

uv− bu,

∂v

∂x

−

∂v

∂y

= −g

uv+ bv. (9.27)

Here x is the distance inside the ﬂuid medium, y the time, g

> 0andg

> 0

are the “gain” coeﬃcients, and b ≥ 0 the absorption coeﬃcient. The u (v)beam

travels in the positive (negative) x direction.

Our goal is to numerically solve this set of nonlinear PDEs for the intensities,

using the method of characteristics, which is based on ﬁnding “characteristic”

directions along which the PDEs can be reduced to ODEs. Since u=u(x, y)and

v = v(x, y), then du=(∂u/∂x) dx+(∂u/∂y) dy and dv =(∂v/∂x) dx+(∂v/∂y) dy.

366 CHAPTER 9. NUMERICAL METHODS

Making use of these results, equations (9.27) can be rewritten as

+(1−

)

∂u

∂y

= f

− (1 +

)

∂v

∂y

= f

, (9.28)

with f

≡−g

uv− bu and f

≡−g

uv+ bv. If we work along a line whose

slope is dy/dx =1, then du/dx = f

, while along a line of slope dy/dx = −1

one has dv/dx = f

. The characteristic directions dy/dx=±1 form a diamond-

shaped grid as in Figure 9.14, the grid spacing ∆y =∆x =h. To solve for u(x, y)

∆x=h

L R

∆y=h

j=0

j=1

j=2

i=0 i=1 i=2

(1,1) (3,1) (5,1)

(2,0) (4,0) (6,0)

(0,2) (2,2) (4,2)

x=0 x=L=1

Figure 9.14: Diamond-shaped grid for solving (9.27).

and v(x, y), we impose the following initial and boundary conditions:

• u(x, 0)=v(x, 0)=0 for 0 <x<L=1, i.e., no pulses initially inside ﬂuid,

• u(0,y)=v(1,y)=sin(2πy) for 0 ≤ y ≤ y

=1/2andzerofory>y

, i.e.,

positive half-sine wave pulses fed in at x=0 and x= L=1.

Starting on the bottom time row, we move along the characteristic directions

dy/dx = 1 (e.g., from the point L to P )anddy/dx = −1(fromR to P )in

calculating the changes in u and v, respectively. Using, say, the Euler approxi-

mation, the values of u and v at P are given by

= u

+ h (f

)

= v

− h (f

)

. (9.29)

Taking, say, g

=0.4, g

= 20, b =0.5, we implement (9.29) in the following

recipe and animate the numerical result. The ﬂuid sample length L =1 is

divided into M = 120 steps, so that the x (or y) stepsize is h = L/M =1/120.

9.2. PARTIAL DIFFERENTIAL EQUATIONS 367

The leading edge of each pulse will just make contact in N = M/2 = 60 time

steps. If we let i and j index the x and y steps, as in Figure 9.14, the grid

coordinates for j =0, 2, ... are (0,j), (2,j), (4,j), ..., while for j =1, 3, ... they

are (1,j), (3,j), (5,j), ... For later convenience, we set m=2i and n =2j.

restart: with(plots): begin:=time():

L:=1: M:=120: h:=L/M; N:=M/2; m:=2*i: n:=2*j:

h :=

120

N := 60

Thevaluesofg

, g

,andb are entered. A functional operator S is introduced

to calculate the input laser shape on the boundaries.

g1:=0.4; g2:=20; b:=0.5; S:=j->sin(n*evalf(Pi)/N);

g1 := 0.4 g2 := 20 b := 0.5 S := j → sin(

n evalf(π)

)

Using S, the boundary conditions are applied at i=0 and M in bc1 and bc2.

bc1:=seq(u[0,n]=S(j),j=0..N/2),seq(u[0,n]=0,j=N/2+1..3*N):

bc2:=seq(v[M,n]=S(j),j=0..N/2),seq(v[M,n]=0,j=N/2+1..3*N):

The initial condition is entered and assigned along with bc1 and bc2.

ic:=seq(u[m,0]=0,i=0..M/2),seq(v[m,0]=0,i=0..M/2):

assign(bc1,bc2,ic):

Functional operators are formed to evaluate f

and f

for a given i, j.

f1:=(i,j)->-g1*u[i,j]*v[i,j]-b*u[i,j]:

f2:=(i,j)->-g2*u[i,j]*v[i,j]+b*v[i,j]:

A double do loop is used to iterate (9.29). The conditional if...then statement

is included so that i starts at i0=0 for j =0, 2, 4, ... and i0=1 for j =1, 3, ....

for j from 0 to 3*N do

if j mod 2 = 0 then i0: = 0 else i0:=1 end if;

for i from i0 to M by 2 do

u[i+1,j+1]: = u[i,j]+h*f1(i,j);

v[i-1,j+1]: = v[i,j]-h*f2(i,j);

end do: end do:

An operator gr is formed to plot u and v on the jth time step.

gr:=j->plot([[seq([m*h,u[m,n]],i=0..M/2)],[seq([m*h,v[m,n]],

i=0..M/2)]],color=[red,blue],labels=["x","u,v"],thickness=2):

The numerically obtained proﬁles are animated with the display command

with the insequence=true option. Click on the computer plot and the start

arrow to initiate the animation.

display([seq(gr(j),j=0..2*N)],insequence=true);

The cpu time for the entire recipe is given below.

cpu:=time()-begin;

cpu := 7.752

368 CHAPTER 9. NUMERICAL METHODS

9.2.7 KdV Solitons

In this recipe, we shall numerically investigate the collision of a taller faster

solitary wave overtaking a shorter slower one. The standard algorithm for nu-

merically integrating the KdV equation (8.25) is obtained as follows. Returning

to the general Taylor expansion (9.3) of y(x ± h), adding the plus and minus

results yields the central diﬀerence approximation to the ﬁrst derivative, viz.,



y(x + h) − y(x − h)

2 h

+ O(h

). (9.30)

This approximation is used for the ﬁrst derivatives in (8.25), viz.,

∂U

∂t

⇒

i,j+1

− U

i,j−1

)

2 k

∂U

∂x

⇒

i+1,j

− U

i−1,j

)

2 h

where k and h are the time and spatial stepsizes, respectively.

A central diﬀerence approximation to the third derivative can be obtained

by modifying the ﬁrst part of Recipe 02-1-2 as follows. A functional operator

t is formed to Taylor expand y(x + ah) to 5th order in h for arbitrary a.The

order of term is removed with the convert( ,polynom) command.

restart:

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

Then, using the operator, the sum y(x+2h)−2 y(x+h)+2y(x−h)−y(x−2 h)

is calculated in eq1 , yielding 2 h



on the rhs. The error here is O(h

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

eq1 := y(x +2h) − 2 y(x + h)+2y(x − h) −y(x − 2 h)=2(D

(3)

)(y)(x) h

A ﬁnite diﬀerence approximation with error O(h

) follows for the third deriva-

tive on using the isolate command.

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

eq2 := (D

(3)

)(y)(x)=−

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

So, we will make the central diﬀerence approximation,

∂

∂x

⇒

i+2,j

− 2 U

i+1,j

+2U

i−1,j

− U

i−2,j

)

2 h

Finally, let’s use the central diﬀerence approximation (U

i+1,j

+ U

i,j

+ U

i−1,j

)/3,

which also has an error O(h

), for the U term in the KdV equation. Putting

all these approximations together, the algorithm for the KdV equation is

i,j+1

= U

i,j−1

−

i+1,j

+ U

i,j

+ U

i−1,j

)

i+1,j

− U

i−1,j

)

−

i+2,j

− 2U

i+1,j

+2U

i−1,j

− U

i−2,j

). (9.31)

Using the Von Neumann stability analysis, it can be shown that this algorithm

is numerically stable for k/h

< 2/(3

√

3)= 0.3849.

Note that the algorithm connects the j + 1 time step to steps j and j −1, so

to iterate (9.31) the U values must be known on the j = 0 and j =1 time steps.