在 Python 中模仿 ode45() 函数

Aditya Raj 2022年7月18日

Python Python Function

MatLab 中使用常微分方程来解决许多科学问题。ode45() 在 MatLab 中用于求解微分方程。

本文将展示我们如何在 Python 中模仿 ode45() 函数。

在 Python 中模仿 `ode45()` 函数

为了模仿 python 中的 ode45() 函数，我们可以使用 scipy 模块中定义的 solve_ivp() 方法。solve_ivp() 方法集成了一个常微分方程 (ODE) 系统。

solve_ivp() 方法将函数作为其第一个输入参数。输入参数中给出的函数必须返回一个包含微分方程系数的数组。
在第二个输入参数中，solve_ivp() 方法采用包含两个数值的元组或列表。这些值表示积分区间，其中元组中的第一个值表示区间的开始，元组的第二个值表示区间中的最大值。
在第三个输入参数中，solve_ivp() 方法采用一个表示初始值的数组。
执行后，solve_ivp() 方法返回具有各种属性的一堆对象。
1. t 属性包含一个包含时间点的 numpy 数组。
2. y 属性包含一个 numpy 数组，其值和时间点在 t 中。
3. sol 属性包含一个 Odesolution 对象，其中包含微分方程的解。如果在 solve_ivp() 方法中 dense_output 参数设置为 false，则 sol 属性包含 None。

为了更好地理解这一点，请参阅以下示例。

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
solution = solve_ivp(function, interval, initial_values)
print("Time:", solution.t)
print("Y:", solution.y)

输出：

Time: [ 0.          0.07578687  0.56581063  1.18741382  1.85887096  2.55035821
  3.25007544  3.95320486  4.65775424  5.36289544  6.06828346  6.77377445
  7.47930839  8.18486026  8.89041961  9.59598208 10.        ]
Y: [[1.00000000e+01 1.16366412e+01 3.10073783e+01 1.07492109e+02
  4.11689241e+02 1.64114780e+03 6.65071446e+03 2.71362627e+04
  1.11036049e+05 4.54874443e+05 1.86437495e+06 7.64300835e+06
  3.13352156e+07 1.28474398e+08 5.26752964e+08 2.15973314e+09
  4.84541488e+09]
 [1.50000000e+01 1.74549617e+01 4.65110674e+01 1.61238163e+02
  6.17533861e+02 2.46172171e+03 9.97607169e+03 4.07043941e+04
  1.66554074e+05 6.82311665e+05 2.79656243e+06 1.14645125e+07
  4.70028233e+07 1.92711598e+08 7.90129446e+08 3.23959970e+09
  7.26812231e+09]
 [2.50000000e+01 2.90916029e+01 7.75184457e+01 2.68730272e+02
  1.02922310e+03 4.10286951e+03 1.66267862e+04 6.78406569e+04
  2.77590123e+05 1.13718611e+06 4.66093739e+06 1.91075209e+07
  7.83380389e+07 3.21185996e+08 1.31688241e+09 5.39933284e+09
  1.21135372e+10]]

在上面的例子中，我们首先定义了一个名为 function 的函数，它以 t 和 y 作为其输入参数，并根据 y 返回一个值。

然后，我们分别使用变量 interval 和 initial_values 为 ODE 定义了一个区间和初始值。我们将 function、interval 和 initial_values 作为输入参数传递给 solve_ivp() 函数，最后，我们在变量解决方案中获得输出。

在输出中，你可以观察到时间值分布在 0 到 10 的区间内。同样，输出包含与每个时间值对应的 y 值。

我们还可以在解决方案的属性 t 中明确指定时间点。为此，我们需要将包含所需时间值的数组传递给 solve_ivp() 方法的 t_eval 参数，如下所示。

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
time_values = [1, 2, 3, 6, 7, 8]
solution = solve_ivp(function, interval, initial_values,t_eval=time_values)
print("Time:", solution.t)
print("Y:", solution.y)

输出：

Time: [1 2 3 6 7 8]
Y: [[7.38683416e+01 5.46053271e+02 4.03089733e+03 1.62618365e+06
  1.20160156e+07 8.87210156e+07]
 [1.10802512e+02 8.19079906e+02 6.04634600e+03 2.43927547e+06
  1.80240234e+07 1.33081523e+08]
 [1.84670854e+02 1.36513318e+03 1.00772433e+04 4.06545912e+06
  3.00400390e+07 2.21802539e+08]]

你可以看到时间值仅包含那些作为输入参数传递给 t_eval 参数的值。类似地，属性 y 仅包含指定 t 值的值。

这种方法可以帮助你获取区间中某些点的值。

在 Python 中模仿 ode45() 函数

相关文章 - Python Function

在 Python 中模仿 `ode45()` 函数