{VERSION 5 0 "IBM INTEL NT" "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 }
{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 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 64 "This worksheet determines Maclaurin polyn
omial solutions to the " }}{PARA 0 "" 0 "" {TEXT -1 52 "initial value \+
ordinary differential equation (IVODE)" }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 19 "          y'(t) =  " }{XPPEDIT 18 0 "
Diff(y(t),t) = G(y);" "6#/-%%DiffG6$-%\"yG6#%\"tGF*-%\"GG6#F(" }{TEXT 
-1 3 " ; " }{XPPEDIT 18 0 "y(0) = u;" "6#/-%\"yG6#\"\"!%\"uG" }{TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 "for " }{XPPEDIT 18 0 "y*epsilon*
R^n;" "6#*(%\"yG\"\"\"%(epsilonGF%)%\"RG%\"nGF%" }{TEXT -1 5 " and " }
{XPPEDIT 18 0 "G;" "6#%\"GG" }{TEXT -1 4 " = (" }{XPPEDIT 18 0 "g[1];
" "6#&%\"gG6#\"\"\"" }{TEXT -1 5 ",...," }{XPPEDIT 18 0 "g[n];" "6#&%
\"gG6#%\"nG" }{TEXT -1 4 ") : " }{XPPEDIT 18 0 "R^n;" "6#)%\"RG%\"nG" 
}{TEXT -1 4 " -> " }{XPPEDIT 18 0 "R^n;" "6#)%\"RG%\"nG" }{TEXT -1 33 
" a second degree polynomial in ( " }{XPPEDIT 18 0 "y[1];" "6#&%\"yG6#
\"\"\"" }{TEXT -1 5 ",...," }{XPPEDIT 18 0 "y[n];" "6#&%\"yG6#%\"nG" }
{TEXT -1 10 "), that is" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 112 
"  g[j] := a[j]+Sum(b[j,m]*y[m],m=1..n)+Sum(c[j,m]*y[m]^2,m=1..n)+Sum(
Sum(d[j,m,i]*y[m]*y[i],i=m+1..n),m=1..n-1);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 18 "   for j = 1,..,n." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 13 "   The terms \+
" }{XPPEDIT 18 0 "a[j];" "6#&%\"aG6#%\"jG" }{TEXT -1 26 " are the cons
tant term in " }{XPPEDIT 18 0 "g[j];" "6#&%\"gG6#%\"jG" }{TEXT -1 12 "
. The terms " }{XPPEDIT 18 0 "b[j,m];" "6#&%\"bG6$%\"jG%\"mG" }{TEXT 
-1 45 " are the coefficients of the linear terms in " }{XPPEDIT 18 0 "
g[j];" "6#&%\"gG6#%\"jG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 
13 "   The terms " }{XPPEDIT 18 0 "c[j,m];" "6#&%\"cG6$%\"jG%\"mG" }
{TEXT -1 48 " are the coefficients of the quadratic terms in " }
{XPPEDIT 18 0 "g[j];" "6#&%\"gG6#%\"jG" }{TEXT -1 12 ". The terms " }
{XPPEDIT 18 0 "d[j,m,i];" "6#&%\"dG6%%\"jG%\"mG%\"iG" }{TEXT -1 24 " a
re the coefficients of" }}{PARA 0 "" 0 "" {TEXT -1 28 "    the mixed p
roduct terms " }{XPPEDIT 18 0 "y[m]*y[i];" "6#*&&%\"yG6#%\"mG\"\"\"&F%
6#%\"iGF(" }{TEXT -1 4 " in " }{XPPEDIT 18 0 "g[j];" "6#&%\"gG6#%\"jG
" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "The solution to the IVODE is dete
rmined using Picard iteration and truncation on the Picard iterates to
 the " }}{PARA 0 "" 0 "" {TEXT -1 92 "Maclaurin polynomial. The comput
ation is all algebraic since we are integrating polynomials." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 42 "Here is an example of the IVODE for n = 4." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 26 "for j from 1 to 4  by 1 do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 117 "  Diff(y[j],t) = a[j]+sum(b[j,m]*y[m],m=1..4)+sum(c[
j,m]*y[m]^2,m=1..4)+sum(sum(d[j,m,i]*y[m]*y[i],i=m+1..4),m=1..3);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 " Give the num
ber of equations in your IVODE" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "number_of_equati
ons := 3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "n := number_of
_equations;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 107 " Give the degree of the Maclaurin polyno
mials you want for the approximation to the solution to your IVODE." }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 36 "degree_of_maclaurin_polynomial := 8;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 33 "  All the coefficients composing " }{XPPEDIT 18 0 "G;" "6#%\"GG
" }{TEXT -1 15 " are set to 0. " }}{PARA 0 "" 0 "" {TEXT -1 64 "  That
 way you only have to set the non-zero ones for your IVODE" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 43 "for j from 1 by 1 to number_of_equations do" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 13 "  a[j] := 0; " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for \+
i from 1 by 1 to number_of_equations do" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 43 "for j from 1 by 1 to number_of_equations do" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 15 "  b[i,j] := 0; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
15 "  c[i,j] := 0; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 43 "for i from 1 by 1 to number_of_equations do" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 43 "for j from 1 by 1 to number_of_equations do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for k from 1 by 1 to number_of_equa
tions do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "  d[i,j,k] := 0; " }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "  Th
is is where you give the values for your coefficients in " }{XPPEDIT 
18 0 "G;" "6#%\"GG" }{TEXT -1 44 ".  (I have set up for the Lorenz equ
ations.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 " b[1,1] := -sigma; b[1,2] := sigma;
 " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 " b[2,1] := r; b[2,2] :
= -1; d[2,1,3] := -1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 " b
[3,3] := -beta;  d[3,1,2] := 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 12 " #a[1] := 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 " b[
1,1] := 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " b[2,1] := -1
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " #c[1,1] := 1;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 42 "  You can now check to see if this is the " }{XPPEDIT 18 
0 "G;" "6#%\"GG" }{TEXT -1 11 " you wanted" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "unass
ign('m,i');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for j from 1
 to number_of_equations by 1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 184 
"  g[j] := a[j]+sum(b[j,m]*y[m],m=1..number_of_equations)+sum(c[j,m]*y
[m]^2,m=1..number_of_equations)+sum(sum(d[j,m,i]*y[m]*y[i],i=m+1..numb
er_of_equations),m=1..number_of_equations-1);" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "  G[j] := g[j];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 
"  for k from 1 to number_of_equations by 1 do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "    G[j] := subs(\{y[k]=u[k]\},G[j]);" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 5 "  od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "  p
rint(g[j],G[j]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 37 "  Set up the initial conditions in u." }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "u[
1] := 0; u[2] := 1; u[3] := 0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "  We now find the Macl
aurin polynomials algebraically!" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "  The 0th and
 1st Maclaurin and Picard polynomials are the same! The 1st " }}{PARA 
0 "" 0 "" {TEXT -1 55 "  Maclaurin polynomial for each component is ou
tputted." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for i from 1 to number_of_equations
 by 1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "  U[i,0] := u[i];   " }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "  alpha[i,0] := u[i];" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 28 "  U[i,1] := U[i,0] + G[i]*t;" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 21 "  alpha[i,1] := G[i];" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 16 "  print(U[i,1]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 
"od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 77 " Now we generate the rest of the story! Since we \+
are squaring polynomials and" }}{PARA 0 "" 0 "" {TEXT -1 86 "then inte
grating polynomials in the Picard iteration it is easy to determine wh
at the " }}{PARA 0 "" 0 "" {TEXT -1 84 "next term in the Maclaurin pol
ynomial is from the Picard iterate. The coefficient of" }}{PARA 0 "" 
0 "" {TEXT -1 59 "the next term in the Maclaurin polynomial for the so
lution " }{XPPEDIT 18 0 "y[j];" "6#&%\"yG6#%\"jG" }{TEXT -1 13 " is gi
ven by " }{XPPEDIT 18 0 "alpha[j,k];" "6#&%&alphaG6$%\"jG%\"kG" }}
{PARA 0 "" 0 "" {TEXT -1 81 "where k is the degree of the next term in
 the Maclaurin polynomial. The Maclaurin" }}{PARA 0 "" 0 "" {TEXT -1 
23 "polynomial solution to " }{XPPEDIT 18 0 "y[j];" "6#&%\"yG6#%\"jG" 
}{TEXT -1 13 " is given by " }{XPPEDIT 18 0 "U[j,k];" "6#&%\"UG6$%\"jG
%\"kG" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "unassign('m','j','i','k');
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 186 "alpha[N,k] := (Sum(b[N
,j]*alpha[j,k-1],j=1..J)+Sum(c[N,j]*Sum(alpha[j,m]*alpha[j,k-m-1],m=0.
.k-1),j=1..J)+Sum(Sum(d[N,m,j]*Sum(alpha[m,i]*alpha[j,k-i-1],i=0..k-1)
,j=m+1..J),m=1..J-1))/k;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "    U[N
,k] := U[N,k-1] + alpha[N,k]*t^k;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "unassig
n('m','j','i');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "for k fr
om 2 by 1 to degree_of_maclaurin_polynomial do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 47 "  for num from 1 by 1 to number_of_equations do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 206 "    alpha[num,k] := expand((sum(b[
num,j]*alpha[j,k-1],j=1..n)+sum(c[num,j]*sum(alpha[j,m]*alpha[j,k-m-1]
,m=0..k-1),j=1..n)+sum(sum(d[num,m,j]*sum(alpha[m,i]*alpha[j,k-i-1],i=
0..k-1),j=m+1..n),m=1..n-1)))/k;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 
"    U[num,k] := U[num,k-1] + alpha[num,k]*t^k;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 5 "  od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "for num from 1 by 1 to numbe
r_of_equations do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "  print(U[num,
degree_of_maclaurin_polynomial]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 
"end;" }}}{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 }