数值积分法求解常微分方程

import numpy as np
from math import sqrt
import sympy
import scipy
from scipy import integrate
from matplotlib import pyplot as plt

# 上篇的向量场绘图函数
def plot_directtion_field(x, y_x, f_xy, x_lim=(-1,1), y_lim=(-1,1), n=50, color='b', lw=0.5, ax=None):
    f_np = sympy.lambdify((x,y_x), f_xy, 'numpy')
    x_vec = np.linspace(x_lim[0], x_lim[1], n)
    y_vec = np.linspace(y_lim[0], y_lim[1], n)
    dx = x_vec[1] - x_vec[0]
    dy = y_vec[1] - y_vec[0]
    if ax is None:
        _, ax = plt.subplots(figsiz=(4,4))

    for xx in x_vec:
        for yy in y_vec:
            Dy = f_np(xx,yy) * dx
            ds = sqrt(dx*dx + Dy*Dy)
            Dx = 0.8 * dx * dx/ds
            Dy = 0.8 * Dy * dy/ds
            ax.plot([xx - Dx/2.0, xx + Dx/2.0], [yy - Dy/2, yy + Dy/2], color=color, lw=lw)

    ax.axis('tight')
    ax.set_title(r"$%s$" % (sympy.latex(sympy.Eq(y(x).diff(x), f_xy))),fontsize=16)
    return ax





if __name__ == '__main__':
    x = sympy.symbols('x')
    y = sympy.Function('y')
    f = y(x)**2 + x
    f_np = sympy.lambdify((y(x), x), f)

    x0, y0 = 0, -0.2
    xn = np.linspace(x0, x0-5, 100) # 初值处向x轴负方向延伸
    xp = np.linspace(x0, x0+2, 100) # 初值处向x轴正方向延伸
    yn = integrate.odeint(f_np, y0, xn) # 数值积分法求解常微分方程，负方向积分
    yp = integrate.odeint(f_np, y0, xp) # 数值积分法求解常微分方程，正方向积分




    fig, ax = plt.subplots(1, 1, figsize=(24, 20))
    # 绘出向量场以作对比
    plot_directtion_field(x, y(x), f, x_lim=(-5,5), y_lim=(-2,10), n=80, ax=ax) # 向量场图

    ax.plot(xn,yn,"r")
    ax.plot(xp,yp,"r") # 两段拼一起

    plt.show()