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

k_1 = f(x_n)\Delta t
k_2 = f(x_n+\frac{1}{2}k_1)\Delta t
k_3 = f(x_n+\frac{1}{2}k_2\Delta t
k_4 = f(x_n+\frac{1}{2}k_3)\Delta t

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
  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):
    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(
  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.plot(tvalues, xvalues)

rk4th(2, 10, 0.1)

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

Recommended Posts

Let's look at a differential equation that cannot be solved normally in Python
A solution to the problem that the Python version in Conda cannot be changed
Let's create a script that registers with Ideone.com in Python.
[Python] It seems that global variables cannot be referenced in Multiprocessing
A record that GAMEBOY could not be done in Python. (PYBOY)
Take a look at the built-in exception tree structure in Python 3.8.2
Let's create a customer database that automatically issues a QR code in Python
Let's make a combination calculation in Python
I took a quick look at the fractions package that handles Python built-in fractions.
I want to create a priority queue that can be updated in Python (2.7)
I made a familiar function that can be used in statistics with Python
A memo that I wrote a quicksort in Python
A program that removes duplicate statements in Python
A mechanism to call a Ruby method from Python that can be done in 200 lines
I want to write in Python! (2) Let's write a test
[Memorandum] Japanese keys cannot be used in python string.Template.substitute
What's in that variable (when running a Python script)
list comprehension because operator.methodcaller cannot be used in python 2.5
In Python, create a decorator that dynamically accepts arguments Create a decorator
Operators ++,-cannot be used in python (difference from php)
Non-linear simultaneous equations can be easily solved in Python.
[Redash] Standard library cannot be used in python function
MALSS, a tool that supports machine learning in Python
Let's take a quick look at CornerNet, an object detector that does not use anchors.
A Python program in "A book that gently teaches difficult programming"
A general-purpose program that formats Linux command strings in python
A function that divides iterable into N pieces in Python
Published a library that hides character data in Python images
Loop through a generator that returns a date iterator in Python
Processing of python3 that seems to be usable in paiza
Let's make a spot sale service 4 (in Python mini Hack-a-thon)
Video cannot be loaded with Spyder in Python development environment
I tried "a program that removes duplicate statements in Python"
Create code that outputs "A and pretending B" in python
[MQTT / Python] Implemented a class that does MQTT Pub / Sub in Python
Error that Qt plugin "cocoa" cannot be found in python-opencv
Take a look at the Python built-in exception tree structure
Scripts that can be used when using bottle in Python
Precautions that must be understood when building a PYTHON environment
A set of script files that do wordcloud in Python3
Let's make a diagram that can be clicked with IPython
AtCoder ABC168 A case expression solved in Ruby and Python
Here's a summary of things that might be useful when dealing with complex numbers in Python
・ <Slack> Write a function to notify Slack so that it can be quoted at any time (Python)