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 enter image description here

  • 1
    $\begingroup$ Some details and/or code snippets of what you have tried may help a lot. $\endgroup$
    – dearN
    Dec 3, 2017 at 18:58
  • $\begingroup$ @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$. $\endgroup$
    – AzJ
    Dec 3, 2017 at 19:07
  • $\begingroup$ @drN What other details do you require? $\endgroup$
    – AzJ
    Dec 3, 2017 at 19:08
  • 1
    $\begingroup$ An important detail I'm missing is the information w'[t]=f[t,w[t]] ! $\endgroup$ Dec 4, 2017 at 8:23
  • 1
    $\begingroup$ 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 $\endgroup$ Dec 5, 2017 at 7:54

1 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]

  • $\begingroup$ 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. $\endgroup$
    – AzJ
    Dec 5, 2017 at 20:58
  • $\begingroup$ Sorry... ! The essential part is NestList[], I have no idea how to simplify this step. $\endgroup$ Dec 5, 2017 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.