[전산물리학] 이차 미분방정식의 수치적분1

(출처: 전산물리학(김범준, 박해진, 손승우, 이미진, 정우성))

 

고2 때(4년 전이라니 충격..) 배운 물리학I을 떠올려보자.  모든 건 F=ma에서 시작했다.

 

 

수정-오일러 방법( The Error-corrected Euler's Method  인듯)

import numpy as np, matplotlib.pyplot as plt
g = 9.8; v_init = 10; theta = np.pi/4.0
dt = 0.01; maxt = 10.0
vx = v_init*np.cos(theta)
vy = v_init*np.sin(theta)
x, y = 0.0, 0.0
t_arr = []; x_arr = []; y_arr = []
for t in np.arange(0, maxt, dt):
    t_arr.append(t); x_arr.append(x); y_arr.append(y)
    x += dt*vx
    vytem = vy + dt* (-g)
    y += 0.5*dt * (vy + vytem)
    vy = vytem
    if (y<0.0) : break #물체가 지면에 닿으면 실행 중단
plt.plot(x_arr,y_arr)
plt.xlabel("x")
plt.ylabel("y")
plt.show()

 

 

뛰어넘기 방법(leap-frog)

import numpy as np, matplotlib.pyplot as plt
g = 9.8; v_init = 10; theta = np.pi/4.0
dt = 0.01; maxt = 10.0
vx = v_init*np.cos(theta)
vy = v_init*np.sin(theta)
x, y = 0.0, 0.0
t_arr = []; x_arr = []; y_arr = []
for t in np.arange(0, maxt, dt):
    t_arr.append(t); x_arr.append(x); y_arr.append(y)
    #여기서부터 다름
    x += dt*vx
    y += dt*vy
    vy += dt*(-g)
    if (y<0.0) : break
plt.plot(x_arr,y_arr)
plt.xlabel("x")
plt.ylabel("y")
plt.show()