{VERSION 3 0 "APPLE_PPC_MAC" "3.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 
0 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 
0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 
{CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 
0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 
256 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 
-1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 257 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 
0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 
258 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 
-1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT -1 11 "Math 2250-3" }}{PARA 
257 "" 0 "" {TEXT -1 22 "Numerical Computations" }}{PARA 258 "" 0 "" 
{TEXT -1 25 "Wednesday, Sept 13,  2000" }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
843 "     In this handout we will study numerical methods for approxim
ating solutions to first order differential equations.  (The fact is, \+
most differential equations do NOT have simple formulas for their solu
tions, despite all the examples you've seen in which they do.  In the \+
case that no nice formula exists for the solution one must approximate
 it numerically.)  A month from now we will see how higher order diffe
rential equations can be converted into first order systems of differe
ntial equations, and that there is a natural way to generalize what we
 are doing now in the context of a single first order differential equ
ation to this more complicated setting.  So understanding today's mate
rial will be an important step in understanding numerical solutions to
 higher order differential equations and to systems of differential eq
uations." }}{PARA 0 "" 0 "" {TEXT -1 86 "     We will be working throu
gh selected material from sections 2.4-2.6 of the text.  " }}{PARA 0 "
" 0 "" {TEXT -1 343 "     The most basic method of approximating solut
ions to differential equations is called Euler's method, after the 170
0's mathematician who first formulated it.  If you want to approximate
 the solution to the initial value problem dy/dx = f(x,y), y(x0)=y0, f
irst pick a step size ``h''.  Then for x between x0 and x0+h, use the \+
constant slope" }}{PARA 0 "" 0 "" {TEXT -1 616 "f(x0,y0).  At x-value \+
x1:=x0+h your y-value will therefore be y1:=y0 + f(x0,y0)h.  Then for \+
x between x1 and x1+h you use the constant slope f(x1,y1), so that at \+
x2:=x1+h your y-value is y2:=y1+f(x1,y1)h.  You continue in this manne
r.  It is easy to visualize if you understand the slope field concept \+
we've been talking about; you just use the slope field where you are a
t the end of each step to get a slope for the next step.  It is straig
htforward to have a programmable calculator or computer software do th
is sort of tedious computation for you.  In Euler's time such computat
ions would have been done by hand!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 467 "     A good first example to illustrate \+
Euler's method is our favorite DE from the time of Calculus, namely dy
/dx =y, say with initial value y(0)=1, so that y=exp(x) is the solutio
n.  Let's take h=0.2 and try to approximate the solution of the x-inte
rval [0,1].  Since the approximate solution will be piecewise affine, \+
we only need to know the approximations at the discrete x values x=0,0
.2,0.4,0.6,0.8,1.   Here's a simple ``do loop'' to make these computat
ions.  " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 46 "restart:  #clear any memory from earlier work " }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 100 "x0:=0.0; xn:=1.0; y0:=1.0; n:=5; h
:=(xn-x0)/n;  \n #specify initial values, number of steps, and size" }
}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x0G\"\"!" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#xnG$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#
y0G$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\"&" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++?!#5" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 100 "f:=(x,y)->y;    #this is the slope function f(x
,y) \n   #in dy/dx = f(x,y), in our example dy/dx = y." }}{PARA 11 "" 
1 "" {XPPMATH 20 "6#>%\"fGR6$%\"xG%\"yG6\"6$%)operatorG%&arrowGF)9%F)F
)F)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 47 "x:=x0; y:=y0;   #initialize x,y for the do loop" }}{PARA 11 "" 
1 "" {XPPMATH 20 "6#>%\"xG\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%
\"yG$\"#5!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for i fro
m 1 to n do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "          k:= f(x,y)
:  #current slope, use : to suppress output" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 45 "          y:= y + h*k: #new y value via Euler" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "          x:= x + h:   #updated x-v
alue:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "          print(x,y,exp(x)
);  #display current values, " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "  \+
                     #and compare to exact solution." }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 47 "    od:                #``od'' ends a do loop  \+
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++?!#5$\"+++++7!\"*$\"+eFS@7
F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"++++S9!\"*$\"+)pC
=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$\"++++G<!\"*$
\"++)=@#=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++!)!#5$\"+++gt?!
\"*$\"+G4aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++5!\"*$\"+++K
)[#F%$\"+G=G=FF%" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 318 "Notice your \+
approximations are all a little too small, in particular your final ap
proximation 2.488... is short of the exact value of exp(1)=e=2.71828..
  The reason for this is that because of the form of our f(x,y) our ap
proximate slope is always less than the actual slope.  We can see this
 graphically using plots: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(
plots):with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "xva
l:=vector(n+1);yval:=vector(n+1);  #to collect all our points" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%xvalG-%&arrayG6$;\"\"\"\"\"'7\"" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%yvalG-%&arrayG6$;\"\"\"\"\"'7\"" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "xval[1]:=x0; yval[1]:=y0;  \+
      #initial values" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%xvalG6#\"
\"\"\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%yvalG6#\"\"\"$\"#5!\"
\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "           #paste in \+
the previous work, and modify for plotting:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 20 "for i from 1 to n do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 33 "          x:=xval[i]:  #current x" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 33 "          y:=yval[i]:  #current y" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 37 "          k:= f(x,y):  #current slope" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 55 "          yval[i+1]:= y + h*k:   #new y v
alue via Euler" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "          xval[i+
1]:= x + h:     #updated x-value:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
57 "              od:                #``od'' ends a do loop  " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "approxsol:=pointplot(\{seq([
xval[i],yval[i]], i=1..n+1)\}):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 49 "exactsol:=plot(exp(t),t=0..1,`color`=`black`):   " }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "     #used t because x was already \+
used above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "display(\{app
roxsol,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 349 262 262 
{PLOTDATA 2 "6&-%'POINTSG6(7$$\"+++++!)!#5$\"+++gt?!\"*7$$\"+++++5F,$
\"+++K)[#F,7$$\"+++++gF)$\"++++G<F,7$$\"+++++?F)$\"+++++7F,7$$\"+++++S
F)$\"++++S9F,7$\"\"!$\"#5!\"\"-%'CURVESG6$7S7$FB$\"\"\"FB7$$\"1nmm;arz
@!#<$\"1GvgZk.A5!#:7$$\"1LL$e9ui2%FP$\"1OdbY\\gT5FS7$$\"1nmm\"z_\"4iFP
$\"1#fT9tfS1\"FS7$$\"1mmmT&phN)FP$\"1Hr`&G_r3\"FS7$$\"1LLe*=)H\\5!#;$
\"1^LEiEj56FS7$$\"1nm\"z/3uC\"F`o$\"1^yZ%yaG8\"FS7$$\"1++DJ$RDX\"F`o$
\"1pPGkJLc6FS7$$\"1nm\"zR'ok;F`o$\"1Am\"\\\\E6=\"FS7$$\"1++D1J:w=F`o$
\"1)e/'[$pj?\"FS7$$\"1LLL3En$4#F`o$\"1SSV5x*GB\"FS7$$\"1nm;/RE&G#F`o$
\"1#3j2pYnD\"FS7$$\"1+++D.&4]#F`o$\"1m9jYu9%G\"FS7$$\"1+++vB_<FF`o$\"1
q#e<$=E78FS7$$\"1+++v'Hi#HF`o$\"1Zw>&\\P*R8FS7$$\"1nm\"z*ev:JF`o$\"1KG
01]dl8FS7$$\"1LLL347TLF`o$\"19H.#p*p'R\"FS7$$\"1LLLLY.KNF`o$\"1Lbfo2iB
9FS7$$\"1++D\"o7Tv$F`o$\"1ZdFH**eb9FS7$$\"1LLL$Q*o]RF`o$\"12/UDl[%[\"F
S7$$\"1++D\"=lj;%F`o$\"1;8&[1^o^\"FS7$$\"1++vV&R<P%F`o$\"1Ksg\"RD$[:FS
7$$\"1LL$e9Ege%F`o$\"1J$))*o>'=e\"FS7$$\"1LLeR\"3Gy%F`o$\"1.c]c%)H8;FS
7$$\"1nm;/T1&*\\F`o$\"1!fh)zw!zk\"FS7$$\"1mm\"zRQb@&F`o$\"1,G!GGVYo\"F
S7$$\"1***\\(=>Y2aF`o$\"16+U5yG<<FS7$$\"1mm;zXu9cF`o$\"1%>jrpbKv\"FS7$
$\"1+++]y))GeF`o$\"1!=v7P07z\"FS7$$\"1****\\i_QQgF`o$\"1OlB#\\E\"H=FS7
$$\"1***\\7y%3TiF`o$\"1fT_<6em=FS7$$\"1****\\P![hY'F`o$\"1V)H1Jn!4>FS7
$$\"1LLL$Qx$omF`o$\"1[cmrs1[>FS7$$\"1+++v.I%)oF`o$\"1RAPIze!*>FS7$$\"1
mm\"zpe*zqF`o$\"1MeRd*=*H?FS7$$\"1+++D\\'QH(F`o$\"1i'*f?z!Q2#FS7$$\"1K
Le9S8&\\(F`o$\"1+43Q,(f6#FS7$$\"1***\\i?=bq(F`o$\"1?.C'Qe4;#FS7$$\"1LL
L3s?6zF`o$\"1(p,G?ne?#FS7$$\"1++DJXaE\")F`o$\"1vYezG)QD#FS7$$\"1nmmm*R
RL)F`o$\"19U**za6,BFS7$$\"1mm;a<.Y&)F`o$\"1^n_#[T/N#FS7$$\"1LLe9tOc()F
`o$\"1g*=(>KS+CFS7$$\"1+++]Qk\\*)F`o$\"1'G]\"H'[sW#FS7$$\"1LL$3dg6<*F`
o$\"1i$4Z;k?]#FS7$$\"1mmmmxGp$*F`o$\"1[Xr/78_DFS7$$\"1++D\"oK0e*F`o$\"
15')HYrh1EFS7$$\"1++v=5s#y*F`o$\"1>xB3j&)fEFS7$FK$\"1X!f%G=G=FFS-%'COL
OURG6&%$RGBGFBFBFB-%+AXESLABELSG6$Q\"t6\"%!G-%%VIEWG6$;FBFK%(DEFAULTG
" 1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 443 "If you connect the Euler dots in your mind, th
e picture above is like the one in figure 2.4.3, on page 109 of the te
xt.  The upper graph is of the exponential function, the lower graph i
s of your Euler approximation.  The reason that the dots lie below the
 true graph is that as y increases the slope f(x,y)=y should also be i
ncreasing, but in the Euler approximation you use the slope at each (l
ower ) point to get to the next (higher) point." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 424 "    It should be that as
 your step size ``h'' gets smaller, your approximations get better to \+
the actual solution.  This is true if your computer can do exact math \+
(which it can't), but in practice you don't want to make the computer \+
do too many computations because of problems with round-off error and \+
computation time, so for example, choosing h=0.0000001 would not be pr
actical. But, trying h=0.01 should be instructive:" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 248 "     Since the width of \+
our x-interval is 1, we stepsize h=0.01 by taking an n-value of subint
ervals to be100.  We keep the other data the same as in our first exam
ple.  The following code only prints out approximations when h is a mu
ltiple of 0.1:" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 48 "x0:=0.0; xn:=1.0; y0:=1.0; n:=100; h:=(xn-x0)/n;" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x0G\"\"!" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#xnG$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#
y0G$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"$+\"" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++5!#6" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 13 "x:=x0; y:=y0;" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%\"xG\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"yG$\"#5!\"\"
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for i from 1 to n do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "          k:= f(x,y):  #current slo
pe" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "          y:= y + h*k: #new y
 value via Euler" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "          x:= x
 + h:   #updated x-value:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "      \+
    if frac(i/10)=0" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "            \+
   then print(x,y,exp(x));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "     \+
     fi;  #use the ``if'' test to decide when to print;" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 60 "               #the command ``frac'' computes \+
the remainder " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "               #o
f a quotient, it will be zero for us if i " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 37 "               #is a multiple of 10. " }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 23 "    od:                " }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6%$\"+++++5!#5$\"+F@i/6!\"*$\"+=4<06F(" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6%$\"+++++?!#5$\"+T+>?7!\"*$\"+eFS@7F(" }}{PARA 11 "" 
1 "" {XPPMATH 20 "6%$\"+++++I!#5$\"+;*[yM\"!\"*$\"+3)e)\\8F(" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"+NP'))[\"!\"*$\"+)pC=\\\"F(" 
}}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++]!#5$\"+B=jW;!\"*$\"+r7s[;F(
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$\"+*p'p;=!\"*$\"++)=@
#=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++q!#5$\"+qLw1?!\"*$\"+2
Fv8?F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++!)!#5$\"+?_r;A!\"*$
\"+G4aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++!*!#5$\"+yEj[C!
\"*$\"+6JgfCF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++5!\"*$\"+MQ
\"[q#F%$\"+G=G=FF%" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "We can also
 see this graphically:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "x
val:=vector(n+1);yval:=vector(n+1);  #to collect all \n    #our points
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%xvalG-%&arrayG6$;\"\"\"\"$,\"7
\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%yvalG-%&arrayG6$;\"\"\"\"$,\"
7\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "xval[1]:=x0; yval[1]
:=y0;        #initial values" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%xv
alG6#\"\"\"\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%yvalG6#\"\"\"$
\"#5!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "           #pa
ste in the previous work, and modify for plotting:" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 20 "for i from 1 to n do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 33 "          x:=xval[i]:  #current x" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 33 "          y:=yval[i]:  #current y" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 37 "          k:= f(x,y):  #current slope" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 55 "          yval[i+1]:= y + h*k:   #new y v
alue via Euler" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "          xval[i+
1]:= x + h:     #updated x-value:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
57 "              od:                #``od'' ends a do loop  " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "approxsol:=pointplot(\{seq([
xval[i],yval[i]], i=1..n+1)\}):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 49 "exactsol:=plot(exp(t),t=0..1,`color`=`black`):   " }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "     #used t because x was already \+
used above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 30 "display(\{approxsol,exactsol\});" }}{PARA 13 
"" 1 "" {GLPLOT2D 349 262 262 {PLOTDATA 2 "6&-%'CURVESG6$7S7$\"\"!$\"
\"\"F(7$$\"1nmm;arz@!#<$\"1GvgZk.A5!#:7$$\"1LL$e9ui2%F.$\"1OdbY\\gT5F1
7$$\"1nmm\"z_\"4iF.$\"1#fT9tfS1\"F17$$\"1mmmT&phN)F.$\"1Hr`&G_r3\"F17$
$\"1LLe*=)H\\5!#;$\"1^LEiEj56F17$$\"1nm\"z/3uC\"FD$\"1^yZ%yaG8\"F17$$
\"1++DJ$RDX\"FD$\"1pPGkJLc6F17$$\"1nm\"zR'ok;FD$\"1Am\"\\\\E6=\"F17$$
\"1++D1J:w=FD$\"1)e/'[$pj?\"F17$$\"1LLL3En$4#FD$\"1SSV5x*GB\"F17$$\"1n
m;/RE&G#FD$\"1#3j2pYnD\"F17$$\"1+++D.&4]#FD$\"1m9jYu9%G\"F17$$\"1+++vB
_<FFD$\"1q#e<$=E78F17$$\"1+++v'Hi#HFD$\"1Zw>&\\P*R8F17$$\"1nm\"z*ev:JF
D$\"1KG01]dl8F17$$\"1LLL347TLFD$\"19H.#p*p'R\"F17$$\"1LLLLY.KNFD$\"1Lb
fo2iB9F17$$\"1++D\"o7Tv$FD$\"1ZdFH**eb9F17$$\"1LLL$Q*o]RFD$\"12/UDl[%[
\"F17$$\"1++D\"=lj;%FD$\"1;8&[1^o^\"F17$$\"1++vV&R<P%FD$\"1Ksg\"RD$[:F
17$$\"1LL$e9Ege%FD$\"1J$))*o>'=e\"F17$$\"1LLeR\"3Gy%FD$\"1.c]c%)H8;F17
$$\"1nm;/T1&*\\FD$\"1!fh)zw!zk\"F17$$\"1mm\"zRQb@&FD$\"1,G!GGVYo\"F17$
$\"1***\\(=>Y2aFD$\"16+U5yG<<F17$$\"1mm;zXu9cFD$\"1%>jrpbKv\"F17$$\"1+
++]y))GeFD$\"1!=v7P07z\"F17$$\"1****\\i_QQgFD$\"1OlB#\\E\"H=F17$$\"1**
*\\7y%3TiFD$\"1fT_<6em=F17$$\"1****\\P![hY'FD$\"1V)H1Jn!4>F17$$\"1LLL$
Qx$omFD$\"1[cmrs1[>F17$$\"1+++v.I%)oFD$\"1RAPIze!*>F17$$\"1mm\"zpe*zqF
D$\"1MeRd*=*H?F17$$\"1+++D\\'QH(FD$\"1i'*f?z!Q2#F17$$\"1KLe9S8&\\(FD$
\"1+43Q,(f6#F17$$\"1***\\i?=bq(FD$\"1?.C'Qe4;#F17$$\"1LLL3s?6zFD$\"1(p
,G?ne?#F17$$\"1++DJXaE\")FD$\"1vYezG)QD#F17$$\"1nmmm*RRL)FD$\"19U**za6
,BF17$$\"1mm;a<.Y&)FD$\"1^n_#[T/N#F17$$\"1LLe9tOc()FD$\"1g*=(>KS+CF17$
$\"1+++]Qk\\*)FD$\"1'G]\"H'[sW#F17$$\"1LL$3dg6<*FD$\"1i$4Z;k?]#F17$$\"
1mmmmxGp$*FD$\"1[Xr/78_DF17$$\"1++D\"oK0e*FD$\"15')HYrh1EF17$$\"1++v=5
s#y*FD$\"1>xB3j&)fEF17$F)$\"1X!f%G=G=FF1-%'COLOURG6&%$RGBGF(F(F(-%'POI
NTSG6aq7$F($\"#5!\"\"7$$\"+++++Q!#5$\"+Qs_f9!\"*7$$\"+++++RFd[l$\"+5D7
u9Fg[l7$$\"+++++OFd[l$\"+&yo2V\"Fg[l7$$\"+++++PFd[l$\"+tk2X9Fg[l7$$\"+
++++NFd[l$\"+dFg;9Fg[l7$$\"+++++MFd[l$\"+()pd-9Fg[l7$$\"+++++KFd[l$\"+
z1%\\P\"Fg[l7$$\"+++++LFd[l$\"+'3!p)Q\"Fg[l7$$\"+++++IFd[l$\"+;*[yM\"F
g[l7$$\"+++++JFd[l$\"+0uKh8Fg[l7$$\"+++++HFd[l$\"+xQ]M8Fg[l7$$\"+++++F
Fd[l$\"+y)3#38Fg[l7$$\"+++++GFd[l$\"+n4H@8Fg[l7$$\"+++++EFd[l$\"+:jD&H
\"Fg[l7$$\"+++++CFd[l$\"+\\Ytp7Fg[l7$$\"+++++DFd[l$\"+&*>V#G\"Fg[l7$$
\"+++++BFd[l$\"+>I;d7Fg[l7$$\"+++++@Fd[l$\"+T>RK7Fg[l7$$\"+++++AFd[l$
\"+gerW7Fg[l7$$\"+++++?Fd[l$\"+T+>?7Fg[l7$$\"+++++>Fd[l$\"+^*3\"37Fg[l
7$$\"+++++<Fd[l$\"+KWI%=\"Fg[l7$$\"+++++=Fd[l$\"+wu9'>\"Fg[l7$$\"+++++
;Fd[l$\"+Y'yD<\"Fg[l7$$\"+++++:Fd[l$\"+c*o4;\"Fg[l7$$\"+++++8Fd[l$\"+
\"G$4Q6Fg[l7$$\"+++++9Fd[l$\"+9UZ\\6Fg[l7$$\"+++++7Fd[l$\"+J]#o7\"Fg[l
7$$\"+++++5Fd[l$\"+F@i/6Fg[l7$$\"+++++6Fd[l$\"+[$oc6\"Fg[l7$$\"+++++!*
!#6$\"+u_o$4\"Fg[l7$$\"+++++!)F\\el$\"+2n&G3\"Fg[l7$$\"+++++gF\\el$\"+
^,_h5Fg[l7$$\"+++++qF\\el$\"+``8s5Fg[l7$$\"+++++]F\\el$\"+]+,^5Fg[l7$$
\"+++++SF\\el$\"+5SgS5Fg[l7$$FdalF\\el$\"+++5?5Fg[l7$$F]^lF\\el$\"++5I
I5Fg[l7$$FadlF\\el$\"++++55Fg[l7$$\"+++++ZFd[l$\"+XME'f\"Fg[l7$$\"++++
+XFd[l$\"+\\2\"[c\"Fg[l7$$\"+++++YFd[l$\"+c)e/e\"Fg[l7$$\"+++++VFd[l$
\"+&zxR`\"Fg[l7$$\"+++++WFd[l$\"+tvJ\\:Fg[l7$$FeflFd[l$\"+NP'))[\"Fg[l
7$$\"+++++TFd[l$\"+sBv.:Fg[l7$$\"+++++UFd[l$\"+'*)*y=:Fg[l7$$\"+++++iF
d[l$\"+.B@`=Fg[l7$$FfelFd[l$\"+*p'p;=Fg[l7$$\"+++++hFd[l$\"+mO'[$=Fg[l
7$$\"+++++fFd[l$\"+.'4()z\"Fg[l7$$\"+++++eFd[l$\"+(f+4y\"Fg[l7$$\"++++
+cFd[l$\"+?)4eu\"Fg[l7$$\"+++++dFd[l$\"+=zEj<Fg[l7$$\"+++++bFd[l$\"+uX
_G<Fg[l7$$\"+++++`Fd[l$\"+6eY%p\"Fg[l7$$\"+++++aFd[l$\"+p/T6<Fg[l7$$\"
+++++[Fd[l$\"+zgA7;Fg[l7$$\"+++++^Fd[l$\"+T\"y5m\"Fg[l7$$\"+++++\\Fd[l
$\"+S$[$G;Fg[l7$$F`flFd[l$\"+B=jW;Fg[l7$$\"+++++vFd[l$\"+q%G\"4@Fg[l7$
$\"+++++wFd[l$\"+b(>-8#Fg[l7$$\"+++++tFd[l$\"+2.dn?Fg[l7$$\"+++++uFd[l
$\"+5gC)3#Fg[l7$$\"+++++sFd[l$\"+9$*4Z?Fg[l7$$\"+++++nFd[l$\"+cZuZ>Fg[
l7$$\"+++++oFd[l$\"+/AAn>Fg[l7$$\"+++++_Fd[l$\"+A*)ox;Fg[l7$$\"+++++mF
d[l$\"+a,YG>Fg[l7$$\"+++++kFd[l$\"+q=Y!*=Fg[l7$$\"+++++lFd[l$\"+*[m$4>
Fg[l7$$\"+++++jFd[l$\"+EWur=Fg[l7$$\"+++++%)Fd[l$\"+ZFs1BFg[l7$$\"++++
+&)Fd[l$\"+u**yHBFg[l7$$\"+++++#)Fd[l$\"+'>r7E#Fg[l7$$\"+++++$)Fd[l$\"
+3R)QG#Fg[l7$$\"+++++\")Fd[l$\"+sB))QAFg[l7$$FaelFd[l$\"+?_r;AFg[l7$$
\"+++++yFd[l$\"+tr.t@Fg[l7$$\"+++++zFd[l$\"+Xvw%>#Fg[l7$$\"+++++pFd[l$
\"+EW*o)>Fg[l7$$\"+++++xFd[l$\"+`>_^@Fg[l7$$F[elFd[l$\"+yEj[CFg[l7$$\"
+++++\"*Fd[l$\"+0!>JZ#Fg[l7$$\"+++++*)Fd[l$\"+!z)QCCFg[l7$$\"+++++')Fd
[l$\"+uy3`BFg[l7$$\"+++++()Fd[l$\"+`(=mP#Fg[l7$$\"+++++))Fd[l$\"+T\\Q+
CFg[l7$$FadlFg[l$\"+MQ\"[q#Fg[l7$$F[flFd[l$\"+qLw1?Fg[l7$$\"+++++rFd[l
$\"+/5$o-#Fg[l7$$\"+++++)*Fd[l$\"+;$=:l#Fg[l7$$\"+++++**Fd[l$\"+*\\L!y
EFg[l7$$\"+++++(*Fd[l$\"+fcEDEFg[l7$$\"+++++&*Fd[l$\"+av`tDFg[l7$$\"++
+++'*Fd[l$\"+IHF*f#Fg[l7$$\"+++++$*Fd[l$\"+(pGG_#Fg[l7$$\"+++++%*Fd[l$
\"+%)p0[DFg[l7$$\"+++++#*Fd[l$\"+&>]y\\#Fg[l-%+AXESLABELSG6$Q\"t6\"%!G
-%%VIEWG6$;F(F)%(DEFAULTG" 1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 
45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 566 "     Actually, considering how many computations you d
id with n=100 you are not so close to the exact solution. In more comp
licated problems it is a very serious issue to find relatively efficie
nt ways of approximating solutions.  An entire field of mathematics, `
`numerical analysis'' deals with such issues for a variety of mathemat
ical problems.   The book talks about some improvements to Euler  in s
ections 2.5 and 2.6.    If you are interested in this important field \+
of mathematics you should read 2.5 and 2.6 carefully. Let's summarize \+
some highlights below." }}{PARA 0 "" 0 "" {TEXT -1 98 "      For any t
ime step h the fundamental theorem of calculus asserts that, since dy/
dx = f(x,y(x)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "y(x+h)=y(x) + Int(dy/dt,t=
x..x+h);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%\"yG6#,&%\"xG\"\"\"%\"h
GF),&-F%6#F(F)-%$IntG6$*&%#dyG\"\"\"%#dtG!\"\"/%\"tG;F(F'F)" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "y(x+h)=y(x) + Int(f(t,y(t)),
t=x..x+h);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%\"yG6#,&%\"xG\"\"\"%
\"hGF),&-F%6#F(F)-%$IntG6$-%\"fG6$%\"tG-F%6#F4/F4;F(F'F)" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 743 "The problem with Euler is that we always
 approximated this integral by h*f(x,y(x)), i.e. we used the left-hand
 endpoint as our approximation of the ``average height''.  The improve
ments to Euler depend on better approximations to that integral.  ``Im
proved Euler'' uses an approximation to the Trapezoid Rule to approxim
ate the integral.  The trapezoid rule for the integral approximation w
ould be (1/2)*h*(f(x,y(x))+f((x+h),y(x+h)).  Since we don't know y(x+h
) we approximate it using unimproved Euler, and then feed that into th
e trapezoid rule.  This leads to the improved Euler do loop below.  Of
 course before you use it you must make sure you initialize everything
 correctly. We'll compare it when n=5, to our first (unimproved) attem
pt." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 69 "x0:=0; y0:=1;xn:=1.0;n:=5;\nh:=(xn-x0)/n; \nx:=x0; y:=y0; \nf:=(
x,y)->y;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x0G\"\"!" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#>%#y0G\"\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#
>%#xnG$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\"&" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++?!#5" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#>%\"xG\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"yG
\"\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGR6$%\"xG%\"yG6\"6$%)op
eratorG%&arrowGF)9%F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
20 "for i from 1 to n do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "    k1:
=f(x,y):          #left-hand slope" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
59 "    k2:=f(x+h,y+h*k1):   #approximation to right-hand slope" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "    k:= (k1+k2)/2:       #approxima
tion to average slope" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "    y:= y+
h*k:            #improved Euler update" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 34 "    x:= x+h:             #update x" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "    print(x,y,exp(x));" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 5 "od:  " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++?!#5$\"++++?7!
\"*$\"+eFS@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"+++S)
[\"!\"*$\"+)pC=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$
\"++![e\"=!\"*$\"++)=@#=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++
!)!#5$\"+gXL:A!\"*$\"+G4aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"++
+++5!\"*$\"+j\"3Fq#F%$\"+G=G=FF%" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "Notice we al
most did as well with improved Euler when n=5 as we did with n=100 in \+
unimproved Euler. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 337 "     One can also use Taylor approximation methods to im
prove upon Euler; by differentiating the equation dy/dx = f(x,y) one c
an solve for higher order derivatives of y in terms of the lower order
 ones, and then use the Taylor approximation for y(x+h) in terms of y(
x).  See the book for more details of this method, we won't do it here
." }}{PARA 0 "" 0 "" {TEXT -1 207 "     In the same vein as ``improved
 Euler'' we can use the Simpson approximation for the integral instead
 of the Trapezoid one, and this leads to the Runge-Kutta method which \+
is widely used in real software." }}{PARA 0 "" 0 "" {TEXT -1 470 " (Yo
u may or may not have talked about Simpson's Rule in Calculus, it is b
ased on a quadratic approximation to the function f, whereas the Trape
zoid rule is based on a first order approximation.)  Here is the code \+
for the Runge-Kutta method.  The text explains it in section 2.6. Simp
son's rule approximates an integral in terms of the integrand values a
t each endpoint and at the interval midpoint.  Runge-Kutta uses two di
fferent approximations for the midpoint value." }}{PARA 0 "" 0 "" 
{TEXT -1 2 "  " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 
"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "x:=x0; y:=y0; n:=5; h:
=(xn-x0)/n;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"xG\"\"!" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#>%\"yG\"\"\"" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%\"nG\"\"&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++?!#
5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for i from 1 to n do" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "    k1:=f(x,y):              #lef
t-hand slope" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "    k2:=f(x+h/2,y+h
*k1/2):   #1st guess at midpoint slope" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 60 "    k3:=f(x+h/2,y+h*k2/2):   #second guess at midpoint slope" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "    k4:=f(x+h,y+h*k3):       #gue
ss at right-hand slope" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "    k:=(k
1+2*k2+2*k3+k4)/6:  #Simpson's approximation for the integral" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "    x:=x+h:                  #x upd
ate" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "    y:=y+h*k:               \+
 #y update" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "    print(x,y,exp(x))
;              #display current values" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 7 "  od:  " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++?!#5$\"+++S@
7!\"*$\"+eFS@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"+gz
\"=\\\"!\"*$\"+)pC=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g
!#5$\"+ck5A=!\"*$\"++)=@#=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++
++!)!#5$\"+D3_DA!\"*$\"+G4aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"
+++++5!\"*$\"+O6D=FF%$\"+G=G=FF%" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "Notice how cl
ose Runge-Kutta gets you to the correct value of e, with n=5." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 207 "Even wit
h code like Runge-Kutta there are subtleties and problems which partic
ular problems will cause.  We will not go into those here; there are g
ood examples in sections 2.4-2.6 and the homework problems." }}}}
{MARK "36 10 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 }