solve_ivp卡住了,我搞不清楚为什么。

huangapple go评论68阅读模式
英文:

solve_ivp gets stuck and I can't figure out why

问题

I am using SciPy's solve_ivp to simulate a model consisting of 5 equations:

Upsilon = 0.305
Alpha = 0.00000055
r = 0.176
b = 0.0000000000005

def model(t, y0, Beta, Eta, Rmin, p1, p2, p3, A, a, Xi, Epsilon, Lambda, Theta, Mu, Delta, Gamma):

    CD, CT, CM, CE, T, TOTAL = y0

    def k_t(t):
        return Rmin + (p1 / (1 + (p2 * t) ** p3))

    def F_t(T):
        return T / (A + T)

    def f_t(CF, T):
        return (CF / T) / (Upsilon + ((a + CF) / T))

    dCD_dt = -(Beta + Eta) * CD
    dCT_dt = (Eta * CD) + (k_t(t) * F_t(T) * CT) - ((Xi + Epsilon + Lambda) * CT) + (Theta * T * CM) - (Alpha * T * CT)
    dCM_dt = (Epsilon * CT) - (Theta * T * CM) - (Mu * CM)
    dCE_dt = (Lambda * CT) - (Delta * CE)
    dT_dt = ((r * T) * (1 - (b * T))) - (Gamma * f_t((CD + CT), T) * T)

    TOTAL = dCD_dt + dCT_dt + dCM_dt + dCE_dt

    print(t)
    return [dCD_dt, dCT_dt, dCM_dt, dCE_dt, dT_dt, TOTAL]

I am performing sensitivity analysis, so I am running the model about 37,000 times with variations in the input parameters. When I run the model, I do so with the following code:

#p_values is a numpy array

final = 100
t_range = (0, final)
t_eval = np.arange(0, final, 0.1)

# run model for each parameter combination
for i in range(len(p_values)):
    CD = p_values[i][15]
    CT, CM, CE = 0, 0, 0
    T = 10 ** 7
    TOTAL = CD
    y0 = [CD, CT, CM, CE, T, TOTAL]

    solution = solve_ivp(model, t_range, y0, t_eval=t_eval, args=(p_values[i][0], p_values[i][1], p_values[i][2], p_values[i][3], p_values[i][4], p_values[i][5], p_values[i][6], p_values[i][7], p_values[i][8], p_values[i][9], p_values[i][10], p_values[i][11], p_values[i][12], p_values[i][13], p_values[i][14])).y

But at some point while it's running it will just... get stuck. Sometimes it happens early (within the first thousand solutions) and sometimes it can get to 10,000 solutions or more before it gets stuck.

To try and solve the problem, I started printing t for each loop of the model, and I noticed something very strange - when it gets stuck, it seems that t stays stagnant, or at least grows extremely slowly. This doesn't make sense to me, as I assume it should always increase by 0.1 at each step. Whenever this happens t is quite high (above 95), but I don't know if this is a coincidence or not. If not, I suppose it could be an issue with the function k(t) in the model producing a number that is far too small, but I don't see why this should be an issue.

If anyone with more experience than me has a suggestion to try and fix it, that would obviously be much appreciated. In lieu of that, I would also be okay with having some way to skip over that run of the model and move onto the next when the running time gets too long, but I don't know how to implement something like this without making the runs very slow.

英文:

I am using SciPy's solve_ivp to simulate a model consisting of 5 equations:

Upsilon = 0.305
Alpha = 0.00000055
r = 0.176
b = 0.0000000000005

def model(t, y0, Beta, Eta, Rmin, p1, p2, p3, A, a, Xi, Epsilon, Lambda, Theta, Mu, Delta, Gamma):

    CD, CT, CM, CE, T, TOTAL = y0

    def k_t(t):
        return Rmin + (p1 / (1 + (p2 * t) ** p3))

    def F_t(T):
        return T / (A + T)

    def f_t(CF, T):
        return (CF / T) / (Upsilon + ((a + CF) / T))

    dCD_dt = -(Beta + Eta) * CD
    dCT_dt = (Eta * CD) + (k_t(t) * F_t(T) * CT) - ((Xi + Epsilon + Lambda) * CT) + (Theta * T * CM) - (Alpha * T * CT)
    dCM_dt = (Epsilon * CT) - (Theta * T * CM) - (Mu * CM)
    dCE_dt = (Lambda * CT) - (Delta * CE)
    dT_dt = ((r * T) * (1 - (b * T))) - (Gamma * f_t((CD + CT), T) * T)

    TOTAL = dCD_dt + dCT_dt + dCM_dt + dCE_dt

    print(t)
    return [dCD_dt, dCT_dt, dCM_dt, dCE_dt, dT_dt, TOTAL]

I am performing sensitivity analysis, so I am running the model about 37,000 times with variations in the input parameters. When I run the model, I do so with the following code:

#p_values is a numpy array

final = 100
t_range = (0, final)
t_eval = np.arange(0, final, 0.1)

# run model for each parameter combination
for i in range(len(p_values)):
    CD = p_values[i][15]
    CT, CM, CE = 0, 0, 0
    T = 10 ** 7
    TOTAL = CD
    y0 = [CD, CT, CM, CE, T, TOTAL]

    solution = solve_ivp(model, t_range, y0, t_eval=t_eval, args=(p_values[i][0], p_values[i][1], p_values[i][2], p_values[i][3], p_values[i][4], p_values[i][5], p_values[i][6], p_values[i][7], p_values[i][8], p_values[i][9], p_values[i][10], p_values[i][11], p_values[i][12], p_values[i][13], p_values[i][14])).y

But at some point while it's running it will just... get stuck. Sometimes it happens early (within the first thousand solutions) and sometimes it can get to 10,000 solutions or more before it gets stuck.

To try and solve the problem, I started printing t for each loop of the model, and I noticed something very strange - when it gets stuck, it seems that t stays stagnant, or at least grows extremely slowly. This doesn't make sense to me, as I assume it should always increase by 0.1 at each step. Whenever this happens t is quite high (above 95), but I don't know if this is a coincidence or not. If not, I suppose it could be an issue with the function k(t) in the model producing a number that is far too small, but I don't see why this should be an issue.

If anyone with more experience than me has a suggestion to try and fix it, that would obviously be much appreciated. In lieu of that, I would also be okay with having some way to skip over that run of the model and move onto the next when the running time gets too long, but I don't know how to implement something like this without making the runs very slow.

答案1

得分: 0

你可以尝试在你的solve_ivp中使用可选的method参数,并将其设置为method = "BDF"。这可能会增加一些准确性。

根据
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

‘BDF’:基于向后差分公式的隐式多步可变阶(1到5)方法,用于导数近似[5]。实现遵循[6]中描述的方法。使用准恒定步长方案,并使用NDF修改来增强准确性。可以在复杂域中应用。

英文:

You could possibly try using the optional method parameter in your solve_ivp, and set it to method = "BDF". It might have a little more accuracy.

According to
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html :

> ‘BDF’: Implicit multi-step variable-order (1 to 5) method based on a
> backward differentiation formula for the derivative approximation [5].
> The implementation follows the one described in [6]. A quasi-constant
> step scheme is used and accuracy is enhanced using the NDF
> modification. Can be applied in the complex domain.

huangapple
  • 本文由 发表于 2023年5月10日 20:24:06
  • 转载请务必保留本文链接:https://go.coder-hub.com/76218384.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定