1
$\begingroup$

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

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

1 Answer 1

Reset to default
3
$\begingroup$

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]

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.