0

我目前正在尝试实现 RK4 例程。在此过程中，我正在努力匹配精度顺序，发现scipy.signal.lsim2()即使使用来自 ODEPACK 的 LSODA（它应该甚至比 RK4 更好）也会遇到同样的问题。

import numpy as np
import scipy as sp
from scipy import signal
from sympy import meijerg, exp, lambdify, vectorize
from sympy.abc import t

# time and input
t0 = np.linspace(0, 0.5, 200+1)  # time vector
u0 = np.ones_like(t0)  # input vector (step)

# transfer function
sys = signal.TransferFunction([147.4, 2948.], [1., 162.14, 2948.])

# exact solution
expr_t = 147.4*meijerg(((-140.272532338765, -20.0, -19.8674676612354, 1), ()),
                     ((), (-141.272532338765, -20.8674676612354, -19.0, 0)),
                     exp(t))
expr_ft = lambdify(t, expr_t)
expr_vft = vectorize(0)(lambda t: float(expr_ft(t)))
y_exact = np.array(expr_vft(t0))

h = t0[1] - t0[0]
print "Time step: h =", h
print "h^4 =", h**4
print "h^5 =", h**5
tout, y1, x1 = signal.lsim(sys, u0, t0)
print "Max global error lsim :", max(abs(y_exact - y1))
print "Max local  error lsim :", max(abs((y_exact[1:] - y_exact[:-1]) - (y1[1:] - y1[:-1])))
tout, y2, x2 = signal.lsim2(sys, u0, t0)
print "Max global error lsim2:", max(abs(y_exact - y2))
print "Max local  error lsim2:", max(abs((y_exact[1:] - y_exact[:-1]) - (y2[1:] - y2[:-1])))

这将打印：

Time step: h = 0.0025
h^4 = 3.90625e-11
h^5 = 9.765625e-14
Max global error lsim : 1.02140518266e-14
Max local  error lsim : 1.80966353014e-14
Max global error lsim2: 2.62325227174e-06
Max local  error lsim2: 5.13960676496e-06

由于 RK4 是 4 阶求解器，并且 LSODA 正在运行更多步骤，我不应该获得更好的精度吗？

为了完整起见，我也用 完成了所有scipy.integrate.ode()操作，将积分器设置为“dopri5”（显式 (4)5 阶 Runge-Kutta Dormand-Prince）：

def dy(t, y0, u0=1.):
    return np.dot(sys.A, y0)+np.dot(sys.B, np.atleast_1d(u0))
def fy(x, u0=1.):
    return np.dot(sys.C, x)+np.dot(sys.D, np.atleast_1d(u0))
ode3 = integrate.ode(dy)
ode3.set_integrator('dopri5').set_initial_value(np.r_[0.,0.])
x3 = np.zeros((len(t0), 2))
for i in range(1, len(t0)):
    x3[i] = ode3.integrate(t0[i])
y3 = np.squeeze(np.apply_along_axis(fy, 1, x3))
print "Max global error dopri5:", max(abs(y_exact - y3))
print "Max local  error dopri5:", max(abs((y_exact[1:] - y_exact[:-1]) - (y3[1:] - y3[:-1])))

哪个打印：

Max global error dopri5: 2.36753650018e-08
Max local  error dopri5: 6.88440715546e-09

最后，问题是：实施是否有任何问题，lsim2()或者那些错误级别是预期的？错误不应该低于h^4 = 3.90625e-11吗？我检查了大部分代码，没有发现任何问题。我还与我手工制作的 RK4 函数进行了比较，它也出现了 Xe-06 错误。我对 O(h^5) 和 o(h^4) 错误符号的解释可能有问题吗？

4

0 回答 0