使用 scipy 的 solve_ivp 在时间步之间运行代码

Question

我正在将我的代码从使用 scipy 的 odeint 转换为 scipy 的 solve_ivp。当使用 odeint 时，我会使用如下的 while 循环：

while solver.successful() : 
    solver.integrate(t_final, step=True)
    # do other operations

这种方法允许我在每个时间步之后存储取决于解决方案的值。

我现在切换到使用solve_ivp，但不确定如何使用solve_ivp 求解器完成此功能。有人用solve_ivp 完成了这个功能吗？

谢谢！

Answer 1

我想我知道你想问什么。我有一个程序，它使用 solve_ivp 在每个时间步之间单独积分，然后使用这些值来计算下一次迭代的值。 （即传热系数、传输系数等）我使用了两个嵌套的 for 循环。内部 for 循环计算或完成您在每个步骤中需要执行的操作。然后将每个值保存在列表或数组中，然后内部循环应该终止。外层循环应该只用于输入时间值并可能重新加载必要的常量。

例如：

for i in range(start_value, end_value, time_step):
start_time = i
end_time = i + time_step
# load initial values and used most recent values
    for j in range(0, 1, 1):


    answer = solve_ivp(function,(start_time,end_time), [initial_values])
    # Save new values at the end of a list storing all calculated values

假设你有一个系统比如

1. d(Y1)/dt = a1*Y2 + Y1

2. d(Y2)/dt = a2*Y1 + Y2

并且您想从 t = 0, 10 解决它。时间步长为 0.1。其中 a1 和 a2 是在别处计算或确定的值。这段代码可以工作。

from scipy.integrate import solve_ivp
import sympy as sp
import numpy as np
import math
import matplotlib.pyplot as plt



def a1(n):
       return 1E-10*math.exp(n)

def a2(n):
       return 2E-10*math.exp(n)

def rhs(t,y, *args):
       a1, a2 = args
       return [a1*y[1] + y[0],a2*y[0] + y[1]]

Y1 = [0.02]
Y2 = [0.01]
A1 = []
A2 = []
endtime = 10 
time_step = 0.1
times = np.linspace(0,endtime, int(endtime/time_step)+1)
tsymb = sp.symbols('t')
ysymb = sp.symbols('y')
for i in range(0,endtime,1):

       for j in range(0,int(1/time_step),1):
              tstart = i + j*time_step
              tend = i + j*time_step + time_step
              A1.append(a1(tstart/100))
              A2.append(a2(tstart/100))
              Y0 = [Y1[-1],Y2[-1]]
              args = [A1[-1],A2[-1]]
              answer = solve_ivp(lambda tsymb, ysymb : rhs(tsymb,ysymb, *args), (tstart,tend), Y0)
              Y1.append(answer.y[0][-1])
              Y2.append(answer.y[1][-1])

fig = plt.figure()
plt1 = plt.plot(times,Y1, label = "Y1")
plt2 = plt.plot(times,Y2, label = "Y2")
plt.xlabel('Time')
plt.ylabel('Y Values')
plt.legend()
plt.grid()
plt.show()