6.5. Runge–Kutta methods#
In this section, we discuss (explicit) Runge–Kutta methods which can take the form
\[\begin{align*}
y_{n+1} &= y_{n}+h\sum _{i=1}^{s}b_{i}kW_{i},\\
k_{i} &= f\left(t_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{i-1}a_{ij}k_{j}\right).
\end{align*}\]
This process is carried forward iteratively until point \(x_{n-1}\) is reached.
The Butcher tableau is often utilized to portray the coefficients of this method. This table puts the coefficients of the method in a table as follows.
\[\begin{align*}
\begin{array}{c|cccc}
c_1 & a_{11} & a_{12}& \dots & a_{1s}\\
c_2 & a_{21} & a_{22}& \dots & a_{2s}\\
\vdots & \vdots & \vdots& \ddots& \vdots\\
c_s & a_{s1} & a_{s2}& \dots & a_{ss} \\
\hline
& b_1 & b_2 & \dots & b_s\\
\end{array}
\end{align*}\]
Notes:
A Runge–Kutta method is consistent if and only if \begin{align*}\sum {i=1}^{s}b{i}=1.\end{align*}
Condition for determining coefficients:\begin{align*}\sum {j=1}^{i-1}a{ij}=c_{i},\quad i=2,\ldots ,s.\end{align*}
Remark
Most of the previous explicit methods can be obtained using the Butcher tableau.
6.5.1. First-order Runge–Kutta method#
For example, when \(s= 1\), we have,
(6.44)#\[\begin{align}
\begin{array}{c|c}
0 & 0 \\
\hline
& 1 \\
\end{array},
\end{align}\]
and
(6.45)#\[\begin{align}
y_{n+1} &= y_{n}+h k_{1},\\
k_{1} &= f\left(t_{n},y_{n}\right).
\end{align}\]
This gives the Euler method, which is first order. The Butcher tableau can be also generated using nodepy python package. To get the Butcher tableau for the Euler method, we have,
from nodepy import rk
def BT(Method):
print(rk.loadRKM()[Method].dj_reduce())
BT('FE')
Forward Euler
0 |
___|___
| 1
6.5.2. Second-order Runge–Kutta methods#
Similarly, the Explicit midpoint method, Heun’s method, and Ralston’s method can be obtained as follows.
BT('Heun22')
BT('Mid22')
BT('MTE22')
SSPRK 22
The optimal 2-stage, 2nd order SSP Runge-Kutta method, also known as Heun's 2nd order method
0 |
1 | 1
_____|__________
| 1/2 1/2
Midpoint Runge-Kutta
0 |
1/2 | 1/2
_____|__________
| 0 1
Minimal Truncation Error 22
0 |
2/3 | 2/3
_____|__________
| 1/4 3/4
Moreover, a generic second-order method can be defined as follows [7]
(6.46)#\[\begin{align}
\begin{array}{c|ccc}0&0&0\\\alpha &\alpha &0\\\hline &1-{\frac {1}{2\alpha }}&{\frac {1}{2\alpha }}\\\end{array}
\end{align}\]
where \(\alpha \neq 0\). Here,
(6.47)#\[\begin{align}
\begin{cases}
\alpha = 1/2, & \text{Explicit midpoint method},\\
\alpha = 2/3, & \text{Ralston's method},\\
\alpha = 1, & \text{Heun's method}.\\
\end{cases}
\end{align}\]
import pandas as pd
import numpy as np
def Runge_Kutta_2nd(f, y0, a, b, alpha = 1/2, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
alpha : float, optional
DESCRIPTION. The default is 1/2. The alpha from generic second-order method.
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
if alpha == 1/2:
Name = 'Explicit midpoint method'
elif alpha == 2/3:
Name = """Ralston's method"""
elif alpha == 1:
Name = """Heun's method"""
else:
Name = """Generic second-order method with $\\alpha = %.2f$""" % alpha
for i in range(N):
k1 = f(t[i], y[i]);
k2 = f(t[i]+ alpha*h, y[i]+ h*alpha*k1);
y[i+1] = y[i] + h*(1 - (1/(2*alpha)))*k1 + h*(1/(2*alpha))*k2
Table = pd.DataFrame({'t':t, 'y':y})
return Table, Name
function [Table,Name] = Runge_Kutta_2nd(f, y0, a, b, N, alpha)
%{
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
alpha : float, optional
DESCRIPTION. The default is 1/2. The alpha from generic second-order method.
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
switch nargin
case 5
alpha = 1/2;
end
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
switch alpha
case 1/2
Name = "Explicit midpoint method";
case 2/3
Name = "Ralston's method";
case 1
Name = "Heun's method";
otherwise
Name = sprintf("Generic second-order method with $\\alpha = %.2f$", alpha);
% loop
for i=1:N
k1 = f(t(i), y(i));
k2 = f(t(i)+ alpha*h, y(i)+ h*alpha*k1);
y(i+1) = y(i) + h*(1 - (1/(2*alpha)))*k1 + h*(1/(2*alpha))*k2;
end
Table = table(t,y);
end
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use second-order Runge–Kutta methods for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
import numpy as np
import pandas as pd
from hd_IVP_Algorithms import Runge_Kutta_2nd
# f(x,y(x)):
f = lambda t, y: t*np.exp(-t**(2))-2*t*y
(a, b) = (0, 1)
# the eact solution y(x)
y_exact = lambda t: (1+(t**2)/2)*np.exp(-t**2)
# Initial Value
y0 = 1
# Table
alph = 1/2
Table, Name = Runge_Kutta_2nd(alpha = alph, f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
y = Table.y
Table = Table.rename(columns = {'y':'y (%s)' % Name})
Table['Error (%s)' % Name] = np.abs(Table['Exact'] - Table['y (%s)' % Name])
for alph in [1, 2/3]:
TB, Name = Runge_Kutta_2nd(alpha = alph, f = f, y0 = y0, a = a, b = b, N = 10)
TB['Exact'] = y_exact(TB['t'])
TB = TB.rename(columns = {'y':'y (%s)' % Name})
TB['Error (%s)' % Name] = np.abs(TB['Exact'] - TB['y (%s)' % Name])
Table = Table.merge(TB, on = ['t', 'Exact'])
Table = Table.reindex(sorted(Table.columns), axis=1)
Cols = [x for x in Table.columns if 'Error' in x]
display(Table[1:].style.set_properties(subset= Cols,**{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'})\
.format(dict(zip(Cols, len(Cols)*['{:.4e}']))))
| Error (Explicit midpoint method) | Error (Heun's method) | Error (Ralston's method) | Exact | t | y (Explicit midpoint method) | y (Heun's method) | y (Ralston's method) | |
|---|---|---|---|---|---|---|---|---|
| 1 | 1.2567e-05 | 4.9834e-05 | 2.2256e-05 | 0.995000 | 0.100000 | 0.994988 | 0.994950 | 0.994978 |
| 2 | 5.1080e-05 | 1.4537e-04 | 7.7404e-05 | 0.980005 | 0.200000 | 0.979954 | 0.979860 | 0.979928 |
| 3 | 1.1676e-04 | 2.7501e-04 | 1.6271e-04 | 0.955058 | 0.300000 | 0.954941 | 0.954783 | 0.954895 |
| 4 | 2.0940e-04 | 4.2042e-04 | 2.7212e-04 | 0.920315 | 0.400000 | 0.920106 | 0.919895 | 0.920043 |
| 5 | 3.2511e-04 | 5.5836e-04 | 3.9535e-04 | 0.876151 | 0.500000 | 0.875826 | 0.875593 | 0.875756 |
| 6 | 4.5452e-04 | 6.6326e-04 | 5.1756e-04 | 0.823258 | 0.600000 | 0.822804 | 0.822595 | 0.822741 |
| 7 | 5.8244e-04 | 7.1061e-04 | 6.2023e-04 | 0.762720 | 0.700000 | 0.762137 | 0.762009 | 0.762100 |
| 8 | 6.8933e-04 | 6.8072e-04 | 6.8352e-04 | 0.696026 | 0.800000 | 0.695337 | 0.695345 | 0.695342 |
| 9 | 7.5452e-04 | 5.6222e-04 | 6.8956e-04 | 0.625026 | 0.900000 | 0.624271 | 0.624463 | 0.624336 |
| 10 | 7.6042e-04 | 3.5464e-04 | 6.2620e-04 | 0.551819 | 1.000000 | 0.551059 | 0.551465 | 0.551193 |
Cols = ['h', 'N']
h = [2**(-i) for i in range(3, 10)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
Methods = []
for alph in [1/2, 1, 2/3]:
for n in range(Table.shape[0]):
TB, Name = Runge_Kutta_2nd(alpha = alph, f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh (%s)' % Name] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
Methods.append(Name)
Cols = [x for x in Table.columns if 'Eh' in x]
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Cols, len(Cols)*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs =[Table[x].values for x in Table[Cols]],
labels = Methods,
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy',
legend_orientation = 'horizontal', ylim = [1.9, 2.25])
| h | N | Eh (Explicit midpoint method) | Eh (Heun's method) | Eh (Ralston's method) | |
|---|---|---|---|---|---|
| 0 | 0.125000 | 8 | 1.2295e-03 | 1.1751e-03 | 1.1300e-03 |
| 1 | 0.062500 | 16 | 2.8534e-04 | 2.5226e-04 | 2.5496e-04 |
| 2 | 0.031250 | 32 | 6.8695e-05 | 5.8166e-05 | 6.0823e-05 |
| 3 | 0.015625 | 64 | 1.6854e-05 | 1.3962e-05 | 1.4834e-05 |
| 4 | 0.007812 | 128 | 4.1743e-06 | 3.4164e-06 | 3.6632e-06 |
| 5 | 0.003906 | 256 | 1.0387e-06 | 8.4501e-07 | 9.1021e-07 |
| 6 | 0.001953 | 512 | 2.5908e-07 | 2.1011e-07 | 2.2686e-07 |
6.5.3. Third-order Runge–Kutta method#
This method is defined as follows,
\[\begin{align*}
y_{n+1} &= y_{n}+\frac{h}{6}\left( k_{1} + 4k_{2} + k_{3} \right),\quad \text{for }n = 0, 1, 2, 3, \ldots,
\end{align*}\]
with
\[\begin{align*}
\begin{cases}
k_{1} & = f\left(t_{n},~y_{n} \right),\\
k_{2} & = f\left(t_{n}+\frac{h}{2},~y_{n}+\frac{h}{2}k_{1} \right),\\
k_{3} & = f\left(t_{n}+h,~y_{n}+ h(2k_{2} -k_{1}) \right)
\end{cases}
\end{align*}\]
or using the Butcher tableau
\[\begin{align*}\begin{array}{c|ccc}0&0&0&0\\1/2&1/2&0&0\\1&-1&2&0\\\hline &1/6&2/3&1/6\\\end{array}\end{align*}\]
import pandas as pd
import numpy as np
def Runge_Kutta_3rd(f, y0, a, b, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/2, y[i]+ h*k1/2)
k3 = f(t[i]+ h, y[i] + h*(- k1 + 2*k2))
y[i+1] = y[i] + h*(k1 + 4*k2 + k3)/6
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
function [Table] = Runge_Kutta_3rd(f, y0, a, b, N)
%{
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
% loop
for i=1:N
k1 = f(t(i), y(i));
k2 = f(t(i)+ h/2, y(i)+ h*k1/2);
k3 = f(t(i)+ h, y(i) + h*(- k1 + 2*k2));
y(i+1) = y(i) + h*(k1 + 4*k2 + k3)/6;
end
Table = table(t,y);
end
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use third-order Runge–Kutta method for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
import numpy as np
import pandas as pd
from hd_IVP_Algorithms import Runge_Kutta_3rd
# f(t, y(t)):
f = lambda t, y: t*np.exp(-t**(2))-2*t*y
(a, b) = (0, 1)
# the eact solution y(t)
y_exact = lambda t: (1+(t**2)/2)*np.exp(-t**2)
#
y0 = 1
Table = Runge_Kutta_3rd(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.100000 | 0.995009 | 0.995000 | 8.4271e-06 |
| 2 | 0.200000 | 0.980022 | 0.980005 | 1.7039e-05 |
| 3 | 0.300000 | 0.955084 | 0.955058 | 2.5794e-05 |
| 4 | 0.400000 | 0.920350 | 0.920315 | 3.4545e-05 |
| 5 | 0.500000 | 0.876194 | 0.876151 | 4.2967e-05 |
| 6 | 0.600000 | 0.823309 | 0.823258 | 5.0491e-05 |
| 7 | 0.700000 | 0.762776 | 0.762720 | 5.6286e-05 |
| 8 | 0.800000 | 0.696085 | 0.696026 | 5.9311e-05 |
| 9 | 0.900000 | 0.625084 | 0.625026 | 5.8454e-05 |
| 10 | 1.000000 | 0.551872 | 0.551819 | 5.2752e-05 |
Cols = ['h', 'N', 'Eh']
h = [2**(-i) for i in range(3, 12)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = Runge_Kutta_3rd(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-1:], 3*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
labels = ['Third-order Runge–Kutta method'],
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy: %s' % 'Third-order Runge–Kutta method',
legend_orientation = 'horizontal', ylim = [2.9, 3.2])
| h | N | Eh | |
|---|---|---|---|
| 0 | 0.125000 | 8 | 1.1968e-04 |
| 1 | 0.062500 | 16 | 1.3832e-05 |
| 2 | 0.031250 | 32 | 1.6616e-06 |
| 3 | 0.015625 | 64 | 2.0371e-07 |
| 4 | 0.007812 | 128 | 2.5221e-08 |
| 5 | 0.003906 | 256 | 3.1375e-09 |
| 6 | 0.001953 | 512 | 3.9125e-10 |
| 7 | 0.000977 | 1024 | 4.8849e-11 |
| 8 | 0.000488 | 2048 | 6.1038e-12 |
6.5.4. Heun’s third-order method#
This method is defined as follows,
\[\begin{align*}
y_{n+1} &= y_{n}+\frac{h}{4}\left( k_{1} + 3k_{3} \right)
\end{align*}\]
with
\[\begin{align*}
\begin{cases}
k_{1} & = f\left(t_{n},~y_{n} \right),\\
k_{2} & = f\left(t_{n}+\frac{h}{3},~y_{n}+\frac{h}{3}k_{1} \right),\\
k_{3} & = f\left(t_{n}+\frac{2}{3}h,~y_{n}+ \frac{2h}{3}k_{2}\right)
\end{cases}
\end{align*}\]
or using the Butcher tableau
\[\begin{align*}
\begin{array}{c|ccc}0&0&0&0\\1/3&1/3&0&0\\2/3&0&2/3&0\\\hline &1/4&0&3/4\\\end{array}
\end{align*}\]
The Butcher tableau for the third-order Heun’s method:
BT('Heun33')
Heun RK 33
Heun's 3-stage, 3rd order
0 |
1/3 | 1/3
2/3 | 2/3
_____|_______________
| 1/4 0 3/4
import pandas as pd
import numpy as np
def Heun_Method_3rd(f, y0, a, b, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/3, y[i]+ h*k1/3)
k3 = f(t[i]+ 2*h/3, y[i] + h*2*k2/3)
y[i+1] = y[i] + h*(k1 + 3*k3)/4
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
function [Table] = Heun_Method_3rd(f, y0, a, b, N)
%{
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
% loop
for i=1:N
k1 = f(t(i), y(i));
k2 = f(t(i)+ h/3, y(i)+ h*k1/3);
k3 = f(t(i)+ 2*h/3, y(i) + h*2*k2/3);
y(i+1) = y(i) + h*(k1 + 3*k3)/4;
end
Table = table(t,y);
end
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use the Heun’s third-order method for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
from hd_IVP_Algorithms import Heun_Method_3rd
Table = Heun_Method_3rd(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.100000 | 0.995000 | 0.995000 | 8.9306e-09 |
| 2 | 0.200000 | 0.980005 | 0.980005 | 4.1287e-08 |
| 3 | 0.300000 | 0.955058 | 0.955058 | 2.5240e-08 |
| 4 | 0.400000 | 0.920316 | 0.920315 | 2.7315e-07 |
| 5 | 0.500000 | 0.876152 | 0.876151 | 1.2114e-06 |
| 6 | 0.600000 | 0.823261 | 0.823258 | 3.1339e-06 |
| 7 | 0.700000 | 0.762726 | 0.762720 | 6.1831e-06 |
| 8 | 0.800000 | 0.696036 | 0.696026 | 1.0131e-05 |
| 9 | 0.900000 | 0.625040 | 0.625026 | 1.4309e-05 |
| 10 | 1.000000 | 0.551837 | 0.551819 | 1.7677e-05 |
Cols = ['h', 'N', 'Eh']
h = [2**(-i) for i in range(3, 12)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = Heun_Method_3rd(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-1:], 3*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
labels = ["""Third-order Heun's method"""],
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy: %s' % """Third-order Heun's method""",
legend_orientation = 'horizontal', ylim = [2.9, 3.2])
| h | N | Eh | |
|---|---|---|---|
| 0 | 0.125000 | 8 | 3.6177e-05 |
| 1 | 0.062500 | 16 | 4.0146e-06 |
| 2 | 0.031250 | 32 | 4.7161e-07 |
| 3 | 0.015625 | 64 | 5.7116e-08 |
| 4 | 0.007812 | 128 | 7.0265e-09 |
| 5 | 0.003906 | 256 | 8.7130e-10 |
| 6 | 0.001953 | 512 | 1.0848e-10 |
| 7 | 0.000977 | 1024 | 1.3533e-11 |
| 8 | 0.000488 | 2048 | 1.6908e-12 |
6.5.5. Ralston’s third-order method#
This method is defined as follows,
\[\begin{align*}
y_{n+1} &= y_{n}+\frac{h}{9}\left( 2k_{1} + 3k_{2} + 4k_{3} \right)
\end{align*}\]
with
\[\begin{align*}
\begin{cases}
k_{1} & = f\left(t_{n},~y_{n} \right),\\
k_{2} & = f\left(t_{n}+\frac{h}{2},~y_{n}+\frac{h}{2}k_{1} \right),\\
k_{3} & = f\left(t_{n}+\frac{3}{4}h,~y_{0} + \frac{3h}{4}k_{2}\right)
\end{cases}
\end{align*}\]
or using the Butcher tableau
\[\begin{align*}
\begin{array}{c|ccc}0&0&0&0\\1/2&1/2&0&0\\3/4&0&3/4&0\\\hline &2/9&1/3&4/9\\\end{array}
\end{align*}\]
import pandas as pd
import numpy as np
def Ralston_Method_3rd(f, y0, a, b, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/2, y[i]+ h*k1/2)
k3 = f(t[i]+ 3*h/4, y[i] + h*3*k2/4)
y[i+1] = y[i] + h*(2*k1 + 3*k2 + 4*k3)/9
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
function [Table] = Ralston_Method_3rd(f, y0, a, b, N)
%{
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
% loop
for i=1:N
k1 = f(t(i), y(i));
k2 = f(t(i)+ h/2, y(i)+ h*k1/2);
k3 = f(t(i)+ 3*h/4, y(i) + h*3*k2/4);
y(i+1) = y(i) + h*(2*k1 + 3*k2 + 4*k3)/9;
end
Table = table(t,y);
end
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use the Ralston’s third-order method for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
from hd_IVP_Algorithms import Ralston_Method_3rd
Table = Ralston_Method_3rd(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.100000 | 0.995002 | 0.995000 | 2.1207e-06 |
| 2 | 0.200000 | 0.980010 | 0.980005 | 4.3623e-06 |
| 3 | 0.300000 | 0.955065 | 0.955058 | 6.7984e-06 |
| 4 | 0.400000 | 0.920325 | 0.920315 | 9.5125e-06 |
| 5 | 0.500000 | 0.876163 | 0.876151 | 1.2538e-05 |
| 6 | 0.600000 | 0.823274 | 0.823258 | 1.5762e-05 |
| 7 | 0.700000 | 0.762739 | 0.762720 | 1.8830e-05 |
| 8 | 0.800000 | 0.696047 | 0.696026 | 2.1097e-05 |
| 9 | 0.900000 | 0.625047 | 0.625026 | 2.1679e-05 |
| 10 | 1.000000 | 0.551839 | 0.551819 | 1.9596e-05 |
Cols = ['h', 'N', 'Eh']
h = [2**(-i) for i in range(3, 12)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = Ralston_Method_3rd(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-1:], 3*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
labels = ["""Third-order Ralston's method"""],
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy: %s' % """Third-order Ralston's method""",
legend_orientation = 'horizontal', ylim = [2.9, 3.2])
| h | N | Eh | |
|---|---|---|---|
| 0 | 0.125000 | 8 | 4.4643e-05 |
| 1 | 0.062500 | 16 | 4.9206e-06 |
| 2 | 0.031250 | 32 | 5.7903e-07 |
| 3 | 0.015625 | 64 | 7.0228e-08 |
| 4 | 0.007812 | 128 | 8.6469e-09 |
| 5 | 0.003906 | 256 | 1.0728e-09 |
| 6 | 0.001953 | 512 | 1.3361e-10 |
| 7 | 0.000977 | 1024 | 1.6670e-11 |
| 8 | 0.000488 | 2048 | 2.0803e-12 |
6.5.6. Strong-Stability preserving Runge-Kutta time-steppers (SSPRK3)#
This method is defined as follows,
\[\begin{align*}
y_{n+1} &= y_{n}+\frac{h}{9}\left( 2k_{1} + 3k_{2} + 4k_{3} \right)
\end{align*}\]
with
\[\begin{align*}
\begin{cases}
k_{1} & = f\left(t_{n},~y_{n} \right),\\
k_{2} & = f\left(t_{n}+h,~y_{n}+hk_{1} \right),\\
k_{3} & = f\left(t_{n}+\frac{h}{2},~y_{0} + \frac{h}{4}k_{1} + \frac{h}{4}k_{2}\right)
\end{cases}
\end{align*}\]
or using the Butcher tableau
\[\begin{align*}
\begin{array}{c|ccc}0&0&0&0\\1&1&0&0\\1/2&1/4&1/4&0\\\hline &1/6&1/6&2/3\\\end{array}
\end{align*}\]
The Butcher tableau for the SSPRK3:
BT('SSP33')
SSPRK 33
The optimal 3-stage, 3rd order SSP Runge-Kutta method
0 |
1 | 1
1/2 | 1/4 1/4
_____|_______________
| 1/6 1/6 2/3
| 0.291 0.291 0.417
import pandas as pd
import numpy as np
def SSPRK3(f, y0, a, b, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h, y[i]+ h*k1)
k3 = f(t[i]+ h/2, y[i] + h*k1/4 + h*k2/4)
y[i+1] = y[i] + h*(k1 + k2 + 4*k3)/6
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use the SSPRK3 method for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
from hd_IVP_Algorithms import SSPRK3
Table = SSPRK3(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.100000 | 0.994992 | 0.995000 | 8.1570e-06 |
| 2 | 0.200000 | 0.979990 | 0.980005 | 1.5254e-05 |
| 3 | 0.300000 | 0.955038 | 0.955058 | 2.0552e-05 |
| 4 | 0.400000 | 0.920292 | 0.920315 | 2.3713e-05 |
| 5 | 0.500000 | 0.876126 | 0.876151 | 2.4931e-05 |
| 6 | 0.600000 | 0.823233 | 0.823258 | 2.4954e-05 |
| 7 | 0.700000 | 0.762695 | 0.762720 | 2.5002e-05 |
| 8 | 0.800000 | 0.695999 | 0.696026 | 2.6573e-05 |
| 9 | 0.900000 | 0.624994 | 0.625026 | 3.1160e-05 |
| 10 | 1.000000 | 0.551779 | 0.551819 | 3.9928e-05 |
Cols = ['h', 'N', 'Eh']
h = [2**(-i) for i in range(3, 12)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = SSPRK3(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-1:], 3*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
labels = ['SSPRK3'],
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy: %s' % 'SSPRK3',
legend_orientation = 'horizontal', ylim = [2.9, 3.1])
| h | N | Eh | |
|---|---|---|---|
| 0 | 0.125000 | 8 | 7.6429e-05 |
| 1 | 0.062500 | 16 | 9.9880e-06 |
| 2 | 0.031250 | 32 | 1.2680e-06 |
| 3 | 0.015625 | 64 | 1.5951e-07 |
| 4 | 0.007812 | 128 | 1.9995e-08 |
| 5 | 0.003906 | 256 | 2.5028e-09 |
| 6 | 0.001953 | 512 | 3.1305e-10 |
| 7 | 0.000977 | 1024 | 3.9144e-11 |
| 8 | 0.000488 | 2048 | 4.8939e-12 |
6.5.7. Fourth-order Runge–Kutta method#
This method is defined as follows,
\[\begin{align*}
y_{n+1}=y_{n}+{\dfrac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\quad \text{for }n = 0, 1, 2, 3, \ldots,
\end{align*}\]
with
\[\begin{align*}
\begin{cases}
k_{1}&=h\ f(t_{n},y_{n}),\\k_{2}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{1}}{2}}\right),\\k_{3}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{2}}{2}}\right),\\k_{4}&=h\ f\left(t_{n}+h,y_{n}+k_{3}\right).
\end{cases}
\end{align*}\]
or using the Butcher tableau
\[\begin{align*}
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0\\
1/2 & 1/2 & 0 & 0 & 0\\
1/2 & 0 & 1/2 & 0 & 0\\
1 & 0 & 0 & 1 & 0\\
\hline
& 1/6 & 1/3 & 1/3 & 1/6\\
\end{array}
\end{align*}\]
import pandas as pd
import numpy as np
def Runge_Kutta_4rd(f, y0, a, b, h= False, N=False):
'''
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
h : float, optional
DESCRIPTION. The default is False. stepsize
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i]);
k2 = f(t[i]+(h/2), y[i]+(h/2)*k1);
k3 = f(t[i]+(h/2), y[i]+(h/2)*k2);
k4 = f(t[i+1], y[i]+h*k3);
y[i+1] = y[i] + h*(k1 + 2*k2 + 2*k3 + k4)/6
del k1, k2, k3, k4
Table = pd.DataFrame({'t':t, 'y':y})
return Table
function [Table] = Runge_Kutta_4rd(f, y0, a, b, N)
%{
Parameters
----------
f : function
DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
the initial value.
a : float
DESCRIPTION. a is the left side of interval [a, b]
b : float
DESCRIPTION. b is the right side of interval [a, b]
N : int, optional
DESCRIPTION. The default is False. number of points.
Returns
-------
Table : dataframe
DESCRIPTION. a summary of the algorithm output
Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
% loop
for i=1:N
k1 = f(t(i), y(i));
k2 = f(t(i)+(h/2), y(i)+(h/2)*k1);
k3 = f(t(i)+(h/2), y(i)+(h/2)*k2);
k4 = f(t(i+1), y(i)+h*k3);
y(i+1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end
Table = table(t,y);
end
Example: Consider the following IVP,
\[\begin{align*}
\begin{cases}
y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\
y(0) = 1,
\end{cases}
\end{align*}\]
with exact solution
\[\begin{align*}
y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}.
\end{align*}\]
Use the fourth-order Runge–Kutta method method for solving this IVP. Also, investigate the order of convergence numerically.
Solution:
from hd_IVP_Algorithms import Runge_Kutta_4rd
Table = Runge_Kutta_4rd(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.100000 | 0.995000 | 0.995000 | 1.0573e-08 |
| 2 | 0.200000 | 0.980005 | 0.980005 | 4.3131e-08 |
| 3 | 0.300000 | 0.955058 | 0.955058 | 9.7938e-08 |
| 4 | 0.400000 | 0.920315 | 0.920315 | 1.7346e-07 |
| 5 | 0.500000 | 0.876151 | 0.876151 | 2.6447e-07 |
| 6 | 0.600000 | 0.823258 | 0.823258 | 3.5890e-07 |
| 7 | 0.700000 | 0.762719 | 0.762720 | 4.3404e-07 |
| 8 | 0.800000 | 0.696026 | 0.696026 | 4.5380e-07 |
| 9 | 0.900000 | 0.625025 | 0.625026 | 3.6887e-07 |
| 10 | 1.000000 | 0.551819 | 0.551819 | 1.2183e-07 |
Cols = ['h', 'N', 'Eh']
h = [2**(-i) for i in range(3, 12)]
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = Runge_Kutta_4rd(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-1:], 3*["{:.4e}"]))))
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
labels = ['Fourth-order Runge–Kutta method'],
xlabel = r"$$i$$",
ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of Accuracy: %s' % 'Fourth-order Runge–Kutta method',
legend_orientation = 'horizontal', ylim = [3.5, 4.5])
| h | N | Eh | |
|---|---|---|---|
| 0 | 0.125000 | 8 | 1.1869e-06 |
| 1 | 0.062500 | 16 | 6.2114e-08 |
| 2 | 0.031250 | 32 | 3.5325e-09 |
| 3 | 0.015625 | 64 | 2.1097e-10 |
| 4 | 0.007812 | 128 | 1.2885e-11 |
| 5 | 0.003906 | 256 | 7.9614e-13 |
| 6 | 0.001953 | 512 | 4.9294e-14 |
| 7 | 0.000977 | 1024 | 3.1086e-15 |
| 8 | 0.000488 | 2048 | 1.4433e-15 |