Python解ODE

Question

如何解决一个ODE，我尝试过使用scipy.integrate.odeint，但是如果我不知道初始值，我该怎么办，

​

​

我就是这么定义的：

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))

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

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。但是上面的代码展示了如何进行符号集成。

另一种方法是数值求解微分方程(对于特定初始条件)，并计算一系列初始条件的解，以探索解是如何依赖于初始条件的。使用package scipy的函数scipy以及packages 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和函数scipy.integrate.solve_ivp之间的参数差异。

上面的代码产生了以下情节：

​

​