一阶常微分方程方向场图的绘制

用户6021899

发布于 2023-03-03 09:07:20

9470

发布于 2023-03-03 09:07:20

文章被收录于专栏：Python编程 pyqt matplotlibPython编程 pyqt matplotlib

方向场图可用于可视化一阶常微分方程的可能解。方向场图由XY平面网格中未知函数斜率的短线组成。y(x) 在XY平面上任意一点的斜率由微分方程

定义给出。在方向场图中，与斜率相切的连续平滑曲线都是该微分方程的可能解。

Python 代码如下:

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

def plot_directtion_field(x, y_x, f_xy, x_lim=(-5,5), y_lim=(-5,5), n=50, color='r', 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')
    fig, axes = plt.subplots(1, 3, figsize=(24, 8))
    plot_directtion_field(x, y(x), -y(x)**2 + y(x) + x, ax=axes[0])
    plot_directtion_field(x, y(x), -x/y(x), ax=axes[1])
    plot_directtion_field(x, y(x), y(x)**3 / x, ax=axes[2])
    plt.show()