{VERSION 4 0 "IBM INTEL LINUX22" "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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" 18 256 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" 18 257 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Head ing 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 " " 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 } } {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 54 "Numerical Treatment of Or dinary Differential Equations" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 36 "Numerical Solution by E uler's Method" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }} {PARA 0 "" 0 "" {TEXT -1 51 "We use some packages: we are already fami liar with " }{TEXT 0 7 "DEtools" }{TEXT -1 5 " and " }{TEXT 0 5 "plots " }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "with(DEtools): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "" 0 " " {TEXT -1 104 "The package plottools contains more functions for plot ting, especially a function to draw a single line:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 16 "with(plottools):" }}{PARA 0 "" 0 "" {TEXT -1 16 "An d finally the " }{TEXT 0 13 "LinearAlgebra" }{TEXT -1 100 "-package - \+ you may guess, what kind of functions it contains. In this worksheet, \+ we make use of the " }{TEXT 0 6 "Vector" }{TEXT -1 5 " and " }{TEXT 0 6 "Matrix" }{TEXT -1 176 " data structures defined there. Unfortunatel y, this package is available only for newer versions of Maple - if you use an old version, you can use lists or the data structures " } {TEXT 0 6 "vector" }{TEXT -1 5 " and " }{TEXT 0 6 "matrix" }{TEXT -1 10 " from the " }{TEXT 0 6 "linalg" }{TEXT -1 17 "-package instead." } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "with(LinearAlgebra):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 26 "The population model again" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "In the previous worksheet, we solved an ordinary differential equa tion for a population model:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "pop eq := diff(p(t),t) = (0.03 - 0.0005*p(t))*p(t);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 8 "p0 := 2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "dsolve(\{popeq, p(0)=p0\},p(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "sol := rhs(%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "DEplot(\{popeq\},[p(t)],t=0..300,[[p(0)=p0]],linecolor=GREEN); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 166 "As we have a formula for the solution, we can easily compute the solution for some value of t - sa y, TEnd=100. We store this value and the corresponding DEplot-graph:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "TEnd := 100;" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 33 "pend := evalf(subs(t=TEnd, sol));" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "graph_exact := DEplot(\{pope q\},[p(t)],t=0..TEnd,[[p(0)=p0]],linecolor=GREEN):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "display(graph_exact);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 38 "Preparation for the numeric al solution" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 159 "Although we get an exact solution for this problem, we will solve it again numerically ( so we can compare the resulting approximation with the exact solution) ." }}{PARA 0 "" 0 "" {TEXT -1 57 "We will compute an approximate solut ion for the value of " }{XPPEDIT 18 0 "pend = p(TEnd);" "6#/%%pendG-% \"pG6#%%TEndG" }{TEXT -1 12 " from above." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 158 "First, for s ome integer number n, we divide the Interval [0,Tend] into n subinterv als - here, we choose all subintervals of the same length delta_t = Te nd/n: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "n := 5;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "delta_t := TEnd/n;" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "The numerical integration will com pute values " }{XPPEDIT 18 0 "pnum[k];" "6#&%%pnumG6#%\"kG" }{TEXT -1 15 ", k=1..n+1 for " }{XPPEDIT 18 0 "p(tp[k]);" "6#-%\"pG6#&%#tpG6#%\" kG" }{TEXT -1 6 " with " }{XPPEDIT 18 0 "tp[k] = (k-1)*delta_t;" "6#/& %#tpG6#%\"kG*&,&F'\"\"\"F*!\"\"F*%(delta_tGF*" }{TEXT -1 27 ". We stor e the grid points " }{XPPEDIT 18 0 "tp;" "6#%#tpG" }{TEXT -1 15 " into a vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "tp := ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "We can draw th e intervals with respect to t in the direction field:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "grid := seq(line([tp[k],0],[tp[k],pend]),k=1..n+ 1):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "display([graph_exact,grid]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 174 "Create an vector with n+1 elements p[1]...p[n+1] (t he method itself will not require this storage, but we save all interm ediate results, so that we can plot them afterwards):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "p := Vector(1..n+1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "We know the first entry of p from the initial condition: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "p[1] := p0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 14 "Euler's method" }}{PARA 0 "" 0 "" {TEXT -1 28 "To compute the other values " }{XPPEDIT 18 0 "p [k];" "6#&%\"pG6#%\"kG" }{TEXT -1 53 ", k=2..n+1, we use Euler's metho d, that approximates " }{XPPEDIT 18 0 "p[k+1];" "6#&%\"pG6#,&%\"kG\"\" \"F(F(" }{TEXT -1 4 " by " }{XPPEDIT 18 0 "p[k]+delta_t*f(t[k],p[k]); " "6#,&&%\"pG6#%\"kG\"\"\"*&%(delta_tGF(-%\"fG6$&%\"tG6#F'&F%6#F'F(F( " }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "f(t,p);" "6#-%\"fG6$%\"tG%\"p G" }{TEXT -1 53 " is the right hand side of the differential equation: " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "for k from 1 to n do\n \+ ft := subs(p(t)=p[k],rhs(popeq));\n p[k+1] := p[k]+delta_t*ft\nod;" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "points := [seq([tp[k],p[k] ],k=1..n+1)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "pointplot( points,style=LINE,thickness=2,color=BLUE);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 30 "display([%,graph_exact,grid]);" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 1 "\n" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 36 " Exercise: Accuracy of Euler's method" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "To study the influence of n to the accuracy of the solution, wr ite a procedure " }{TEXT 0 9 "eulerp(n)" }{TEXT -1 36 " that computes \+ the approximation to " }{TEXT 0 4 "pend" }{TEXT -1 53 " for different \+ values of n. Therefore, the result of " }{TEXT 0 9 "eulerp(4)" }{TEXT -1 11 " should be " }{XPPEDIT 256 0 "16.05748735;" "6#$\"+N([dg\"!\") " }{TEXT -1 16 ", the result of " }{TEXT 0 10 "eulerp(20)" }{TEXT -1 11 " should be " }{XPPEDIT 257 0 "22.44722956;" "6#$\"+cHsWA!\")" } {TEXT -1 90 ". You don't need to store the whole array p, it is suffic ient to keep the actual value of " }{XPPEDIT 18 0 "p[k];" "6#&%\"pG6#% \"kG" }{TEXT -1 31 " in a simple (scalar) variable." }{MPLTEXT 1 0 0 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "For the fol lowing computations, tell Maple to use 10 digits for floating point co mputations:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 65 "Then, compute a list of twelve approximations with n=2,4,...2^12:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "Plot the result (" }{TEXT 0 11 "pointplot ()" }{TEXT -1 35 ", x-axis: 1..12, y-axis: result of " }{TEXT 0 8 "eul erp()" }{TEXT -1 2 "):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "Next, compute for each of the 12 \+ values the (absolute) error (the difference to the exact solution):" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "... and plot it in the same way as the previous plot:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "To get a better impression of what is happening, plot it again with a logarithmic scale for the vertical axis (" }{TEXT 0 9 "logplot ()" }{TEXT -1 4 ")..." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 45 ".. and compute and plot the ratio \+ (error for " }{XPPEDIT 18 0 "n = 2^k;" "6#/%\"nG)\"\"#%\"kG" }{TEXT -1 13 ")/(error for " }{XPPEDIT 18 0 "n = 2^(k+1);" "6#/%\"nG)\"\"#,&% \"kG\"\"\"F)F)" }{TEXT -1 14 ") for k=1..11:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "Explain the l ast result with respect to the convergence rate of Euler's method!" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 139 "Now, do the same computations (except for the last two plots) \+ again for an floating point accuracy of only 4 digits and explain the \+ result!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{PARA 4 "" 0 "" {TEXT -1 0 "" }}{PARA 3 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 13 "Other Metho ds" }}{PARA 0 "" 0 "" {TEXT -1 78 "In this section we compare the resu lts obtained by the Euler method to results" }}{PARA 0 "" 0 "" {TEXT -1 83 "computed by Runge-Kutta 2nd order method, Adams-Bashforth metho d and Heun's method." }}{PARA 0 "" 0 "" {TEXT -1 82 "As an example ser ve the Lotka-Volterra equations (see lecture slides of Wednesday," }} {PARA 0 "" 0 "" {TEXT -1 13 "November 19)." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 21 "restart;\nwith(plots):" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 159 "First we define the equations. For conve nience we do so here without applying\nsymbolical notation. Next we de fine directly the right hand side of our equations" }}{PARA 0 "" 0 "" {TEXT -1 87 "as functions of p and q. This will turn out to be more su itable since it allows shorter" }}{PARA 0 "" 0 "" {TEXT -1 9 "notation ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 161 "L1 := diff(p(t),t) = \+ (-0.5+0.005*q(t))*p(t);\nL2 := diff(q(t),t) = (0.2-0.02*p(t))*q(t);\nr _L1:=(p,q)->(-0.5+0.005*q)*p;\nr_L2:=(p,q)->(0.2-0.02*p)*q;\nTEnd := 1 00;\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "For the purpose of compa rison we again implement Euler's method such as it" }}{PARA 0 "" 0 "" {TEXT -1 78 "computes us a numerical solution for the coupled differen tial equation system." }}{PARA 0 "" 0 "" {TEXT -1 91 "Define a vector \+ for p, q and t of length n. As initial values take p=6 and q=80. Updat e all" }}{PARA 0 "" 0 "" {TEXT -1 91 "variables within a single loop. \+ The vector t is not needed for the computation of the other" }}{PARA 0 "" 0 "" {TEXT -1 97 "variables, however, you could use it to plot p \+ and q vs. time. For rather high n (e.g. n=500) you" }}{PARA 0 "" 0 "" {TEXT -1 96 "will find that the solution does not stay on closed circl es, but moves away from the theoretical" }}{PARA 0 "" 0 "" {TEXT -1 12 "equilibrium." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 105 "n1 := 1000;\np v := vector(n1);\nqv := vector(n1);\ndelta_t := TEnd/n1;\npv[1] := 6; \nqv[1] := 80;\ntv[1] := 0;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "Implement now in the sam e way Heun's 2nd order scheme. For n=100 steps" }}{PARA 0 "" 0 "" {TEXT -1 94 "you will see that the approximation is much better, i.e. \+ the solution is kept more efficiently" }}{PARA 0 "" 0 "" {TEXT -1 18 " on closed circles." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 104 "n2 : = 100;\npv := vector(n2);\nqv := vector(n2);\ndelta_t := TEnd/n2;\npv[ 1] := 6;\nqv[1] := 80;\ntv[1] := 0;\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "Apply the Adams-Bashforth scheme to the problem. The firs t step should be" }}{PARA 0 "" 0 "" {TEXT -1 29 "a Runge-Kutta 2nd ord er step." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 104 "n3 := 100;\npv := vect or(n3);\nqv := vector(n3);\ndelta_t := TEnd/n3;\npv[1] := 6;\nqv[1] := 80;\ntv[1] := 0;\n" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "Finally implement the 2nd order Runge Kutta sch eme." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 104 "n4 := 100;\npv := \+ vector(n4);\nqv := vector(n4);\ndelta_t := TEnd/n4;\npv[1] := 6;\nqv[1 ] := 80;\ntv[1] := 0;\n" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{MARK "2" 0 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }