[Wissenschaftlich-technische Berechnung nach Python] Lösen der eindimensionalen Newton-Gleichung nach der Runge-Kutta-Methode 4. Ordnung
Einführung
** Die [Lunge-Kutta-Methode] vierter Ordnung (https://ja.wikipedia.org/wiki/%E3%83%AB%E3%83%B3%E3%82), die eine der Lösungen für gewöhnliche Differentialgleichungen darstellt. Ein Beispiel für die numerische Lösung der Newtonschen Gleichung durch% B2% EF% BC% 9D% E3% 82% AF% E3% 83% 83% E3% 82% BF% E6% B3% 95) wird gegeben. ** ** **
Inhalt
(1) [Aufwärmen] Lösen der gewöhnlichen Differentialgleichung erster Ordnung nach der Runge-Kutta-Methode 4. Ordnung. (2) Oberschwingungsoszillator (3) Gedämpfte Vibration (4) Erzwungene Vibration
Inhalt (1): Lösen der gewöhnlichen Differentialgleichung erster Ordnung
import numpy as np
import matplotlib.pyplot as plt
"""
4. Runge-Gewöhnliche Differentialgleichung erster Ordnung nach der Kutta-Methode
The fourth-order Runge-Kutta method for the 1st order ordinary differencial equation
dx/dt = -x^3+sin(t)Lösen
"""
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()
Inhalt (2): Harmonischer Oszillator
import numpy as np
import matplotlib.pyplot as plt
"""
4. Runge-Lösen der eindimensionalen Newtonschen Gleichung nach der Kutta-Methode
The fourth-order Runge-Kutta method for the 1D Newton's equation
#Beispiel:Lineare Belastbarkeit(harmonic oscilation)
"""
def f(x, v, t):
return -k*x #Strom wiederherstellen
M = 1.0 # mass: 1 Kg
k = 1.0 #Federkonstante
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()
Inhalt (3) Dämpfungsvibration
import numpy as np
import matplotlib.pyplot as plt
"""
Gedämpfte 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()
Inhalt (4) Erzwungene Vibration
import numpy as np
import matplotlib.pyplot as plt
"""
Erzwungene 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()
Nachtrag
Bei der hier durchgeführten Lunge-Kutta-Methode wird die Gesamtenergie, die eine der konservierten Mengen des Systems ist, mit zunehmendem Zeitschritt nicht numerisch konserviert. Eine der Lösungen, um dies zu verbessern (symmetrisch zu machen), ist die Berley-Methode (siehe Qiita-Artikel).
Recommended Posts