Enns R.H. Computer Algebra Recipes for Mathematical Physics

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

8.3. NONLINEAR ODES: APPROXIMATE METHODS 319

E =



(t) dt, (8.22)

over a range a to b of interest. One then tries to choose a Φ which minimizes the

total error. As in the Rayleigh–Ritz method, one or more adjustable parameters

may be included and the error is minimized with respect to these parameters.

As an example, let’s apply the Ritz method to the hard spring equation,

¨x + x + x

=0, (8.23)

with x(0) ≡ A =2, ˙x(0) = 0. Although we have already solved the hard spring

equation analytically in terms of Jacobian elliptic functions, the Ritz solution

will be compared to the numerical solution. As a trial function, let’s take

Φ=C

cos(ωt)+(A −C

) cos(3 ωt), (8.24)

with C

and ω adjustable parameters. At t=0, Φ(0)= A,

Φ(0)= 0, so the initial

conditions are satisﬁed.

The value of A is entered along with the upper limit T = n 2 π/|ω| in the

time range t =0 to T “of interest” that will be used for plotting purposes. For

sake of deﬁniteness, let’s take n=4, sothat T =8π/|ω|,withω to be determined.

restart: with(plots): A:=2: n:=4: T:=n*2*Pi/abs(omega):

The trial function Φ is entered.

Phi:=C1*cos(omega*t)+(A-C1)*cos(3*omega*t);

Φ:=C1 cos(ωt)+(2−C1) cos(3 ωt)

A functional operator e for generating the residual of the NLODE (8.23) is

formed. Entering e(Phi) then will produce the residual.

e:=x->diff(x,t,t)+x+xˆ3;

e := x → (

x)+x + x

The total error is calculated, the assumption that ω>0 being provided to

facilitate the integration. A lengthy expression in terms of C1 and ω results.

Error:=int(e(Phi)ˆ2,t=0..T) assuming omega>0;

Error := π(1256 C1

+ 4032 ω

C1 − 1376 C1 + 2416 C1

+ 2592 ω

− 390 C1

− 2304 ω

+55C1

+ 544 − 2592 C1 ω

+ 656 C1

+1728 C1

− 3712 C1

− 2256 C1

− 312 C1

)/(2 ω)

The total error must be minimized with respect to ω and C1 . To this end, the

derivatives of the Error with respect to ω and C1 are set equal to zero in cond1

and cond2 . The very lengthy algebraic equations are not shown here.

cond1:=diff(Error,omega)=0; cond2:=diff(Error,C1)=0;

The two conditions are numerically solved for C1, ω and the solution assigned.

sol:=fsolve({cond1,cond2},{C1,omega}); assign(sol):

sol := {ω = −1.978539105, C1 =1.937051656}

320 CHAPTER 8. NLODES & PDES OF PHYSICS

The fully evaluated total error, trial function Φ, and time T , are now displayed.

Error:=evalf(abs(Error)); Phi:=Phi; T:=evalf(T);

Error := 0.2146672855

Φ:=1.937051656 cos(1.978539105 t)+0.062948344 cos(5.935617315 t)

T := 12.70267601

The trial function Φ is plotted over the time interval t=0 to T .SinceΦisour

“best” approximation to the exact solution x, the ordinate is labeled as x.

gr1:=plot(Phi,t=0..T,labels=["t","x"]):

The NLODE is solved numerically, subject to the same initial condition,

Y:=dsolve({e(y(t))=0,y(0)=A,D(y)(0)=0},y(t),numeric,

output=listprocedure):

and plotted over the same time range. A point style is chosen, the numerical

curve being represented by 150 size 12 black circles.

gr2:=odeplot(Y,[t,y(t)],0..T,style=point,symbol=circle,

symbolsize=12,numpoints=150,color=black):

The graphs of the Ritz and numerical solutions are superimposed,

display([gr1,gr2]);

the resulting picture being shown in Figure 8.11. The Ritz solution is in excel-

–2

Figure 8.11: Circles: numerical solution. Curve: Ritz solution.

lent agreement with the numerical points over the time range of the plot. It

should be mentioned that there are several variations on the Ritz approximation

method, which are discussed in Zwillinger’s book.

8.4. NONLINEAR PDES 321

8.4 Nonlinear PDEs

Most of the nonlinear PDEs that my graduate students and I have had to solve

in our research in nonlinear optics and ﬂuid dynamics have involved either

seeking special solutions, or more often using numerical methods. Here, a few

examples of the former are presented, the latter being covered in Chapter 9.

8.4.1 John Scott Russell’s Chance Interview

The boisterous sea of liberty is never without a wave.

Thomas Jeﬀerson, American president, (1743–1826)

Shortly after Jeﬀerson’s death, the great Scottish naval architect, John Scott

Russell, made the following important scientiﬁc observation:

I was observing the motion of a boat which was rapidly drawn along a

narrow channel by a pair of horses, when the boat suddenly stopped –

not so the mass of water in the channel which it had put in motion;

it accumulated round the prow of the vessel in a state of violent

agitation, then suddenly leaving it behind, rolled forward with great

velocity, assuming the form of a large solitary elevation, a rounded

smooth and well-deﬁned heap of water, which continued its course

along the channel apparently without change of form or diminution

of speed. I followed it on horseback, and overtook it still rolling on

at a rate of some eight or nine miles an hour, preserving its original

ﬁgure some thirty feet long and a foot to a foot and a half in height.

Its height gradually diminished, and after a chase of one or two

miles I lost it in the windings of the channel. Such, in the month of

August 1834, was my ﬁrst chance interview with that singular and

beautiful phenomenon ...

The relevant nonlinear PDE describing the water waves for the above situa-

tion was derived several decades later by Korteweg and deVries, the Korteweg-

deVries (KdV) equation taking the form [KdV95]

∂ψ

∂t

+ ψ

∂ψ

∂x

∂

∂x

=0, (8.25)

with x the distance, t the time, and ψ the transverse displacement of the water

from equilibrium. Damping of the water waves is completely neglected here.

What Scott Russell observed was a special localized solution of the KdV equa-

tion, referred to as a solitary wave solution. The form of the solitary wave may

be obtained by assuming that ψ(x, t)=U(z = x−ct), subject to the asymptotic

boundary conditions that U(z) and all its derivatives vanish as |z|→∞.The

parameter c will be the speed of the solitary wave. In this recipe, U will be

derived and the solitary wave animated for two diﬀerent values of c.

322 CHAPTER 8. NLODES & PDES OF PHYSICS

The plots and PDEtools packages are loaded, the former containing the

animation command, while the PDEtools package has the dchange command

which will be used to transform the variables in the KdV equation.

restart: with(plots): with(PDEtools):

The KdV equation (8.25) is entered.

pde:=diff(psi(x,t),t)+psi(x,t)*diff(psi(x,t),x)

+diff(psi(x,t),x,x,x)=0;

pde := (

∂

∂t

ψ(x, t)) + ψ(x, t)(

∂

∂x

ψ(x, t)) + (

∂

∂x

ψ(x, t)) = 0

The variable transformation x = z + cτ, t = τ,andψ(x, t)=U(z) is entered,

x, t and ψ(x, t) being the “old” variables, while z, τ,andU(z)arethe“new”

variables. This transformation will have the eﬀect of reducing the nonlinear

pde to a nonlinear ODE.

tr:={x=z+c*tau,t=tau,psi(x,t)=U(z)};

tr := {x = z + cτ, t= τ, ψ(x, t)=U (z)}

The variable transformation tr is applied to pde using the dchange command.

ode1:=dchange(tr,pde,[z,tau,U(z)]);

ode1 := −(

U (z)) c + U (z)(

U (z)) + (

U (z)) = 0

The lhs of ode1 is now integrated. Since U(z) and its second derivative must

vanish at z =∞, the integration constant is equal to zero. The integrated result

is set equal to zero, yielding the second order NLODE ode2 .

ode2:=int(lhs(ode1),z)=0;

ode2 := −c U (z)+

U (z)

U (z)) = 0

Then ode2 is analytically solved for U(z), yielding two implicit solutions.

sol:=dsolve(ode2,U(z));

sol :=



U (z)

√

−3 a

+9c a

+9 C1

a − z − C2 =0,



U (z)

−

√

−3 a

+9c a

+9 C1

a − z − C2 =0

The positive square root solution (the ﬁrst one here) is now diﬀerentiated with

respect to z,anddU/dz is isolated to the lhs of the slope equation.

slope:=isolate(diff(sol[1],z),diff(U(z),z));

slope :=

U (z)=



−3 U (z)

+9c U (z)

+9 C1

Since U(z)anddU(z)/dz vanish at z =∞, the integration constant

C1canbe

set equal to zero. The radical simpliﬁcation command is also applied in ode3 .

ode3:=radsimp(eval(slope,_C1=0));

8.4. NONLINEAR PDES 323

ode3 :=

U (z)=

√

3 U (z)



−U (z)+3c

The ﬁrst order ODE ode3 is analytically solved.

sol2:=dsolve(ode3,U(z));

sol2 := z +

2 arctanh(



−3 U (z)+9c

√

)

√

C1 =0

Then sol2 is solved for U(z). The integration constant

C1 determines the

location of the peak of the solitary wave and can be set equal to zero without

loss of generality.

U:=eval(solve(sol2,U(z)),_C1=0);

U := 3 c − 3 tanh(

√

)

Converting U to a sine/cosine form and simplifying with the trig option gives

us the ﬁnal form of the solitary wave solution in terms of the variable z.

U:=simplify(convert(U,sincos),trig);

U :=

3 c

cosh(

√

)

U has a maximum height at z = 0 which is proportional to c, a width which

decreases with increasing c, and goes to zero as |z|→∞. We can convert back

to the original variables by substituting z = x−ctinto U. ψ, given in the output

of the following command line, is the solitary wave observed so long ago by John

Scott Russell. This is clearly a localized pulse which travels unchanged in shape

in the positive x direction with velocity c. Since the height is proportional to

c, this implies that taller KdV solitary waves travel faster than shorter ones.

We can conﬁrm this by animating the solitary wave proﬁle for two diﬀerent

velocities. First, let’s apply the unapply command to ψ, turning it into a

functional operator f depending on the value of c.

psi:=subs(z=x-c*t,U); f:=unapply(psi,c):

ψ :=

3 c

cosh(

(x − ct)

√

)

Taking c = 1 and c = 3, two solitary waves are animated by entering f(1) and

f(3) in the animate command. On executing the following command line on

the computer, and clicking on the plot, the animation will begin.

animate({f(1),f(3)},x=-10..40,t=0..10,numpoints=200,

frames=50,thickness=2,axes=framed,labels=["x","psi"]);

The initial proﬁle of the animation is shown in Figure 8.12. As time progresses,

the shorter solitary wave lags behind the taller one.

324 CHAPTER 8. NLODES & PDES OF PHYSICS

psi

02040

Figure 8.12: Initial frame of the solitary wave animation.

This diﬀerence in speeds between solitary waves of diﬀerent amplitudes will be

used in the next chapter to numerically test their stability against collisions

with one another. If they survive unchanged, except perhaps for a phase shift,

they are called solitons.

8.4.2 There is a Similarity

Ought not education to bring out and fortify the diﬀerences

rather than the similarities?

Virginia Woolf, British novelist, referring to the gender issue, (1882–1941).

In the previous recipe, the number of independent variables was reduced from

two to one, thus converting a nonlinear PDE to an ODE for which a physically

important exact solution could be obtained. It was a simple application of the

similarity method. Now I will show you a slightly more complicated similarity

solution, solving a nonlinear diﬀusion equation by taking a diﬀerent algebraic

combination of the independent variables.

Consider the one-dimensional nonlinear diﬀusion equation,

∂C

∂t

∂

∂x



D(C)

∂C

∂x



,D(C)=C

, (8.26)

with the diﬀusion coeﬃcient D no longer constant, but a function of the con-

centration C. Several special cases of (8.26) have been studied in the litera-

ture, viz., n = 3 to model the spreading of thin liquid ﬁlms under the action

of gravity [Buc77], n ≥ 1 to describe the percolation of gas through porous

media [Mus37], and n=6 to model certain radiative heat transfer [LP80].

After loading the PDEtools and plots library packages,

restart: with(PDEtools): with(plots):

8.4. NONLINEAR PDES 325

the nonlinear diﬀusion equation (8.26) is entered. The inert form of the deriva-

tive operator is used to prevent the rhs from being explicitly diﬀerentiated.

pde:=Diff(C(x,t),t)-Diff(C(x,t)ˆn*diff(C(x,t),x),x)=0;

pde := (

∂

∂t

C (x, t)) − (

∂

∂x

(C (x, t)

(

∂

∂x

C (x, t)))) = 0

New variables z, τ ,andU(z) are introduced, related to the old variables x, t,

and C(x, t) by the transformation x = zτ

, t = τ ,andC(x, t)=U(z)/τ

,with

the parameter m to be determined.

tr:={x=z*tauˆm,t=tau,C(x,t)=U(z)/tauˆm};

tr := {x = zτ

,t= τ, C (x, t)=

U (z)

}

The transformation is applied to pde using the dchange command, the result

then being multiplied by −τ

m+1

and expanded.

ode:=expand(-tauˆ(m+1)*dchange(tr,pde,[z,tau,U(z)]));

ode := (

U (z)) zm+ U (z) m +

τ (

U (z)

)

n (

U (z))

(τ

)

U (z)

τ (

U (z)

)

(

U (z))

(τ

)

To remove all τ factors from ode and thus really produce an ODE, we must

have (τ

)

=τ (1/τ

)

which is now entered.

eq:=(tauˆm)ˆ2=tau*(1/tauˆm)ˆn;

eq := (τ

)

= τ (

)

Then eq is solved for m, yielding m =1/(n + 2). This m value will be automat-

ically substituted into ode .

m:=solve(eq,m);

m :=

2+n

Dividing ode by m, simplifying symbolically with the power option, and then

applying the general simpliﬁcation command, yields the formidable looking

NLODE given in ode2

ode2:=simplify(simplify(ode/m,power,symbolic));

ode2 := (

U (z)) z + U (z)+2U (z)

(n−1)

n (

U (z))

+ U (z)

(

U (z)) n

+ U (z)

(n−1)

(

U (z))

+2U (z)

(

U (z)) = 0

A general analytic solution to ode2 is sought using the dsolve command, the

result then being simpliﬁed symbolically with the power option.

sol:=simplify(dsolve(ode2,U(z)),power,symbolic);

326 CHAPTER 8. NLODES & PDES OF PHYSICS

sol := U (z)=

b( a)&where

[{

b( a)

(

b( a)) +

a b( a)+2 C1 + C1 n

2+n

=0},

{

a = z, b( a)=U (z)}, {z = a, U (z)= b( a)}]

A general solution for U(z) is not produced in sol, but U(z)=

b( a)where b( a)

satisﬁes a ﬁrst order NLODE. A solution can be generated if the integration

constant

C1 is set equal to zero. The following operand command is used to

extract the resulting ODE from sol.

_C1:=0: ode3:=op([2,2,1],sol)[1];

ode3 :=

(

d a

b( a)) +

a b( a)

2+n

Then, ode3 is analytically solved for

b( a).

sol2:=dsolve(ode3,_b(_a));

sol2 :=

b( a)=

(

−

n +4 C2 +2 C2 n

2+n

)

(

)

(

)

The integration constant

C2 controls the height and width of the solution.

Here, I will take

C2 =1. Substituting

a =z on the rhs of sol2 yields U(z).

U:=subs({_a=z,_C2=1},rhs(sol2));

U :=

(

−z

n +4+2n

2+n

)

(

)

(

)

The above solution only exists for z values such that (−z

n+4+2 n)/(2+n) ≥ 0.

Outside this range, we can take U = 0, which clearly satisﬁes the NLODE.

Substituting z = x/t

into U/t

gives us a solution, labeled c, in terms of the

original variables.

c:=subs(z=x/tˆm,U/tˆm);

c :=

⎛

⎜

⎝

−

(

2+n

)

+4+2n

2+n

⎞

⎟

⎠

(

)

(

)

(

2+n

)

The x range for which c is valid is found by setting c =0 and solving for x.The

result is simpliﬁed with the radical simpliﬁcation command and labeled X.

X:=radsimp(solve(c=0,x));

X :=

√



n (2 + n) t

(

2+n

)

As a speciﬁc example, let’s take n = 3, so the solution can be used to describe

the spreading of a thin liquid ﬁlm under the action of gravity. Other n values

8.4. NONLINEAR PDES 327

can be chosen, if so desired. The total time range is taken to be T =1500, and

X is evaluated at this time and assigned the name X0 .

n:=3: T:=1500: X0:=eval(X,t=T):

The complete concentration proﬁle C is then described by the following piece-

wise function.

C:=piecewise(abs(x)<X,c,abs(x)>X,0);

C :=

⎧

⎪

⎨

⎪

⎩

(2/3)

(−

3 x

5 t

(2/5)

+2)

(1/3)

2 t

(1/5)

|x| <

√

15 t

(1/5)

√

15 t

(1/5)

< |x|

A 3-dimensional plot of C is now produced,

plot3d(C,x=-X-5..X+5,t=1..T,numpoints=1500,axes=box,labels

=["x","t","C"],orientation=[40,70],tickmarks=[3,3,3]);

the resulting picture being shown in Figure 8.13.

–10

500

1000

0.5

Figure 8.13: Similarity solution of the nonlinear diﬀusion equation for n=3.

The solution captures some of the experimentally observed features of spreading

thin ﬁlms. Unlike the situation for the linear (D constant) diﬀusion equation,

there is now a sharp interface separating the regions of nonzero and zero con-

centration. Further, the interface propagates with ﬁnite speed, in contrast with

the inﬁnite speed of linear diﬀusion. The propagation of this interface can be

observed by executing the following animate command, clicking on the resulting

computer plot, and on the start arrow.

animate(C,x=-X0..X0,t=1..T,frames=50,numpoints=500,

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

For more on similarity methods, you are referred to Bluman and Cole [BC74].

328 CHAPTER 8. NLODES & PDES OF PHYSICS

8.4.3 Creating Something Out Of Nothing

Say nothing good of yourself, you will be distrusted;

say nothing bad of yourself, you will be taken at your word.

Joseph Roux, French priest, writer, (1834–86)

The sine-Gordon equation (SGE),

∂

∂x

−

∂

∂t

=sinu, (8.27)

is a model equation for describing the motion of a Bloch domain wall between

two ferromagnetic domains. A solitary wave solution to the SGE, describing the

Bloch wall motion, could be obtained in a similar manner to that for the KdV

equation. An alternate way is to make use of an auto-B¨acklund transformation.

Given a solution of a nonlinear PDE, such a transformation allows us to ﬁnd a

diﬀerent solution of the same PDE. More, generally a B¨acklund transformation

may enable one to use the solution of one nonlinear PDE to determine the

solution of another nonlinear PDE. B¨acklund transformations are diﬃcult to

ﬁnd, so in this recipe I will merely conﬁrm the auto-B¨acklund transformation

for the SGE and use it to create a non-trivial solitary wave solution of (8.27),

starting with the trivial null (u=0) solution.

After loading the PDEtools and plots library packages,

restart: with(PDEtools): with(plots):

the sine-Gordon equation is entered.

sge:=diff(u(x,t),x,x)-diff(u(x,t),t,t)=sin(u(x,t));

sge := (

∂

∂x

u(x, t)) − (

∂

∂t

u(x, t)) = sin(u(x, t))

Introducing the variable transformation x=X +T , t= X −T , u(x, t)=U(X, T ),

tr:={x=X+T,t=X-T,u(x,t)=U(X,T)};

tr := {u(x, t)=U (X, T ),t= X −T, x = X + T }

and using the dchange command, the SGE takes the form shown in sge2 .

sge2:=dchange(tr,sge,[X,T,U(X,T)]);

sge2 :=

∂

∂X ∂T

U (X, T)=sin(U (X, T ))

From Zwillinger [Zwi89], the auto-B¨acklund transformation for the SGE is given

by ab1 and ab2 ,wherea is an arbitrary parameter.

ab1:=diff(V(X,T),X)=diff(U(X,T),X)+2*a*sin((V(X,T)+U(X,T))/2);

ab1 :=

∂

∂X

V (X,T)=(

∂

∂X

U (X,T)) + 2 a sin(

V (X,T)+

U (X,T))

ab2:=diff(V(X,T),T)=-diff(U(X,T),T)+(2/a)*sin((V(X,T)-U(X,T))/2);

ab2 :=

∂

∂T

V (X,T)=−(

∂

∂T

U (X,T)) +

2sin(

V (X,T) −

U (X,T))