In[1]:=

(* Preliminary Setup : Define Modules *)Off[General :: "spell1"]  ... ng the method *)ClearAll[yhalf1, yhalf2, yhalf3, yhalf4, EX1, EX2, EX3, ye1, ye2, ye3] ;

Out[16]=

y0 + (f1[y0] + f2[y0] + f3[y0]) h + (1/2 f1[y0] f1^′[y0] + 1/2 f2[y0] f1^′[y0] + ... ]/2) f3^′[y0]^2 + (f1[y0] + f2[y0] + f3[y0]/2) f3[y0] f3^′′[y0])) h^3 + O[h]^4

In[18]:=

(* Define the Commutator for two vector fields *)Com[fa_, fb_] := fb * D[fa, y[s]] - ... rence between method and exact *)DY = Expand[Collect[y1[h] - ye[h], h]] + O[h]^(i + 1) ;

Out[22]=

y0 + (f1[y0] + f2[y0] + f3[y0]) h + 1/2 (f1[y0] + f2[y0] + f3[y0]) (f1^′[y0] + f2^R ... 8242;[y0] - f3[y0] (f2^′[y0] f3^′[y0] + f2[y0] f3^′′[y0]))) h^3 + O[h]^5

In[24]:=

(* Solve for the residual *)T = Simplify[Flatten[Solve[DY  O[h]^(i + 1), gres[y0]] //. {y0  y}]]

Out[24]=

{gres[y] f3[y]^2 (-(-1 + a133) f1^′′[y] - (-1 + a233) f2^′′[ ... ′′[y] + (2 + a123 - a231) f2^′′[y] + (-1 + a133) f3^′′[y]))}

In[25]:=

(* By inspection we can determine many of the coefficients *)a121 = 2 ; a122 = 1 ; a ... l *)T = Simplify[Flatten[Solve[DY  O[h]^(i + 1), gres[y0]] //. {y0  y}]]

Out[26]=

{gres[y]  -(-2 + a123 + a132) f1[y] f2^′[y] f3^′[y] + f3[y] ((-4 + a132 ... ((2 + a123 - a231) f1^′[y] f3^′[y] + (-4 + a132 + a231) f1[y] f3^′′[y])}

In[27]:=

(* Due to cyclic permutation symmetry this is underdetermined *)(* Hence, any soluti ... l *)T = Simplify[Flatten[Solve[DY  O[h]^(i + 1), gres[y0]] //. {y0  y}]]

Out[28]=

{gres[y] 0}

In[29]:=

(* Hence the first modified vector field can be written as : f[y] = f1[y] + ... ere a is an arbitrary constant, and {f1, f2} is the commutator of f1 and f2 *)   

Created by Mathematica  (March 23, 2007) Valid XHTML 1.1!