{VERSION 6 0 "IBM INTEL NT" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 2 1 2 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 2 0 2 0 2 2 0 1 }{PSTYLE "_psty le4" -1 203 1 {CSTYLE "" -1 -1 "Courier" 1 12 255 0 0 1 2 1 2 2 1 2 1 1 1 1 }1 1 0 0 0 0 2 0 2 0 2 2 0 1 }{PSTYLE "_pstyle5" -1 204 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 2 0 2 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 8 "restart;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 204 "" 0 " " {TEXT -1 2 " " }}{PARA 204 "" 0 "" {TEXT -1 70 " This work sheet l ooks at numerical solutions to an initial value ODE" }}{PARA 204 "" 0 "" {TEXT -1 73 " using the modified Picard method through integration . It will generate " }}{PARA 0 "" 0 "" {TEXT -1 63 " a Maclaurin seri es for your solution or a numerical solution." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 28 " We set up the example f or " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 53 " \+ x'' + sin(x) = 0 ; x(0) = pi/4 x'(0) = 0" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 " Let v = x', \+ y = sin(x) and w = cos(x) we have" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 34 " x' = v ; x(0) = pi/4" }} {PARA 0 "" 0 "" {TEXT -1 32 " v' = -y ; v(0) = 0;" }} {PARA 0 "" 0 "" {TEXT -1 40 " y' = wv ; y(0) = sin(pi/4); " }}{PARA 0 "" 0 "" {TEXT -1 18 " w' = " }{TEXT -1 23 "-yv ; w(0) = cos(pi/4);" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 16 "with(numapprox):" }}{PARA 203 "> " 0 " " {MPLTEXT 1 -1 12 "with(plots):" }}}{EXCHG {PARA 204 "" 0 "" {TEXT -1 1 " " }}{PARA 204 "" 0 "" {TEXT -1 66 " m is the number of eqns, n \+ is the order of the Maclaurin solution" }}{PARA 204 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 7 "m := 4;" }}} {EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 7 "n := 8;" }}}{EXCHG {PARA 204 "" 0 "" {TEXT -1 1 " " }}{PARA 204 "" 0 "" {TEXT -1 21 " h is th e time step" }}{PARA 204 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 203 "> \+ " 0 "" {MPLTEXT 1 -1 12 "h := 1./2^3;" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 17 "time_steps := 60;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 18 "num \+ := time_steps;" }}}{EXCHG {PARA 204 "" 0 "" {TEXT -1 2 " " }}{PARA 204 "" 0 "" {TEXT -1 78 " Set up the memory. ic is for the initial co nditions, y is for the solutions," }}{PARA 0 "" 0 "" {TEXT -1 89 " g \+ is for the RHS of the ODEs and out is for output of solutions - numeri cal or symbolic" }}{PARA 0 "" 0 "" {TEXT -1 89 " If you want the symb olic solution, set num := 1. If you do numerical then you can plot " } }{PARA 0 "" 0 "" {TEXT -1 30 " points using the vector out." }}{PARA 204 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 18 "ic := array(1..m):" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 17 "y := array(1..m):" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 17 "g := array(1..m):" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 26 "ou t := array(1..m,1..num):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 21 "ic[1] := evalf(Pi/4) ;" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 11 "ic[2] := 0;" }}} {EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 26 "ic[3] := evalf(sin(Pi/4)) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "ic[4] := evalf(cos(Pi/ 4));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 204 "" 0 "" {TEXT -1 0 "" }}{PARA 204 "" 0 "" {TEXT -1 29 " Set up t he RHS for the ODE" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 14 "g[1] := y[2];" }}}{EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 14 "g[2] := -y[3];" }}}{EXCHG {PARA 203 "> " 0 " " {MPLTEXT 1 -1 19 "g[3] := y[2]*y[4];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "g[4] := -y[2]*y[3];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 32 " The integration is done here." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 203 "> " 0 "" {MPLTEXT 1 -1 27 "for j from 1 by 1 to num \+ do" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 15 " print(j,j *h);" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 27 " for i from 1 by 1 to m do" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 18 " y[i] : = ic[i]:" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 5 " od: " }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 27 " for k from 2 by 1 to n do " }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 14 " #print(k) ;" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 30 " for i fr om 1 by 1 to m do " }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 31 " y[i] := subs(\{t=s\},y[i]):" }{TEXT -1 0 "" }}{PARA 203 " > " 0 "" {MPLTEXT 1 -1 7 " od:" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 29 " for i from 1 by 1 to m do" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 45 " v[i] := ic[i]+expand(int(g[i],s=0..t)):" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 7 " od:" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 29 " for i from 1 by 1 \+ to m do" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 19 " \+ y[i] := v[i];" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 7 " \+ od:" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 29 " for l from 1 by 1 to m do" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 40 " w := collect(sort(y[l],t),t): " }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 46 " for jj from degree(y[l],t) by -1 to k do" } {TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 65 " if (degr ee(w,t)=jj) then w := w - lcoeff(w,t)*t^(jj): fi:" }{TEXT -1 0 "" }} {PARA 203 "> " 0 "" {MPLTEXT 1 -1 9 " od:" }{TEXT -1 0 "" }} {PARA 203 "> " 0 "" {MPLTEXT 1 -1 16 " y[l] := w:" }{TEXT -1 0 " " }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 7 " od:" }{TEXT -1 0 "" }} {PARA 203 "> " 0 "" {MPLTEXT 1 -1 5 " od:" }{TEXT -1 0 "" }}{PARA 203 "> " 0 "" {MPLTEXT 1 0 27 " for i from 1 by 1 to m do" }}{PARA 203 "> " 0 "" {MPLTEXT 1 0 22 " out[i,j] := ic[i];" }}{PARA 203 "> \+ " 0 "" {MPLTEXT 1 0 25 " f := unapply(y[i],t);" }}{PARA 203 "> " 0 "" {MPLTEXT 1 0 18 " ic[i] := f(h);" }}{PARA 203 "> " 0 "" {MPLTEXT 1 -1 5 " od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 " In this example, if w e plot (out[4,j],out[3,j]) for j = 1 to num then we will get the plot \+ of the motion" }}{PARA 0 "" 0 "" {TEXT -1 45 " of a pendulum which i s an arc of a circle." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "plot([seq([out[4,j],out[3,j]],j=1..num)],scaling=cons trained);" }}}{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 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 8 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }