{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 "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 -1 11 "Math 2250-3" }}{PARA 257 "" 0 " " {TEXT -1 22 "Numerical Computations" }}{PARA 258 "" 0 "" {TEXT -1 24 "Wednesday Sept 10, 2003" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 843 " In this handout we will study numerical m ethods for approximating solutions to first order differential equatio ns. (The fact is, most differential equations do NOT have simple form ulas for their solutions, despite all the examples you've seen in whic h they do. In the case that no nice formula exists for the solution o ne must approximate it numerically.) A month from now we will see how higher order differential equations can be converted into first 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 first ord er differential equation to this more complicated setting. So underst anding today's material will be an important step in understanding num erical solutions to higher order differential equations and to systems of differential equations." }}{PARA 0 "" 0 "" {TEXT -1 86 " We wi ll be working through selected material from sections 2.4-2.6 of the t ext. " }}{PARA 0 "" 0 "" {TEXT -1 960 " The most basic method of \+ approximating solutions to differential equations is called Euler's me thod, 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 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 between x1 a nd 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 been tal king about; you just use the slope field where you are at the end of e ach step to get a slope for the next step. It is straightforward to h ave a programmable calculator or computer software do this sort of ted ious computation for you. In Euler's time such computations would hav e been done by hand!" }}{PARA 0 "" 0 "" {TEXT -1 467 " A good firs t example to illustrate Euler's method is our favorite DE from the tim e 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 th e solution of the x-interval [0,1]. Since the approximate solution wi ll 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 computations. " }}{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 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 100 "f:=(x,y)->y; #this is the slope function f(x,y) \n #in d y/dx = f(x,y), in our example dy/dx = y." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "x:=x0; y:=y0; #initialize x,y for the do loop" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 322 "for i from 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-va lue:\n print(x,y,exp(x)); #display current values, \n \+ #and compare to exact solution.\n od: #` `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 64 "xval:=vector(n+1);yval:=vector(n+1) ; #to collect all our points" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "xval[1]:=x0; yval[1]:=y0; #initial values" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 355 " #paste in the previous \+ work, and modify for plotting:\nfor i from 1 to n do\n x:=xva l[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 v ia 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 because x was a lready used above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "displa y(\{approxsol,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 349 262 262 {PLOTDATA 2 "6&-%'POINTSG6(7$$\"\"!F($\"#5!\"\"7$$\"+++++?!#5$\"+++++7 !\"*7$$\"+++++SF/$\"++++S9F27$$\"+++++gF/$\"++++G$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%FQ$\"3aNdbY\\gT5FT7$$\"3nmmm \"z_\"4iFQ$\"3#=fT9tfS1\"FT7$$\"3[mmmT&phN)FQ$\"3;Hr`&G_r3\"FT7$$\"3BL Le*=)H\\5!#=$\"3L^LEiEj56FT7$$\"3fmm\"z/3uC\"Fao$\"3q]yZ%yaG8\"FT7$$\" 3%)***\\7LRDX\"Fao$\"3')oPGkJLc6FT7$$\"3]mm\"zR'ok;Fao$\"3-Am\"\\\\E6= \"FT7$$\"3v***\\i5`h(=Fao$\"36)e/'[$pj?\"FT7$$\"3WLLL3En$4#Fao$\"3]RSV 5x*GB\"FT7$$\"3pmm;/RE&G#Fao$\"3)>3j2pYnD\"FT7$$\"3\")*****\\K]4]#Fao$ \"3Sm9jYu9%G\"FT7$$\"3$******\\PAvr#Fao$\"3%*p#e<$=E78FT7$$\"3)****** \\nHi#HFao$\"3PYw>&\\P*R8FT7$$\"3jmm\"z*ev:JFao$\"3@KG01]dl8FT7$$\"3?L LL347TLFao$\"3?9H.#p*p'R\"FT7$$\"3+LLLLY.KNFao$\"35Lbfo2iB9FT7$$\"3v** *\\7o7Tv$Fao$\"3\"ouv#H**eb9FT7$$\"3&GLLLQ*o]RFao$\"3r1/UDl[%[\"FT7$$ \"3@++D\"=lj;%Fao$\"3\\;8&[1^o^\"FT7$$\"31++vV&R'=e\"FT7$$\"3GLLeR\"3Gy%Fao$\"3+.c]c%) H8;FT7$$\"3cmm;/T1&*\\Fao$\"3#**eh)zw!zk\"FT7$$\"3%em;zRQb@&Fao$\"3G,G !GGVYo\"FT7$$\"3[***\\(=>Y2aFao$\"3L6+U5yG<jrpbKv\"FT7$$\"3k******\\y))GeFao$\"33!=v7P07z\"FT7$$\"3'*)***\\i_QQg Fao$\"3XOlB#\\E\"H=FT7$$\"3@***\\7y%3TiFao$\"3*)eT_<6em=FT7$$\"35**** \\P![hY'Fao$\"3DV)H1Jn!4>FT7$$\"3jKLL$Qx$omFao$\"3uZcmrs1[>FT7$$\"3!)* ****\\P+V)oFao$\"39RAPIze!*>FT7$$\"3?mm\"zpe*zqFao$\"3hLeRd*=*H?FT7$$ \"3%)*****\\#\\'QH(Fao$\"3zh'*f?z!Q2#FT7$$\"3GKLe9S8&\\(Fao$\"3^**33Q, (f6#FT7$$\"3R***\\i?=bq(Fao$\"3m>.C'Qe4;#FT7$$\"3\"HLL$3s?6zFao$\"3C(p ,G?ne?#FT7$$\"3a***\\7`Wl7)Fao$\"3EvYezG)QD#FT7$$\"3#pmmm'*RRL)Fao$\"3 ^8U**za6,BFT7$$\"3Qmm;a<.Y&)Fao$\"3+^n_#[T/N#FT7$$\"3KS+CFT7$$\"3t******\\Qk\\*)Fao$\"3t&G]\"H'[sW#FT7$$\"3CLL$3dg6<* Fao$\"3bh$4Z;k?]#FT7$$\"3HmmmmxGp$*Fao$\"3d[Xr/78_DFT7$$\"3A++D\"oK0e* Fao$\"3j4')HYrh1EFT7$$\"3A++v=5s#y*Fao$\"3!*=xB3j&)fEFT7$FL$\"3`X!f%G= G=FFT-%'COLOURG6&%$RGBGF(F(F(-%+AXESLABELSG6$Q\"t6\"Q!Fa]l-%%VIEWG6$;F 'FL%(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 443 "If you connect the Euler dots in your mind, the picture above is like the one in fig ure 2.4.3, on page 109 of the text. The upper graph is of the exponen tial function, the lower graph is of your Euler approximation. The re ason that the dots lie below the true graph is that as y increases the slope f(x,y)=y should also be increasing, but in the Euler approximat ion you use the slope at each (lower ) point to get to the next (highe r) point." }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 424 " It should be t hat as your step size ``h'' gets smaller, your approximations get bett er 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 com puter do too many computations because of problems with round-off erro r and computation time, so for example, choosing h=0.0000001 would not be practical. But, trying h=0.01 should be instructive:" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 248 " Since the wi dth of our x-interval is 1, we stepsize h=0.01 by taking an n-value of subintervals to be100. We keep the other data the same as in our fir st example. The following code only prints out approximations when h \+ is a multiple of 0.1:" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "x0:=0.0; xn:=1.0; y0:=1.0; n:=100; \+ h:=(xn-x0)/n;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "f:=(x,y)-> y;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "x:=x0; y:=y0;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 438 "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 print(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 wil l be zero for us if i \n #is a multiple of 10. \n od: " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++5!#5$\"+E@i/6!\"*$\"+=4<0 6F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"+++++?!#5$\"+S+>?7!\"*$\"+eF S@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?!\"*$\"+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%" }}}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 33 "We can also see this grap hically:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "xval:=vector(n+1 );yval:=vector(n+1); #to collect all \n #our points" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "xval[1]:=x0; yval[1]:=y0; #i nitial values" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 355 " \+ #paste in the previous work, and modify for plotting:\nfor i from 1 \+ to n do\n x:=xval[i]: #current x\n y:=yval[i]: #cu rrent 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: #u pdated x-value:\n od: #``od'' ends a do lo op " }}}{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 91 "exactsol:=plot(exp(t),t=0..1,`color`=`black`): \n # used t because x was already used above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "display(\{approxsol,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 413 322 322 {PLOTDATA 2 "6&-%'CURVESG6$7S7$$\"\"!F)$\"\"\"F) 7$$\"3dmmm;arz@!#>$\"3)y_2wWO?-\"!#<7$$\"3[LL$e9ui2%F/$\"3aNdbY\\gT5F2 7$$\"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]dl8F27$$\"3?LL L347TLFE$\"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\"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>F27$$\"3jKLL$Qx$omFE$\"3uZcmrs1[>F27$$\"3!)*****\\P+V)oFE$\"39R APIze!*>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$$\"3Hmmmmx Gp$*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)-%'PO INTSG6aq7$$\"+++++H!#5$\"+wQ]M8!\"*7$$\"+++++GFa[l$\"+m4H@8Fd[l7$F($\" #5!\"\"7$$\"+++++5!#6$\"++++55Fd[l7$$\"+++++?Fa\\l$\"+++5?5Fd[l7$$\"++ +++IFa\\l$\"++5II5Fd[l7$$\"+++++SFa\\l$\"+5SgS5Fd[l7$$\"+++++]Fa\\l$\" +]+,^5Fd[l7$$\"+++++gFa\\l$\"+],_h5Fd[l7$$\"+++++qFa\\l$\"+_`8s5Fd[l7$ $\"+++++!)Fa\\l$\"+1n&G3\"Fd[l7$$\"+++++!*Fa\\l$\"+t_o$4\"Fd[l7$$F`\\l Fa[l$\"+E@i/6Fd[l7$$\"+++++6Fa[l$\"+Z$oc6\"Fd[l7$$\"+++++7Fa[l$\"+I]#o 7\"Fd[l7$$\"+++++8Fa[l$\"+!G$4Q6Fd[l7$$\"+++++9Fa[l$\"+8UZ\\6Fd[l7$$\" +++++:Fa[l$\"+b*o4;\"Fd[l7$$\"+++++;Fa[l$\"+X'yD<\"Fd[l7$$\"+++++\"Fd[l7$$\"+++++>Fa[l$\"+]*3\"3 7Fd[l7$$Ff\\lFa[l$\"+S+>?7Fd[l7$$\"+++++@Fa[l$\"+S>RK7Fd[l7$$\"+++++AF a[l$\"+ferW7Fd[l7$$\"+++++BFa[l$\"+=I;d7Fd[l7$$\"+++++CFa[l$\"+[Ytp7Fd [l7$$\"+++++DFa[l$\"+%*>V#G\"Fd[l7$$\"+++++EFa[l$\"+9jD&H\"Fd[l7$$\"++ +++FFa[l$\"+x)3#38Fd[l7$$F[]lFa[l$\"+:*[yM\"Fd[l7$$\"+++++JFa[l$\"+/uK h8Fd[l7$$\"+++++KFa[l$\"+y1%\\P\"Fd[l7$$\"+++++LFa[l$\"+&3!p)Q\"Fd[l7$ $\"+++++MFa[l$\"+')pd-9Fd[l7$$\"+++++NFa[l$\"+cFg;9Fd[l7$$\"+++++OFa[l $\"+%yo2V\"Fd[l7$$\"+++++PFa[l$\"+sk2X9Fd[l7$$\"+++++QFa[l$\"+Ps_f9Fd[ l7$$\"+++++RFa[l$\"+4D7u9Fd[l7$$F`]lFa[l$\"+MP'))[\"Fd[l7$$\"+++++TFa[ l$\"+rBv.:Fd[l7$$\"+++++UFa[l$\"+&*)*y=:Fd[l7$$\"+++++VFa[l$\"+%zxR`\" Fd[l7$$\"+++++WFa[l$\"+svJ\\:Fd[l7$$\"+++++XFa[l$\"+[2\"[c\"Fd[l7$$\"+ ++++YFa[l$\"+b)e/e\"Fd[l7$$\"+++++ZFa[l$\"+WME'f\"Fd[l7$$\"+++++[Fa[l$ \"+ygA7;Fd[l7$$\"+++++\\Fa[l$\"+R$[$G;Fd[l7$$Fe]lFa[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$$\"+++++dFa[l$\"+Fd[l7$$\"+++++mFa[l$\"+`,YG>Fd[l7$$\"+++++nFa[l$\"+bZuZ>Fd[l7$$ \"+++++oFa[l$\"+.AAn>Fd[l7$$\"+++++pFa[l$\"+DW*o)>Fd[l7$$F_^lFa[l$\"+p Lw1?Fd[l7$$\"+++++rFa[l$\"+.5$o-#Fd[l7$$\"+++++sFa[l$\"+8$*4Z?Fd[l7$$ \"+++++tFa[l$\"+1.dn?Fd[l7$$\"+++++uFa[l$\"+4gC)3#Fd[l7$$\"+++++vFa[l$ \"+p%G\"4@Fd[l7$$\"+++++wFa[l$\"+a(>-8#Fd[l7$$\"+++++xFa[l$\"+_>_^@Fd[ l7$$\"+++++yFa[l$\"+sr.t@Fd[l7$$\"+++++zFa[l$\"+Wvw%>#Fd[l7$$Fd^lFa[l$ \"+>_r;AFd[l7$$\"+++++\")Fa[l$\"+rB))QAFd[l7$$\"+++++#)Fa[l$\"+&>r7E#F d[l7$$\"+++++$)Fa[l$\"+2R)QG#Fd[l7$$\"+++++%)Fa[l$\"+YFs1BFd[l7$$\"+++ ++&)Fa[l$\"+t**yHBFd[l7$$\"+++++')Fa[l$\"+ty3`BFd[l7$$\"+++++()Fa[l$\" +_(=mP#Fd[l7$$\"+++++))Fa[l$\"+S\\Q+CFd[l7$$\"+++++*)Fa[l$\"+*y)QCCFd[ l7$$Fi^lFa[l$\"+xEj[CFd[l7$$\"+++++\"*Fa[l$\"+/!>JZ#Fd[l7$$\"+++++#*Fa [l$\"+%>]y\\#Fd[l7$$\"+++++$*Fa[l$\"+'pGG_#Fd[l7$$\"+++++%*Fa[l$\"+$)p 0[DFd[l7$$\"+++++&*Fa[l$\"+`v`tDFd[l7$$\"+++++'*Fa[l$\"+HHF*f#Fd[l7$$ \"+++++(*Fa[l$\"+ecEDEFd[l7$$\"+++++)*Fa[l$\"+:$=:l#Fd[l7$$\"+++++**Fa [l$\"+)\\L!yEFd[l7$$F`\\lFd[l$\"+LQ\"[q#Fd[l-%+AXESLABELSG6$Q\"t6\"Q!F cjm-%%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" }}}}{PAGEBK }{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 11 "" 1 "" {XPPMATH 20 "6#/-%\"yG6#,&%\"xG\"\"\"%\"hGF),&-F%6#F(F)-%$IntG6$* &%#dyGF)%#dtG!\"\"/%\"tG;F(F'F)" }}}{EXCHG {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 Eul er is that we always approximated this integral by h*f(x,y(x)), i.e. w e used the left-hand endpoint as our approximation of the ``average he ight''. The improvements to Euler depend on better approximations to \+ that integral. ``Improved Euler'' uses an approximation to the Trapez oid Rule to approximate the integral. The trapezoid rule for the inte gral 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 approximate it using unimproved Euler, and th en feed that into the trapezoid rule. This leads to the improved Eule r do loop below. Of course before you use it you must make sure you i nitialize everything correctly. We'll compare it when n=5, to our firs t (unimproved) attempt." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}{EXCHG {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;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 289 "f or i from 1 to n do\n k1:=f(x,y): #left-hand slope\n k2 :=f(x+h,y+h*k1): #approximation to right-hand slope\n k:= (k1+k2) /2: #approximation to average slope\n y:= y+h*k: #i mproved Euler update\n x:= x+h: #update x\n print(x, y,exp(x));\nod:" }}{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 with improved Euler when n=5 as we did w ith n=100 in unimproved Euler. " }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " One can also use Taylor appr oximation methods to improve upon Euler; by differentiating the equati on dy/dx = f(x,y) one can solve for higher order derivatives of y in t erms of the lower order ones, and then use the Taylor approximation fo r y(x+h) in terms of y(x). See the book for more details of this meth od, 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 f or the integral instead of the Trapezoid one, and this leads to the Ru nge-Kutta method which is widely used in real software." }}{PARA 0 "" 0 "" {TEXT -1 470 " (You may or may not have talked about Simpson's Ru le in Calculus, it is based on a quadratic approximation to the functi on f, whereas the Trapezoid rule is based on a first order approximati on.) Here is the code for the Runge-Kutta method. The text explains \+ 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 value. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "x:=x0; y:=y0; n:=5; h:=(xn-x0)/n;\nf:=(x,y)->y;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 456 "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 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%" }}}{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 w ith code like Runge-Kutta there are subtleties and problems which part icular problems will cause. We will not go into those here; there are good examples in sections 2.4-2.6 and the homework problems." }}} {MARK "51 3 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }