为什么 LSODA 无法整合逻辑功能?

计算科学 Python scipy
2021-12-18 14:34:36

scipy.integrate.ode在解决逻辑函数时比较了一些不同的 ODE 积分器:

x(t)=11+ert
x˙=rx(1x)

我听说 LSODA 应该很好,所以我有点惊讶地发现它完全失败了r>4. dopri5另一方面,对于非常大的值,积分器似乎没有问题r(见下图)。

为什么 LSODA 在这里失败了?我是否需要调整一些参数,或者对于这个特定问题来说它只是一个糟糕的选择?

这是我的代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

# Logistic function
r = 5
def func(t):
    return 1. / (1 + np.exp(-r * t))

def dfunc(t, X):
    x, = X
    return [r * x * (1 - x)]

t0 = -10
t1 = 10
dt = 0.01
X0 = func(t0)

integrator = 'lsoda'

t = [t0]
Xi = [func(t0)]
ode = integrate.ode(dfunc)
ode.set_integrator(integrator)
ode.set_initial_value(X0, t0)

while ode.successful() and ode.t < t1:
    t.append(ode.t + dt)
    Xi.append(ode.integrate(ode.t + dt))

t = np.array(t)     # Time
X = func(t)         # Solution
Xi = np.array(Xi)   # Numerical

# Plot analytical and numerical solutions
X = func(t)

plt.subplot(211)
plt.title("Solution")
plt.plot(t, X, label="analytical", color='g')
plt.plot(t, Xi, label=integrator)
plt.xlabel("t")
plt.ylabel("x")
plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.05))

# Plot numerical integration errors
err = np.abs(X - Xi.flatten())
print("{}: mean error: {}".format(integrator, np.mean(err)))

plt.subplot(212)
plt.title("Error")
plt.plot(t, err, label=integrator)
plt.xlabel("t")
plt.ylabel("error")
plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.05))

plt.tight_layout()
plt.show()

结果integrator = 'lsoda'

lsoda:平均误差:0.500249750249742

结果与 lsoda

结果integrator = 'dopri5'

dopri5:平均错误:3.7564128626655536e-11

dopri5 的结果

1个回答

当你使用r=5,初始条件为

x(10)e501.9×1022.
这比机器 epsilon 小得多,2×1016,并且 LSODA 很可能只是得出解同样为零的结论。

为了检查这个想法,我更换了t05代替10, 以便x(5)1.4×1011,然后它似乎工作。解决方案不是很准确,错误上升到0.01, 但这很可能是因为数值误差将凸块的位置放置在稍微错误的位置,而解的形状是正确的。

LSODA 似乎决定解决方案相同为零的方式是它使用其绝对容差,并且可能x(10)比那个小,导致问题。所以要解决这个问题,而不是改变t0,一个也可以换线

ode.set_integrator(integrator)

ode.set_integrator(integrator, atol=0.0)

然后它似乎也可以工作,但错误是105. 要将错误降低到与 DOPRI5 相同的级别,我可以设置rtol=1e-10,但我不太确定为什么 DOPRI5 使用其默认参数进行管理,或者这些默认值是否不同。无论如何,我认为这是这里工作的默认公差。