{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 }{CSTYLE "" -1 258 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 255 0 0 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 "Hea ding 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 "Heading 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 0 {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 0 {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 0 {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 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 36 "Exercise: Accuracy of Eul er's method" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "To study the influe nce of n to the accuracy of the solution, write 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 do n't need to store the whole array p, it is sufficient to keep the actu al value of " }{XPPEDIT 18 0 "p[k];" "6#&%\"pG6#%\"kG" }{TEXT -1 31 " \+ in a simple (scalar) variable." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "n \+ := 'n';" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 158 "eulerp := proc(n)\n lo cal delta_t,i,pnum;\n delta_t := TEnd/n;\n pnum := p0;\n for i to n do \n pnum := pnum + delta_t*subs(p(t)=pnum,rhs(popeq));\n od\nen d;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "eulerp(4);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "eulerp(20);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "pend;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "For the following computations, tell Maple to use 10 digi ts for floating point computations:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Digits := 10;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "Then, comput e a list of twelve approximations with n=2,4,...2^12:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "res := [seq(eulerp(2^i),i=1..12)];" }}}{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 37 "pointplot ([seq([i,res[i]],i=1..12)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 " Next, compute for each of the 12 values the (absolute) error (the diff erence to the exact solution):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "e rr := [seq(pend-res[i],i=1..12)];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "... and plot it in the same way as the previous plot:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "pointplot([seq([i,err[i]],i=1..12)]);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "To get a better impression of wha t is happening, plot it again with a logarithmic scale for the vertica l axis (" }{TEXT 0 9 "logplot()" }{TEXT -1 4 ")..." }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 47 "logplot([seq([i,err[i]],i=1..12)],style=POINT);" }} }{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 46 "pointplot([seq([i,err[i]/err[i+1]],i=1..11)]);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "Explain the last result with respe ct to the convergence rate of Euler's method!" }}{PARA 0 "" 0 "" {TEXT 258 307 "Computing err[i+1] - compared with the computation of e rr[i] - we use twice as many grid points, i.e. divide delta_t by two a nd receive an error that is (for small values of delta_t) about half t he size as err[i] (err[i]/err[i+1] converges towards 2): the error is \+ of the same order of magnitude as delta_t!" }}{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 floatin g point accuracy of only 4 digits and explain the result!" }}{PARA 0 " " 0 "" {TEXT 259 182 "Here, the error increases again, if delta_t beco mes too small: the accumulated round-off-error form more and more floa ting point operations dominates the overall error of the method!" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {PARA 3 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 13 " Other Methods" }}{PARA 0 "" 0 "" {TEXT -1 78 "In this section we compa re the results obtained by the Euler method to results" }}{PARA 0 "" 0 "" {TEXT -1 83 "computed by Runge-Kutta 2nd order method, Adams-Bash forth method and Heun's method." }}{PARA 0 "" 0 "" {TEXT -1 82 "As an \+ example serve the Lotka-Volterra equations (see lecture slides of Wedn esday," }}{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. F or convenience we do so here without applying\nsymbolical notation. Ne xt we define 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 suitable 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;\nT End := 100;\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "For the purpose \+ of comparison we again implement Euler's method such as it" }}{PARA 0 "" 0 "" {TEXT -1 78 "computes us a numerical solution for the coupled \+ differential 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=8 0. Update all" }}{PARA 0 "" 0 "" {TEXT -1 91 "variables within a singl e 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 close d circles, but moves away from the theoretical" }}{PARA 0 "" 0 "" {TEXT -1 12 "equilibrium." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 493 "n1 := 1000;\npv := vector(n1);\nqv := vector(n1);\ndelta_t := TEnd/n1;\npv[ 1] := 6;\nqv[1] := 80;\ntv[1] := 0;\nfor i from 1 to n1-1 do\n pv[i+1 ] := pv[i] + delta_t*subs(p(t)=pv[i],q(t)=qv[i],rhs(L1)):\n qv[i+ 1] := qv[i] + delta_t*subs(p(t)=pv[i],q(t)=qv[i],rhs(L2)):\n tv[i +1] := tv[i] + delta_t:\nod:\npoints_eu :=[seq([pv[i],qv[i]],i=1..n1)] :\nplot([points_eu],p=5..20,q=50..150);\np_vs_t := [seq([tv[i],pv[i]], i=1..n1)]:\nq_vs_t := [seq([tv[i],qv[i]],i=1..n1)]:\nplot([p_vs_t,q_vs _t],t=0..100);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "Implement now in the same way Heun 's 2nd order scheme. For n=100 steps" }}{PARA 0 "" 0 "" {TEXT -1 94 "y ou will see that the approximation is much better, i.e. the solution i s kept more efficiently" }}{PARA 0 "" 0 "" {TEXT -1 18 "on closed circ les." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 571 "n2 := 100;\npv := \+ vector(n2);\nqv := vector(n2);\ndelta_t := TEnd/n2;\npv[1] := 6;\nqv[1 ] := 80;\ntv[1] := 0;\nfor i from 1 to n2-1 do\n i;\n p_tpq := r_L1( pv[i],qv[i]):\n q_tpq := r_L2(pv[i],qv[i]):\n p_tphpqhq := r_L1(pv[i ]+delta_t*p_tpq,qv[i]+delta_t*q_tpq):\n q_tphpqhq := r_L2(pv[i]+delta _t*p_tpq,qv[i]+delta_t*q_tpq):\n Fp_1 := p_tpq + p_tphpqhq;\n Fq_1 : = q_tpq + q_tphpqhq;\n pv[i+1] := pv[i] + delta_t/2*Fp_1:\n qv[i+1] \+ := qv[i] + delta_t/2*Fq_1:\n tv[i+1] := tv[i] + delta_t:\nod:\npoints _he :=[seq([pv[i],qv[i]],i=1..n2)]:\nplot([points_eu,points_he],p=5..2 0,q=50..150);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "Apply the Adams- Bashforth scheme to the problem. The first step should be" }}{PARA 0 " " 0 "" {TEXT -1 29 "a Runge-Kutta 2nd order step." }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 791 "n3 := 100;\npv := vector(n3);\nqv := vector(n3);\n delta_t := TEnd/n3;\npv[1] := 6;\nqv[1] := 80;\ntv[1] := 0;\ni := 1;\n p_tpq := r_L1(pv[i],qv[i]):\nq_tpq := r_L2(pv[i],qv[i]):\np_tphpqhq := r_L1(pv[i]+delta_t*p_tpq,qv[i]+delta_t*q_tpq):\nq_tphpqhq := \+ r_L2(pv[i]+delta_t*p_tpq,qv[i]+delta_t*q_tpq):\nFp_1 := p_tpq + p_tphpqhq:\nFq_1 := q_tpq + q_tphpqhq:\npv[i+1] := pv[i] + delta_t/2* Fp_1:\nqv[i+1] := qv[i] + delta_t/2*Fq_1:\ntv[i+1] := tv[i] + delta_t: \nfor i from 2 to n3-1 do\n Fp_1 := 3*r_L1(pv[i],qv[i])-r_L1(pv[i-1], qv[i-1]):\n Fq_1 := 3*r_L2(pv[i],qv[i])-r_L2(pv[i-1],qv[i-1]):\n pv[ i+1] := pv[i] + delta_t/2*Fp_1:\n qv[i+1] := qv[i] + delta_t/2*Fq_1: \n tv[i+1] := tv[i] + delta_t:\nod:\npoints_ab :=[seq([pv[i],qv[i]],i =1..n3)]:\nplot([points_eu,points_he,points_ab],p=5..20,q=50..150);" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 " Finally implement the 2nd order Runge Kutta scheme." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 587 "n4 := 100;\npv := vector(n4);\nqv := vec tor(n4);\ndelta_t := TEnd/n4;\npv[1] := 6;\nqv[1] := 80;\ntv[1] := 0; \nfor i from 1 to n4-1 do\n i;\n p_tpq := r_L1(pv[i],qv[i]):\n q_tp q := r_L2(pv[i],qv[i]):\n p_tphpqhq := r_L1(pv[i]+delta_t/2*p_tpq ,qv[i]+delta_t/2*q_tpq):\n q_tphpqhq := r_L2(pv[i]+delta_t/2*p_tp q,qv[i]+delta_t/2*q_tpq):\n Fp_1 := p_tphpqhq;\n Fq_1 := q_tphpqhq; \n pv[i+1] := pv[i] + delta_t*Fp_1:\n qv[i+1] := qv[i] + delta_t*Fq_ 1:\n tv[i+1] := tv[i] + delta_t:\nod:\npoints_rk :=[seq([pv[i],qv[i]] ,i=1..n4)]:\nplot([points_eu,points_he,points_ab,points_rk],p=5..20,q= 50..150);" }{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 }