Enns R.H. Computer Algebra Recipes for Mathematical Physics
Подождите немного. Документ загружается.
9.2. PARTIAL DIFFERENTIAL EQUATIONS 369
The former are given by U
i,0
=f(x
i
), where f(x) is the input profile. For the U
values on the first step, we shall use the result (9.25), i.e., U
i,1
=f(x
i
)+kg(x
i
),
where g(x) is the initial transverse velocity.
We begin the main part of our recipe by considering M = 150 spatial steps
and N =300 time steps. The starting time of the recipe is also recorded.
>
restart: with(plots): M:=150: N:=300: begin:=time():
A functional operator F is formed to calculate a time-dependent solitary wave
profile centered at time t=0 at x=X and having speed c.
>
F:=(X,c)->3*c *(sech((sqrt(c)/2)*((x-X)-c*t)))ˆ2:
The spatial stepsize is taken to be h = 1, the time stepsize k =0.25, so k/h
3
=
0.25 and the scheme will be stable. For f(x), we will consider two separated
solitary waves, a taller one having speed c1 =0.8 and a shorter one with speed
c2 =0.1. The taller, faster, solitary wave will be initially centered at X1 =
M/3= 50, the shorter, slower, wave at X2 =M/2= 75. As time progresses, the
taller solitary wave will overtake the shorter one and a collision take place.
>
h:=1: k:=.25: c1:=.8; c2:=.1; X1:=M/3; X2:=M/2;
c1 := 0.8 c2 := 0.1 X1 := 50 X2 := 75
Then F(X1,c1)+F(X2,c2) is entered and the time derivative taken in g.
>
f:= F(X1,c1)+F(X2,c2): g:=diff(f,t):
On setting t = 0, the initial profile is given by f =3c1 sech(
√
c1
2
(x−X1 ))
2
+
3 c2 sech(
√
c2
2
(x−X2 ))
2
which is now plotted in Figure 9.15.
>
t:=0: plot(f,x=0..M);
0
0.5
1
1.5
2
20 40
60
x
100 120 140
Figure 9.15: Initial profile.
To evaluate f at the spatial mesh points on the zeroth time row and g at the
mesh points on the first time row, they are turned into functional operators in
terms of x by using the unapply command.
>
f2:=unapply(f,x): g2:=unapply(g,x):
Using f2 and g2, U
i,0
and U
i,1
are calculated in ic1 and ic2 for i=0 to M.
370 CHAPTER 9. NUMERICAL METHODS
>
ic1:=seq(U(i,0)=evalf(f2(i)),i=0..M):
>
ic2:=seq(U(i,1)=evalf(f2(i))+k*g2(i),i=0..M):
To avoid unknown U values creeping into the do loop, we will set all U values
equal to zero for i =0 to M and j =2 to N. These “initialized” zeros will be
overwritten as the do loop progresses.
>
init:=seq(seq(U(i,j)=0,i=0..M),j=2..N):
Then ic1 , ic2 ,andinit are assigned.
>
assign(ic1,ic2,init):
The algorithm (9.31) is then iterated from i=2 to M − 2andj =1 to N.We
must start i at 2 and end at M − 2 to avoid unknown U values entering the
iterative calculation.
>
for j from 1 to N do;
>
for i from 2 to M-2 do
>
U(i,j+1):=U(i,j-1)-(k/(3*h))*(U(i+1,j)+U(i,j)+U(i-1,j))
*(U(i+1,j)-U(i-1,j))-(k/hˆ 3)*(U(i+2,j)-2*U(i+1,j)
+2*U(i-1,j)-U(i-2,j));
>
end do: end do:
An operator gr is created to plot U
i,j
for i=2 to M −2onthejth time step.
>
gr:=j->plot([seq([i,U(i,j)],i=2..M-2)]):
Every second graph is displayed and the sequence of pictures animated with the
insequence=true option. Click on the computer plot and on the start arrow
to initiate the animation.
>
display([seq(gr(2*j),j=0..N/2)],insequence=true);
It is observed that the tall solitary wave passes through the shorter one and both
solitary waves appear to be unchanged after the collision, aside from a very small
numerical ripple which arises because of the relatively coarse grid which was
used and the fact that the initial solitary wave profiles overlapped slightly and
their tails were artificially truncated at the edges of the numerical grid. As the
tall pulse passes through the shorter one, you should notice that the amplitudes
do not add linearly, which would be expected from linear superposition. Finally,
the cpu time for the recipe is now recorded.
>
cpu:=time()-begin;
cpu := 29.913
9.3 Supplementary Recipes
09-S01: White Dwarf Equation
Chandresekhar’s theory of white dwarf stars produces the NLODE [Cha39]
xy
(x)+2y
(x)+x (y
2
− C)
3/2
=0, with y(0) = 1,y
(0) = 0.
Making use of the first principles and dial-up Euler methods with h = .01 and
18 - digits accuracy, numerically compute y(x) over the range 0 ≤ x ≤ 4with
9.3 SUPPLEMENTARY RECIPES 371
C =0.1 and plot both numerical results together. Show that the numerical
value of y(x = 1) is exactly the same for both methods. Hint: Start at x =0.01
to avoid any problem at the origin.
09-S02: Spruce Budworm Infestation
The sudden outbreak of the spruce budworm which can rapidly defoliate a for-
est and kill the trees can be described [LJH78] by the dimensionless ODE
˙x(τ )=rx(1 −
x
K
) −
x
2
1+x
2
.
Here x(τ) is proportional to the budworm population number at time τ and
the growth coefficient r and carrying capacity parameter K are positive con-
stants. The first term describes the growth of the budworm population with a
saturation effect included due to the finite forest available, while the last term
models the decrease in population due to bird predation. Taking K =300, use
the first-principles Euler method with h =0.01 and x(0) = 0.5todetermine
the time evolution of the budworm population over the time interval t= 0 to 50
for r =0.1, r =0.5, and r =1.0. Plot and discuss the results.
09-S03: A Math Example
Consider the first-order nonlinear ODE
dy/dx =
√
x (y − x)(y −2)
with y(0) = 1. Show that an analytic closed form of solution is not obtain-
able with the dsolve command. Taking h =0.01 and 10-digit accuracy, solve
for y(x)outtox = 5 using the modified Euler method and plot every 10th
numerical point.
09-S04: Hermione Hippo
Hermione Hippo swims across the muddy Mombopo river by steadily aiming
at a tree directly across the river on the water’s edge. The Mombopo is 1 km
wide and has a speed of 1.5 km/hour while Hermione’s speed is 2 km/hour.
In Cartesian coordinates, Hermione is initially at (x =1,y = 0) while the tree
is at (0, 0). Derive Hermione’s equations of motion in Cartesian coordinates.
Using the first-principles 4th-order RK method with h =0.001, numerically
solve these equations to determine how long it takes Hortense to reach the tree.
Determine the analytic solution y(x) for Hortense’s path across the river and
plot the analytic and numerical solutions together in the same figure.
09-S05: The Oregonator
A kinetic model ([FKN72],[FN74],[EM00]) describing the oscillations in a cer-
tain chemical system is the Oregonator system of nonlinear ODEs, viz.,
˙x = x + y − qx
2
− xy, ˙y = −y +2hz− xy, p ˙z = x − z,
where x, y,andz are the dimensionless concentrations of three chemical species
and , p, q,andh are positive parameters. Taking =0.003, p =2, q =0.006,
h =0.75, x(0) = 30, y(0) = 10, and z(0) = 30, numerically solve the system of
equations using the RKF45 method with the option maxfun=0. Then use the
372 CHAPTER 9. NUMERICAL METHODS
odeplot command with numpoints=3000 to plot (a) x vs. y vs. z for t=0..50,
(b) x vs. t, y vs. t,andz vs. t in the same figure, for t=0..30. Numerically show
that a closed loop results in the x-y-z space no matter what initial conditions
are chosen. The Oregonator is an example of a 3-dimensional limit cycle.
09-S06: Lorenz’s Butterfly
The Lorenz system [Lor63] of nonlinear ODES are given by
˙x = σ (y − x), ˙y = rx− y −xz, ˙z = xy− bz,
with x, y, z real and σ, r,andb real, positive constants. Derive a second-order-
accurate semi-implicit numerical scheme for the Lorenz system. Taking σ =10,
b =8/3, r = 28, x(0) = 2, y(0) = 5, z(0) = 5, a stepsize h =0.01, and n = 4000
iterations, iterate the algorithm. Plot x vs. y vs. z,usingazhue shading.
The resulting figure should resemble the wings of a butterfly. After a transient
interval, the trajectory remains confined to the region of the wings, but never
retraces the same path. This is another example of a strange chaotic attractor.
09-S07: A Stiff Harmonic Oscillator
Consider the following heavily overdamped forced harmonic oscillator equation,
¨x +50˙x + x =2sint,
with the initial condition x(0) = 0, ˙x(0) = 10. Solve the ODE analytically and
determine the two widely different characteristic times in the solution. Create
operators to numerically solve the ODE using the 4th-order Rung´e–Kutta dial-
up method and plot the numerical solution for arbitrary stepsize h. By plotting
the analytic and numerical solution together, determine the approximate h value
at which numerical instability sets in. This critical h value should be roughly
comparable to the shortest characteristic time.
09-S08: Courant Stability Condition
Prove the Courant stability condition for the linear wave equation.
09-S09: Poisson’s Equation
Consider Poisson’s equation
∂
2
V
∂x
2
+
∂
2
V
∂y
2
= xe
y
, 0 ≤ x ≤ 2, 0 ≤ y ≤ 1,
with V (0,y)=0, V (2,y)=2e
y
, V (x, 0) =x,andV (x, 1) = ex, as the boundary
conditions. Using a central difference approximation for each second derivative
at a mesh point (i, j) and evaluating the rhs of Poisson’s equation at this point,
derive a finite difference scheme for solving the equation. Dividing the x-interval
into 30 steps and the y interval into 15 steps, determine V at each mesh point
and plot the numerical solution.
09-S10: Crank–Nicolson Method
The Crank–Nicolson (CN) method results on averaging the forward and back-
ward difference methods. For the normalized linear temperature equation, the
9.3 SUPPLEMENTARY RECIPES 373
forward and backward schemes connecting the j and j + 1 time steps are,
(T
i,j+1
− T
i,j
)
k
=
(T
i+1,j
− 2 T
i,j
+ T
i−1,j
)
h
2
,
(T
i,j+1
− T
i,j
)
k
=
(T
i+1,j+1
− 2 T
i,j+1
+ T
i−1,j+1
)
h
2
.
Adding the two equations, dividing by 2, and rearranging, yields,
−rT
i−1,j+1
+ 2 (1 + r) T
i,j+1
− rT
i+1,j+1
= rT
i−1,j
+2(1−r) T
i,j
+ rT
i+1,j
with r = k/h
2
. This finite difference CN formula relates three unknown values
of T on the j + 1 time row to three known values on the j time row. Modify
recipe 09-2-2 to implement the CN formula, taking all parameters the same,
except now with r =0.6. As you may verify the explicit scheme in 09-2-2 is
numerically unstable for this r value,buttheCNschemeisstable.
09-S11: Klein–Gordon Equation
The transverse vibrations of a horizontal string embedded in a vertical elastic
membrane which exerts a Hooke’s law restoring force on the string can be de-
scribed by the Klein–Gordon equation,
∂
2
U
∂x
2
−
∂
2
U
∂t
2
= aU,
with a a positive constant. If such a string of length L = 1 is fixed at both
ends and has the initial shape U(x, 0) = x
4
(L − x)
2
and transverse velocity
˙
U(x, 0) = x (L − x), numerically determine the transverse displacement of the
string for t ≥ 0 using a matrix approach and animate the motion. Take a=100,
spatial stepsize h=0.025 and time stepsize k =0.01.
Bibliography
[AS72] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Func-
tions with Formulas, Graphs, and Mathematical Tables. National
Bureau of Standards, Washington, DC, 1972.
[AW00] G. B. Arfken and H. J. Weber. Mathematical Methods for Physi-
cists, 5th ed. Academic Press, New York, 2000.
[Bay99] W. E. Baylis. Electrodynamics: A Modern Geometric Approach.
Birkh¨auser, Boston, MA, 1999.
[BC74] G. W. Bluman and J. D. Cole. Similarity Methods for Differential
Equations. Springer-Verlag, New York, 1974.
[BF89] R. L. Burden and J. D. Faires. Numerical Analysis, 4th ed. PWS-
Kent, Boston, 1989.
[Boa83] M. L. Boas. Mathematical Methods in the Physical Sciences, 2nd
ed. John Wiley, New York, 1983.
[Bra68] T. C. Bradbury. Theoretical Mechanics. Wiley, New York, 1968.
[Buc77] J. Buckmaster. Viscous sheets advancing over dry beds. Journal
of Fluid Mechanics, 81:735, 1977.
[Cha39] S. Chandresekhar. An Introduction to the Study of Stellar Struc-
ture. Dover Reprint, Chicago, 1939.
[Cha03] B. W. Char. Maple 9 Learning Guide. Waterloo Maple, Waterloo,
Canada, 2003.
[Dav62] H. T. Davis. Introduction to Nonlinear Differential and Integral
Equations. Dover, New York, 1962.
[DL03] C. Doran and A. Lasenby. Geometric Algebra for Physicists. Cam-
bridge University Press, Cambridge, UK, 2003.
[EM00] R. H. Enns and G. C. McGuire. Nonlinear Physics with Maple for
Scientists and Engineers, 2nd ed. Birkh¨auser, Boston, MA, 2000.
376 BIBLIOGRAPHY
[FC99] G. R. Fowles and G. L Cassiday. Analytical Mechanics, 6th ed.
Saunders College, Orlando, FL, 1999.
[Feh70] E. Fehlberg. Klassische runge-kutta formeln vierter und
niedrigerer ordnung mit schrittweiten-kontrolle und ihre andwen-
dung auf w¨armeleitungsprobleme. Computing, 6:61, 1970.
[FKN72] R. J. Field, E. K¨or¨os, and R. M. Noyes. Oscillations in chemical
systems, Part 2. Thorough analysis of temporal oscillations in the
bromate–cerium–malonic acid system. Journal of the American
Chemical Society, 94:8649, 1972.
[FN74] R. J. Field and R. M. Noyes. Oscillations in chemical systems,
IV. Limit cycle behavior in a model of a real chemical reaction.
Journal of Chemical Physics, 60:1877, 1974.
[GPS02] H. Goldstein, C. Poole, and J. Safco. Classical Mechanics, 3rd ed.
Addison Wesley, New York, 2002.
[Gri95] D. J. Griffiths. Introduction to Quantum Mechanics. Prentice-Hall,
Englewood Cliffs, N.J., 1995.
[Gri99] D. J. Griffiths. Introduction to Electrodynamics, 3rd ed. Prentice-
Hall, Upper Saddle River, N.J., 1999.
[Hel50] H. Helmholtz. Messungen ¨uber den zeitlichen verlauf der zuchung
animalischer muskeln und die fortplanzungsgeschwindigkeit der
reizung in der nerven. Arch. Anat. Physiol., 276, 1850.
[Hes99] D. Hestenes. New Foundations for Classical Mechanics, 2nd ed.
Kluwer Academic, Dordrecht,The Netherlands, 1999.
[Hil57] F. B. Hildebrand. Advanced Calculus for Engineers. Prentice-Hall,
Englewood Cliffs, N. J., 1957.
[IM69] J. Irving and N. Mullineux. Mathematics in Physics and Engi-
neering. Academic Press, New York, 1969.
[KB43] N. Krylov and N. Bogoliubov. Introduction to Nonlinear Mechan-
ics. Princeton University Press, Princeton, 1943.
[KdV95] D. J. Korteweg and G. de Vries. On the change of form of long
waves advancing in a rectangular canal, and a new type of long
stationary wave. Philosophical Magazine, 39:422, 1895.
[LJH78] D. Ludwig, D. D. Jones, and C. S. Holling. Qualitative analysis
of insect outbreak systems: the spruce budworm and forest. J.
Anim. Ecol., 47:315, 1978.
BIBLIOGRAPHY 377
[Lor63] E. N. Lorenz. Deterministic nonperiodic flow. J. Atmospheric Sci.,
20:130, 1963.
[LP80] E. W. Larsen and G. C. Pomraning. Asymptotic analysis of non-
linear marshak waves. SIAM Journal of Applied Mathematics,
39:201, 1980.
[LS71] L. Lapidus and J. H. Seinfeld. Numerical Solution of Ordinary
Differential Equations. Academic Press, New York, 1971.
[MF53] P. M. Morse and H. Feshbach. Methods of Theoretical Physics.
McGraw Hill, New York, 1953.
[MGH
+
03a] M. B. Monagan, K. O. Geddes, K. M. Heal, G. Labahn, S. M.
Vorkoetter, J. McCarron, and P. DeMarco. Maple 9 Advanced
Programming Guide. Waterloo Maple, Waterloo, Canada, 2003.
[MGH
+
03b] M. B. Monagan, K. O. Geddes, K. M. Heal, G. Labahn, S. M.
Vorkoetter, J. McCarrona, and P. DeMarco. Maple 9 Introductory
Programming Guide. Waterloo Maple, Waterloo, Canada, 2003.
[MM57] H. Margenau and G. M. Murphy. The Mathematics of Physics
and Chemistry, 2nd ed. D. Van Nostrand, New York, 1957.
[Mor48] P. M. Morse. Vibration and Sound. McGraw-Hill, New York, 1948.
[MT95] J. B. Marion and S. T. Thornton. Classical Dynamics of Particles
and Systems, 4th ed. Saunders College, Orlando, FL, 1995.
[Mus37] M. Muskat. The Flow of Homogeneous Fluids Through Porous
Media. McGraw Hill, New York, 1937.
[MW71] J. Mathews and R. L. Walker. Mathematical Methods of Physics,
2nd ed. Addison-Wesley, New York, 1971.
[NYA65] J. Nagumo, S. Yoshizawa, and S. Arimoto. Bistable transmission
lines. Trans. IEEE on Circuit Theory, CT-12:400, 1965.
[PFTV90] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling.
Numerical Recipes. Cambridge University Press, Cambridge, 1990.
[PLA92] M. Peastrel, R. Lynch, and A. Armenti. Terminal velocity of a
shuttlecock in vertical fall. In Angelo Armenti (Jr.), editor, The
Physics of Sports. American Institute of Physics, New York, 1992.
[RE76] S. S. Rangnekar and R. H. Enns. Numerical solution of the tran-
sient gain equations for stimulated backward scattering in absorb-
ing fluids. Canadian Journal of Physics, 54:1564, 1976.
[Ros76] O. E. Rossler. An equation for continous chaos. Physics Letters,
57A:397, 1976.
378 BIBLIOGRAPHY
[Sco87] D. E. Scott. An Introduction to Circuit Analysis. McGraw-Hill,
New York, 1987.
[Spi71] M. Spiegel. Advanced Mathematics for Engineers and Scientists.
McGraw Hill, New York, 1971.
[SR66] I. S. Sokolnikoff and R. M. Redheffer. Mathematics of Physics and
Modern Engineering. McGraw-Hill, New York, 1966.
[Ste87] J. Stewart. Calculus. Brooks/Cole, Pacific Grove, CA, 1987.
[Ste94] B. Stearn. Simply Heartsmart Cooking. Random House, Toronto,
Ont., Canada, 1994.
[Wal84] P. R. Wallace. Mathematical Analysis of Physical Problems.Dover,
New York, 1984.
[Wie73] S. Wieder. The Foundations of Quantum Theory. Academic Press,
New York, 1973.
[Zwi89] D. Zwillinger. Handbook of Differential Equations.Academic
Press, San Diego, 1989.
Index
Abramowitz & Stegun, 82
acceleration, 96
air resistance
Newton’s law, 18
Stokes’s law, 18, 28
Airy wave function, 43
algebraic equations, nonlinear, 284
alpha decay, 55
analytic, 65
analytic function, 185
angular integral, 201
anode, 334
anti-kink soliton, 329
approximation
backward difference, 338, 362
central difference, 62, 355, 362
finite difference, 62, 85
forward difference, 338
fourth derivative, 63
third derivative, 368
arbitrary function, 129
Argand plane, 185
argument, 185
aurora, 278
auto-B¨acklund transformation, 328
autonomous ODE, 302
B¨acklund transformation, 328
back emf, 141
badminton, 18
balalaika, 179
band matrix, 358
bandwidth theorem, 255
bar, 153
beam, 33, 232
bending moment, 21
bending of starlight, 286
bending,criticallength,43
Bernoulli number, 84
Bernoulli ODE, 292, 333
Bessel function, 43, 139, 167, 250
first kind, 37
order of, 37
properties, 38
recurrence relations, 39
second kind, 37
spherical, 170
zeros, 38, 281
Bessel ODE, 37
beta decay, 55
bispherical, 158
Bloch wall, 328
Bouguer, Pierre, 294
boundary condition
asymptotic, 146, 163
clamped end, 21, 153
fixed edge, 149, 167
fixed-end, 130, 133, 139
free-end, 21, 139
split-boundary, 183
Sturm–Liouville, 41
boundary value problem, 16, 146
brachistochrone, 258
branch cut, 203, 209, 244
branch point, 185, 256
bright soliton, 336
Bromwich contour, 244
Bromwich integral, 244, 247, 256
cable, 270
calculus of variations, 257
subsidiary condition, 265