{VERSION 6 0 "IBM INTEL NT" "6.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
267 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 271 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
272 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 0 1 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 276 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
277 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 283 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
287 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 288 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 291 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
292 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 293 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Tim
es" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }
{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 
0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Normal" -1 256 1 
{CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 
0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT 256 31 "Differential Equations \+
in Maple" }}{PARA 256 "" 0 "" {TEXT -1 10 "April 2009" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT 
-1 6 "Basics" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "The principal command to solve dif
ferential equations is  " }{TEXT 257 9 "dsolve.  " }{TEXT -1 95 "It re
quires that we communicate the differential equation, and possibly the
 inital conditions.\n" }{TEXT 273 8 "Example:" }{TEXT -1 19 "  Solve  \+
y' = x - y" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "dsolve( D(y)(x) = x-y
(x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 316 "Notice that we had to e
nter  D(y)(x) for the derivative, and we could not just write y'.  We \+
also had to spell out y(x) explicitly, rather than just entering x-y f
or the right hand side of the equation.  Maple returned the general so
lution, although the appearance of the arbitrary constant looks somewh
at awkward.  " }}{PARA 0 "" 0 "" {TEXT -1 79 "If we add the  initial c
ondition  y(0) = 1 to our example, we need to enter    " }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "dsolve(\{D(y)(x)=x - y(x) , y(0)=1
\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 263 "and maple returns the so
lution to the initial value problem  y'=x-y, y(0)=1.\nIt is advisable \+
to define differential equation and initial conditions separately, and
 then enter these definitions into the dsolve command.  Reworking our \+
example in this way results in" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 25 "de := D(y)(x) = x - y(x);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 13 "ic := y(0)=1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 63 
"The general solution of the differential equation is found from" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "dsolve( de );" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 41 "and the inital value problem has solution" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "dsolve(\{de,ic\});" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 8 "Here is " }{TEXT 274 15 "another example" 
}{TEXT -1 41 ":  Solve the initial value problem  y' = " }{XPPEDIT 18 
0 "1-y^2;" "6#,&\"\"\"F$*$%\"yG\"\"#!\"\"" }{TEXT -1 12 "  , y(0)=0  \+
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "de2 := D(y)(x) = 1 - y(x)^2;\ni
c2 := y(0) = 0;\ndsolve(\{de2, ic2\});" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{SECT 1 {PARA 3 "" 0 "" {TEXT -1 17 "Graphical Methods" }}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 160 "Now we turn to graphical methods.  We lo
ok for ways to display the direction field, with or without solution c
urves.  Graphical methods require that we use the " }{TEXT 258 7 "DEto
ols" }{TEXT -1 18 " package.  \nAs an " }{TEXT 275 7 "example" }{TEXT 
-1 56 ", let us again use the first order equation  y' = x - y." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "restart; with(DEtools):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "de := D(y)(x) = x - y(x);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT 259 19 "Direction field:  \n" }{TEXT -1 
29 "The command for this task is " }{TEXT 260 10 "dfieldplot" }{TEXT 
-1 145 " and we need to communicate the differential equation, the nam
e of the solution curve and the viewing area, i.e. the range of the x \+
and y values." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "dfieldplot( de , y
(x) , x=-2..2,y=-2..2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "Let's \+
add some extras (title, more colorful)" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 78 "dfieldplot( de , y(x) , x=-2..2,y=-2..2, color=(x-y),title=`Di
rection Field`);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 261 41 "Direction Fi
eld with Selected Solutions:\n" }{TEXT -1 57 "For this task we need to
 use a different command, namely " }{TEXT 267 6 "DEplot" }{TEXT -1 67 
", which uses a numerical technique to approximate the solutions.   " 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "DEplot( de , y(x) , x=-2..2,y=-0.
5..3, [[y(0)=0],[y(0)=1]],linecolor=[blue,black]);" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 8 "Here is " }{TEXT 276 15 "another example" }{TEXT 
-1 9 ":  y' =  " }{XPPEDIT 18 0 "x^2+y^2-1;" "6#,(*$%\"xG\"\"#\"\"\"*$
%\"yGF&F'F'!\"\"" }{TEXT -1 4 "  . " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
32 "de := D(y)(x) = x^2 +y(x)^2 - 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 91 "The general solution is a mess involving special functions (Whi
ttaker) and complex numbers." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "dso
lve(de,y(x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "Using graphs we \+
get some idea of what goes on." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "D
Eplot(de,y(x),x=-2..2,y=-2..2,[[y(0)=0],[y(0)=1],[y(0)=-1]],linecolor=
blue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 
{PARA 3 "" 0 "" {TEXT -1 40 "Elastic Spring (A Second Order Equation)
" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 251 "In this part we look at second order dif
ferential equations, and we use the equation of motion for a mass-spri
ng system as an example.  It does not make sense to consider direction
 fields, because second derivatives are involved.  With m as the mass,
 " }{TEXT 279 1 "a" }{TEXT -1 111 " as the friction coefficient and k \+
reflecting the strength of the spring, we consider the second order eq
uation" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "ODE := m*(D@@2)(x)(t) + a
*D(x)(t) + k*x(t) = 0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 137 "where \+
x(t) describes the position of the mass at time t, with x=0 being the \+
equilibrium position. Its general solution can be found with " }{TEXT 
262 7 "dsolve:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "dsolve(ODE);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT 277 8 "Example:" }{TEXT -1 36 " We consid
er the special case where " }{TEXT 280 8 "m=1, a=0" }{TEXT -1 27 " (fr
ictionless motion) and " }{TEXT 281 3 "k=4" }{TEXT -1 69 ".  We substi
tute these values into the differential equation via the " }{TEXT 263 
4 "subs" }{TEXT -1 68 " command and give a new name to the resulting d
ifferential equation:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "ode := sub
s(\{m=1, a=0, k=4\},ODE);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "The \+
general solution of the new equation is" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 12 "dsolve(ode);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 137 "The two \+
constants can be determined if we add initial conditions to our proble
m.  Suppose we require  x(0)=2 and x'(0) = 1.  Then we have" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 25 "IC := x(0)= 2, D(x)(0)=1;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 53 "and the solution to the initial value pro
blem becomes" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "dsolve(\{ode,IC\});
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 268 62 "Example (Dependence of Solut
ions on the Friction Parameter):  " }{TEXT -1 220 "Suppose that the ma
ss is held fixed at m = 4 and the spring constant is fixed at k=1.  Ho
w does varying the friction parameter affect the solution curves?  In \+
all cases we will use the initial conditions x(0)=2, x'(0)=1." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "m := 4; k := 1;\nic := x(0)=2, D(x)
(0)=1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Frictionless motion:" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "x1 := rhs(dsolve(\{subs(a=0,ODE),i
c\}));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 6 "Cases " }{TEXT 282 11 "a
=0.5, a=4 " }{TEXT -1 4 "and " }{TEXT 283 3 "a=8" }{TEXT -1 1 ";" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 118 "x2 := rhs(dsolve(\{subs(a=1/2,ODE)
,IC\}));\nx3 := rhs(dsolve(\{subs(a=4,ODE),IC\}));\nx4 := rhs(dsolve(
\{subs(a=8,ODE),IC\}));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "Graphs
 of the solution curves" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plot([x1
,x2,x3,x4], t=0..40, color=[blue, black, red, maroon]);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 87 "We see that in the frictionless case (blu
e curve) the oscillations go on forever, with " }{TEXT 284 5 "a=0.5" }
{TEXT -1 98 " (black curve) there is a noticable decay of the amplitud
e.  An important transition happens when " }{TEXT 285 3 "a=4" }{TEXT 
-1 141 " (purple curve) called critical damping:  The friction is so s
trong that the solution approaches the x-axis without ever dropping be
low it.  " }{TEXT 286 3 "a=4" }{TEXT -1 43 " is the first value for wh
ich this occurs, " }{TEXT 287 6 "a=3.99" }{TEXT -1 41 " is worked out \+
below for comparison. For " }{TEXT 288 3 "a=8" }{TEXT -1 56 " the solu
tion (red curve) appoaches zero rather quickly." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 41 "x5 := rhs(dsolve(\{subs(a=3.99,ODE),IC\}));" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "and a further investigation shows \+
that the curve crosses the x-axis for the first time at " }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 81 "T := (Pi-arctan(1/sqrt(799)))*800/sqrt(799);
\nevalf(%);\nsubs(t=T,x5);\nsimplify(%);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 35 "Predator/Pr
ey (System of Equations)" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "
restart; with(DEtools):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 192 "A cou
pled system of differential equations is usually more difficult to sol
ve than a single equation.  We illustrate solving systems of different
ial equations using a predator/prey example.   " }{TEXT 289 4 "R(t)" }
{TEXT -1 50 " represents the numbers of rabbits (prey) at time " }
{TEXT 290 1 "t" }{TEXT -1 2 ", " }{TEXT 291 4 "W(t)" }{TEXT -1 34 " ar
e the number of wolves at time " }{TEXT 292 1 "t" }{TEXT -1 237 ".  Th
e rabbit population by itself grows at 8%, wolves by themselves decrea
se at 2%.  The rabbit-wolf interactions decrease the rabbit population
 a and increase the wolf population.  The resulting system of differen
tial equations becomes" }}{PARA 0 "" 0 "" {TEXT -1 5 "     " }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 103 "sys :=  [D(R)(t) =   0.08*R(t) - 0.001*R
(t)*W(t),\n         D(W)(t) = - 0.02*W(t) + 0.00002*R(t)*W(t)];\n" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "Notice that the system does not co
ntain the variable t explicitly.  Such systems are called " }{TEXT 
264 11 "autonomous," }{TEXT -1 44 " and we can obtain a direction fiel
d in the " }{TEXT 265 12 "phase plane " }{TEXT -1 9 "with the " }
{TEXT 293 10 "dfieldplot" }{TEXT -1 9 " command." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 60 "dfieldplot(sys, \{R(t),W(t)\}, t=0..1, R = 0..3000, W
=0..150);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 110 "The equivalent of t
he DEplot command, which sketches selected solution curves in the sing
le equation case, is " }{TEXT 266 13 "phaseportrait" }{TEXT -1 1 "." }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 143 "phaseportrait(sys,[R(t),W(t)],t=0
..180,[ [R(0)=1000,W(0)=40], [R(0)=1000,W(0)=60], [R(0)=1200,W(0)=80]]
,linecolor=blue,thickness=2,stepsize=1);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 70 "The stepsize stipulation was added to obtain a smoother l
ooking curve." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
49 "The closed form solution of the system is a mess " }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 12 "dsolve(sys);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 127 "We resort to numerical solutions.  First we define the system \+
again, and include the initial conditions  R(0)=1000 and W(0)=40." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 131 "SYS :=  [D(R)(t) =   0.08*R(t) - 0
.001*R(t)*W(t),\n         D(W)(t) = - 0.02*W(t) + 0.00002*R(t)*W(t),\n
         R(0)=1000, W(0)=40];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "
Next we define sln to be the resulting numerical solution." }}{PARA 0 
"> " 0 "" {MPLTEXT 1 0 44 "sln := dsolve(SYS,[R(t),W(t)],type=numeric)
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "We can take snapshots at cer
tain instants.  For example,  the states at t=25 or at t=100 are  " }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "sln(25);\nsln(100);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 228 "We see that R grows rather quickly from \+
1000 initially to 2423 at time t=25, and at time t=100 it has already \+
been reduced to 224.   The graphs of the solution tell the dynamics of
 the system.  The plots package is required and " }{TEXT 270 7 "odeplo
t" }{TEXT -1 37 " will display the numerical solution." }{MPLTEXT 1 0 
0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "odeplot(sln,[[t,R(t)],[t,W(t
)]],t=0..500);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 40 "Here is the sol
ution in the phase plane:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "odeplo
t(sln,[R(t),W(t)],0..200);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 14 "Euler's Method" }}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 74 "We begin with a general purpose program which impl
ements Euler's method.  " }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT 278 49 "Euler's Method for a Sngle Differential Equation." }
{TEXT -1 17 "\nInput parameters" }}{PARA 0 "" 0 "" {TEXT -1 29 "f     \+
  differential equation" }}{PARA 0 "" 0 "" {TEXT -1 25 "a,b    interva
l endpoints" }}{PARA 0 "" 0 "" {TEXT -1 29 "y0     initial value, y(a)
=y0" }}{PARA 0 "" 0 "" {TEXT -1 83 "h      step size\nOutput\npts    p
oints containing the apporximation of the solution\n" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 223 "Euler := proc(f,a,b,y0,h)\nlocal N,x,y,k,pt,dy,
pts;\nN := round(b/h)+1;\nx := a; y := y0;\npt(1) := [x,y];\nfor k fro
m 2 to N do\ndy := h*f(x,y);\nx := x+h; y := y + dy;\npt(k) := [x,y];
\nend do;\npts:= [seq(pt(k),k=1..N)];\nend proc:" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT 271 8 "Example:" }{TEXT -1 59 "  Approximate the solution
 of differential equation  y' =  " }{XPPEDIT 18 0 "1-xy^2;" "6#,&\"\"
\"F$*$%#xyG\"\"#!\"\"" }{TEXT -1 167 "  on the interval [0,5] with Eul
er's method.  Use the initial value y(0)=1, and the step size h=0.1. \+
\nFirst we define the right hand side of the differential equation:" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "F := (x,y) -> 1-x*y^2;  " }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 182 "Next we call up Euler's method, w
hich was defined above.  The input arguments (equation, solution inter
val, initial value, step size) have to be entered to match the problem
 at hand." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "pts := Euler(F,0,5,1,0
.1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 159 "The output represents an
 approximation of the solution.  For example we see that y(1) = 1.1231
 apporximately.  A plot of all points sketches this curve nicely." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "plot(pts,style=point);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 159 "As a follow up we display the exact solu
tion along with its approximation in a common graph.  First we invoke \+
the plots package and save the above plot in G1. " }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 28 "G1 := plot(pts,style=point):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 45 "The exact solution involves special functions" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 53 "Y := rhs( dsolve( \{D(y)(x) = F(x,y(x)), y(0)=1 \+
\}) ); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "and we save its graph \+
in G2.  Finally, we display both graphs simultaneously." }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 69 "G2 := plot(Y,x=0..5,color=blue,view=[0..5, 0..
 1.5]):\ndisplay(G1,G2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 251 "This
 was a nice example.  The original approximations move away from the e
xact solution, but instead of amplifying the error in subsequent steps
, the approximations close in on the true solution eventually.  This b
ehaviour is associated with the term " }{TEXT 269 9 "stability" }
{TEXT -1 1 "." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
8 "Here is " }{TEXT 272 15 "another example" }{TEXT -1 460 ".  This ti
me the gap between the exact solution and the approximation widens as \+
t increases.  We also look at the effect of changing the step size.   \+
The problem at hand is the initial value problem  y' = x + y,  y(0) = \+
1.  We test the performance for the step sizes h=0.5 (red) and h=0.1 (
black).  The exact solution is graphed in blue.  Note that with(plots)
 and Euler program are still in effect, since we didn't restart at the
 beginning of the new example." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "F
 := (x,y) -> y+x;  " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 286 "a :
= 0: b:= 2: y0 := 1:  \npts := Euler(F,a,b,y0,0.5):\nG1 := plot(pts,st
yle=point,color=red):\npts := Euler(F,a,b,y0,0.1):\nG2 := plot(pts,sty
le=point,color=black):\nY := rhs( dsolve( \{D(y)(x) = F(x,y(x)), y(a)=
y0 \}) ):\nG3 := plot(Y,x=a..b, color=blue,view=[a..b,0..12]):\ndispla
y(G1,G2,G3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 1 0" 10 }
{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }