时滞微分方程求数值解的solve_ivp代码

以下是使用Python中的`scipy.integrate.solve_ivp`函数求解时滞微分方程的示范代码：

import numpy as np
from scipy.integrate import solve_ivp

# 定义时滞微分方程
def delay_system(t, y, tau):
    y_tau = np.interp(t-tau, t, y)
    dydt = -0.1*y[0] + 0.2*y[1] + 0.1*y_tau
    dzdt = -0.2*y[1] + 0.2*y_tau
    return [dydt, dzdt]

# 定义初始条件和时间范围
y0 = [0.5, 0.5]
t_span = [0, 100]

# 定义时滞
tau = 3

# 使用solve_ivp函数求解时滞微分方程
sol = solve_ivp(delay_system, t_span, y0, args=(tau,), dense_output=True)

# 绘制结果
t_eval = np.linspace(0, 100, 1000)
y_eval = sol.sol(t_eval)
import matplotlib.pyplot as plt
plt.plot(t_eval, y_eval[0], label='y')
plt.plot(t_eval, y_eval[1], label='z')
plt.legend()
plt.show()

在上述代码中，我们首先定义了一个时滞微分方程`delay_system`，其中`t`和`y`分别表示时间和系统状态向量，`tau`表示时滞。然后，我们定义了初始条件`y0`和时间范围`t_span`。接下来，我们使用`solve_ivp`函数求解时滞微分方程，其中`args=(tau,)`表示将时滞`tau`作为额外参数传递给`delay_system`函数，`dense_output=True`表示返回一个可用于插值的函数对象。最后，我们使用`sol.sol`函数计算时间范围内的状态，并绘制结果。