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

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.

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

• Some details and/or code snippets of what you have tried may help a lot. Dec 3, 2017 at 18:58
• @drN I added code examples. But I don't think they will be very helpful. I have not even been able to get an expansion for $k_2$.
– AzJ
Dec 3, 2017 at 19:07
• @drN What other details do you require?
– AzJ
Dec 3, 2017 at 19:08
• An important detail I'm missing is the information w'[t]=f[t,w[t]] ! Dec 4, 2017 at 8:23
• Is the answer of further interest? If so you have to substitute w(i)=w[ti],w(i+1)=w[ti+h] in the taylor expansion too. The Series-Expansion of this expression shows the consistence of the considered rungekutta method Dec 5, 2017 at 7:54

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]

• Thank you for your answer. Though I understand what you did but combing so many commands in 1 or 2 lines makes it difficult to understand what steps mathematica is doing.
– AzJ
Dec 5, 2017 at 20:58
• Sorry... ! The essential part is NestList[], I have no idea how to simplify this step. Dec 5, 2017 at 21:55