python - 使用 scipysolve_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()