使用scipy.integrate进行常微分方程组的数值求解

scipy.integrate是Python的科学计算库scipy中的一个子包，它提供了常微分方程组的数值求解功能。常微分方程组是由多个微分方程组成的方程组，其中每个方程都涉及一个未知函数及其导数。例如，考虑以下的常微分方程组：

dy1/dt = f1(y1, y2, t)

dy2/dt = f2(y1, y2, t)

其中y1(t)和y2(t)是未知函数，t是独立变量，f1和f2是给定的函数。我们的目标是找到y1(t)和y2(t)的数值解。

首先，我们需要导入必要的模块和函数：

import numpy as np

from scipy.integrate import solve_ivp

接下来，我们定义f1和f2函数：

def f1(y1, y2, t):

    return 2*y1 + y2

def f2(y1, y2, t):

    return -y1 + 2*y2 + np.sin(t)

这些函数定义了常微分方程组的右侧。现在，我们可以使用solve_ivp函数求解常微分方程组的数值解：

# 创建时间间隔

t_span = (0, 10)

# 创建初始条件

y0 = [1, 0]

# 求解常微分方程组

sol = solve_ivp(fun = lambda t, y: [f1(y[0], y[1], t), f2(y[0], y[1], t)], t_span = t_span, y0 = y0)

solve_ivp函数采用几个参数来控制求解过程。fun参数指定常微分方程组的右侧函数，t_span参数指定求解的时间间隔，y0参数指定初始条件。解sol包含解的所有信息。

我们可以使用sol.t和sol.y属性获取时间和解的数组：

# 获取时间和解的数组

t = sol.t

y = sol.y

下面是一个完整的例子，将上述步骤整合在一起：

import numpy as np

from scipy.integrate import solve_ivp

# 定义常微分方程组的右侧函数

def f1(y1, y2, t):

    return 2*y1 + y2

def f2(y1, y2, t):

    return -y1 + 2*y2 + np.sin(t)

# 创建时间间隔

t_span = (0, 10)

# 创建初始条件

y0 = [1, 0]

# 求解常微分方程组

sol = solve_ivp(fun = lambda t, y: [f1(y[0], y[1], t), f2(y[0], y[1], t)], t_span = t_span, y0 = y0)

# 获取时间和解的数组

t = sol.t

y = sol.y

# 打印结果

print("Time:", t)

print("Solution:", y)

在这个例子中，我们使用solve_ivp函数求解了常微分方程组，并打印出了时间和解的数组。注意，sol.y是解矩阵，其中每一行对应一个未知函数的解。