{VERSION 6 0 "IBM INTEL NT" "6.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 "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{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 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 134 " This worksheet generates the Maclauri n power series solution for J quadratic polynomial (degree 2) differen tial equations. That is," }}{PARA 0 "" 0 "" {TEXT -1 73 " we generat e the Maclaurin polynomial approximation for the solution to" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 33 " y' = P(y) ; y(0) = a" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 128 " where P is a polynomial of degree 2 to the degree cho sen by the user. If the user wants a 10th degree Maclaurin polynomial \+ we" }}{PARA 0 "" 0 "" {TEXT -1 125 " generate it as easily as the 10 0th degree Maclaurin polynomial. We call this method the Power Series (of Parker-Sochacki)" }}{PARA 0 "" 0 "" {TEXT -1 10 " method." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 " We write the polynomial P as" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "P[j] = a[j]+Sum(b [j,m]*y[m],m = 1 .. J)+Sum(Sum(d[j,m,i]*y[i]*y[m],i = m ..J),m = 1 .. \+ J);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 18 " for j = 1,..,J." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 30 " For example if J = 4 we have" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "J := 4;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 90 "P[j] = a[j]+sum(b[j,m]*y[m],m = 1 .. J)+sum(sum(d[j ,m,i]*y[i]*y[m],i = m ..J),m = 1 .. J);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "unassign('J');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{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 3 " " }}{PARA 0 "" 0 "" {TEXT -1 13 " The terms " } {XPPEDIT 18 0 "a[j];" "6#&%\"aG6#%\"jG" }{TEXT -1 26 " are the constan t term in " }{XPPEDIT 18 0 "P[j];" "6#&%\"PG6#%\"jG" }{TEXT -1 12 ". T he 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 "P[j ];" "6#&%\"PG6#%\"jG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 13 " \+ The terms " }{XPPEDIT 18 0 "d[j,m,i];" "6#&%\"dG6%%\"jG%\"mG%\"iG" } {TEXT -1 49 " are the coefficients of the mixed product terms " } {XPPEDIT 18 0 "y[m]*y[i];" "6#*&&%\"yG6#%\"mG\"\"\"&F%6#%\"iGF(" } {TEXT -1 4 " in " }{XPPEDIT 18 0 "P[j];" "6#&%\"PG6#%\"jG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 101 " The solution to the IVODE is determined using C auchy products to determine the coefficients of the " }}{PARA 0 "" 0 " " {TEXT -1 93 " Maclaurin polynomial. The computation is all algebrai c since we are multiplying and adding." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 27 " We assume each so lution " }{XPPEDIT 18 0 "y[j];" "6#&%\"yG6#%\"jG" }{TEXT -1 13 " for j = 1..J" }{TEXT 258 1 " " }{TEXT -1 26 "has the power series form " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "y[j] := sum(alpha[ j,i],i=0..nd);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 76 " where nd is the degree of the Maclaurin polynomial desired by the \+ user. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 59 " The coefficients for j = 1,..,J are given by the \+ following" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 " The firs t coefficient is given by the initial conditions of the solutions " } {TEXT 256 1 "y" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "alpha[j,0] = y[j ](0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 " The second coefficient is given by the diffe rential equation y' = P(y)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "alpha[j,1] = a[j]+su m(b[j,m]*y[m],m=1..J)+sum(sum(d[j,m,i]*y[m]*y[i],i=m..J),m=1..J);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 " The rest of the coeffi cients are given by Cauchy Products as detailed in " }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 84 " SOME PROPERTIES OF SOLUTIONS TO POLYNOMIAL SYSTEMS OF DIFFERENTIAL EQUATIONS" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 44 " by Car others, Parker, Sochacki and Warne." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 138 "alpha[j,k] = (Sum(b[j,n]*alpha[n,k-1],n = 1 .. J) +Sum(Sum(d[j,m,n]*Sum(alpha[m,i]*alpha[n,k-i-1],i = 0 .. k-1),n = m .. J),m = 1 .. J))/k;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 35 " for k = 2,3,... and j = 1,..., J" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 10 " We let " } {XPPEDIT 18 0 "U[j];" "6#&%\"UG6#%\"jG" }{TEXT -1 33 " be the Maclauri n polynomial for " }{XPPEDIT 18 0 "y[j];" "6#&%\"yG6#%\"jG" }{TEXT -1 6 ". Then" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "U[j,k ] = U[j,k-1] + alpha[j,k]*t^k;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 " We now give an example of the ODE for n = 4." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for j from 1 to 4 by 1 d o" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 " Diff(y[j],t) = a[j]+sum(b[j, m]*y[m],m=1..4)+sum(sum(d[j,m,i]*y[m]*y[i],i=1..4),m=1..4);" }}{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 74 " You now choose the number of equations you have in y our quadratic ODE" }}{PARA 0 "" 0 "" {TEXT -1 92 " and the degree you want for your Maclaurin polynomial approximate solution to this O DE" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "number_ of_equations := 4;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "degre e_of_maclaurin_polynomial := 8;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 34 " All \+ the coefficients composing " }{XPPEDIT 18 0 "P;" "6#%\"PG" }{TEXT -1 15 " are set to 0. " }}{PARA 0 "" 0 "" {TEXT -1 64 " That way you onl y have to set the non-zero ones for your IVODE" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{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_equa tions 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 16 " b[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_equations 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 0 "" }}{PARA 0 "" 0 "" {TEXT -1 96 " Give the values for the coefficients a,b and c that make up the RHS of the ODE which is P.." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "a[1] := 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "b[2,3] := 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "d[3,3,4] := -2;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "b[4,3] := 1; d[4,4,4] := -2; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 42 "#b[1,1] := 1; b[1,2] := 1; d[1,1,2] := -1;" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "#d[2,2,2] := -1;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 37 " Set up the initial co nditions in u." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 " u[1] := 0; u[2] := 0; u[3] := 1; u[4] := 0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 42 " You can now check to see if this is the " }{XPPEDIT 18 0 "P;" "6#%\" PG" }{TEXT -1 12 " you wanted." }}{PARA 0 "" 0 "" {TEXT -1 45 " Your \+ ODE is printed below in terms of y[j] " }}{PARA 0 "" 0 "" {TEXT -1 70 " following your ODE is the value of the ODE at the initial conditio n" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "unassig n('m,i');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for j from 1 t o number_of_equations by 1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 136 " \+ g[j] := a[j]+sum(b[j,m]*y[m],m=1..number_of_equations)+sum(sum(d[j,m, i]*y[m]*y[i],i=m..number_of_equations),m=1..number_of_equations);" }} {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 " print(g[j],G[j]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 " unassign('m','i');" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 87 " We now find the Mac laurin polynomials algebraically up to the degree you chose above." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 73 " The 0th and 1st Maclaurin and Picard polynomials are th e same! The 1st " }}{PARA 0 "" 0 "" {TEXT -1 55 " Maclaurin polynomia l for each component is outputted." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} }{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 45 " print(sort(collect(simplify(U[i,1]),t),t));" }}{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 80 " Now we generate the rest of the story! Since we are squaring p olynomials and" }}{PARA 0 "" 0 "" {TEXT -1 88 " then integrating poly nomials in the Picard iteration it is easy to determine what the " }} {PARA 0 "" 0 "" {TEXT -1 86 " next term in the Maclaurin polynomial i s from the Picard iterate. The coefficient of" }}{PARA 0 "" 0 "" {TEXT -1 61 " the next term in the Maclaurin polynomial for the solut ion " }{XPPEDIT 18 0 "y[j];" "6#&%\"yG6#%\"jG" }{TEXT -1 13 " is given by " }{XPPEDIT 18 0 "alpha[j,k];" "6#&%&alphaG6$%\"jG%\"kG" }}{PARA 0 "" 0 "" {TEXT -1 83 " where k is the degree of the next term in the Maclaurin polynomial. The Maclaurin" }}{PARA 0 "" 0 "" {TEXT -1 25 " \+ 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 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "unassign('m','j',' i');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "J := number_of_equa tions;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "for k from 2 by 1 to degree_of_maclaurin_polynomial do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " print(degree=k);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 " for n um from 1 by 1 to number_of_equations do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 141 " alpha[num,k] := expand((sum(b[num,j]*alpha[j,k-1 ],j=1..J)+sum(sum(d[num,m,j]*sum(alpha[m,i]*alpha[j,k-i-1],i=0..k-1),j =m..J),m=1..J)))/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 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 63 " Here is the Maclaurin polynomial for ea ch of your components." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "for j from 1 by 1 to J do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " print(U[j,k-1]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 " \+ If you want a numerical solution to your ODE you can now choose the t ime step" }}{PARA 0 "" 0 "" {TEXT -1 31 " and the number of time step s." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "h := 1/ 2^5.;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "num_time_steps := \+ 500;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "h*num_time_steps;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 " Restart for the nume rical solution. The vector out[] contains the solution for j = 1..J." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 "\015" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "unassign('m' ,'j','i');\015" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "ne := num ber_of_equations;\015" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "nd := degree_of_maclaurin_polynomial;\015" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "for n from 1 by 1 to num_time_steps do\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " print(step,n); \015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 " for num from 1 by 1 to number_of_equations do \015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 129 " alpha[num,1] := a[num] +sum(b[num,mm]*alpha[mm,0],mm=1..ne)+sum(sum(d[num,mn,jj]*alpha[mn,0]* alpha[jj,0],jj=1..ne),mn=1..ne);\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od;\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 " for k from 2 by 1 to degree_of_maclaurin_polynomial do\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 " for num from 1 by 1 to number_of_equations do\015 " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " #print(eqn,num);\015" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 139 " alpha[num,k] := (sum(b[num,j ]*alpha[j,k-1],j=1..ne)+sum(sum(d[num,m,j]*sum(alpha[m,i]*alpha[j,k-i- 1],i=0..k-1),j=1..ne),m=1..ne))/k;\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 " od:\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od:\015" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 " for num from 1 by 1 to number_of_ equations do\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 " temp := al pha[num,nd-1]+alpha[num,nd]*h;\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 " for jm from 1 by 1 to nd-1 do\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 " temp := alpha[num,nd-jm-1]+temp*h;\015" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 " od:\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 " alpha[num,0] := temp;\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od:\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 " for \+ num from 1 by 1 to number_of_equations do\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 " out[num,n] := alpha[num,0];\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od:\015" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "od :\015" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "for j from 1 by 1 to J do listplot([seq([ i*h,out[j,i]],i=1..num_time_steps)]); end;" }}}{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 }