Using Mathematica to expand a taylor series to prove the order of a Runge-Kutta method

Question

Given the Runge-Kutta method, $$w_{i+1}=w_i+\frac{1}{4}k_1 +\frac{3}{8}k_2+\frac{3}{8}k_3 $$ where \begin{align} k_1 &= hf(t_i,w_i)\\ k_2 &= hf(t_i+\frac{2}{3}h,w_i+\frac{2}{3}k_1)\\ k_3 &= hf(t_i+\frac{2}{3}h,w_i+\frac{2}{3}k_2) \end{align} and $w_i$ is the approximation of $y$ the true solution at step $i$.

Show by matching Taylor series expansion of the true solution $y$, that the above method has an error of order $O(h^4)$.

Notes: This is not a question on how to solve the above question but how to do the necessary multivariate Taylor series expansion in Mathematica. Whatever I try I gives strange answers. I would like to have a general procedure for this type of problem as I am having trouble proving the order of a 5th order Runge-Kutta methods by hand.

If you need any clarification or explanation feel free to leave a comment.

Code Examples

Clear[k1, k2]
n = 2;
k1 = h f[t, w]
k2 = Series[h f[t + (2 h)/3, w + (2 k1)/3], {t, 0, n}, {w, 0, n}]
Clear[k1, k2]
k1 = h f[t, w]
k2 = Series[h f[t, w], {t, (2 h)/3, n}, {w, (2 k1)/3, n}]

Output

Answer 1

The procedure(solve ode y'[t]== f[t,y[t]] ) you want to analyse is of the form

k1 = f[t, y[t]];
k2 = f[t + 2/3 h, y[t] + 2/3 h k1];
k3 = f[t + 2/3 h, y[t] + 2/3 h k2];
proc = y[t + h] - y[t] - h (1/4 k1 + 3/8 k2 + 3/8 k3 ) 

In the last equation the series expansion of y[t+h]-y[t] consists of the derivatives y'[t],y''[t],...which have to be replaced by the knowledge of the ode y'[t]== f[t,y[t]], y''[t]=...,..:

n = 3; (* maximum order *)
subst = MapThread[#1 -> #2 &, {Table[D[y[t], {t, k}], {k, 1, n}],NestList[(D[#, t] /. y'[t] -> f[t, y[t]]) &, f[t, y[t]], n - 1]}];   

With this substitution you can do the seriesexpansion of your procedure

Collect[ Normal[Series[proc/.subst, {h, 0, 4}] ] /. subst , h, Simplify]

which contains h^4 as smallest part. So your procedure is of order O[h^3]