常微分方程数值求解【python】

简述

这里只考虑最为简单的一种常微分方程

dydx=f(x,y)\frac{dy}{dx} = f(x,y)dxdy=f(x,y)

然后这里的实例都是以下面这个方程来做展示的。

dydx=y∗(1−y2)\frac{dy}{dx} = y*(1-y^2)dxdy=y∗(1−y2)
初值给定

y(0)=2y(0) = 2y(0)=2

这个方程的精确解结果是下面这个方程

y(x)=4∗e2x4∗e2x−3y(x) = \sqrt{\frac{4*e^{2x}}{4*e^{2x}-3}}y(x)=4∗e2x−34∗e2x

文章目录

简述
欧拉公式求解
- 简单的理论推理
- 代码实现
- - 实现后的效果
  - 代码
  - 误差画图
  - 误差画图代码
改进版欧拉公式
- 理解这个公式
- 改进版本的画图
- 欧拉算法和改进版欧拉算法的比较
- - 加上绝对值再来看
累积误差和分步的误差
- 图像
- 代码

欧拉公式求解

欧拉公式非常简洁。(欧拉果然大佬！！！)

yn+1=yn+h∗f(xn,yn)xn=x0+n∗hy_{n+1} = y_n + h*f(x_n, y_n)\\ x_n=x_0+n*hyn+1=yn+h∗f(xn,yn)xn=x0+n∗h

h是步长

简单的理论推理

其实非常直观，将上面的第一个式子变形一下，有

dydx=f(x,y)=yn+1−ynh\frac{dy}{dx} = f(x,y)=\frac{y_{n+1}- y_n}{h}dxdy=f(x,y)=hyn+1−yn

是不是非常熟悉，

微商的推导公式右边的加上一个取极限（h→0h\to0h→0）
之前的第一个式子，就是一阶泰勒展开。

代码实现

实现后的效果

x取值	精确值	数值逼近值	误差
0	2.0	2	0.0
0.1	1.6096571705090292	1.4	0.20965717050902932
0.2	1.4181045558702239	1.2656	0.15250455587022382
0.3	1.303667649645571	1.1894433603584	0.11422428928717099
0.4	1.2281238419433613	1.1401081630148027	0.08801567892855866
0.5	1.175177856120688	1.1059224047188059	0.06925545140188216
0.6	1.1365806251915203	1.0812532167953721	0.05532740839614814
0.7	1.1076620270914186	1.0629683037980915	0.0446937232933271
0.8	1.085560957285936	1.0491601738751934	0.03640078341074249
0.9	1.0684188538447474	1.0385912416407315	0.029827612204015974
1	1.0549729219451955	1.0304204608015664	0.02455246114362919

代码

import numpy as np# dy/dx = f(x,y)
# Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1
while x <= 1:loss = F(x) - yprint(x, '|', F(x), '|', y, '|', loss)y = y + h * f(y)x += 0.1

误差画图

误差画图代码

import numpy as np
import matplotlib.pyplot as plt# dy/dx = f(x,y)
# Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1
losses = []
while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y = y + h * f(y)x += 0.1plt.plot(np.arange(0, 1.1, 0.1), losses)
plt.plot(np.arange(0, 1.1, 0.1), losses, 'r*')
plt.show()

改进版欧拉公式

之前提到的欧拉公式简单，容易实现，但是效率却有点低。
所以有欧拉公式的改进版本。

yˉn+1=yn+h∗f(xn,yn)yn+1=yn+h2∗(f(xn,yn)+f(xn+1,yˉn+1))\bar y_{n+1} = y_n + h * f(x_n, y_n)\\ y_{n+1} = y_n + \frac{h}{2} * (f(x_n, y_n) + f(x_{n+1}, \bar y_{n+1}))yˉn+1=yn+h∗f(xn,yn)yn+1=yn+2h∗(f(xn,yn)+f(xn+1,yˉn+1))

x的话，任然是按照步长移动的，所以，这里就不作赘述了。

理解这个公式

书上的话，是用定积分的方式推理出来的。
但是这个具体含义，也没有给出什么东西来。
但是实际上，这个公式的逻辑含义是非常漂亮的！

将上面的两个公式联立起来。
yn+1−ynh=12∗(yˉn+1−ynh+f(xn+1,yˉn+1))\frac{y_{n+1}-y_n}{h} = \frac{1}{2}* (\frac{\bar y_{n+1} -y_n}{h} + f(x_{n+1}, \bar y_{n+1}))hyn+1−yn=21∗(hyˉn+1−yn+f(xn+1,yˉn+1))

之前的欧拉方法是用这里的 yˉn+1−ynh\frac{\bar y_{n+1} - y_n}{h}hyˉn+1−yn来近似dydx\frac{dy}{dx}dxdy
而现在我们逼近的其实是，用欧拉公式得到的导数，和实际上的在该点的导数值的均值。
从这个角度上来看，其实会是更贴近于正确的结果。

注意到： 我们这里用了yn+1y_{n+1}yn+1或者说是yny_nyn是正确解。
用这个其实是合理的，因为，这里的算法是基于之前的欧拉算法的，而欧拉算法本来就是会近似逼近正确的结果。所以，yny_nyn本身就是会逐渐接近于正确的结果。所以之前的假设是合理的。
而且，留心到，这是在原来的基础上，做了新的加速的过程。

欧拉公式的预估值，和代入其的校正值的平均值，就是改进欧拉方法的计算结果

改进版本的画图

import numpy as np
import matplotlib.pyplot as plt# dy/dx = f(x,y)
# Advanced Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1
losses = []
while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1+y2)/2x += 0.1plt.plot(np.arange(0, 1.1, 0.1), losses)
plt.plot(np.arange(0, 1.1, 0.1), losses, 'r*')
plt.show()

欧拉算法和改进版欧拉算法的比较

import numpy as np
import matplotlib.pyplot as plt# dy/dx = f(x,y)
# Advanced Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1def AdvancedEular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1 + y2) / 2x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(loss)y = y + h * f(y)x += 0.1return lossesadEl = AdvancedEular()
El = Eular()
plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Advance Eular')
plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular')
plt.plot(np.arange(0, 1.1, 0.1), El, 'b*')
plt.legend()
plt.show()

加上绝对值再来看

import numpy as np
import matplotlib.pyplot as plt# dy/dx = f(x,y)
# Advanced Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1def AdvancedEular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y1 = y + h * f(y)y2 = y + h * f(y1)y = (y1 + y2) / 2x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y = y + h * f(y)x += 0.1return lossesadEl = AdvancedEular()
El = Eular()
plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Advance Eular')
plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular')
plt.plot(np.arange(0, 1.1, 0.1), El, 'b*')
plt.legend()
plt.show()

累积误差和分步的误差

前面的话，我们是用一个不那么正确的结果来放到公式中来推导出正确的结果。下面我们尝试每一步都用正确的结果来推导出对应y值。

x取值	精确值	数值逼近值(每步使用对应的精确值推出)	误差
0	2.0	2	0.0
0.1	1.6096571705090292	1.4	0.20965717050902932
0.2	1.4181045558702239	1.35356132529304	0.06454323057718381
0.30000000000000004	1.303667649645571	1.274731273707409	0.028936375938162007
0.4	1.2281238419433613	1.2124696651611984	0.015654176782162965
0.5	1.175177856120688	1.1656997597866852	0.009478096334002872
0.6	1.1365806251915203	1.1303985272996449	0.006182097891875404
0.7	1.1076620270914186	1.103413438852542	0.004248588238876527
0.7999999999999999	1.085560957285936	1.082527495787655	0.003033461498280987
0.8999999999999999	1.0684188538447474	1.0661899261885104	0.0022289276562370564
0.9999999999999999	1.0549729219451955	1.0532987133870213	0.0016742085581742394

图像

代码

import numpy as np
import matplotlib.pyplot as plt# dy/dx = f(x,y)
# Euler formula
f = lambda x: x * (1 - x ** 2)
F = lambda x: ((4 * np.exp(2 * x) / (4 * np.exp(2 * x) - 3))) ** 0.5
x = 0
y = 2
h = 0.1def stepEular():x = 0y = 2h = 0.1losses = []while x <= 1:ty = F(x)loss = ty - yprint(x, '|', ty, '|', y, '|', loss)losses.append(abs(loss))y = ty + h * f(ty)x += 0.1return lossesdef Eular():x = 0y = 2h = 0.1losses = []while x <= 1:loss = F(x) - y# print(x, '|', F(x), '|', y, '|', loss)losses.append(abs(loss))y = y + h * f(y)x += 0.1return lossesadEl = stepEular()
El = Eular()
plt.plot(np.arange(0, 1.1, 0.1), adEl, label='Step Eular')
plt.plot(np.arange(0, 1.1, 0.1), adEl, 'r*')plt.plot(np.arange(0, 1.1, 0.1), El, label='Eular')
plt.plot(np.arange(0, 1.1, 0.1), El, 'b*')
plt.legend()
plt.show()