如何用Python求解微分方程组

要用Python求解微分方程组，需要使用一些数值求解工具库，例如Scipy库。以下是一个使用Scipy库解决微分方程组的简单示例：

 

首先，安装Scipy库：

pip install scipy

然后，导入必要的库：

import numpy as np

from scipy.integrate import solve_ivp

接下来，定义微分方程组。例如，假设要求解以下的 Lorenz 方程：

dx/dt = sigma * (y - x)

dy/dt = rho * x - y - x * z

dz/dt = -beta * z + x * y

其中 sigma，rho 和 beta 是常数。

 

你可以将这个微分方程组转换为一个函数，让 solve_ivp() 函数对它进行求解。例如：

def lorenz(t, xyz, sigma, rho, beta):

    x, y, z = xyz

    dxdt = sigma * (y - x)

    dydt = rho * x - y - x * z

    dzdt = -beta * z + x * y

    return [dxdt, dydt, dzdt]

接下来，设置初值和参数，然后调用 solve_ivp() 函数进行求解：

# 初值和参数

t0 = 0

t_end = 50

t_eval = np.linspace(t0, t_end, 1000)

xyz0 = [1.0, 1.0, 1.0]

sigma, rho, beta = 10, 28, 8/3

 

# 求解微分方程组

sol = solve_ivp(lorenz, [t0, t_end], xyz0, t_eval=t_eval, args=(sigma, rho, beta))

最后，你可以将求解结果可视化。例如，你可以绘制每个变量随时间的变化曲线：

import matplotlib.pyplot as plt

 

fig, axs = plt.subplots(3, 1, figsize=(8, 8))

axs[0].plot(sol.t, sol.y[0], 'b', lw=2)

axs[0].set_ylabel('x')

axs[1].plot(sol.t, sol.y[1], 'r', lw=2)

axs[1].set_ylabel('y')

axs[2].plot(sol.t, sol.y[2], 'g', lw=2)

axs[2].set_ylabel('z')

axs[2].set_xlabel('t')

plt.show()

这就是使用Python求解微分方程组的简单示例。你可以根据自己的需要修改微分方程组的定义、初值和参数等。