### 使用SciPy求解微分方程数值解 #### 单个常微分方程的求解方法 对于单个常微分方程，可以利用`scipy.integrate.solve_ivp`函数来获得其数值解。此函数提供了多种算法用于解决初值问题。 下面是一个简单的例子，展示如何使用该函数求解一阶线性ODE \( \frac{dy}{dt} = -2y, y(0)=1\) 的数值解： ```python from scipy.integrate import solve_ivp import numpy as np import matplotlib.pyplot as plt def model(t, y): # 定义要解的微分方程模型 dydt = -2 * y return dydt t_span = (0, 5) # 时间范围 y_init = [1] # 初始条件 sol = solve_ivp(model, t_span=t_span, y0=y_init) plt.plot(sol.t, sol.y[0], 'o', label='Numerical solution') plt.xlabel('Time $t$', fontsize=14) plt.ylabel('$y(t)$', fontsize=14) plt.legend() plt.show() ``` 这段代码定义了一个简单的一阶线性ODE并调用了`solve_ivp()`来进行求解[^1]。 #### 微分方程组的求解方式 当面对由多个相互关联的微分方程组成的系统时，同样可以通过传递一个返回向量形式导数列表给定的函数至`solve_ivp`实现求解。这里提供了一种处理两个耦合在一起的第一级非线性ODE系统的实例: 考虑洛伦兹吸引子(Lorenz attractor)，它是一套三个互连的一阶非线性ODEs，在混沌理论中有重要地位。这些方程式如下所示： \[ \begin{aligned} &\dot{x}=σ(y-x),\\ &\dot{y}=ρx-y-xz, \\ &\dot{z}=xy-βz. \end{aligned} \] 其中\( σ \)(sigma)、\( ρ \)(rho) 和 \( β \)(beta) 是参数。下面是具体的Python代码片段用来模拟这个过程： ```python def lorenz(t, state, sigma, rho, beta): x, y, z = state dx_dt = sigma*(y - x) dy_dt = x*(rho - z)-y dz_dt = x*y-beta*z return [dx_dt, dy_dt, dz_dt] initial_state = [-8., 7., 27.] params = {'sigma': 10., 'rho': 28., 'beta': 8./3} solution = solve_ivp(lorenz, t_span=(0, 50), y0=initial_state, args=(params['sigma'], params['rho'], params['beta'])) fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(projection="3d") ax.plot(solution.y[0], solution.y[1], solution.y[2]) ax.set_xlabel("X Axis") ax.set_ylabel("Y Axis") ax.set_zlabel("Z Axis") ax.set_title("Lorenz Attractor Trajectory") plt.show() ``` 上述代码展示了怎样设置以及可视化三维空间内的轨迹图。