Lecture Notes

2.6 Runge - Kutta Method

Processed February 25, 2026 at 09:48 AM
Original Handwritten Notes
PAGE 1

2.6 Runge - Kutta Method

solve \( y' = f(t, y) \)

basic idea : replacing the slope part in Euler

\[ y_{n+1} = y_n + \underbrace{f(t_n, y_n)}_{\text{slope}} h \]

Runge - Kutta order 4 (RK4) :

\[ y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)h \]
  • \( k_1 = f(t_n, y_n) \) slope at beginning (Euler)
  • \( k_2 = f(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_1 h) \) slope at mid pt, using Euler (\(k_1\)) to get there
  • \( k_3 = f(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_2 h) \) slope at mid pt, using \(k_2\) to get there
  • \( k_4 = f(t_n + h, y_n + k_3 h) \) slope at end, using \(k_3\) to get there

\( \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \) is weighted avg. of slopes

PAGE 2
A coordinate graph illustrating the Runge-Kutta 4th order (RK4) numerical method. The horizontal axis represents time \( t \) with marked points at \( t_0 \), \( t_0 + h/2 \), and \( t_0 + h \). The vertical axis represents \( y \) with marked points at \( y_0 \), \( y_0 + hk_1/2 \), \( y_0 + hk_2/2 \), and \( y_0 + hk_3 \). A blue curve labeled \( y(t) \) represents the solution. Four red vectors represent the slope estimates: \( k_1 \) is the initial slope at \( (t_0, y_0) \); \( k_2 \) is the slope at the midpoint \( t_0 + h/2 \) estimated using \( k_1 \); \( k_3 \) is a second slope estimate at the midpoint using \( k_2 \); and \( k_4 \) is the slope at the end of the interval \( t_0 + h \) estimated using \( k_3 \). The final computed point is labeled \( (t_1, y_1) \) and is located at the end of the interval.
Figure: A coordinate graph illustrating the Runge-Kutta 4th order (RK4) numerical method. The horizontal axis represents time \( t \) with marked points at \( t_0 \), \( t_0 + h/2 \), and \( t_0 + h \). The vertical axis represents \( y \) with marked points at \( y_0 \), \( y_0 + hk_1/2 \), \( y_0 + hk_2/2 \), and \( y_0 + hk_3 \). A blue curve labeled \( y(t) \) represents the solution. Four red vectors represent the slope estimates: \( k_1 \) is the initial slope at \( (t_0, y_0) \); \( k_2 \) is the slope at the midpoint \( t_0 + h/2 \) estimated using \( k_1 \); \( k_3 \) is a second slope estimate at the midpoint using \( k_2 \); and \( k_4 \) is the slope at the end of the interval \( t_0 + h \) estimated using \( k_3 \). The final computed point is labeled \( (t_1, y_1) \) and is located at the end of the interval.
PAGE 3

advantages over Euler :

Stability

A coordinate graph showing a vertical y-axis and a horizontal x-axis. A solid green curve starts high on the y-axis and decays exponentially towards the x-axis. A dashed black line starts at the same point on the y-axis but has a steeper negative slope, crossing the x-axis and continuing into the negative y-region, ending at a red point. This illustrates how a large step size in a numerical method can lead to instability.
Figure: A coordinate graph showing a vertical y-axis and a horizontal x-axis. A solid green curve starts high on the y-axis and decays exponentially towards the x-axis. A dashed black line starts at the same point on the y-axis but has a steeper negative slope, crossing the x-axis and continuing into the negative y-region, ending at a red point. This illustrates how a large step size in a numerical method can lead to instability.

\( y' = -10y \)

could blow up if step is too big

RK4 uses two mid pt slopes and the "corrects" the error w/ large h, so the situation above is less likely to happen

accuracy:

  • error in Euler: global error is proportional to \( h \)
  • in Imp. Euler (RK2): global error is proportional to \( h^2 \)
  • in RK4: global error is proportional to \( h^4 \)

half the step size \( \rightarrow \) reduce error by factor of 16

PAGE 4

example

\( y' = 1 - t + 4y \)
\( y(0) = 1 \)
\( h = 0.1 \) to estimate \( y(0.1) \)

\( t_0 = 0 \)
\( y_0 = 1 \) given

\( t_1 = t_0 + h = 0.1 \) (target)

\( k_1 = f(t_0, y_0) = 1 - (0) + 4(1) = 5 \)

\( k_2 = f(t_0 + \frac{1}{2}h, y_0 + \frac{1}{2}k_1 h) \)
\( = 1 - (0 + 0.05) + 4(1 + \frac{1}{2} \cdot 5 \cdot 0.1) = 5.95 \)

\( k_3 = f(t_0 + \frac{1}{2}h, y_0 + \frac{1}{2}k_2 h) \)
\( = 1 - (0 + 0.05) + 4(1 + \frac{1}{2} \cdot 5.95 \cdot 0.1) = 6.14 \)

\( k_4 = f(t_0 + h, y_0 + k_3 h) \)
\( = 1 - (0 + 0.1) + 4(1 + 6.14 \cdot 0.1) = 7.356 \)

\( y_1^* = y_0 + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)h \)
\( = 1.60893 \)

PAGE 5

How good?

  • RK4 w/ \( h=0.1 \rightarrow 1.60893 \) error \( 1.085 \times 10^{-4} \)
  • true value \( \rightarrow 1.60904 \)
  • RK4 w/ \( h=0.05 \rightarrow 1.609034 \) error \( 8 \times 10^{-6} \) (reduction by factor of 14)
  • Euler w/ \( h=0.1 \rightarrow 1.5 \)
  • \( h=0.01 \rightarrow 1.59529 \) error \( 1 \times 10^{-2} \)
  • \( h=0.005 \rightarrow 1.60206 \) error \( 7 \times 10^{-3} \) (50 steps)
  • RK2 w/ \( h=0.1 \rightarrow 1.595 \)
  • w/ \( h=0.01 \rightarrow 1.60886 \)

in this situation, 1 RK4 step beats all w/ same \( h \)

beyond RK4 (4th order), the math gets messy

common ones: RK 4/5

PAGE 6

Runge-Kutta-Fehlberg (RKF45)

The RKF45 method uses six stages to calculate both a 4th-order and a 5th-order approximation, allowing for adaptive step-size control.

\[k_1 = f(t_n, y_n)\]\[k_2 = f(t_n + \frac{1}{4}h, y_n + \frac{1}{4}hk_1)\]\[k_3 = f(t_n + \frac{3}{8}h, y_n + \frac{3}{32}hk_1 + \frac{9}{32}hk_2)\]\[k_4 = f(t_n + \frac{12}{13}h, y_n + \frac{1932}{2197}hk_1 - \frac{7200}{2197}hk_2 + \frac{7296}{2197}hk_3)\]\[k_5 = f(t_n + h, y_n + \frac{439}{216}hk_1 - 8hk_2 + \frac{3680}{513}hk_3 - \frac{845}{4104}hk_4)\]\[k_6 = f(t_n + \frac{1}{2}h, y_n - \frac{8}{27}hk_1 + 2hk_2 - \frac{3544}{2565}hk_3 + \frac{1859}{4104}hk_4 - \frac{11}{40}hk_5)\]

Final Estimates

The 4th-order estimate and the 5th-order estimate are calculated as follows:

\[y_{n+1} = y_n + h \left( \frac{25}{216}k_1 + \frac{1408}{2565}k_3 + \frac{2197}{4104}k_4 - \frac{1}{5}k_5 \right)\]\[\hat{y}_{n+1} = y_n + h \left( \frac{16}{135}k_1 + \frac{6656}{12825}k_3 + \frac{28561}{56430}k_4 - \frac{9}{50}k_5 + \frac{2}{55}k_6 \right)\]
PAGE 7

Matlab ode45 uses Dormand-Prince 4/5-order

compares \( y_{n+1} \) from 4th vs 5th

if \( |\text{difference}| \) is small \( \rightarrow \) \( h \) is ok \( \rightarrow \) moves on

if \( |\text{difference}| > \text{tolerance} \), shrink \( h \), repeat

"adaptive step size"

  • straight portion: big \( h \) ok
  • windy portion: need small \( h \)
PAGE 8

Adaptive Stepping in ode45: \( y' = 1 - t + 4y \)

A line graph showing the solution y against time t for the differential equation y' = 1 - t + 4y. The vertical axis is labeled 'Solution (y)' and ranges from 0 to 4000 in increments of 500. The horizontal axis is labeled 'Time (t)' and ranges from 0 to 2 in increments of 0.2. The plot features a blue line connecting discrete points marked with red-filled circles, identified in the legend as 'Steps chosen by ode45'. The curve starts near (0,0) and exhibits exponential growth, reaching a value of approximately 3500 at t=2. The markers are most densely clustered at the start (t=0) and become more widely spaced as the time increases, illustrating the adaptive step size selection of the solver.
Figure: A line graph showing the solution y against time t for the differential equation y' = 1 - t + 4y. The vertical axis is labeled 'Solution (y)' and ranges from 0 to 4000 in increments of 500. The horizontal axis is labeled 'Time (t)' and ranges from 0 to 2 in increments of 0.2. The plot features a blue line connecting discrete points marked with red-filled circles, identified in the legend as 'Steps chosen by ode45'. The curve starts near (0,0) and exhibits exponential growth, reaching a value of approximately 3500 at t=2. The markers are most densely clustered at the start (t=0) and become more widely spaced as the time increases, illustrating the adaptive step size selection of the solver.
PAGE 9

ode45 reacting to a 'Sharp Event'

A line graph titled 'ode45 reacting to a Sharp Event' showing a numerical solution over time. The horizontal axis is labeled 'Time (t)' and ranges from 0 to 4. The vertical axis is labeled 'y' and ranges from -2 to 12. The plot shows a series of red data points connected by a blue line. From \( t=0 \) to approximately \( t=1.8 \), the value of \( y \) remains constant at 0. Near \( t=1.8 \), the density of data points increases dramatically, indicating a reduction in step size by the solver. The curve then rises sharply to a peak of approximately \( y=10.1 \) at \( t=2.1 \). After the peak, the curve follows a smooth exponential-like decay back towards \( y=0 \) as \( t \) approaches 4, with the data points becoming more widely spaced again.
Figure: A line graph titled 'ode45 reacting to a Sharp Event' showing a numerical solution over time. The horizontal axis is labeled 'Time (t)' and ranges from 0 to 4. The vertical axis is labeled 'y' and ranges from -2 to 12. The plot shows a series of red data points connected by a blue line. From \( t=0 \) to approximately \( t=1.8 \), the value of \( y \) remains constant at 0. Near \( t=1.8 \), the density of data points increases dramatically, indicating a reduction in step size by the solver. The curve then rises sharply to a peak of approximately \( y=10.1 \) at \( t=2.1 \). After the peak, the curve follows a smooth exponential-like decay back towards \( y=0 \) as \( t \) approaches 4, with the data points becoming more widely spaced again.