** The 4th order Runge-Kutta method, which is one of the solutions to ordinary differential equations. An example of the numerical solution of the Newton equation by% B2% EF% BC% 9D% E3% 82% AF% E3% 83% 83% E3% 82% BF% E6% B3% 95) is given. ** **
(1) [Warming up] Solving the first-order ordinary differential equation by the 4th-order Runge-Kutta method. (2) Harmonic oscillator (3) Damped vibration (4) Forced vibration
import numpy as np
import matplotlib.pyplot as plt
"""
4th runge-First-order ordinary differential equation by Kutta method
The fourth-order Runge-Kutta method for the 1st order ordinary differencial equation
dx/dt = -x^3+sin(t)Solve
"""
def f(x,t):
return -x**3 +np.sin(t)
a = 0.0
b = 10.0
N = 100
h = (b-a)/N
tpoints = np.arange(a,b,h)
xpoints = []
x = 0.0
for t in tpoints:
xpoints.append(x)
k1 = h*f(x,t)
k2 = h*f(x+0.5*k1, t+0.5*h)
k3 = h*f (x+0.5*k2, t+0.5*h)
k4 = h*f(x+k3, t+h)
x += (k1+2*k2+2*k3+k4)/6
plt.plot (tpoints, xpoints)
plt.xlabel("t")
plt.ylabel("x(t)")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
"""
4th runge-Solving a one-dimensional Newton equation by the Kutta method
The fourth-order Runge-Kutta method for the 1D Newton's equation
#Example:Linear restoring force(harmonic oscilation)
"""
def f(x, v, t):
return -k*x #Restoring force
M = 1.0 # mass: 1 Kg
k = 1.0 #Spring constant
t0 = 0.0
t1 = 10.0
N = 50
del_t = (t1-t0)/N # time grid
tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []
# initial condition
x0 = 0.0
v0 = 1.0
x, v = x0, v0
for t in tpoints:
xpoints.append(x)
vpoints.append(v)
k1v =f(x, v, t)*del_t /M
k1x = v * del_t
k2v = f(x+k1x/2, v+k1v/2, t+del_t/2 )*del_t /M
k2x =(v+k1v/2)*del_t
k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
k3x =(v+k2v/2 ) *del_t
k4v = f(x+k3x, v+k3v, t+del_t )*del_t /M
k4x = (v+k3v )*del_t
v += (k1v + 2 * k2v + 2* k3v + k4v)/6
x += (k1x + 2 * k2x + 2* k3x + k4x)/6
plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("x(t)", fontsize=24)
tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper left')
plt.show()
plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)
tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper left')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
"""
Damped vibration
Damped oscillation
"""
def f(x, v, t):
k=1.0
a=1.0
return -k*x-a*v
M = 1.0 # mass: 1 Kg
t0 = 0.0
t1 = 10.0
N = 50
del_t = (t1-t0)/N # time grid
tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []
# initial condition
x0 = 0.0
v0 = 1.0
x, v = x0, v0
for t in tpoints:
xpoints.append(x)
vpoints.append(v)
k1v =f(x, v, t)*del_t /M
k1x = v * del_t
k2v = f(x+k1/2, v+k1v/2, t+del_t/2 )*del_t /M
k2x =(v+k1v/2)*del_t
k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
k3x =(v+k2v/2 ) *del_t
k4v = f(x+k3, v+k3v, t+del_t )*del_t /M
k4x = (v+k3v )*del_t
v += (k1v + 2 * k2v + 2* k3v + k4v)/6
x += (k1x + 2 * k2x + 2* k3x + k4x)/6
plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("x(t)", fontsize=24)
#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper right')
plt.show()
plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)
#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper right')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
"""
Forced vibration
Forced vibration
"""
def f(x, v, t):
k=1.0
gamma=0.1
return -k*x-2* gamma *v + 2*sin(t)
M = 1.0 # mass: 1 Kg
t0 = 0.0
t1 = 50.0
N = 1000
del_t = (t1-t0)/N # time grid
tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []
# initial condition
x0 = 0.0
v0 = 1.0
x, v = x0, v0
for t in tpoints:
xpoints.append(x)
vpoints.append(v)
k1v =f(x, v, t)*del_t /M
k1x = v * del_t
k2v = f(x+k1/2, v+k1v/2, t+del_t/2 )*del_t /M
k2x =(v+k1v/2)*del_t
k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
k3x =(v+k2v/2 ) *del_t
k4v = f(x+k3, v+k3v, t+del_t )*del_t /M
k4x = (v+k3v )*del_t
v += (k1v + 2 * k2v + 2* k3v + k4v)/6
x += (k1x + 2 * k2x + 2* k3x + k4x)/6
plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("x(t)", fontsize=24)
#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper right')
plt.show()
plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)
#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper right')
plt.show()
In the Runge-Kutta method performed here, the total energy, which is one of the conserved quantities of the system, is not numerically conserved as the time step increases. One of the solutions to improve this (make it symplectic) is the Berley method (see Qiita article).