python微分求解

发表于 2021-10-15 更新于 2024-08-31 分类于 Python 阅读次数：本文字数： 1.6k 阅读时长 ≈ 6 分钟

高阶微分方程，可以用python“数值求解”或者“解析式求解”，很方便，解析式求解有些可能无解，数值求解画图方便一些

这里用简单函数

假设我们现在知道原函数 \(f(x)=x^2\)

那么二阶微分方程为： \[f''(x)+f'(x)+(-2-2\times x) = 0\] 初始条件 \[f(0)=0 \quad f'(0)=0\]

先看代码

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


def eq_ode(x, ys):
    return -ys[1] + 2 + 2 * x


def fvdp(x, ys):
    return [
        ys[1],
        eq_ode(x, ys),
    ]


def eq_sim(x, f):
    return sy.diff(f, x, 2) + sy.diff(f, x, 1) + (-2 - 2 * x)


# 数值求解
def solve_ode():
    x = np.linspace(0, 5, 1000)
    y0 = [0, 0]
    ys = odeint(fvdp, y0, x, tfirst=True)

    y, = plt.plot(x, ys[:, 0], label='y')
    y1, = plt.plot(x, ys[:, 1], label='y\'')
    y2, = plt.plot(x, eq_ode(x, ys[:, 0]), label='y\'\'')
    plt.legend(handles=[y, y1, y2])
    plt.savefig("test")


# 解析式求解
def slove_sim_eq():
    x = sy.symbols('x')
    f = sy.Function('f')(x)
    s = sy.dsolve(eq_sim(x, f), f)
    sy.pprint(s)
    return s


slove_sim_eq()

solve_ode()

解析式求解

使用 sympy 库 dsolve 函数

微分方程在 eq_test 函数里，其中 3 个注释掉了，分别是

\[f''(x)+sin(x)=0 \quad f''(x)+x^2=0 \quad f''(x)+f(x)=0\]

import sympy as sy


# 二阶常系数齐次微分方程
def eq_test(x, f):
    # return sy.diff(f(x), x, 2) + sy.sin(x)  # f''+sin(x)=0
    # return sy.diff(f, x, 2) + x**2  # f''+x^2=0
    # return sy.diff(f, x, 2) + f  # f(x)''+f(x)=0
    return sy.diff(f, x, 2) + sy.diff(f, x, 1) + (-2 - 2 * x)


# 化简方程
def slove_eq():
    x = sy.symbols('x')  # 约定变量
    f = sy.Function('f')(x)  # 约定函数
    s = sy.dsolve(eq_test(x, f), f)
    # 美化输出，在vscode-ipynb里，甚至可以输出为latex
    sy.pprint(s)
    print(sy.printing.latex(s))
    return s


slove_eq()

使用函数 sympy，eq_test 函数即为微分方程，sy.pprint可以美化输出，print(sy.printing.latex(s))可以输出为 latex 代码

解析解：

\[\displaystyle f{\left(x \right)} = C_{1} + C_{2} e^{- x} + x^{2}\]

数值求解

scipy 里内置函数 odeint, solve_ivp 两个都可以，区别如下

函数	作用	速度	区别
odeint	微分方程数值求解	快	参数简洁，推荐
solve_ivp	如果使用 `method='LSODA'` 与 `odeint` 一样	相对慢一些	通用，含多个计算方式

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


def eq(t, ys):
    return -ys[1] + 2 + 2 * t


def fvdp(t, ys):
    return [
        ys[1],
        eq(t, ys),
    ]


def solve_second_order_ode():
    '''
    求解二阶ODE
    '''
    t = np.linspace(0, 5, 1000)
    tspan = (0, 5)
    # ! 注意啊，这个初值一定要对应上面t的初值
    y0 = [0, 0]  # 初值条件
    # 初值[1.0, 2.0]表示y(0)=1,y'(0)=2

    # 返回y，其中y[:,0]是y[0]的值，就是最终解，y[:,1]是y'(x)的值，tfirst为了与下面的solve_ivp定义的fvdp2中（y，t）顺序保持一致
    y_ = odeint(fvdp, y0, t, tfirst=True)
    ys = np.transpose(y_)
    title = "ode"

    # y_ = solve_ivp(fvdp, t_span=tspan, method="LSODA", y0=y0, t_eval=t)
    # ys = y_.y
    # title = "ivp"

    print(ys)

    y = ys[0]
    y_1 = ys[1]
    y_2 = eq(t, ys)

    y2, = plt.plot(t, y, label='y(0)')
    y2_1, = plt.plot(t, y_1, label='y(1)')
    y2_2, = plt.plot(t, y_2, label='y(2)')
    plt.legend(handles=[y2, y2_1, y2_2])
    plt.title(title)
    plt.savefig("test_" + title)


solve_second_order_ode()

odeint 函数

odeint 函数一般只需要传递 3 个参数

微分方程矩阵函数、自变量、初始值

1	odeint(fvdp, y0, t, tfirst=True)

微分方程矩阵函数

对于高阶，需要构造微分方程函数 fvdp，即最开始提到的方程为

\[ f''(x)+f'(x)+(-2-2x) = 0 \quad f(0)=0 \quad f'(0)=0 \]

要把这个二阶方程全部转化为矩阵

\[[y', y'']\]

即 fvdp 函数的返回值，这个函数的参数为 t 自变量，y 因变量的导数矩阵 \([y, y']\)

odeint 函数的参数顺序 (y, t)，如果要与 solve_ivp 一致参数顺序变成 (t, y)，odeint 函数多传递一个参数 tfirst=True

所以，返回值 \(y'=y[1]\)，y'' 用微分方程 (\(f''(x)+f'(x)+(-2-2x) = 0\)) 表示

\(y''=-y[1]+2+2t\)

最终 fvdp 函数返回值

\[[y', y'']=[y[1], -y[1]+2+2t]\]

自变量

用 numpy 数组表示即可 t = np.linspace(0, 5, 1000)

初始值

注意与自变量对应，二阶微分方程初始值为

\[[f(t_0), f'(t_0)]\]

这里的\(t_0\)必须是自变量的第一个值，不然会有很大的误差

返回值

odeint 函数返回值只有 ys=\([[y], [y']]\) ，但是注意转置一下，行列不太方便，如果不转置

\(y=ys[:, 0] \qquad y'=ys[:, 1]\)

solve_ivp 函数

1	y_ = solve_ivp(fvdp, t_span=tspan, method="LSODA", y0=y0, t_eval=t)

微分方程矩阵函数

solve_ivp 不一样的仅仅只有参数顺序，这个函数参数应该是 (t, y) 与 matlab 是一致的，

自变量

这个也有区别，自变量数组在 t_eval=t，但是还需要一个 t_span=(min, max)，其中 max min 为自变量数组的最大值最小值，如果只有t_span，取得值太少，不行

初始值

与 odeint 一致

返回值

ys = y_.y

\(y=ys[0] \qquad y'=ys[1]\)

总结

还是用 odeint吧，参数少一些，方便

其他

solve_ivp 积分：

详细看官网文档，包含各种方法的具体参考文献

'RK45'（默认）：显式 Runge-Kutta 5(4)阶。假设四阶方法的精度控制误差，但使用五阶精确公式（完成局部外推）采取步骤。四次插值多项式用于密集输出[2]。可应用于复杂领域。
'RK23'：显式 Runge-Kutta (2)3 阶。假设二阶方法的精度，控制误差，但使用三阶精确公式（完成局部外推）采取步骤。三次 Hermite 多项式用于密集输出。可应用于复杂领域。
'DOP853'：8 阶显式 Runge-Kutta。最初用 Fortran 编写的“DOP853”算法的 Python 实现。精确到 7 阶的 7 阶插值多项式用于密集输出。可应用于复杂领域。
'Radau'：5 阶 Radau IIA 的隐式 Runge-Kutta。误差由三阶精确嵌入公式控制。满足搭配条件的三次多项式用于密集输出。
'BDF'：基于导数逼近的后向微分公式的隐式多步变阶（1 到 5）方法。使用准恒定步长方案，并使用 NDF 修改提高精度。可应用于复杂领域。
'LSODA'：具有自动刚度检测和切换功能的 Adams/BDF 方法。这是来自 ODEPACK 的 Fortran 求解器的封装。

显式 Runge-Kutta 方法（'RK23'、'RK45'、'DOP853'）应用于非刚性问题，隐式方法（'Radau'、'BDF'）应用于刚性问题。在 Runge-Kutta 方法中，建议使用“DOP853”进行高精度求解（rtol 和 atol 的值较低）。

如果不确定，请先尝试运行“RK45”。如果它进行了异常多次的迭代、发散或失败，则您的问题可能很僵硬，您应该使用“Radau”或“BDF”。'LSODA' 也可以是一个很好的通用选择，但它可能不太方便，因为它包装了旧的 Fortran 代码。

本文作者：yuhldr
本文地址： https://yuhldr.github.io/posts/13877.html
版权声明：转载请注明出处！