{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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 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 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 267 "" 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 "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 }{PSTYLE "" 0 259 1 {CSTYLE "" -1 -1 "" 0 1 0 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 260 1 {CSTYLE "" -1 -1 "" 0 1 0 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 261 1 {CSTYLE "" -1 -1 "" 0 1 0 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 262 1 {CSTYLE "" -1 -1 "" 0 1 0 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 263 1 {CSTYLE "" -1 -1 "" 0 1 0 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 "" 13 264 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 265 1 {CSTYLE "" -1 -1 "" 0 1 0 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 266 1 {CSTYLE "" -1 -1 "" 0 1 0 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 "" 13 268 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 269 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 271 1 {CSTYLE "" -1 -1 "" 0 1 0 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 272 1 {CSTYLE "" -1 -1 "" 0 1 0 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 {PARA 256 "" 0 "" {TEXT 256 11 "Math 2250-1" }}{PARA 257 "" 0 "" {TEXT 257 57 "Numerical Solutions to First Order Differential Equat ions" }}{PARA 258 "" 0 "" {TEXT -1 18 "September 15, 2008" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 116 " You may wish to download this file from our Math 2250 homework or lecture page, or by directly opening the URL " }}{PARA 263 "" 0 "" {TEXT -1 61 " http: //www.math.utah.edu/~korevaar/2250fall08/numerical1.mws" }}{PARA 0 "" 0 "" {TEXT -1 127 "from Maple. It contains discussion and Maple comma nds which will help you answer your homework problems from sections 2. 4-2.6." }}{PARA 0 "" 0 "" {TEXT -1 617 " In this handout we will s tudy numerical methods for approximating solutions to first order diff erential equations. We will soon see how higher order differential eq uations can be converted into first order systems of differential equa tions. It turns out that there is a natural way to generalize what w e do here in the context of a single first order differential equation s, to systems of first order differential equations. So understanding this project 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 75 " We wi ll be working through material from sections 2.4-2.6 of the text." }} {PARA 0 "" 0 "" {TEXT -1 232 " The most basic method of approximat ing solutions to differential equations is called Euler's method, afte r the 1700's mathematician who first formulated it. If you want to ap proximate the solution to the initial value problem " }}{PARA 265 "" 0 "" {XPPEDIT 18 0 "dy/dx = f(x,y);" "6#/*&%#dyG\"\"\"%#dxG!\"\"-%\"fG 6$%\"xG%\"yG" }{TEXT -1 0 "" }}{PARA 266 "" 0 "" {XPPEDIT 18 0 "y(x[0] ) = y[0];" "6#/-%\"yG6#&%\"xG6#\"\"!&F%6#F*" }{TEXT -1 2 " ," }}{PARA 0 "" 0 "" {TEXT -1 50 "first pick a step size ``h''. Then for x betwe en " }{XPPEDIT 18 0 "x[0];" "6#&%\"xG6#\"\"!" }{TEXT -1 5 " and " } {XPPEDIT 18 0 "x[0]+h;" "6#,&&%\"xG6#\"\"!\"\"\"%\"hGF(" }{TEXT -1 25 ", use the constant slope " }{XPPEDIT 18 0 "f(x[0],y[0]);" "6#-%\"fG6$ &%\"xG6#\"\"!&%\"yG6#F)" }{TEXT -1 14 ". At x-value " }{XPPEDIT 18 0 "x[1] = x[0]+h;" "6#/&%\"xG6#\"\"\",&&F%6#\"\"!F'%\"hGF'" }{TEXT -1 32 " your y-value will therefore be " }{XPPEDIT 18 0 "y[1] := y[0]+f(x [0],y[0])*h;" "6#>&%\"yG6#\"\"\",&&F%6#\"\"!F'*&-%\"fG6$&%\"xG6#F+&F%6 #F+F'%\"hGF'F'" }{TEXT -1 22 ". Then for x between " }{XPPEDIT 18 0 " x[1];" "6#&%\"xG6#\"\"\"" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "x[1]+h; " "6#,&&%\"xG6#\"\"\"F'%\"hGF'" }{TEXT -1 28 " you use the constant sl ope " }{XPPEDIT 18 0 "f(x[1],y[1]);" "6#-%\"fG6$&%\"xG6#\"\"\"&%\"yG6# F)" }{TEXT -1 13 ", so that at " }{XPPEDIT 18 0 "x[2] := x[1]+h;" "6#> &%\"xG6#\"\"#,&&F%6#\"\"\"F+%\"hGF+" }{TEXT -1 17 " your y-value is " }{XPPEDIT 18 0 "y[2] := y[1]+f(x[1],y[1])*h;" "6#>&%\"yG6#\"\"#,&&F%6# \"\"\"F+*&-%\"fG6$&%\"xG6#F+&F%6#F+F+%\"hGF+F+" }{TEXT -1 494 ". You \+ continue in this manner. It is easy to visualize if you understand th e slope field concept we've been talking about; you just use the slope field with finite rather than infinitesimal stepping in the x-variabl e. You use the value of the slope field at your current point to get \+ a slope which you then use to move to the next point. It is straightf orward to have the computer do this sort of tedious computation for yo u. In Euler's time such computations would have been done by hand!" } }{PARA 0 "" 0 "" {TEXT -1 1 " " }{TEXT 258 1 " " }{TEXT -1 131 " A g ood first example to illustrate Euler's method is our favorite DE from the time of Calculus, namely the initial value problem" }}{PARA 261 " " 0 "" {XPPEDIT 18 0 "dy/dx = y;" "6#/*&%#dyG\"\"\"%#dxG!\"\"%\"yG" } {TEXT -1 0 "" }}{PARA 262 "" 0 "" {XPPEDIT 18 0 "y(0) = 1;" "6#/-%\"yG 6#\"\"!\"\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 14 "We know \+ that " }{XPPEDIT 18 0 "y = exp(x);" "6#/%\"yG-%$expG6#%\"xG" }{TEXT -1 501 " is the solution. Let's take h=0.2 and try to approximate the \+ solution on the x-interval [0,1]. Since the approximate solution will be piecewise affine, we only need to know the approximations at the d iscrete x values x=0,0.2,0.4,0.6,0.8,1. I've drawn a picture on the n ext page with appoximated (x,y) points and the exact the solution grap h. Use the empty space to fill in the hand-computations which produce these y values and points. Then compare to the ``do loop'' computat ion directly below." }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 259 "" 0 "" {GLPLOT2D 298 257 257 {PLOTDATA 2 "6'-%'POINTSG6(7$$\"\"!F ($\"#5!\"\"7$$\"+++++?!#5$\"+++++7!\"*7$$\"+++++SF/$\"++++S9F27$$\"+++ ++gF/$\"++++G$\"3)y_2wWO?-\"!#<7$$\"3[LL$e 9ui2%FQ$\"3aNdbY\\gT5FT7$$\"3nmmm\"z_\"4iFQ$\"3#=fT9tfS1\"FT7$$\"3[mmm T&phN)FQ$\"3;Hr`&G_r3\"FT7$$\"3BLLe*=)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]RSV5x*GB\"FT7$$\"3pmm;/RE&G#Fao$\"3)>3j 2pYnD\"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?LLL347TLFao$\"3?9H.#p*p'R\"FT7$$\"3+L LLLY.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_QQgFao$\"3XOlB#\\E\"H=FT7$$\"3@***\\7y%3 TiFao$\"3*)eT_<6em=FT7$$\"35****\\P![hY'Fao$\"3DV)H1Jn!4>FT7$$\"3jKLL$ Qx$omFao$\"3uZcmrs1[>FT7$$\"3!)*****\\P+V)oFao$\"39RAPIze!*>FT7$$\"3?m m\"zpe*zqFao$\"3hLeRd*=*H?FT7$$\"3%)*****\\#\\'QH(Fao$\"3zh'*f?z!Q2#FT 7$$\"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$\"3EvYez G)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(-%+AXESLAB ELSG6$Q\"t6\"Q!Fa]l-%&TITLEG6#%Fapproximate|+and~exact~solution~graphs G-%%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" }}{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 261 17 "Class exercise 1:" }{TEXT -1 44 " Hand work for the app roximate solution to " }}{PARA 271 "" 0 "" {XPPEDIT 18 0 "dy/dx = y" " 6#/*&%#dyG\"\"\"%#dxG!\"\"%\"yG" }{TEXT -1 0 "" }}{PARA 272 "" 0 "" {XPPEDIT 18 0 "y(0) = 1" "6#/-%\"yG6#\"\"!\"\"\"" }{TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 67 "on the interval [0,1], with n=5 subdivisi ons, using Euler's method." }}{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 "" }}{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 262 17 "Clas s exercise 2:" }{TEXT -1 53 " Why are your approximations too small i n this case?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 34 "Here is the automated \+ computation:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "restart: #c lear any memory from earlier work " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "Digits:=6: #number of significant digits you desire \n #default is 10." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 150 "x 0:=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$\"'++?!\"'" }}} {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#>%\"fGf*6$%\"xG%\"yG6\"6$%) operatorG%&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 338 "for i fr om 1 to n do\n k:= f(x,y): #current slope,use: to suppress o utput\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 solution. \n od: \n #``od'' ends a do loop " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++?!\"'$\"'++7!\"&$\"'S@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++S!\"'$\"'+S9!\"&$\"'#=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++g!\"'$\"'+G " 0 "" {MPLTEXT 1 0 28 "resta rt: #clear all memory" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 " Digits:=6:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "with(plots): \+ #load plotting library\nwith(linalg): #linear algebra library" }}} {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 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "n:=5:h:=1/n:x0:=0:y0: =1: #initialize" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "xval:=ve ctor(n+1);yval:=vector(n+1); \n #to collect all our points" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "xval[1]:=x0; yval[1]:=y0; \+ \n #initial values" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 362 " #paste in the previous work, and modify to store\n #all values in an array:\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-va lue:\n end do: #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 76 "display(\{approxsol,exactsol\},title=`approximate\nand exact solut ion graphs`);" }}{PARA 13 "" 1 "" {TEXT -1 0 "" }}}{PARA 264 "" 1 "" {TEXT -1 1 " " }}{PAGEBK }{PARA 268 "" 0 "" {TEXT -1 560 " 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 th e computer do too many computations because of problems with round-off error and computation time, so for example, choosing h=0.0000001 woul d not be practical. But, trying h=0.01 in our previous initial value p roblem should be instructive.\n If we change the n-value to 100 an d keep the other data the same we can rerun our experiment:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "x0:=0.0; xn:=1.0; y0:=1.0; n:=100; \+ h:=(xn-x0)/n;\nx:=x0; y:=y0;" }{TEXT -1 0 "" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 612 "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 \+ end if; #use the ``if'' test to decide when to print;\n \+ #the command ``frac'' computes the remainder \n #o f a quotient, it will be zero for us if i \n #is a multi ple of 10. This way we only print\n #the approximations \+ every 0.1 increments, even\n #though our actual time ste p is 0.01.\n end do: #end the do loop " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++5!\"'$\"'j/6!\"&$\"'<06F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"'++?!\"'$\"'??7!\"&$\"'S@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++I!\"'$\"'%yM\"!\"&$\"'')\\8F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"'++S!\"'$\"'())[\"!\"&$\"'#=\\\"F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++]!\"'$\"'jW;!\"&$\"'s[;F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"'++g!\"'$\"'q;=!\"&$\"'7A=F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++q!\"'$\"'w1?!\"&$\"'v8?F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++!)!\"'$\"'s;A!\"&$\"'aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++!*!\"'$\"'k[C!\"&$\"'gfCF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++5!\"&$\"'$[q#F%$\"'G=FF%" }}}{PARA 0 "" 0 "" {TEXT -1 134 " So you can see we got closer to the actual value of e, \+ but really, considering how much work Maple did this was not a great r esult. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 267 18 "Class exercise 3: " }{TEXT -1 156 " For this very special initial \+ value problem which has exp(x) as the solution, set up Euler on the x- interval [0,1], with \"n\" subdivisions, and step size " }{XPPEDIT 18 0 "h = 1/n;" "6#/%\"hG*&\"\"\"F&%\"nG!\"\"" }{TEXT -1 95 ". Write \+ down the resulting Euler estimate for exp(1). What is the limit of th is estimate as " }{XPPEDIT 18 0 "proc (n) options operator, arrow; in finity end proc;" "6#f*6#%\"nG7\"6$%)operatorG%&arrowG6\"%)infinityGF* F*F*" }{TEXT -1 46 "? You learned this special limit in Calculus!" }} {PAGEBK }{PARA 0 "" 0 "" {TEXT -1 139 "We can make a picture of what w e did as follows, using the mouse to cut and paste previous work, and \+ then editing it for the new situation:" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Digits:=6:with(plots):with(linalg): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 486 "f:=(x,y)->y;\nn:=100; \+ x0:=0.0; y0:=1.0;\nxn:=1.0; h:=(xn-x0)/n;\nxval:=vector(n+1);yval:=vec tor(n+1); \n #to collect all our points. Now n=100\nxval[1]:=x0; yv al[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 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 184 "approxsol2:=pointplot(\{seq([xval[ i],yval[i]], i=1..n+1)\}):\nexactsol:=plot(exp(t),t=0..1,`color`=`red` ): \n #used t because x was already used above\ndisplay(\{approx sol2,exactsol\});" }}{PARA 13 "" 1 "" {GLPLOT2D 304 260 260 {PLOTDATA 2 "6&-%'POINTSG6aq7$$\"+++++U!#5$\"+&*)*y=:!\"*7$$\"+++++VF)$\"+%zxR` \"F,7$$\"+++++SF)$\"+MP'))[\"F,7$$\"+++++TF)$\"+rBv.:F,7$$\"+++++WF)$ \"+svJ\\:F,7$$\"+++++XF)$\"+[2\"[c\"F,7$$\"+++++!*F)$\"+xEj[CF,7$$\"++ +++\"*F)$\"+/!>JZ#F,7$$\"+++++)*F)$\"+:$=:l#F,7$$\"+++++**F)$\"+)\\L!y EF,7$$\"+++++%*F)$\"+$)p0[DF,7$$\"+++++&*F)$\"+`v`tDF,7$$\"+++++#*F)$ \"+%>]y\\#F,7$$\"+++++$*F)$\"+'pGG_#F,7$$\"+++++'*F)$\"+HHF*f#F,7$$\"+ ++++(*F)$\"+ecEDEF,7$$\"\"!Fdp$\"#5!\"\"7$$\"+++++5!#6$\"++++55F,7$$\" +++++?F[q$\"+++5?5F,7$$\"+++++IF[q$\"++5II5F,7$$F4F[q$\"+5SgS5F,7$$\"+ ++++]F[q$\"+]+,^5F,7$$\"+++++gF[q$\"+],_h5F,7$$\"+++++qF[q$\"+_`8s5F,7 $$\"+++++!)F[q$\"+1n&G3\"F,7$$FHF[q$\"+t_o$4\"F,7$$FjpF)$\"+E@i/6F,7$$ \"+++++6F)$\"+Z$oc6\"F,7$$\"+++++7F)$\"+I]#o7\"F,7$$\"+++++8F)$\"+!G$4 Q6F,7$$\"+++++9F)$\"+8UZ\\6F,7$$\"+++++:F)$\"+b*o4;\"F,7$$\"+++++;F)$ \"+X'yD<\"F,7$$\"+++++\"F,7$$\"+ ++++>F)$\"+]*3\"37F,7$$F`qF)$\"+S+>?7F,7$$\"+++++@F)$\"+S>RK7F,7$$\"++ +++AF)$\"+ferW7F,7$$\"+++++BF)$\"+=I;d7F,7$$\"+++++CF)$\"+[Ytp7F,7$$\" +++++DF)$\"+%*>V#G\"F,7$$\"+++++EF)$\"+9jD&H\"F,7$$\"+++++FF)$\"+x)3#3 8F,7$$\"+++++GF)$\"+m4H@8F,7$$\"+++++HF)$\"+wQ]M8F,7$$FeqF)$\"+:*[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$$\"+++++O F)$\"+%yo2V\"F,7$$\"+++++PF)$\"+sk2X9F,7$$\"+++++QF)$\"+Ps_f9F,7$$\"++ +++RF)$\"+4D7u9F,7$$\"+++++YF)$\"+b)e/e\"F,7$$\"+++++ZF)$\"+WME'f\"F,7 $$\"+++++[F)$\"+ygA7;F,7$$\"+++++\\F)$\"+R$[$G;F,7$$F^rF)$\"+A=jW;F,7$ $\"+++++^F)$\"+S\"y5m\"F,7$$\"+++++_F)$\"+@*)ox;F,7$$\"+++++`F)$\"+5eY %p\"F,7$$\"+++++aF)$\"+o/T6)4eu\"F,7$$\"+++++dF)$\"+F,7$$\"+++++mF)$\"+`,YG>F,7$$\"+++++nF)$\"+bZuZ>F,7$ $\"+++++oF)$\"+.AAn>F,7$$\"+++++pF)$\"+DW*o)>F,7$$FhrF)$\"+pLw1?F,7$$ \"+++++rF)$\"+.5$o-#F,7$$\"+++++sF)$\"+8$*4Z?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]sF)$\"+>_r;AF,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$$\"+++++()F)$\"+_(=mP#F ,7$$\"+++++))F)$\"+S\\Q+CF,7$$\"+++++*)F)$\"+*y)QCCF,7$$FjpF,$\"+LQ\"[ q#F,-%'CURVESG6$7S7$Fcp$\"\"\"Fdp7$$\"3dmmm;arz@!#>$\"3)y_2wWO?-\"!#<7 $$\"3[LL$e9ui2%F][m$\"3aNdbY\\gT5F`[m7$$\"3nmmm\"z_\"4iF][m$\"3#=fT9tf S1\"F`[m7$$\"3[mmmT&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%)***\\7 LRDX\"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$\"3Sm9jYu9%G\"F`[m7$$\"3$******\\PAvr#Fc\\m$\"3%*p#e<$=E7 8F`[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$\"35Lbfo2iB9F`[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$$\"31++vV&R'=e\"F`[m7$$\"3GLLeR\"3Gy%Fc\\m$\"3+.c]c%)H8;F`[m7$$\"3cm m;/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\"zp e*zqFc\\m$\"3hLeRd*=*H?F`[m7$$\"3%)*****\\#\\'QH(Fc\\m$\"3zh'*f?z!Q2#F `[m7$$\"3GKLe9S8&\\(Fc\\m$\"3^**33Q,(f6#F`[m7$$\"3R***\\i?=bq(Fc\\m$\" 3m>.C'Qe4;#F`[m7$$\"3\"HLL$3s?6zFc\\m$\"3C(p,G?ne?#F`[m7$$\"3a***\\7`W l7)Fc\\m$\"3EvYezG)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$\"3j4')HYrh1EF`[m7$$\"3A++v=5s#y*Fc\\m$\"3!*=xB3j&)fEF`[m 7$Fhjl$\"3`X!f%G=G=FF`[m-%'COLOURG6&%$RGBG$\"*++++\"!\")FcpFcp-%+AXESL ABELSG6$Q\"t6\"Q!Ffjm-%%VIEWG6$;FcpFhjl%(DEFAULTG" 1 2 0 1 10 0 2 9 1 4 2 1.000000 45.000000 44.000000 0 0 "Curve 1" "Curve 2" }}}}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 470 " In \+ more complicated problems it is a very serious issue to find relativel y efficient ways of approximating solutions. An entire field of mathe matics, ``numerical analysis'' deals with such issues for a variety of mathematical problems. The book talks about some improvements to Eu ler in sections 2.5 and 2.6, in particular it discusses improved Eule r, and Runge Kutta. Runge Kutta-type codes are actually used in comme rical numerical packages, e.g. in Maple. " }}{PARA 0 "" 0 "" {TEXT -1 50 " Let's summarize some highlights from 2.5-2.6." }}{PARA 0 " " 0 "" {TEXT -1 99 " For any time step h the fundamental theorem o f calculus asserts that, since dy/dx = f(x,y(x))," }}{PARA 260 "" 0 " " {XPPEDIT 18 0 "y(x+h) = y(x)+Int(f(t,y(t)),t = x .. `x+h`);" "6#/-% \"yG6#,&%\"xG\"\"\"%\"hGF),&-F%6#F(F)-%$IntG6$-%\"fG6$%\"tG-F%6#F4/F4; F(%$x+hGF)" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 818 "The proble m 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 improvements to Euler depend on better approxim ations to that integral. These are subtle, because we don't yet have \+ an approximation for y(t) when t is greater than x, so also not for th e integrand. ``Improved Euler'' uses an approximation to the Trapezoi d Rule to approximate the integral. Recall, the trapezoid rule for th is 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 approximate it using unimproved Euler, and then feed that into the trapezoid rule. This leads to the improv ed Euler do loop below. Of course before you use it you must make sur e you initialize everything correctly." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "x:=x0; y:=y0; n:=5; h:= (xn-x0)/n;" }}}{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 right-hand slope\n k:= (k1+k2)/2: # approximation to average slope\n y:= y+h*k: #improved Eu ler 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%$\"'++?!\"'$\"'+?7!\"&$\"'S@7F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++S!\"'$\"'S)[\"!\"&$\"'#=\\\"F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"'++g!\"'$\"'&e\"=!\"&$\"'7A=F(" }}{PARA 11 "" 1 " " {XPPMATH 20 "6%$\"'++!)!\"'$\"'M:A!\"&$\"'aDAF(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++5!\"&$\"'r-FF%$\"'G=FF%" }}}{PARA 0 "" 0 "" {TEXT -1 100 "Notice you almost did as well with n=5 in improved Euler as you did with n=100 in unimproved Euler. " }}{PARA 0 "" 0 "" {TEXT -1 337 " One can also use Taylor approximation methods to improve \+ upon Euler; by differentiating the equation dy/dx = f(x,y) one can sol ve 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). S ee the book for more details of this method, we won't do it here." }} {PAGEBK }{PARA 0 "" 0 "" {TEXT -1 728 " In the same vein as ``impr oved Euler'' we can use the Simpson approximation for the integral ins tead of the Trapezoid one, and this leads to the Runge-Kutta method. \+ (You may or may not have talked about Simpson's Parabolic Rule for app roximating definite integrals in Calculus, it is based on a quadratic \+ approximation to the function f, whereas the Trapezoid rule is based o n 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 y(x+h/2) to plug into f(x+h/2,y(x+h/2))." } }{PARA 0 "" 0 "" {TEXT -1 90 " Before you use the loop on the next page you must initialize your values, as before. " }}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 33 "x:=x0; y:=y0; n:=5; h:=(xn-x0)/n;" }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}}{PARA 269 "" 1 "" {TEXT -1 33 "We're still using slope function " }{XPPEDIT 18 0 "f(x,y) = y;" "6#/-%\"fG6$%\"xG %\"yGF(" }{TEXT -1 75 " but the code below will work for whatever you \+ define this function to be! " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 458 "for i from 1 to n do\n k1:=f(x,y): #left-hand slo pe\n k2:=f(x+h/2,y+h*k1/2): #1st guess 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 slope\n k:=(k1+2*k2+2*k3+k4)/6: \+ #Simpson's approximation for the integral\n x:=x+h: \+ #x update\n y:=y+h*k: #y update\n print(x,y,ex p(x)); #display current values\n od: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++?!\"'$ \"'S@7!\"&F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++S!\"'$\"'#=\\\"! \"&F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++g!\"'$\"'6A=!\"&$\"'7A= F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++!)!\"'$\"'_DA!\"&$\"'aDAF( " }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"'++5!\"&$\"'D=FF%$\"'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 303 " As we know, solutions to non-linear D E's can blow up, and there are other interesting pathologies as well, \+ so if one is doing numerical solutions there is a real need for care a nd understanding. One such example is in your homework. Another exce llent example to look at is Example 4, pages 138-139." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "6 0" 82 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }