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)
\begin{eqnarray}
dx&=&cos(t)dt\\
\int{dx}&=&\int{cos(t)}dt\\
x&=&sin(t)+x(t=0)
\end{eqnarray}
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
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()
Δt=0.01
Δt=0.001
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).
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()
\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) $
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()
Recommended Posts