[Numerical calculation method, python] Solving ordinary differential equations by Eular method

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.

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.

_data.png

The characteristics of the system can be seen by changing various parameters in the initial conditions $ t0, x0 $ and f.py.

Recommended Posts

[Numerical calculation method, python] Solving ordinary differential equations by Eular method
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation
[Science / technical calculation by Python] Numerical solution of first-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
Solve ordinary differential equations in Python
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
Numerical calculation of differential equations with TensorFlow 2.0
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Numerical calculation to find the value of derivative (differential)
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
Numerical calculation with Python
Numerical calculation of compressible fluid by finite volume method
Solve simultaneous ordinary differential equations with Python and SymPy.
Solving one-dimensional wave equation by finite difference method (python)
[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
Python --Differential equation numerical solution Euler method & central difference method & Runge-Kutta method
[Python] Calculation method with numpy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
Numerical analysis of ordinary differential equations with Scipy's odeint and ode
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Science and technology calculation by Python] (Partial) differential, mathematical formula, sympy
Solving simultaneous linear equations in Python (sweeping method and fractional representation)
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
Create a new Python numerical calculation project
Alignment algorithm by insertion method in Python
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
I tried it with Wolfram Alpha by referring to "Solving simultaneous linear equations with Python (sweeping method and fractional expression)".
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics