python - Accuracy lost when using scipy.signal.lsim2 -

i'm trying implement rk4 routine. in process, i'm struggling match accuracy order , found out scipy.signal.lsim2() suffers same problem, if using lsoda odepack (which should better rk4).

import numpy np import scipy sp scipy import signal sympy import meijerg, exp, lambdify, vectorize sympy.abc import t  # time , 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]))) 

this print:

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 

since rk4 4th order solver , lsoda running many more steps, shouldn't better accuracy?

for sake of completeness, did scipy.integrate.ode(), setting integrator "dopri5" (explicit (4)5th order 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)) 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]))) 

which prints:

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

finally, question is: is there problem lsim2() implementation or error levels expected?. shouldn't errors lower h^4 = 3.90625e-11? have checked of code , couldn't find problems. have compared hand-made rk4 function , gets xe-06 errors. could problem interpretation of o(h^5) , o(h^4) error notation?