python - Python解决ODE

Question

一个ODE怎么解，我试过用scipy.integrate.odeint，但是不知道初始值怎么办，</p>

这是我定义的：</p>

alpha=0.204
beta=0.196
b=5.853
c=241.38

def Left_equation(x):
   return alpha * x

def Right_equation(x):
   return (b * (x ** 2)) + c

def diff_equation(x, t):
   return (t ** beta) * (Right_equation(x) - Left_equation(x))

我也想得到结果图，不知道是否需要先得到解析解。

标签： python

Answer 1

如果初始条件未知，则需要象征性地进行积分。Python包sympy可用于常微分方程的符号积分，如下（使用函数sympy.dsolve）：

"""How to integrate symbolically an ordinary differential equation."""
import sympy


def main():
    alpha = 0.204
    beta = 0.196
    b = 5.853
    c = 241.38
    t = sympy.symbols('t')
    x = sympy.Function('x')(t)
    dxdt = sympy.Derivative(x, t)
    e = (t**beta) * ((b * (x**2)) + c - alpha * x)
    x_eq = sympy.dsolve(dxdt - e)


if __name__ == '__main__':
    main()

对于此示例，该函数sympy.solvers.ode.ode.dsolve引发异常PolynomialDivisionFailed。但是上面的代码显示了如何进行符号积分。

另一种方法是数值求解微分方程（针对特定初始条件），并计算一系列初始条件的解，以探索解如何依赖于初始条件。使用scipy.integrate.odepack.odeintPython 包scipy和 Python 包matplotlib和的功能numpy，可以如下完成：

"""How to integrate numerically an ordinary differential equation."""
from matplotlib import pyplot as plt
import numpy as np
import scipy.integrate


def main():
    x0 = np.linspace(0, 0.2, 10)  # multiple alternative initial conditions
    t = np.linspace(0, 0.01, 100)  # where to solve
    x = scipy.integrate.odeint(deriv, x0, t)
    # plot results
    plt.plot(t, x)
    plt.xlabel('t')
    plt.ylabel('x')
    plt.grid(True)
    plt.show()


def deriv(x, t):
    alpha = 0.204
    beta = 0.196
    b = 5.853
    c = 241.38
    return (t ** beta) * ((b * (x**2)) + c - alpha * x)


if __name__ == '__main__':
    main()

正如函数文档scipy.integrate.odeint中所述：

对于新代码，用于scipy.integrate.solve_ivp求解微分方程。

该函数scipy.integrate.solve_ivp可以按如下方式使用：

"""How to integrate numerically an ordinary differential equation."""
from matplotlib import pyplot as plt
import numpy as np
import scipy.integrate


def main():
    x0 = np.linspace(0, 0.2, 10)  # multiple alternative initial conditions
    t = (0.0, 0.01)  # where to solve: (start, end)
    res = scipy.integrate.solve_ivp(deriv, t, x0)  # different order of
        # arguments than for the function `scipy.integrate.odeint`
    # plot results
    plt.plot(res.t, res.y.T)
    plt.xlabel('t')
    plt.ylabel('x')
    plt.grid(True)
    plt.show()


def deriv(t, x):  # not x, t as for the function `scipy.integrate.odeint`
    alpha = 0.204
    beta = 0.196
    b = 5.853
    c = 241.38
    return (t ** beta) * ((b * (x**2)) + c - alpha * x)


if __name__ == '__main__':
    main()

scipy.integrate.odeint请注意 function和 function之间的参数差异scipy.integrate.solve_ivp。

上面的代码产生以下图：