Let's look at a differential equation that cannot be solved normally in Python

Euler method

When $ \ dot {x} = f (x) $, I would like to find a function $ g $ that ideally integrates $ x $ over time so that $ x = g (t) $. However, analytical integration is not possible in nonlinear systems. Therefore, let us consider solving it numerically by discretizing the time. It would be good if we could find a numerical value similar to the integral of $ f (x) $ over time for each t. Therefore, the Euler method, which is an application of the segmented quadrature method that combines rectangles, can be used.

If you take a small time width $ \ Delta {t} $ and find the area of ​​a rectangle with a height of $ f (x_0) $ and a width of $ \ Delta {t} $, you can find the integral value in that time width. A value close to the solution $ g (t) $) is obtained. So, if you use this and consider the value $ x_1 $ corresponding to the place where you advanced by $$ Delta {t} $ from the initial value, $ x (t_0 + \ Delta {t}) \ approx x_0 + f (x_0) It can be \ Delta {t} = x_1 $. In generalization, $ x_ {n + 1} = x_n + f (x_n) \ Delta {t} $, and by repeating this, it is possible to obtain something close to $ g (t) $.

However, since this is a method of finding the next approximate value from the approximate value, errors will accumulate more and more. In order to improve the accuracy of the approximation, we will take a closer look at the area of ​​the rectangle.

Modified Euler method

A simple idea is to use the average height of the rectangle with $ f (x_ {n + 1}) $ instead of just $ f (x_n) $. This is called the modified Euler method. The error converges faster than the Euler method, but the theory of computation of the error is not well understood. If you imagine a monotonous function in the $ \ Delta {t} $ interval, it's easy to see that this will make it stick out of the rectangle and reduce scraping.

Runge-Kutta method

The Euler method adopted one end (for the younger one) and the modified Euler method adopted both ends as the height of a rectangle, but the Runge-Kutta method is a more elaborate selection of the numerical values ​​to be adopted. The runge-Kutta method of the fourth order is just right in terms of the balance with the amount of calculation. In the quaternary Runge-Kutta method, the rectangular values ​​are calculated for the four points in the interval, and the integrated values ​​in the interval are calculated by averaging with weights. Specifically as follows

Often used is Runge-Kutta to draw the solution trajectory when given the most interesting initial value many times like a graph of a function, and give various other initial values ​​to Runge- You can grasp the whole atmosphere by trying Kutta only once for each and drawing it as a vector field or a tilt field. I will implement it below.

import numpy as np
import matplotlib.pyplot as plt

#Define the function you want to look up
def dotx(x):
  return x*(1-x) 

#Runge-Kutta method-Plot on x plane
def rk4th(x0, n, dt):#x0:Initial value, n: t range, dt: step width
  t=0
  x=x0
  xvalues=[]
  tvalues=np.arange(0,n,dt)
  scope=len(tvalues)#Keep the number of elements unified with this

  #Apply RungeKutta to dotx and draw a trajectory corresponding to x0
  for i in range(scope):
    xvalues.append(x)
    k1 = dotx(x)*dt
    k2 = dotx(x+k1/2)*dt
    k3 = dotx(x+k2/2)*dt
    k4 = dotx(x+k3)*dt
    x = x + (k1+ 2*k2 + 2*k3 + k4)/6
    t += dt  

  #Draw a hard spot by finding it from other initial values
  T, X = np.meshgrid(
    np.arange(0,n,1),
    np.arange(0,2,0.1)
  )
  k1 = dotx(X)*dt
  k2 = dotx(X+k1/2)*dt
  k3 = dotx(X+k2/2)*dt
  k4 = dotx(X+k3)*dt
  V = (k1+2*k2+2*k3+k4)/6
  plt.quiver(T,X,dt,V)
  plt.plot(tvalues, xvalues)
  plt.show()

rk4th(2, 10, 0.1)

The results are shown below. This is really a slope, not a vector, because the horizontal axis is t.