Editing Runge–Kutta methods (section)

==The Runge–Kutta method==

[[File:Runge-Kutta slopes.svg|thumb|Slopes used by the classical Runge-Kutta method|alt=|300x300px]]
The most widely known member of the Runge–Kutta family is generally referred to as "RK4", the "classic Runge–Kutta method" or simply as "the Runge–Kutta method".

Let an [[initial value problem]] be specified as follows:

:<math> \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0. </math>

Here <math>y</math> is an unknown function (scalar or vector) of time <math>t</math>, which we would like to approximate; we are told that <math>\frac{dy}{dt}</math>, the rate at which <math>y</math> changes, is a function of <math>t</math> and of <math>y</math> itself. At the initial time <math>t_0</math> the corresponding <math>y</math> value is <math>y_0</math>. The function <math>f</math> and the [[initial conditions]] <math>t_0</math>,  <math>y_0</math> are given.

Now we pick a step-size ''h'' &gt; 0 and define:

:<math>\begin{align}
y_{n+1} &= y_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4 \right),\\
t_{n+1} &= t_n + h \\
\end{align}</math>

for ''n'' = 0, 1, 2, 3, ..., using<ref>{{harvnb|Press|Teukolsky|Vetterling|Flannery|2007|p=908}}; {{harvnb|Süli|Mayers|2003|p=328}}</ref>
:<math>
\begin{align}
 k_1 &= \ f(t_n, y_n), \\
 k_2 &= \ f\!\left(t_n + \frac{h}{2}, y_n + h \frac{k_1}{2}\right), \\ 
 k_3 &= \ f\!\left(t_n + \frac{h}{2}, y_n + h \frac{k_2}{2}\right), \\
 k_4 &= \ f\!\left(t_n + h, y_n + h k_3\right).
\end{align}
</math>
(''Note: the above equations have different but equivalent definitions in different texts.''<ref name=notation>{{harvtxt|Atkinson|1989|p=423}}, {{harvtxt|Hairer|Nørsett|Wanner|1993|p=134}}, {{harvtxt|Kaw|Kalu|2008|loc=§8.4}} and {{harvtxt|Stoer|Bulirsch|2002|p=476}} leave out the factor ''h'' in the definition of the stages. {{harvtxt|Ascher|Petzold|1998|p=81}}, {{harvtxt|Butcher|2008|p=93}} and {{harvtxt|Iserles|1996|p=38}} use the ''y'' values as stages.</ref>)

Here <math>y_{n+1}</math> is the RK4 approximation of <math>y(t_{n+1})</math>, and the next value (<math>y_{n+1}</math>) is determined by the present value (<math>y_n</math>) plus the [[weighted average]] of four increments, where each increment is the product of the size of the interval, ''h'', and an estimated slope specified by function ''f'' on the right-hand side of the differential equation.
* <math>k_1</math> is the slope at the beginning of the interval, using <math> y </math> ([[Euler's method]]);
* <math>k_2</math> is the slope at the midpoint of the interval, using <math> y </math> and <math> k_1 </math>;
* <math>k_3</math> is again the slope at the midpoint, but now using <math> y </math> and <math> k_2 </math>;
* <math>k_4</math> is the slope at the end of the interval, using <math> y </math> and <math> k_3 </math>.

In averaging the four slopes, greater weight is given to the slopes at the midpoint. If <math>f</math> is independent of <math>y</math>, so that the differential equation is equivalent to a simple integral, then RK4 is [[Simpson's rule]].<ref name="Süli 2003 328">{{harvnb|Süli|Mayers|2003|p=328}}</ref>

The RK4 method is a fourth-order method, meaning that the [[Truncation error (numerical integration)|local truncation error]] is [[Big O notation|on the order of]] <math>O(h^5)</math>, while the  [[Truncation error (numerical integration)|total accumulated error]] is on the order of <math>O(h^4)</math>.

In many practical applications the function <math>f</math> is independent of <math>t</math> (so called [[Autonomous system (mathematics)|autonomous system]], or time-invariant system, especially in physics), and their increments are not computed at all and not passed to function <math>f</math>, with only the final formula for <math>t_{n+1}</math> used.