EM
import numpy as np
import matplotlib.pyplot as plt
def f(r, t, mkc=[0.5, 50, 10]):
x=r[0]; dx=r[1]
f_x=dx
f_dx=-(mkc[2]*f_x+mkc[1]*x)/mkc[0]
return np.array([f_x, f_dx])
if __name__=="__main__":
mkc=[[0.5, 50, 8], [0.5, 50, 10], [0.5, 50, 12]]
t=np.linspace(0, 3, 1000)
h=t[1]-t[0]
xx=[]
for _mkc in mkc:
x=np.zeros(len(t))
y=np.zeros(len(t))
# INITIAL STATE
x[0]=10.
y[0]=0
for i in range(len(t)-1): # 4-th order RKM
r=[x[i], y[i]]
x[i+1], y[i+1]=r+h*f(r, t[i], _mkc)
xx.append(x)
plt.plot(t, xx[0], label="under damping")
plt.plot(t, xx[1], label="critical damping")
plt.plot(t, xx[2], label="over damping")
plt.title("Euler's Method")
plt.legend()
plt.show()
2nd-RKM
import numpy as np
import matplotlib.pyplot as plt
def f(r, t, mkc=[0.5, 50, 10]):
x=r[0]; dx=r[1]
f_x=dx
f_dx=-(mkc[2]*f_x+mkc[1]*x)/mkc[0]
return np.array([f_x, f_dx])
if __name__=="__main__":
mkc=[[0.5, 50, 8], [0.5, 50, 10], [0.5, 50, 12]]
t=np.linspace(0, 3, 1000)
h=t[1]-t[0]
xx=[]
for _mkc in mkc:
x=np.zeros(len(t))
y=np.zeros(len(t))
# INITIAL STATE
x[0]=10.
y[0]=0
for i in range(len(t)-1): # 4-th order RKM
r=[x[i], y[i]]
k_1=h*f(r, t[i], _mkc)
k_2=h*f(r+k_1/2, t[i]+h/2, _mkc)
x[i+1], y[i+1]=r+k_2
xx.append(x)
plt.plot(t, xx[0], label="under damping")
plt.plot(t, xx[1], label="critical damping")
plt.plot(t, xx[2], label="over damping")
plt.title("2nd-order RK")
plt.legend()
plt.show()
4th-RKM
import numpy as np
import matplotlib.pyplot as plt
def f(r, t, mkc=[0.5, 50, 10]):
x=r[0]; dx=r[1]
f_x=dx
f_dx=-(mkc[2]*f_x+mkc[1]*x)/mkc[0]
return np.array([f_x, f_dx])
if __name__=="__main__":
mkc=[[0.5, 50, 8], [0.5, 50, 10], [0.5, 50, 12]]
t=np.linspace(0, 3, 1000)
h=t[1]-t[0]
xx=[]
for _mkc in mkc:
x=np.zeros(len(t))
y=np.zeros(len(t))
# INITIAL STATE
x[0]=10.
y[0]=0
for i in range(len(t)-1): # 4-th order RKM
r=[x[i], y[i]]
k_1=h*f(r, t[i], _mkc)
k_2=h*f(r+k_1/2, t[i]+h/2, _mkc)
k_3=h*f(r+k_2/2, t[i]+h/2, _mkc)
k_4=h*f(r+k_3, t[i]+h, _mkc)
x[i+1], y[i+1]=r+(k_1+2*k_2+2*k_3+k_4)/6
xx.append(x)
plt.plot(t, xx[0], label="under damping")
plt.plot(t, xx[1], label="critical damping")
plt.plot(t, xx[2], label="over damping")
plt.title("4th-order RK")
plt.legend()
plt.show()