Python --Differential equation numerical solution Euler method & central difference method & Runge-Kutta method

The Euler method, the central difference method, and the Runge-Kutta method are used to calculate the solution of the following ordinary differential equations.

\frac{dx}{dt}=f(t)=cos(t)

General solution

\begin{eqnarray}
dx&=&cos(t)dt\\
\int{dx}&=&\int{cos(t)}dt\\
x&=&sin(t)+x(t=0)
\end{eqnarray}

Euler method

Taylor expansion of $ x (t + Δt) $

x(t+Δt)=x(t)+\frac{dx(t)}{dt}Δt+O(Δt^2)

If $ O (Δt ^ 2) $ is truncated as a truncation error

\begin{eqnarray}
x(t+Δt)&≃&x(t)+\frac{dx(t)}{dt}Δt\\
&=&x(t)+cos(t)Δt...①
\end{eqnarray}

Then, the method of sequentially calculating $ x $ at $ t + Δt $ using $ x $ at $ t $.

The algorithm is

1. Give initial values to $ t and x $.
2. Find $ x $ at $ t + Δt $ by equation (1).
3. Repeat range 2 you want to calculate.

·Source code

import math
import matplotlib.pyplot as plt
import numpy as np

def f(t):
    return math.cos(t)

DELTA_T = 0.001
MAX_T = 100.0

t = 0.0 #t initial value
x = 0.0 # t=X at 0

x_hist = [x]
t_hist = [t]

#Sequential calculation
while t < MAX_T:
    x += f(t)*DELTA_T
    t += DELTA_T
    x_hist.append(x)
    t_hist.append(t)

#Numerical solution plot
plt.plot(t_hist, x_hist)

#Exact solution(sin(t))Plot
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)

plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)

plt.show()

·result

Δt=0.01

Δt=0.001

Central difference method

Taylor expansion of $ x (t + Δt) $ and $ x (t-Δt) $

\begin{eqnarray}
x(t+Δt)&=&x(t)+\frac{dx(t)}{dt}Δt+\frac{1}{2}\frac{d^2x}{dt^2}Δt^2+O(Δt^3)...①\\
x(t-Δt)&=&x(t)-\frac{dx(t)}{dt}Δt+\frac{1}{2}\frac{d^2x}{dt^2}Δt^2+O(Δt^3)...②
\end{eqnarray}

①-② and truncating $ O (Δt ^ 3) $

\begin{eqnarray}
x(t+Δt)-x(t-Δt)&≃&2\frac{dx(t)}{dt}Δt\\
x(t+Δt)&≃&x(t-Δt)+2cos(t)Δt...③\\
\end{eqnarray}

If you replace it with $ t + Δt $ ⇒ $ t $

\begin{eqnarray}
x(t)&≃&x(t-2Δt)+2cos(t-Δt)Δt...③\\
\end{eqnarray}

When $ x $ at $ t-2Δt $ is known, the method of sequentially calculating $ x $ at $ t $ by equation (3).

·Source code

import math
import matplotlib.pyplot as plt
import numpy as np

def f(t):
    return math.cos(t)

DELTA_T = 0.001
MAX_T = 100.0

t = 0.0 #t initial value
x = 0.0 # t=X at 0

x_hist = [x]
t_hist = [t]

while t < MAX_T:
    x += 2*f(t-DELTA_T)*DELTA_T
    t += 2*DELTA_T
    x_hist.append(x)
    t_hist.append(t)

#Numerical solution plot
plt.plot(t_hist, x_hist)

#Exact solution(sin(t))Plot
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)

plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)

plt.show()

·result

Δt=0.01

Δt=0.001

4th order Runge-Kutta method

\begin{eqnarray}
k_1&=&Δtf(t,x)\\
k_2&=&Δtf(t+\frac{Δt}{2}, x(t)+\frac{k_1}{2})\\
k_3&=&Δtf(t+\frac{Δt}{2}, x(t)+\frac{k_2}{2})\\
k_4&=&Δtf(t+Δt, x(t)+k_3)\\
x(t+Δt)&=&x(t)+\frac{1}{6}(k_1+2k_2+2k_3+k_4)
\end{eqnarray}

By, $ x $ at $ t + Δt $ is calculated sequentially using $ x $ at $ t $.

When $ f (t, x) = f (t) $, from $ k_2 = k_3 $

\begin{eqnarray}
x(t+Δt)&=&x(t)+\frac{1}{6}(k_1+4k_2+k_4)
\end{eqnarray}

The error contained in $ x (t + Δt) $ is $ O (Δt ^ 5) $

·Source code

import math
import matplotlib.pyplot as plt
import numpy as np

def f(t):
    return math.cos(t)

DELTA_T = 0.001
MAX_T = 100.0

t = 0.0 #t initial value
x = 0.0 # t=X at 0

x_hist = [x]
t_hist = [t]

#Sequential calculation
while t < MAX_T:
    k1 = DELTA_T*f(t)
    k2 = DELTA_T*f(t+DELTA_T/2)
    k3 = DELTA_T*f(t+DELTA_T/2)
    k4 = DELTA_T*f(t+DELTA_T)
    x += (k1 + 2*k2 + 2*k3 + k4)/6
    t += DELTA_T
    x_hist.append(x)
    t_hist.append(t)
    
#Numerical solution plot
plt.plot(t_hist, x_hist)

#Exact solution(sin(t))Plot
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)

plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)

plt.show()

·result

Δt=0.001

Δt=0.001