我scipy.integrate.ode
在解决逻辑函数时比较了一些不同的 ODE 积分器:
我听说 LSODA 应该很好,所以我有点惊讶地发现它完全失败了. dopri5
另一方面,对于非常大的值,积分器似乎没有问题(见下图)。
为什么 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
结果integrator = 'dopri5'
:
dopri5:平均错误:3.7564128626655536e-11