本文介绍了绘制一阶ODE的解及其导数的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

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

来自数学导入sqrt的

 从 numpy 导入零、linspace、数组从scipy.integrate导入odeint导入 matplotlib.pylab 作为 pltdef f(y,r):f =零(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)返回f米 = 1.15k = 3 *平方(1-1/m)l = 2.0r = 1.0e-10rf = 1.0rspan = linspace(r,rf,1000)y0 =数组([l])溶胶= odeint(f,y0,rspan)plt.plot(rspan,sol,'k:',lw = 1.5)
解决方案

来自 odeint

要解决有关 squeeze 的评论,请参阅下文(摘录自

其中 n 的定义依据:

I have this code to solve a simple first order ODE using odeint. I managed to plot the solution y(r), but I also want to plot the derivative y'= dy/dr. I know y' it is given by f(y,r), but I'm not sure how to call the function with the output of the integration. Thank you in advance.

    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)
解决方案

From odeint documentation:

I have modified your code in the following manner and obtained the figure below.

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

To address the comment about squeeze, please see below (extracted from scipy.integrate.solve_ivp):

where n is defined according to:

这篇关于绘制一阶ODE的解及其导数的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

05-29 04:54
查看更多