{VERSION 5 0 "SUN SPARC SOLARIS" "5.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 Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 " " 0 1 0 0 255 1 0 0 0 0 0 0 1 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 1 }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 1 }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 1 }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 1 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 256 "" 0 "" {TEXT -1 9 "Math 2250" }}{PARA 257 "" 0 "" {TEXT -1 22 "Numerical Computations" }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 843 " In this handout we will study nume rical methods for approximating solutions to first order differential \+ equations. (The fact is, most differential equations do NOT have simp le formulas for their solutions, despite all the examples you've seen \+ in which they do. In the case that no nice formula exists for the sol ution one must approximate it numerically.) A month from now we will \+ see how higher order differential equations can be converted into firs t order systems of differential equations, and that there is a natural way to generalize what we are doing now in the context of a single fi rst order differential equation to this more complicated setting. So \+ understanding today's material will be an important step in understand ing numerical solutions to higher order differential equations and to \+ systems of differential equations." }}{PARA 0 "" 0 "" {TEXT -1 86 " \+ We will be working through selected material from sections 2.4-2.6 o f the text. " }}{PARA 0 "" 0 "" {TEXT -1 960 " The most basic met hod of approximating solutions to differential equations is called Eul er's method, after the 1700's mathematician who first formulated it. \+ If you want to approximate the solution to the initial value problem d y/dx = f(x,y), y(x0)=y0, first pick a step size ``h''. Then for x bet ween x0 and x0+h, use the constant slope f(x0,y0). At x-value x1:=x0+ h your y-value will therefore be y1:=y0 + f(x0,y0)h. Then for x betwe en 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 manner. It \+ is easy to visualize if you understand the slope field concept we've b een talking about; you just use the slope field where you are at the e nd of each step to get a slope for the next step. It is straightforwa rd to have a programmable calculator or computer software do this sort of tedious computation for you. In Euler's time such computations wo uld 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 solution. Le t's take h=0.2 and try to approximate the solution of the x-interval [ 0,1]. Since the approximate solution will be piecewise affine, we onl y 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 computations. \+ " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {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$\"\"!F&" }}{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 functio n f(x,y) \n #in dy/dx = f(x,y), in our example dy/dx = y." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGf*6$%\"xG%\"yG6\"6$%)operatorG%&arrowG F)9%F)F)F)" }}}{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$\"\"!F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"yG$\"#5!\"\"" }} }{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 24 "print('x','y','exp(x)');" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for i from 1 to n do" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 63 " k:= f(x,y): #current slope, use : t o 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-value:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 " print(x,y,exp(x)); #display current values , " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 " #and c ompare to exact solution." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 " od : #``od'' ends a do loop " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%%\"xG%\"yG-%$expG6#F#" }}{PARA 11 "" 1 "" {XPPMATH 20 " 6%$\"+++++A!\"*$\"+1P3IuF%$\"+*\\8]-*F%" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++C!\"*$\"+Z/5;*)F%$\"+QwJ-6!\")" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++E!\"*$\"+a?$*p5!\")$\"+/QPY8F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++G!\"*$\"+l%=RG\"!\")$\"+xYYW;F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++I!\"*$\"+e@qS:!\")$\"+#p`&3?F(" }}} {PARA 0 "" 0 "" {TEXT -1 318 "Notice your approximations are all a lit tle too small, in particular your final approximation 2.488... is shor t of the exact value of exp(1)=e=2.71828.. The reason for this is tha t because of the form of our f(x,y) our approximate slope is always le ss than the actual slope. We can see this graphically using plots: " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(plots):with(linalg): " }}{PARA 7 "" 1 "" {TEXT -1 80 "Warning, the protected names norm and trace have been redefined and unprotected\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "xval:=vector(n+1);yval:=vector(n+1); #to collec t 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#\"\"\"$\"\"!F)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%yvalG6#\"\"\"$\"#5!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 " #paste in the previous work, and modify fo r 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 " yva l[i+1]:= y + h*k: #new y value 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(\{approxsol,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 349 262 262 {PLOTDATA 2 "6&-%'CURVESG6$7S7$$\"\"!F)$\"\" \"F)7$$\"3dmmm;arz@!#>$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%F/$\"3aNdbY\\g T5F27$$\"3nmmm\"z_\"4iF/$\"3#=fT9tfS1\"F27$$\"3[mmmT&phN)F/$\"3;Hr`&G_ r3\"F27$$\"3BLLe*=)H\\5!#=$\"3L^LEiEj56F27$$\"3fmm\"z/3uC\"FE$\"3q]yZ% yaG8\"F27$$\"3%)***\\7LRDX\"FE$\"3')oPGkJLc6F27$$\"3]mm\"zR'ok;FE$\"3- Am\"\\\\E6=\"F27$$\"3v***\\i5`h(=FE$\"36)e/'[$pj?\"F27$$\"3WLLL3En$4#F E$\"3]RSV5x*GB\"F27$$\"3pmm;/RE&G#FE$\"3)>3j2pYnD\"F27$$\"3\")*****\\K ]4]#FE$\"3Sm9jYu9%G\"F27$$\"3$******\\PAvr#FE$\"3%*p#e<$=E78F27$$\"3)* *****\\nHi#HFE$\"3PYw>&\\P*R8F27$$\"3jmm\"z*ev:JFE$\"3@KG01]dl8F27$$\" 3?LLL347TLFE$\"3?9H.#p*p'R\"F27$$\"3+LLLLY.KNFE$\"35Lbfo2iB9F27$$\"3v* **\\7o7Tv$FE$\"3\"ouv#H**eb9F27$$\"3&GLLLQ*o]RFE$\"3r1/UDl[%[\"F27$$\" 3@++D\"=lj;%FE$\"3\\;8&[1^o^\"F27$$\"31++vV&R'=e\"F27$$\"3GLLeR\"3Gy%FE$\"3+.c]c%)H8;F27 $$\"3cmm;/T1&*\\FE$\"3#**eh)zw!zk\"F27$$\"3%em;zRQb@&FE$\"3G,G!GGVYo\" F27$$\"3[***\\(=>Y2aFE$\"3L6+U5yG<jrpbKv\"F 27$$\"3k******\\y))GeFE$\"33!=v7P07z\"F27$$\"3'*)***\\i_QQgFE$\"3XOlB# \\E\"H=F27$$\"3@***\\7y%3TiFE$\"3*)eT_<6em=F27$$\"35****\\P![hY'FE$\"3 DV)H1Jn!4>F27$$\"3jKLL$Qx$omFE$\"3uZcmrs1[>F27$$\"3!)*****\\P+V)oFE$\" 39RAPIze!*>F27$$\"3?mm\"zpe*zqFE$\"3hLeRd*=*H?F27$$\"3%)*****\\#\\'QH( FE$\"3zh'*f?z!Q2#F27$$\"3GKLe9S8&\\(FE$\"3^**33Q,(f6#F27$$\"3R***\\i?= bq(FE$\"3m>.C'Qe4;#F27$$\"3\"HLL$3s?6zFE$\"3C(p,G?ne?#F27$$\"3a***\\7` Wl7)FE$\"3EvYezG)QD#F27$$\"3#pmmm'*RRL)FE$\"3^8U**za6,BF27$$\"3Qmm;a<. Y&)FE$\"3+^n_#[T/N#F27$$\"3KS+CF27$$\"3t******\\ Qk\\*)FE$\"3t&G]\"H'[sW#F27$$\"3CLL$3dg6<*FE$\"3bh$4Z;k?]#F27$$\"3Hmmm mxGp$*FE$\"3d[Xr/78_DF27$$\"3A++D\"oK0e*FE$\"3j4')HYrh1EF27$$\"3A++v=5 s#y*FE$\"3!*=xB3j&)fEF27$F*$\"3`X!f%G=G=FF2-%'COLOURG6&%$RGBGF)F)F)-%' POINTSG6(7$F($\"#5!\"\"7$$\"+++++?!#5$\"+++++7!\"*7$$\"+++++SFe[l$\"++ ++S9Fh[l7$$\"+++++gFe[l$\"++++G " 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$\" \"!F&" }}{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 10 "y0 := 1.0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#y0G$\"#5!\"\"" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 13 "x:=x0; y:=y0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>% \"xG$\"\"!F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"yG$\"#5!\"\"" }}} {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$\"+E@i/6!\"*$\"+=4<06F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"+++++?!#5$\"+S+>?7!\"*$\"+eFS@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++I!#5$\"+:*[yM\"!\"*$\"+3)e)\\8F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"+MP'))[\"!\"*$\"+)pC=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++]!#5$\"+A=jW;!\"*$\"+r7s[;F( " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$\"+)p'p;=!\"*$\"++)=@ #=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++q!#5$\"+pLw1?!\"*$\"+2 Fv8?F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++!)!#5$\"+>_r;A!\"*$ \"+G4aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++!*!#5$\"+xEj[C! \"*$\"+6JgfCF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++5!\"*$\"+LQ \"[q#F%$\"+G=G=FF%" }}}{PARA 0 "" 0 "" {TEXT -1 33 "We can also see th is graphically:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "xval:=vec tor(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#>&%%xvalG6# \"\"\"$\"\"!F)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%yvalG6#\"\"\"$\" #5!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 " #past e 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)$ \"\"\"F)7$$\"3dmmm;arz@!#>$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%F/$\"3aNdb Y\\gT5F27$$\"3nmmm\"z_\"4iF/$\"3#=fT9tfS1\"F27$$\"3[mmmT&phN)F/$\"3;Hr `&G_r3\"F27$$\"3BLLe*=)H\\5!#=$\"3L^LEiEj56F27$$\"3fmm\"z/3uC\"FE$\"3q ]yZ%yaG8\"F27$$\"3%)***\\7LRDX\"FE$\"3')oPGkJLc6F27$$\"3]mm\"zR'ok;FE$ \"3-Am\"\\\\E6=\"F27$$\"3v***\\i5`h(=FE$\"36)e/'[$pj?\"F27$$\"3WLLL3En $4#FE$\"3]RSV5x*GB\"F27$$\"3pmm;/RE&G#FE$\"3)>3j2pYnD\"F27$$\"3\")**** *\\K]4]#FE$\"3Sm9jYu9%G\"F27$$\"3$******\\PAvr#FE$\"3%*p#e<$=E78F27$$ \"3)******\\nHi#HFE$\"3PYw>&\\P*R8F27$$\"3jmm\"z*ev:JFE$\"3@KG01]dl8F2 7$$\"3?LLL347TLFE$\"3?9H.#p*p'R\"F27$$\"3+LLLLY.KNFE$\"35Lbfo2iB9F27$$ \"3v***\\7o7Tv$FE$\"3\"ouv#H**eb9F27$$\"3&GLLLQ*o]RFE$\"3r1/UDl[%[\"F2 7$$\"3@++D\"=lj;%FE$\"3\\;8&[1^o^\"F27$$\"31++vV&R'=e\"F27$$\"3GLLeR\"3Gy%FE$\"3+.c]c%)H 8;F27$$\"3cmm;/T1&*\\FE$\"3#**eh)zw!zk\"F27$$\"3%em;zRQb@&FE$\"3G,G!GG VYo\"F27$$\"3[***\\(=>Y2aFE$\"3L6+U5yG<jrpb Kv\"F27$$\"3k******\\y))GeFE$\"33!=v7P07z\"F27$$\"3'*)***\\i_QQgFE$\"3 XOlB#\\E\"H=F27$$\"3@***\\7y%3TiFE$\"3*)eT_<6em=F27$$\"35****\\P![hY'F E$\"3DV)H1Jn!4>F27$$\"3jKLL$Qx$omFE$\"3uZcmrs1[>F27$$\"3!)*****\\P+V)o FE$\"39RAPIze!*>F27$$\"3?mm\"zpe*zqFE$\"3hLeRd*=*H?F27$$\"3%)*****\\# \\'QH(FE$\"3zh'*f?z!Q2#F27$$\"3GKLe9S8&\\(FE$\"3^**33Q,(f6#F27$$\"3R** *\\i?=bq(FE$\"3m>.C'Qe4;#F27$$\"3\"HLL$3s?6zFE$\"3C(p,G?ne?#F27$$\"3a* **\\7`Wl7)FE$\"3EvYezG)QD#F27$$\"3#pmmm'*RRL)FE$\"3^8U**za6,BF27$$\"3Q mm;a<.Y&)FE$\"3+^n_#[T/N#F27$$\"3KS+CF27$$\"3t** ****\\Qk\\*)FE$\"3t&G]\"H'[sW#F27$$\"3CLL$3dg6<*FE$\"3bh$4Z;k?]#F27$$ \"3HmmmmxGp$*FE$\"3d[Xr/78_DF27$$\"3A++D\"oK0e*FE$\"3j4')HYrh1EF27$$\" 3A++v=5s#y*FE$\"3!*=xB3j&)fEF27$F*$\"3`X!f%G=G=FF2-%'COLOURG6&%$RGBGF) F)F)-%'POINTSG6aq7$$\"+++++$*!#5$\"+'pGG_#!\"*7$$\"+++++%*Fa[l$\"+$)p0 [DFd[l7$$\"+++++5!#6$\"++++55Fd[l7$$\"+++++?F]\\l$\"+++5?5Fd[l7$$\"+++ ++IF]\\l$\"++5II5Fd[l7$$\"+++++CFa[l$\"+[Ytp7Fd[l7$$\"+++++DFa[l$\"+%* >V#G\"Fd[l7$$\"+++++EFa[l$\"+9jD&H\"Fd[l7$$\"+++++FFa[l$\"+x)3#38Fd[l7 $$\"+++++GFa[l$\"+m4H@8Fd[l7$$\"+++++HFa[l$\"+wQ]M8Fd[l7$$Fb\\lFa[l$\" +S+>?7Fd[l7$$\"+++++@Fa[l$\"+S>RK7Fd[l7$$\"+++++AFa[l$\"+ferW7Fd[l7$$ \"+++++BFa[l$\"+=I;d7Fd[l7$$\"+++++\"Fd[l7$$\"+++++>Fa[l$\"+]*3\"37Fd[l7$$F\\\\lFa[l$\"+E@i/6Fd[ l7$$\"+++++6Fa[l$\"+Z$oc6\"Fd[l7$$\"+++++7Fa[l$\"+I]#o7\"Fd[l7$$\"++++ +8Fa[l$\"+!G$4Q6Fd[l7$$\"+++++9Fa[l$\"+8UZ\\6Fd[l7$$\"+++++:Fa[l$\"+b* o4;\"Fd[l7$$\"+++++;Fa[l$\"+X'yD<\"Fd[l7$$\"+++++!*F]\\l$\"+t_o$4\"Fd[ l7$$\"+++++SF]\\l$\"+5SgS5Fd[l7$$\"+++++]F]\\l$\"+]+,^5Fd[l7$$\"+++++g F]\\l$\"+],_h5Fd[l7$$\"+++++qF]\\l$\"+_`8s5Fd[l7$$\"+++++!)F]\\l$\"+1n &G3\"Fd[l7$F($\"#5!\"\"7$$\"+++++lFa[l$\"+)[m$4>Fd[l7$$\"+++++mFa[l$\" +`,YG>Fd[l7$$\"+++++nFa[l$\"+bZuZ>Fd[l7$$\"+++++oFa[l$\"+.AAn>Fd[l7$$ \"+++++pFa[l$\"+DW*o)>Fd[l7$$FbdlFa[l$\"+pLw1?Fd[l7$$\"+++++rFa[l$\"+. 5$o-#Fd[l7$$\"+++++eFa[l$\"+'f+4y\"Fd[l7$$\"+++++fFa[l$\"+-'4()z\"Fd[l 7$$F]dlFa[l$\"+)p'p;=Fd[l7$$\"+++++hFa[l$\"+lO'[$=Fd[l7$$\"+++++iFa[l$ \"+-B@`=Fd[l7$$\"+++++jFa[l$\"+DWur=Fd[l7$$\"+++++kFa[l$\"+p=Y!*=Fd[l7 $$FhclFa[l$\"+A=jW;Fd[l7$$\"+++++^Fa[l$\"+S\"y5m\"Fd[l7$$\"+++++_Fa[l$ \"+@*)ox;Fd[l7$$\"+++++`Fa[l$\"+5eY%p\"Fd[l7$$\"+++++aFa[l$\"+o/T6)4eu\"Fd[l7$$\"+++++d Fa[l$\"+ JZ#Fd[l7$$\"+++++#*Fa[l$\"+%>]y\\#Fd[l7$$\"+++++$)Fa[l$\"+2R)QG#Fd[l7$ $\"+++++%)Fa[l$\"+YFs1BFd[l7$$\"+++++&)Fa[l$\"+t**yHBFd[l7$$\"+++++')F a[l$\"+ty3`BFd[l7$$\"+++++()Fa[l$\"+_(=mP#Fd[l7$$\"+++++))Fa[l$\"+S\\Q +CFd[l7$$\"+++++*)Fa[l$\"+*y)QCCFd[l7$$F^clFa[l$\"+xEj[CFd[l7$$\"+++++ wFa[l$\"+a(>-8#Fd[l7$$\"+++++xFa[l$\"+_>_^@Fd[l7$$\"+++++yFa[l$\"+sr.t @Fd[l7$$\"+++++zFa[l$\"+Wvw%>#Fd[l7$$FgdlFa[l$\"+>_r;AFd[l7$$\"+++++sF a[l$\"+8$*4Z?Fd[l7$$\"+++++tFa[l$\"+1.dn?Fd[l7$$\"+++++uFa[l$\"+4gC)3# Fd[l7$$\"+++++vFa[l$\"+p%G\"4@Fd[l7$$\"+++++\\Fa[l$\"+R$[$G;Fd[l7$$\"+ ++++\")Fa[l$\"+rB))QAFd[l7$$\"+++++#)Fa[l$\"+&>r7E#Fd[l-%+AXESLABELSG6 $Q\"t6\"Q!Fcjm-%%VIEWG6$;F(F*%(DEFAULTG" 1 2 0 1 10 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 0 "Curve 1" "Curve 2" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 566 " Actually, consid ering how many computations you did with n=100 you are not so close to the exact solution. In more complicated problems it is a very serious issue to find relatively efficient ways of approximating solutions. \+ An entire field of mathematics, ``numerical analysis'' deals with such issues for a variety of mathematical problems. The book talks about some improvements to Euler in sections 2.5 and 2.6. If you are in terested in this important field of mathematics you should read 2.5 an d 2.6 carefully. Let's summarize some highlights below." }}{PARA 0 "" 0 "" {TEXT -1 98 " For any time 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\"\"\"%\"hGF),&-F%6#F(F)-%$IntG6$*&%#dyGF)%#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)" }}}{PARA 0 "" 0 "" {TEXT -1 743 "The problem with Euler is th at we always approximated this integral by h*f(x,y(x)), i.e. we used t he left-hand endpoint as our approximation of the ``average height''. \+ The improvements to Euler depend on better approximations to that int egral. ``Improved Euler'' uses an approximation to the Trapezoid Rule to approximate the integral. The trapezoid rule for the integral app roximation would 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 the trapezoid rule. This leads to the improved Euler do loo p below. Of course before you use it you must make sure you initializ e everything correctly. We'll compare it when n=5, to our first (unimp roved) attempt." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "x0:=0; y0:=1;xn:=1.0;n:=5;\nh:=(xn-x0)/n; \nx:=x 0; 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#>%\"fGf *6$%\"xG%\"yG6\"6$%)operatorG%&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: #approximation 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%" }}}{PARA 0 "" 0 "" {TEXT -1 100 "Notice we almost did as well wit h 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 " O ne can also use Taylor approximation methods to improve upon Euler; by differentiating the equation dy/dx = f(x,y) one can 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 fo r 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 o ne, and this leads to the Runge-Kutta method which is widely used in r eal software." }}{PARA 0 "" 0 "" {TEXT -1 470 " (You may or may not ha ve talked about Simpson's Rule in Calculus, it is based on a quadratic approximation to the function f, whereas the Trapezoid 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. Simpson's rule approxim ates an integral in terms of the integrand values at each endpoint and at the interval midpoint. Runge-Kutta uses two different approximati ons 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): #left-hand slope" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 " k2:=f(x+h/2,y+h*k1/2): #1st g uess 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): #guess at right-hand \+ slope" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 " k:=(k1+2*k2+2*k3+k4)/6 : #Simpson's approximation for the integral" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 " x:=x+h: #x update" }}{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@7 F(" }}{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 "" }}} {PARA 0 "" 0 "" {TEXT -1 74 "Notice how close 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 with code like Runge-Kutta there ar e subtleties and problems which particular problems will cause. We wi ll not go into those here; there are good examples in sections 2.4-2.6 and the homework problems." }}}{MARK "2 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }