Implicit-Euler.nb

In[1]:=

$(* Preliminary Setup : Define Modules *) CreateRules(EQ_, N_) := Module[{}, Flatten[Ta ... Expand the method in a series to order i *)y1[h] = Simplify[Series[yb[h], {h, 0, i}] //.R]$

Out[6]=

In[7]:=

$(* The perturbed vector field is defined here *)FP[y_] = f[y] + h^(i - 1) * G[y] ; ... i], y[s] y0}] ; ye[h] = Simplify[Series[y[s + h], {h, 0, i}]//.R]$

Out[10]=

In[11]:=

(* Expand the difference between method and exact *)DY = Expand[Collect[y1[h] - ye[h], ... ;(* Solve for the perturbation *)T1 = Flatten[Solve[DY  0, G[y0]]] //. y0y

Out[12]=

${G[y] 1/2 f[y] f^′[y]}$

In[13]:=

(* ROUND TWO : COMPUTE THE SECOND TERM *)(* Set i to be the top - ... Expand the method in a series to order i *)y1[h] = Simplify[Series[yb[h], {h, 0, i}] //.R]

Out[15]=

y0 + f[y0] h + f[y0] f^′[y0] h^2 + (f[y0] f^′[y0]^2 + 1/2 f[y0]^2 f^′′[y0]) h^3 + O[h]^4

In[16]:=

$(* The perturbed vector field is defined here *)FP[y_] = (FP[y] /. T1) + h^(i - 1) * G ... i], y[s] y0}] ; ye[h] = Simplify[Series[y[s + h], {h, 0, i}]//.R]$

Out[19]=

y0 + f[y0] h + f[y0] f^′[y0] h^2 + (G[y0] + 1/12 f[y0] (8 f^′[y0]^2 + 5 f[y0] f^′′[y0])) h^3 + O[h]^4

In[20]:=

(* Expand the difference between method and exact *)DY = Expand[Collect[y1[h] - ye[h], ... 71;(* Solve for the perturbation *)T1 = Flatten[Solve[DY  0, G[y0]]]//.y0y

Out[21]=

${G[y] 1/12 (4 f[y] f^′[y]^2 + f[y]^2 f^′′[y])}$

In[22]:=

(* ROUND THREE : COMPUTE THE THIRD TERM *)(* Set i to be the top - ... Expand the method in a series to order i *)y1[h] = Simplify[Series[yb[h], {h, 0, i}] //.R]

Out[24]=

y0 + f[y0] h + f[y0] f^′[y0] h^2 + 1/2 f[y0] (2 f^′[y0]^2 + f[y0] f^′′ ... (6 f^′[y0]^3 + 9 f[y0] f^′[y0] f^′′[y0] + f[y0]^2 f^(3)[y0]) h^4 + O[h]^5

In[25]:=

$(* The perturbed vector field is defined here *)FP[y_] = (FP[y] /. T1) + h^(i - 1) * G ... i], y[s] y0}] ; ye[h] = Simplify[Series[y[s + h], {h, 0, i}]//.R]$

Out[28]=

y0 + f[y0] h + f[y0] f^′[y0] h^2 + 1/2 f[y0] (2 f^′[y0]^2 + f[y0] f^′′ ... f^′[y0]^3 + 16 f[y0] f^′[y0] f^′′[y0] + 2 f[y0]^2 f^(3)[y0])) h^4 + O[h]^5

In[29]:=

(* Expand the difference between method and exact *)DY = Expand[Collect[y1[h] - ye[h], ... 71;(* Solve for the perturbation *)T1 = Flatten[Solve[DY  0, G[y0]]]//.y0y

Out[30]=

${G[y] 1/12 (3 f[y] f^′[y]^3 + 2 f[y]^2 f^′[y] f^′′[y])}$

In[31]:=

(* ROUND FOUR : VERIFY THE EXPANSION *)(* Set i to be the top - order retained ... i], y[s] y0}] ; ye[h] = Simplify[Series[y[s + h], {h, 0, i}]//.R]

Out[35]=

y0 + f[y0] h + f[y0] f^′[y0] h^2 + 1/2 f[y0] (2 f^′[y0]^2 + f[y0] f^′′ ... (6 f^′[y0]^3 + 9 f[y0] f^′[y0] f^′′[y0] + f[y0]^2 f^(3)[y0]) h^4 + O[h]^5

In[36]:=

(* Expand the difference between method and exact *)DY = Expand[Collect[y1[h] - ye[h], h]] + O[h]^(i + 1)

Out[36]=

In[37]:=

(* Final Modified Vector Field *)FP[y] + O[h]^i

Out[37]=

f[y] + 1/2 f[y] f^′[y] h + (1/3 f[y] f^′[y]^2 + 1/12 f[y]^2 f^′′[y]) h^2 + (1/4 f[y] f^′[y]^3 + 1/6 f[y]^2 f^′[y] f^′′[y]) h^3 + O[h]^4

Created by Mathematica (June 11, 2004)