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


t.png


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

t1.png

t2.png


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

t1.png t2.png


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

t1.png t2.png


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

Recommended Posts

[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics
Solving one-dimensional wave equation by finite difference method (python)
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Solving the Schrodinger equation in the steady state in the 3D isotropic harmonic oscillator potential by the matrix method, boundary value problem, quantum mechanics
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Inverse matrix calculation, numpy
[Scientific / technical calculation by Python] Histogram, visualization, matplotlib
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
[Scientific / technical calculation by Python] Logarithmic graph, visualization, matplotlib
[Scientific / technical calculation by Python] Polar graph, visualization, matplotlib
[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Scientific / technical calculation by Python] 3rd order spline interpolation, scipy
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Vector field visualization example, electrostatic field, matplotlib
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Calculation of matrix product by @ operator, python3.5 or later, numpy
[Scientific / technical calculation by Python] Drawing animation of parabolic motion with locus, matplotlib
[Scientific / technical calculation by Python] Plot, visualize, matplotlib 2D data with error bars
[Scientific / technical calculation by Python] Drawing of 3D curved surface, surface, wireframe, visualization, matplotlib
Solving the equation of motion in Python (odeint)
[Scientific / technical calculation by Python] Plot, visualization, matplotlib of 2D data read from file
[Scientific / technical calculation by Python] Drawing, visualization, matplotlib of 2D (color) contour lines, etc.
[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy
[Scientific and technical calculation by Python] Drawing of fractal figures [Sierpinski triangle, Bernsley fern, fractal tree]
Python --Differential equation numerical solution Euler method & central difference method & Runge-Kutta method
[Python] LASSO regression with equation constraints using the multiplier method
Scientific / technical calculation by Python] Drawing and visualization of 3D isosurface and its cross section using mayavi
[Science / technical calculation by Python] Comparison of convergence speeds of SOR method, Gauss-Seidel method, and Jacobi method for Laplace equation, partial differential equation, boundary value problem