[PYTHON] Perturbation development and simulation Perturbation method

What is a perturbation method?

The perturbation method is an approximate calculation method of equations that has been highly developed in the field of physics. By using this solution method, it is possible to approximate the differential equations and partial differential equations that are difficult to solve by ordinary calculations. When the perturbation method is classified into two types, it can be divided into two types: "perturbation problem" and "singular perturbation problem". The former can be solved step by step, but there are few established solutions to the singular perturbation problem.

In this article, I will introduce the perturbation problem method. The results of deriving the solution of a simple differential equation and performing numerical simulation are shown.

First-order nonlinear ODE

\dot{y} + y + \epsilon y^2 =0 \\
s.t.   y(0)=1, \epsilon=o(1)

It is well known that solving nonlinear problems is generally very difficult. Here we assume that $ \ epsilon $ is small enough. Simply put, if $ \ epsilon $ is small enough, some might wonder if it should be treated as 0. The solution at that time is easy

y=exp(-x)

Can be derived. Is this approximation a good approximation?

def sim(time,y=1,e=0.01):
    save=[y]
    dt=0.01*e
    step=int(time/dt)
    for i in range(step):
        y=y-dt*y-e*dt*y**2
        save.append(y)
    return save
y1=sim(2)
a=np.linspace(0,2,len(y1))
y=np.exp(-a)
plt.plot(sim(2),a)
plt.plot(y,a)
plt.show()

As you can see, it is not much different from the true distribution. This approximation is called a 0th-order approximation in perturbation theory. To improve this approximation, let's do a perturbation expansion.

y=y_0+\epsilon y_1+\epsilon^2 y_2+...

Create an asymptotic series of $ \ epsilon $ in this way. Substituting this into the differential equation and solving up to the second order,

y_0=exp(-x)\\
\dot{y_1}+y_1+y_o^2=0\\
y_1=-exp(-x)+exp(-2x)\\
\dot{y_2}+y_2+2y_0y_1=0\\
y_2=exp(-x)-2exp(-2x)+exp(-3x)\\

Simulate this.

y1=sim(2)
a=np.linspace(0,2,len(y1))
y=np.exp(-a)
y2=np.exp(-2*a)
y3=np.exp(-3*a)
plt.plot(sim(2),a)
plt.plot(y,a)
plt.plot(y1+(y2-y)*0.01,a)
plt.plot(y1+(y2-y)*0.01+(y-2*y2+y3)*0.0001,a)
plt.savefig("simu")
plt.show()

As you can see from the results, the accuracy is very good.

Conclusion

It can be seen that the accuracy of the approximation is considerably improved by performing the approximation calculation in this way.

Next time, I will introduce the perturbation problem that this approximate calculation does not work, and explain the singular perturbation method.