常微分方程（ODE）的积分曲线和方向场可视化（Python）

​今天老师布置了个一阶线性微分方程的python可视化作业，由于作者本人水平有限（爆哭），之后再把非线性和高阶微分方程学会了再一并补充进来。

一阶微分方程

一阶线性微分方程

基本概念

对于一阶微分方程，方程内只存在一阶导，常见的是如下形式：
$$\frac{dy}{dx}=f(x,y)$$

积分曲线：

该方程的解会得到一个$y$关于$x$的函数（无论是显式还是隐式的），构成一条曲线即为积分曲线。

方向场图：

同时，我们将上述方程的右边$f(x,y)$ 的大小作为斜率绘制在各个点处（即各个点箭头的方向斜率即为$f(x,y)$的值），得到方向场图。

等倾斜线图：

之后，我们将各个箭头方向一致的点连接起来，构成等斜线图，方程$\frac{dy}{dx}=f(x,y)$的等斜线即为$f(x,y)=k$

例子1： $\frac{dy}{dx}=x^2-y$

其积分曲线和方向场的图如下：


此外老师还说了要画等值线，或许这和等高线差不多？

Python（pycharm）代码如下：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义微分方程
def model(x, y):
    return x-y


def draw1(x_values, y_values):  #绘制积分曲线
    # 绘制结果
    plt.plot(x_values, y_values)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title(' $\\frac{dy}{dx} = x^2 - y$ ')
    plt.grid(True)
    plt.show()

def gradient(x,y):
    return x-y

def draw2(x,y):  # 绘制方向场的图
    u=0.5  # 矢量箭头的长度均为1
    alpha = np.arctan(gradient(x, y))
    plt.quiver(x, y, u*np.cos(alpha), u*np.sin(alpha), angles='xy', scale_units='xy', scale=1, color='blue')
    plt.show()


def draw_contour(X, Y):  # 等倾斜线图,注意输入的是网格化的x，y
    k = gradient(X, Y)  # 计算每个点的斜率

    plt.figure(figsize=(8, 6))
    contours = plt.contour(X, Y, k, levels=20, cmap='viridis')
    plt.colorbar(contours, label='k')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid(True)
    plt.show()
# 设置初始条件
y0 = [0]  # 初始条件 y(0) = 0
# 设置时间范围
x_span = (-4, 4)  # 从 x = 0 到 x = 10
# 调用 solve_ivp
solution = solve_ivp(model, x_span, y0, method='RK45', t_eval=np.linspace(-4 ,4, 100))
# 提取结果
x_values = solution.t
y_values = solution.y[0]
draw1(x_values, y_values)

x_vec = np.linspace(-3, 3, 20)
y_vec = np.linspace(-3, 3, 20)
X_vec,Y_vec=np.meshgrid(x_vec, y_vec)
draw2(X_vec, Y_vec)
draw_contour(X_vec,Y_vec)

scipy.integrate.solve_ivp函数常用于求解微分方程，无论是线性或是非线性，低阶或是高阶，该函数都可以获得相当准确的数值解，传入参数解释如下：

solution = solve_ivp(model, x_span, y0, method=‘RK45’, t_eval=np.linspace(-4 ,4, 100))

传入参数	解释
model	传入微分方程的模型，一般来说是$f(x,y)$
x_span	需要以$(a,b)$格式传入两个参数，表示$x$的范围
y0	初值条件
method=‘RK45’	Runge-Kutta-Fehlberg 方法,常见的还有’RK23’、‘DOP853’、‘LSODA’、‘BDF’
t_eval	选取一些点方便后续作图
solution.y	是一个二维数组，其中 solution.y[0] 提取出第一个变量（在一阶微分方程中就是 y）在所有时间点上的解

例子2：$\frac{dy}{dx}=x-y$

积分曲线和方向场如下：

不好意思，这题没代码咯，拿例1代码玩去吧