Introduction

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.

Preparing to use the Eular method

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}

Actual example

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.

It is possible to solve it analytically, but we will consider solving it by a numerical calculation method.

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}

Eular method

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.

Program (python)

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.