An ordinary differential equation is an equation that includes the independent variables t, the function y (t) that depends on t, and its nth derivative (n = 0,1,2, ..., N). In other words
g(t,y,y',y'',y^{(3)},...,y^{(N)}) = 0 \hspace{50pt}(1)
An equation that can be described in the form of.
Such equations can be approximated by a numerical solution called the Eular method (Euler method). I will explain using an actual example.
The differential equation $ (1) $ above is actually in a form that is not suitable for using the Euler method. You have to somehow transform this into something like $ (2) $ below.
\frac{dx(t)}{dt} = f(t,x(t)) \hspace{50pt}(2)
\\However, x(t)=\begin{pmatrix}y\\y'\\y''\\...\\y^{(N-1)}\end{pmatrix}
Consider a system with a mass $ m $ at the tip of a spring with a spring constant of $ k $. An oscillating external force is applied to the system, and the differential equation when considering the air resistance is as follows.
m\frac{d^2y(t)}{dt^2}+2\gamma \frac{dy(t)}{dt}+ky(t)=a\sin(\omega t)\hspace{50pt}(3)
When $ (3) $ is transformed,
\frac{d^2y(t)}{dt^2}=-\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t) \hspace{50pt}(3-1)
Also,
\frac{dy(t)}{dt}=\frac{dy(t)}{dt}\hspace{50pt}(3-2)
Where the vector $ x (t) $
x(t)=\begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} = \begin{pmatrix} y(t) \\ \frac{dy(t)}{dt} \end{pmatrix}
If you put, from $ (3-1) $, $ (3-2) $,
\frac{dx(t)}{dt}
=\frac{d}{dt} \begin{pmatrix} y(t) \\ \frac{dy(t)}{dt} \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ \frac{d^2y(t)}{dt^2} \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ -\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t) \end{pmatrix}
=\begin{pmatrix} x_2(t) \\ -\frac{2\gamma}{m} x_2(t)-\frac{k}{m}x_1(t)+\frac{a}{m} \sin(\omega t) \end{pmatrix}
Therefore, it can be expressed as follows.
\frac{dx(t)}{dt} = f(t,x(t))\hspace{50pt}(4)
\\where : f(t,x)=\begin{pmatrix} x_2 \\ -\frac{2\gamma}{m} x_2-\frac{k}{m}x_1+\frac{a}{m} \sin(\omega t) \end{pmatrix}
With the work up to this point, we were able to create a differential equation with a good-looking shape. Now, let's talk about the Eular method itself.
First, from the definition of differentiation
\frac{dx(t)}{dt} = f(t,x(t))
⇒\lim_{\Delta t \to 0} \frac{x(t+\Delta t)-x(t)}{\Delta t} = f(t,x(t))
Here, if you use $ \ Delta t $, which is small enough,
\frac{x(t+\Delta t)-x(t)}{\Delta t} \fallingdotseq f(t,x(t))
⇒x(t+\Delta t) \fallingdotseq x(t)+f(t,x(t))\Delta t\hspace{50pt}(5)
Using this formula, if you know the first $ x (t) $, you can calculate $ x (t + \ Delta t) $ after $ \ Delta t $ seconds.
The Eular method is implemented in python as follows.
f.py
import numpy as np
a = 10 #a(Amplitude of external force)
(m,gm2,k) = (30,5,1)#m,2Γ,k
omg = np.pi #ω
def f(t,x):
f0 = x[1]
f1 = -(gm2/m)*x[1] -(k/m)*x[0] + (a/m)*np.sin(omg*t)
return np.array([f0,f1])
f_D.py
import numpy as np
from f import f as f
dt = 0.01
def f_D_eular(t,x):
return (t + dt, x + f(t,x)*dt)
main.py
import numpy as np
from f_D import f_D_eular as f_D
def main():
(t0,x0)=(0,np.array([1,0]))
t1 = 100
(t,x)=(t0,x0)
while(t < t1):
print(t,x[0])
(t,x)=f_D(t,x)
main()
It will be like that. By the way, the graph of the execution result of this program is as follows.
The characteristics of the system can be seen by changing various parameters in the initial conditions $ t0, x0 $ and f.py.