수치적으로 방정식을 풀기(Runge Kunta 메소드)
Welcome to the Numerical-Analysis wiki!
How to we solve the complicated equations without boring hand writing calculation? especially differential equation.
Here is some solution!
This wiki based on professor Ha’s lecture note for computational physics class at Chung ang University.
Computer is just a tool.. we are learning physics!
1. Euler method
Sprit of method
Let’s think about Simple Harmonic Oscillations examples for Euler method(also called 1st order Runge Kunta method).
$\ddot{x}+w^2x=0 \to \ \frac{v_f-v_i}{\Delta{t}}+w^2x=0 \to v_f=v_i - x w^2\times\Delta{t}$
$G \equiv -x w^2$
we define G as above. then we can write the equation like this.
$v_f=v_i+ G\times\Delta{t}$
$x_f=x_i+v_i\times \Delta{t}$
This is just 1st order polynomial equation and G is a slope! if we find what the G is, then we can find velocity after $\Delta{t}$.
and this velocity determined position after $\Delta{t}$. and by using this velocity and position we can find next G(in case of us, $G \equiv -x w^2$).
and if we know what the G is, then we can find velocity after $\Delta{t}$ and this velocity determined position after $\Delta{t}$..
and.. and.. and..
Repeat this process and update the value(in this process, store each value in the list). That is the key concept of the Euler method. Eventually we got a numerical solution!
note that Euler method is approximation. more small $\Delta{t}$ makes more accurate result.
Here is the example. We can briefly write down a python code.
import matplotlib.pyplot as plt
import numpy as np
from math import *
#setting
delta_t=0.001
t_min=0
t_max=30
time=np.arange(t_min,t_max,delta_t)
sol_x=[] #array of solution of positon
sol_v=[] #array of solution of velocity
#initial condition
f=0.0
m=2.0
k=1.0
c=0.0
x=1.0
v=0.0
omega=sqrt(k/m)
#update with time interval delta t
for t in time:
#define slope
G=-x*omega*omega
#equation of motion
v=v+G*delta_t
x=x+v*delta_t
#store each value
sol_x.append(x)
#plot graph
plt.plot(time,sol_x,'r',label="position")
plt.legend()
plt.grid(True)
plt.show()
There is another example for driving oscillations.
$m\ddot{x} + c\dot{x} + kx = F(t)$
this equation of motion can be expressed of matrix form
$A\dot{y}+By=F$
$y_f=y_i+A^{-1}(F-By_i)\Delta{t}$
$G\equiv A^{-1}(F-By_i)$
And then, we can write the equations like this.
$y_f=y_i+G\times\Delta{t}$
this is the python code that implement above equations.
import matplotlib.pyplot as plt
import numpy as np
from math import *
#setting
delta_t=0.001
t_min=0
t_max=30
time=np.arange(t_min,t_max,delta_t)
sol_x=[] #array of solution of positon
sol_v=[] #array of solution of velocity
#initial condition
f=0.0
m=2.0
k=1.0
c=0.0
x=1.0
v=0.0
omega=sqrt(k/m)
A=np.array([[m,0],[0,1]])
B=np.array([[c,k],[-1,0]])
F=np.array([f,0])
y=np.array([v,x])
#print(A)
#update with time interval delta t
for t in time:
G=np.matmul(np.linalg.inv(A),(F-np.matmul(B,y)))
y=y+G*delta_t
sol_x.append(y[1])
sol_v.append(y[0])
name="Driven Oscillation"
plt.plot(time, sol_x,'r',label="position")
plt.grid(True)
plt.legend()
plt.title(name)
plt.show()
2. Runge Kunta method
We will see that Runge Kunta method(also called 4th order RK method) is more accurate than 1st RK method.
$y_f=y_i+G\times\Delta{t}$
$G \equiv \frac{1}{6}k1+\frac{2}{6}k2+\frac{2}{6}k3+\frac{1}{6}k4$
k1~k4 is the slope in the time interval.
G can be defined as different coefficient of each k.
import matplotlib.pyplot as plt
import numpy as np
from math import *
#setting
delta_t=0.001
t_min=0
t_max=100
time=np.arange(t_min,t_max,delta_t)
sol_x=[]
sol_v=[]
#initial condition
f=0.0
m=2.0
k=1.0
c=0.0
x=1.0
v=0.0
omega=sqrt(k/m)
A=np.array([[m,0],[0,1]])
B=np.array([[c,k],[-1,0]])
F=np.array([0,0])
y=np.array([v,x])
A_inverse=np.linalg.inv(A)
def G(t,y):
F[0]=f*np.cos(omega*t)
return A_inverse.dot(F-B.dot(y))
def RK4(t,y,delta_t):
k1=G(t,y)
k2=G(t+0.5*delta_t,y+0.5*delta_t*k1)
k3=G(t+0.5*delta_t,y+0.5*delta_t*k2)
k4=G(t+0.5*delta_t,y+0.5*delta_t*k3)
return 1./6.0*k1 + 2./6*k2+2./6*k3+1./6*k4
#return k1
#print(A)
for t in time:
y=y+delta_t*RK4(t,y,delta_t)
sol_x.append(y[1])
sol_v.append(y[0])
#G=RK4(t,y,delta_t)
name="Driven Oscillation"
Force=f
#label=f'Force:{Force}'
plt.plot(time, sol_x,'r',label="position")
#plt.plot(time, sol_v,'--',label="velocitiy")
plt.grid(True)
plt.legend()
plt.title(name)
plt.show()
