matplotlib - 绘制一阶 ODE 的解及其导数

Question

我有这段代码可以使用 odeint 解决简单的一阶 ODE。我设法绘制了解 y(r)，但我还想绘制导数 y'= dy/dr。我知道 y' 它由 f(y,r) 给出，但我不确定如何使用积分输出调用该函数。先感谢您。

    from math import sqrt
    from numpy import zeros,linspace,array
    from scipy.integrate import odeint
    import matplotlib.pylab as plt

    def f(y,r):
        f = zeros(1)
        f[0] = -(y[0]*(y[0]-1.0)/r)-y[0]*(2/r+\
        ((r/m)/(1-r**2/m))*(2*sqrt(1-r**2/m)-k)/(k-sqrt(1-r**2/m)))\
        -(1/(1-r**2/m))*(-l*(l+1)/r+\
         (3*r/m)*(k+2*sqrt(1-r**2/m))/(k-sqrt(1-r**2/m)))\
        +((4*r**3)/((m**2)*(1-r**2/m)))*(1/(k-sqrt(1-r**2/m))**2)
        return f

    m = 1.15            
    k = 3*sqrt(1-1/m)
    l = 2.0
    r = 1.0e-10                         
    rf = 1.0                         

    rspan = linspace(r,rf,1000)
    y0 = array([l])
    sol = odeint(f,y0,rspan)
    plt.plot(rspan,sol,'k:',lw=1.5)

标签： matplotlibodeint

Answer 1

从odeint 文档：

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

我按照以下方式修改了你的代码，得到了下图。

import matplotlib.pyplot as plt
from numpy import gradient, squeeze, sqrt
from scipy.integrate import solve_ivp


def fun(t, y):
    l = 2
    m = 1.15
    k = 3 * sqrt(1 - 1 / m)
    return (-y * (y - 1) / t - y * (2 / t + t / m / (1 - t ** 2 / m)
                                    * (2 * sqrt(1 - t ** 2 / m) - k)
                                    / (k - sqrt(1 - t ** 2 / m)))
            - 1 / (1 - t ** 2 / m) * (-l * (l + 1) / t + 3 * t / m
                                      * (k + 2 * sqrt(1 - t ** 2 / m))
                                      / (k - sqrt(1 - t ** 2 / m)))
            + 4 * t ** 3 / m ** 2 / (1 - t ** 2 / m)
            / (k - sqrt(1 - t ** 2 / m)) ** 2)


sol = solve_ivp(fun, t_span=[1e-10, 1], y0=[2], method='BDF',
                dense_output=True)
if sol.success is True:
    print(sol.t.shape, sol.y.shape)
    plt.plot(sol.t, squeeze(sol.y), color='xkcd:avocado',
             label='Scipy Solution')
    plt.plot(sol.t, fun(sol.t, squeeze(sol.y)), linestyle='dashed',
             color='xkcd:purple', label='Derivative Using the Function')
    plt.plot(sol.t, gradient(squeeze(sol.y), sol.t), linestyle='dotted',
             color='xkcd:bright orange', label='Derivative Using Numpy')
    plt.legend()
    plt.tight_layout()
    plt.savefig('so.png', bbox_inches='tight', dpi=300)
    plt.show()

要解决关于 的评论squeeze，请参见下文（摘自scipy.integrate.solve_ivp）：

其中n根据以下定义：