求解高阶微分方程的数值解

题目

        使用odeint求解这个方程的数值解: 

思路: 分解为方程组进行计算. 也就是带入, 化简为下面的方程: 

 代码实现

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

# 定义微分方程组
# 输入两个function以及自变量, 输出两个function的导数
def dydx(y, x):
    y, u = y
    dydx = [u, 20 * (1 - y**2) * u - y]
    return dydx

"""else
def fvdp(y,t):
    # 也可以不解包参数, 直接使用list, 其中y[0]=y,y[1]=u
    dy1 = y[1]
    dy2 = 20*(1-y[0]**2) * y[1] - y[0]
    return [dy1, dy2] 

"""

# 初始条件
# y(0) = 0, u(0)=y'(0) = 2
y0 = [0, 2]

# 自变量范围
x = np.arange(0, 0.25, 0.01)

# 求解微分方程
# 使用ode_integration, 输入微分方程组&初始值(与x的初始值对应)&自变量
sol = odeint(dydx, y0, x)

# 提取解 第一列是y,第二列是u
y = sol[:, 0]
u = sol[:, 1]

# 打印结果
print("x={}\n对应的数值解y={}".format(x, y))

# 绘制结果
# plot函数, 输入自变量, 因变量. label给出标签
plt.plot(x, y, label="y(x)")
plt.plot(x, u, label="y'(x)")
plt.xlabel("x")
plt.ylabel("y, y'")
plt.legend()
plt.show()

 结果: 

x=[0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13
 0.14 0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24]
对应的数值解y=[0.         0.02213973 0.04917542 0.08217383 0.12241129 0.17139008
 0.23083061 0.30261249 0.38862014 0.49043034 0.6087762  0.7427756
 0.8890802  1.041387   1.19095297 1.32843735 1.44644047 1.54133536
 1.61342302 1.6657038  1.70228635 1.72720565 1.74383888 1.7547616
 1.76182797]

进程已结束，退出代码为 0

参考结果: 

参考

马马老师文章