[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method

Introduction

** 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. ** **

Contents

(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

Contents (1): Solving the first-order ordinary differential equation

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()

Contents (2): Harmonic oscillator

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()

Contents (3) Damped vibration

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()

Contents (4) Forced vibration

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()

Addendum

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).