{VERSION 4 0 "SUN SPARC SOLARIS" "4.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 } {CSTYLE "" -1 256 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 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 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 }{PSTYLE "" 0 258 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 256 11 "Math 2280-2" }}{PARA 257 "" 0 "" {TEXT 257 23 "Maple Project 1, Part 2" }}{PARA 258 "" 0 "" {TEXT -1 16 "January 25, 2001" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 60 " You should download this file from our Maple ho me page " }{TEXT 259 62 "http://www.math.utah.edu/~korevaar/2280spring 01/2280maple.html" }{TEXT -1 190 ", and then open it from Maple. It c ontains discussion and Maple commands which will help you answer the q uestions on the Part 2 solution template, which is also available at t he Maple page." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 629 " In this part of the project we will study numerical methods for approximating solutions to first order differential equat ions. We will see soon how higher order differential equations can be converted into first order systems of differential equations. It tu rns out that there is a natural way to generalize what we do here in t he context of a single first order differential equations, to systems \+ of first order differential equations. So understanding this project \+ material will be an important step in understanding numerical solutio ns to higher order differential equations and to systems of differenti al equations." }}{PARA 0 "" 0 "" {TEXT -1 75 " We will be working \+ through material from sections 2.4-2.6 of the text." }}{PARA 0 "" 0 " " {TEXT -1 1036 " The most basic method of approximating solutions to differential equations is called Euler's method, after the 1700'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, first pick a step size ``h''. Then for x between x0 and x0+h, use the cons tant 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 between x1 and x1+h you use the con stant 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 un derstand the slope field concept we've been talking about; you just us e the slope field with finite rather than infinitesimal stepping in th e x-variable. You use the value of the slope field at your current po int to get a slope which you use to move to the next point. It is str aightforward to have the computer do this sort of tedious computation \+ for you. In Euler's time such computations would have been done by ha nd!" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{TEXT 258 1 " " }{TEXT -1 464 " A good first example to illustrate Euler's method is our favorite D E from the time of Calculus, namely dy/dx =y, say with initial value y (0)=1, so that y=exp(x) is the solution. Let's take h=0.2 and try to \+ approximate the solution of the x-interval [0,1]. Since the approxima te solution will be piecewise affine, we only need to know the approxi mations at the discrete x values x=0,0.2,0.4,0.6,0.8,1. Here's a simp le ``do loop'' to make these computations. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "restart: #clear a ny memory from earlier work " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 150 "x0:=0.0; xn:=1.0; y0:=1.0; n:=5; h:=(xn-x0)/n; \n # specify \+ initial\n # values, number of steps, step 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; \n #this \+ is the \"slope\" function f(x,y) \n #in dy/dx = f(x,y). We want dy/ dx = y." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGR6$%\"xG%\"yG6\"6$%)o peratorG%&arrowGF)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!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 340 "for i fr om 1 to n do\n k:= f(x,y): #current slope, use : to suppress output\n y:= y + h*k: #new y value via Euler\n x:= \+ x + h: #updated x-value:\n print(x,y,exp(x)); \n \+ #display current values, \n #and compare to exact solutio n.\n od: \n #``od'' ends a do loop " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++?!#5$\"+++++7!\"*$\"+eFS@7F(" }} {PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++S!#5$\"++++S9!\"*$\"+)pC=\\\"F (" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$\"++++G " 0 "" {MPLTEXT 1 0 25 "with(plots):with(linalg):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 68 "xval:=vector(n+1);yval:=vector(n+1); \n #to collect all our points" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%xvalG-%& arrayG6$;\"\"\"\"\"'7\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%yvalG-%& arrayG6$;\"\"\"\"\"'7\"" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "xval[1]:=x0; yval[1]:=y0; \n #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 347 " #paste in the pr evious work, and modify for plotting:\nfor i from 1 to n do\n \+ x:=xval[i]: #current x\n y:=yval[i]: #current y\n \+ k:= f(x,y): #current slope\n yval[i+1]:= y + h*k: #new y \+ value via Euler\n xval[i+1]:= x + h: #updated x-value:\n \+ 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 95 "exactsol:=plot(exp(t),t=0..1,`color`=`black`): \n #used t becau se x was already used above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "display(\{approxsol,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 298 257 257 {PLOTDATA 2 "6&-%'CURVESG6$7S7$$\"\"!F)$\"\"\"F)7$$\"3dmmm ;arz@!#>$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%F/$\"3aNdbY\\gT5F27$$\"3nmmm \"z_\"4iF/$\"3#=fT9tfS1\"F27$$\"3[mmmT&phN)F/$\"3;Hr`&G_r3\"F27$$\"3BL Le*=)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=\"F 27$$\"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$\"3Sm9j Yu9%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;%F E$\"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\"F27$$\"3k**** **\\y))GeFE$\"33!=v7P07z\"F27$$\"3'*)***\\i_QQgFE$\"3XOlB#\\E\"H=F27$$ \"3@***\\7y%3TiFE$\"3*)eT_<6em=F27$$\"35****\\P![hY'FE$\"3DV)H1Jn!4>F2 7$$\"3jKLL$Qx$omFE$\"3uZcmrs1[>F27$$\"3!)*****\\P+V)oFE$\"39RAPIze!*>F 27$$\"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$\"3Ev YezG)QD#F27$$\"3#pmmm'*RRL)FE$\"3^8U**za6,BF27$$\"3Qmm;a<.Y&)FE$\"3+^n _#[T/N#F27$$\"3KS+CF27$$\"3t******\\Qk\\*)FE$\"3 t&G]\"H'[sW#F27$$\"3CLL$3dg6<*FE$\"3bh$4Z;k?]#F27$$\"3HmmmmxGp$*FE$\"3 d[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)-%'POINTSG6(7$$ \"+++++S!#5$\"++++S9!\"*7$$\"+++++gFa[l$\"++++G " 0 "" {MPLTEXT 1 0 62 "x0:=0.0; xn:=1.0; y0:=1.0; n:=100; h:=(xn-x0)/n;\nx:=x0; y:=y0;" } }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 591 "for i from 1 to n do\n k:= f(x,y): #current slope\n y:= y + h*k: #new y value via Euler\n x:= x + h: \+ #updated x-value:\n if frac(i/10)=0\n then prin t(x,y,exp(x));\n fi; #use the ``if'' test to decide when to \+ print;\n #the command ``frac'' computes the remainder \n #of a quotient, it will be zero for us if i \n \+ #is a multiple of 10. This way we only print\n #th e approximations every 0.1 increments, even\n #though ou r actual time step is 0.01.\n 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;!\"*$\"+r 7s[;F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++g!#5$\"+)p'p;=!\"*$ \"++)=@#=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++q!#5$\"+pLw1?! \"*$\"+2Fv8?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 1 " " }{TEXT -1 268 "So you can see we got closer to the actual value of e, but really , considering how much work we did this was not a great result. We ca n make a picture of what we did as follows, using the mouse to cut and paste previous work, and then editing it for the new situation" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "restart:with(plots):with(lin alg):" }}{PARA 7 "" 1 "" {TEXT -1 50 "Warning, the name changecoords h as been redefined\n" }}{PARA 7 "" 1 "" {TEXT -1 80 "Warning, the prote cted names norm and trace have been redefined and unprotected\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 487 "\nf:=(x,y)->y;\nn:=100; x0: =0.0; y0:=1.0;\nxn:=1.0; h:=(xn-x0)/n;\nxval:=vector(n+1);yval:=vector (n+1); \n #to collect all our points. Now n=100\nxval[1]:=x0; yval[ 1]:=y0; \n #initial values\nfor i from 1 to n do\n x :=xval[i]: #current x\n y:=yval[i]: #current y\n k := f(x,y): #current slope\n yval[i+1]:= y + h*k: #new y va lue via Euler\n xval[i+1]:= x + h: #updated x-value:\n \+ od: #``od'' ends a do loop " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGR6$%\"xG%\"yG6\"6$%)operatorG%&arrowGF)9%F)F )F)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"$+\"" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%#x0G$\"\"!F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%# y0G$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#xnG$\"#5!\"\"" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++5!#6" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%%xvalG-%&arrayG6$;\"\"\"\"$,\"7\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%yvalG-%&arrayG6$;\"\"\"\"$,\"7\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%xvalG6#\"\"\"$\"\"!F)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%%yvalG6#\"\"\"$\"#5!\"\"" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 184 "approxsol2:=pointplot(\{seq([xval[i],yval[i]], i=1 ..n+1)\}):\nexactsol:=plot(exp(t),t=0..1,`color`=`red`): \n #use d t because x was already used above\ndisplay(\{approxsol2,exactsol\}) ;" }}{PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6&-%'POINTSG6 aq7$$\"+++++@!#5$\"+S>RK7!\"*7$$\"+++++AF)$\"+ferW7F,7$$\"+++++qF)$\"+ pLw1?F,7$$\"+++++pF)$\"+DW*o)>F,7$$\"+++++_F)$\"+@*)ox;F,7$$\"+++++`F) $\"+5eY%p\"F,7$$\"+++++sF)$\"+8$*4Z?F,7$$\"+++++rF)$\"+.5$o-#F,7$$\"++ +++>F)$\"+]*3\"37F,7$$\"+++++?F)$\"+S+>?7F,7$$\"\"!Ffn$\"#5!\"\"7$$\"+ ++++5!#6$\"++++55F,7$$FWF]o$\"+++5?5F,7$$\"+++++IF]o$\"++5II5F,7$$\"++ +++SF]o$\"+5SgS5F,7$$\"+++++]F]o$\"+]+,^5F,7$$\"+++++gF]o$\"+],_h5F,7$ $F4F]o$\"+_`8s5F,7$$\"+++++!)F]o$\"+1n&G3\"F,7$$\"+++++!*F]o$\"+t_o$4 \"F,7$$F\\oF)$\"+E@i/6F,7$$\"+++++6F)$\"+Z$oc6\"F,7$$\"+++++7F)$\"+I]# o7\"F,7$$\"+++++8F)$\"+!G$4Q6F,7$$\"+++++9F)$\"+8UZ\\6F,7$$\"+++++:F)$ \"+b*o4;\"F,7$$\"+++++;F)$\"+X'yD<\"F,7$$\"+++++\"F,7$$\"+++++tF)$\"+1.dn?F,7$$\"+++++uF)$\"+4gC)3#F,7 $$\"+++++vF)$\"+p%G\"4@F,7$$\"+++++wF)$\"+a(>-8#F,7$$\"+++++xF)$\"+_>_ ^@F,7$$\"+++++yF)$\"+sr.t@F,7$$\"+++++zF)$\"+Wvw%>#F,7$$F^qF)$\"+>_r;A F,7$$\"+++++\")F)$\"+rB))QAF,7$$\"+++++#)F)$\"+&>r7E#F,7$$\"+++++$)F)$ \"+2R)QG#F,7$$\"+++++%)F)$\"+YFs1BF,7$$\"+++++&)F)$\"+t**yHBF,7$$\"+++ ++')F)$\"+ty3`BF,7$$\"+++++aF)$\"+o/T6)4eu\"F,7$$\"+++++dF)$\"+F,7$$\"+++++mF)$\"+`,YG>F,7$$\"+++++nF)$ \"+bZuZ>F,7$$\"+++++oF)$\"+.AAn>F,7$$\"+++++HF)$\"+wQ]M8F,7$$FfoF)$\"+ :*[yM\"F,7$$\"+++++JF)$\"+/uKh8F,7$$\"+++++KF)$\"+y1%\\P\"F,7$$\"+++++ LF)$\"+&3!p)Q\"F,7$$\"+++++MF)$\"+')pd-9F,7$$\"+++++NF)$\"+cFg;9F,7$$ \"+++++OF)$\"+%yo2V\"F,7$$\"+++++PF)$\"+sk2X9F,7$$\"+++++QF)$\"+Ps_f9F ,7$$\"+++++RF)$\"+4D7u9F,7$$F[pF)$\"+MP'))[\"F,7$$\"+++++TF)$\"+rBv.:F ,7$$\"+++++UF)$\"+&*)*y=:F,7$$\"+++++VF)$\"+%zxR`\"F,7$$\"+++++WF)$\"+ svJ\\:F,7$$\"+++++XF)$\"+[2\"[c\"F,7$$\"+++++YF)$\"+b)e/e\"F,7$$\"++++ +ZF)$\"+WME'f\"F,7$$\"+++++[F)$\"+ygA7;F,7$$\"+++++\\F)$\"+R$[$G;F,7$$ F`pF)$\"+A=jW;F,7$$\"+++++^F)$\"+S\"y5m\"F,7$$\"+++++BF)$\"+=I;d7F,7$$ \"+++++CF)$\"+[Ytp7F,7$$\"+++++DF)$\"+%*>V#G\"F,7$$\"+++++EF)$\"+9jD&H \"F,7$$\"+++++FF)$\"+x)3#38F,7$$\"+++++GF)$\"+m4H@8F,7$$\"+++++()F)$\" +_(=mP#F,7$$\"+++++))F)$\"+S\\Q+CF,7$$\"+++++*)F)$\"+*y)QCCF,7$$FcqF)$ \"+xEj[CF,7$$\"+++++\"*F)$\"+/!>JZ#F,7$$\"+++++#*F)$\"+%>]y\\#F,7$$\"+ ++++$*F)$\"+'pGG_#F,7$$\"+++++%*F)$\"+$)p0[DF,7$$\"+++++&*F)$\"+`v`tDF ,7$$\"+++++'*F)$\"+HHF*f#F,7$$\"+++++(*F)$\"+ecEDEF,7$$\"+++++)*F)$\"+ :$=:l#F,7$$\"+++++**F)$\"+)\\L!yEF,7$$F\\oF,$\"+LQ\"[q#F,-%'CURVESG6$7 S7$Fen$\"\"\"Ffn7$$\"3dmmm;arz@!#>$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%F] [m$\"3aNdbY\\gT5F`[m7$$\"3nmmm\"z_\"4iF][m$\"3#=fT9tfS1\"F`[m7$$\"3[mm mT&phN)F][m$\"3;Hr`&G_r3\"F`[m7$$\"3BLLe*=)H\\5!#=$\"3L^LEiEj56F`[m7$$ \"3fmm\"z/3uC\"Fc\\m$\"3q]yZ%yaG8\"F`[m7$$\"3%)***\\7LRDX\"Fc\\m$\"3') oPGkJLc6F`[m7$$\"3]mm\"zR'ok;Fc\\m$\"3-Am\"\\\\E6=\"F`[m7$$\"3v***\\i5 `h(=Fc\\m$\"36)e/'[$pj?\"F`[m7$$\"3WLLL3En$4#Fc\\m$\"3]RSV5x*GB\"F`[m7 $$\"3pmm;/RE&G#Fc\\m$\"3)>3j2pYnD\"F`[m7$$\"3\")*****\\K]4]#Fc\\m$\"3S m9jYu9%G\"F`[m7$$\"3$******\\PAvr#Fc\\m$\"3%*p#e<$=E78F`[m7$$\"3)***** *\\nHi#HFc\\m$\"3PYw>&\\P*R8F`[m7$$\"3jmm\"z*ev:JFc\\m$\"3@KG01]dl8F`[ m7$$\"3?LLL347TLFc\\m$\"3?9H.#p*p'R\"F`[m7$$\"3+LLLLY.KNFc\\m$\"35Lbfo 2iB9F`[m7$$\"3v***\\7o7Tv$Fc\\m$\"3\"ouv#H**eb9F`[m7$$\"3&GLLLQ*o]RFc \\m$\"3r1/UDl[%[\"F`[m7$$\"3@++D\"=lj;%Fc\\m$\"3\\;8&[1^o^\"F`[m7$$\"3 1++vV&R'=e \"F`[m7$$\"3GLLeR\"3Gy%Fc\\m$\"3+.c]c%)H8;F`[m7$$\"3cmm;/T1&*\\Fc\\m$ \"3#**eh)zw!zk\"F`[m7$$\"3%em;zRQb@&Fc\\m$\"3G,G!GGVYo\"F`[m7$$\"3[*** \\(=>Y2aFc\\m$\"3L6+U5yG<jrpbKv\"F`[m7 $$\"3k******\\y))GeFc\\m$\"33!=v7P07z\"F`[m7$$\"3'*)***\\i_QQgFc\\m$\" 3XOlB#\\E\"H=F`[m7$$\"3@***\\7y%3TiFc\\m$\"3*)eT_<6em=F`[m7$$\"35**** \\P![hY'Fc\\m$\"3DV)H1Jn!4>F`[m7$$\"3jKLL$Qx$omFc\\m$\"3uZcmrs1[>F`[m7 $$\"3!)*****\\P+V)oFc\\m$\"39RAPIze!*>F`[m7$$\"3?mm\"zpe*zqFc\\m$\"3hL eRd*=*H?F`[m7$$\"3%)*****\\#\\'QH(Fc\\m$\"3zh'*f?z!Q2#F`[m7$$\"3GKLe9S 8&\\(Fc\\m$\"3^**33Q,(f6#F`[m7$$\"3R***\\i?=bq(Fc\\m$\"3m>.C'Qe4;#F`[m 7$$\"3\"HLL$3s?6zFc\\m$\"3C(p,G?ne?#F`[m7$$\"3a***\\7`Wl7)Fc\\m$\"3EvY ezG)QD#F`[m7$$\"3#pmmm'*RRL)Fc\\m$\"3^8U**za6,BF`[m7$$\"3Qmm;a<.Y&)Fc \\m$\"3+^n_#[T/N#F`[m7$$\"3KS+CF`[m7$$\"3t*** ***\\Qk\\*)Fc\\m$\"3t&G]\"H'[sW#F`[m7$$\"3CLL$3dg6<*Fc\\m$\"3bh$4Z;k?] #F`[m7$$\"3HmmmmxGp$*Fc\\m$\"3d[Xr/78_DF`[m7$$\"3A++D\"oK0e*Fc\\m$\"3j 4')HYrh1EF`[m7$$\"3A++v=5s#y*Fc\\m$\"3!*=xB3j&)fEF`[m7$Fhjl$\"3`X!f%G= G=FF`[m-%'COLOURG6&%$RGBG$\"*++++\"!\")FenFen-%+AXESLABELSG6$Q\"t6\"Q! 6\"-%%VIEWG6$;FenFhjl%(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 517 " In more complicated probl ems it is a very serious issue to find relatively efficient ways of ap proximating solutions. An entire field of mathematics, ``numerical an alysis'' deals with such issues for a variety of mathematical problems . The book talks about some improvements to Euler in sections 2.5 a nd 2.6, in particular it discusses improved Euler, and Runge Kutta. R unge Kutta codes are actually used in commerical numerical packages, e .g. in Maple. Let's summarize some highlights from 2.5-2.6 below." } }{PARA 0 "" 0 "" {TEXT -1 100 " For any time step h the fundament al theorem of calculus asserts that, since dy/dx = f(x,y(x))," }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "y(x+h)=y(x) + Int(f(t,y(t)),t=`x`.. `x+h`);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%\"yG6#,&%\"xG\"\"\"%\"hG F),&-F%6#F(F)-%$IntG6$-%\"fG6$%\"tG-F%6#F4/F4;F(%$x+hGF)" }}{PARA 0 " " 0 "" {TEXT -1 681 "The problem with Euler is that we always approxim ated this integral by h*f(x,y(x)), i.e. we used the left-hand endpoint as our approximation of the ``average height''. The improvements to \+ Euler depend on better approximations to that integral. ``Improved Eu ler'' uses an approximation to the Trapezoid Rule to approximate the i ntegral. The trapezoid rule for the integral approximation would be ( 1/2)*h*(f(x,y(x))+f((x+h),y(x+h)). Since we don't know y(x+h) we appr oximate it using unimproved Euler, and then feed that into the trapezo id rule. This leads to the improved Euler do loop below. Of course b efore you use it you must make sure you initialize everything correctl y." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{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$\"\"!F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>% \"yG$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\"&" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++?!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 292 "for i from 1 to n do\n k1:=f(x,y): \+ #left-hand slope\n k2:=f(x+h,y+h*k1): #approximation to righ t-hand slope\n k:= (k1+k2)/2: #approximation to average slope \n y:= y+h*k: #improved Euler update\n x:= x+h: \+ #update x\n print(x,y,exp(x));\nod: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{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 82 " Notice you almost did as well with n=5 as you did with n=100 in unimpr oved Euler. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " One can also use Taylor a pproximation methods to improve upon Euler; by differentiating the equ ation dy/dx = f(x,y) one can solve for higher order derivatives of y i n 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 m ethod, we won't do it here." }}{PARA 0 "" 0 "" {TEXT -1 641 " In t he same vein as ``improved Euler'' we can use the Simpson approximatio n for the integral instead of the Trapezoid one, and this leads to the Runge-Kutta method. (You may or may not have talked about Simpson's \+ Rule in Calculus, it is based on a quadratic approximation to the func tion f, whereas the Trapezoid rule is based on a first order approxima tion.) Here is the code for the Runge-Kutta method. The text explain s it in section 2.6. Simpson's rule approximates an integral in terms of the integrand values at each endpoint and at the interval midpoint . Runge-Kutta uses two different approximations for the midpoint valu e." }}{PARA 0 "" 0 "" {TEXT -1 89 " Before you use the loop below \+ you must initialize your values, like you did above. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "x:=x0; \+ y:=y0; n:=5; h:=(xn-x0)/n;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"xG$ \"\"!F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"yG$\"#5!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"nG\"\"&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"+++++?!#5" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 458 "for i from 1 to n do\n k1:=f(x, y): #left-hand slope\n k2:=f(x+h/2,y+h*k1/2): #1st g uess at midpoint slope\n k3:=f(x+h/2,y+h*k2/2): #second guess at \+ midpoint slope\n k4:=f(x+h,y+h*k3): #guess at right-hand slop e\n k:=(k1+2*k2+2*k3+k4)/6: #Simpson's approximation for the integ ral\n x:=x+h: #x update\n y:=y+h*k: \+ #y update\n print(x,y,exp(x)); #display current va lues\n od: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{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!\"*$\"+G4aD AF(" }}{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 va lue of e, with n=5." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 427 " As we know, solutions to non-linear DE's can blow u p, and there are other interesting pathologies as well, so if one is d oing numerical solutions there is a real need for care. The text has \+ a number of examples of this. In the solution template which accompan ies this handout you are asked to work two problems in this vein. You may also want to use the routines in this handout to do some of your \+ book homework problems." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "4 1" 48 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }