Python scipy.integrate.solve_ivp() Examples

def integrate_TOV(self):
        """
        Evolves TOV+k2 equations and returns final quantities
        """

        # integration settings the same as in lalsimulation
        rel_err = 1e-4
        abs_err = 0.0

        result = solve_ivp(self.__tov_eqns, (self.pseudo_enthalpy, 1e-16), self.y, rtol=rel_err,
                           atol=abs_err)
        m_fin = result.y[0, -1]
        r_fin = result.y[1, -1]
        H_fin = result.y[2, -1]
        B_fin = result.y[3, -1]

        k_2 = self.__calc_k2(r_fin, B_fin, H_fin, m_fin / r_fin)

        return m_fin, r_fin, k_2 

def simulate(self, x_0, t_eval, rtol=1e-6, atol=1e-6):
        """Simulate closed-loop system using Runge-Kutta 4,5.

        Solution is evaluated at N time steps. Outputs times and corresponding
        solutions as numpy array (N,) * numpy array (N, n).

        Inputs:
        Initial condition, x_0: numpy array (n,)
        Solution times, t_eval: numpy array (N,)
        RK45 relative tolerance, rtol: float
        RK45 absolute tolerance, atol: float
        """

        t_span = [t_eval[0], t_eval[-1]]
        sol = solve_ivp(self.dx, t_span, x_0, t_eval=t_eval, rtol=rtol, atol=atol)
        return sol.t, sol.y.T 

def reconstruction():
    """Reconstruction of the dynamics of the Kuramoto model"""

    def approximated_dynamics(_, theta):
        """Construction of the right-hand side of the system from the coefficient tensor"""

        cores = [np.zeros([1, theta.shape[0] + 1, 1, 1])] + [np.zeros([1, theta.shape[0] + 1, 1, 1]) for _ in
                                                             range(1, p)]
        for q in range(p):
            cores[q][0, :, 0, 0] = [1] + [psi[q](theta[r]) for r in range(theta.shape[0])]
        psi_x = TT(cores)
        psi_x = psi_x.full().reshape(np.prod(psi_x.row_dims), 1)

        rhs = psi_x.transpose() @ xi
        rhs = rhs.reshape(rhs.size)
        return rhs

    sol = spint.solve_ivp(approximated_dynamics, [0, time], x_0, method='BDF', t_eval=np.linspace(0, time, m))
    sol = sol.y

    return sol 

def test_empty():
    def fun(t, y):
        return np.zeros((0,))

    y0 = np.zeros((0,))

    for method in ['RK23', 'RK45', 'Radau', 'BDF', 'LSODA']:
        sol = assert_no_warnings(solve_ivp, fun, [0, 10], y0,
                                 method=method, dense_output=True)
        assert_equal(sol.sol(10), np.zeros((0,)))
        assert_equal(sol.sol([1, 2, 3]), np.zeros((0, 3)))

    for method in ['RK23', 'RK45', 'Radau', 'BDF', 'LSODA']:
        sol = assert_no_warnings(solve_ivp, fun, [0, np.inf], y0,
                                 method=method, dense_output=True)
        assert_equal(sol.sol(10), np.zeros((0,)))
        assert_equal(sol.sol([1, 2, 3]), np.zeros((0, 3))) 

def take_step(self, s, u):
        # clip the action
        u = np.clip(u, self.action_range[0], self.action_range[1])

        # concatenate s and u to pass through ds_dt
        s_aug = np.append(s, u)

        # solve the differientable equation to compute next state
        s_next = solve_ivp(self.ds_dt, (0., self.time_interval), s_aug).y[0:2, -1] # last index is the action applied

        # project state to the valid range
        s_next[StateIndex.THETA] = wrap(s_next[StateIndex.THETA], self.angle_range[0],
                     self.angle_range[1])
        s_next[StateIndex.THETA_DOT] = np.clip(s_next[StateIndex.THETA_DOT],
                                         self.angular_velocity_range[0],
                                         self.angular_velocity_range[1])

        return s_next 

def take_step(self, s, u):
        # clip the action
        u = np.clip(u, self.action_range[0], self.action_range[1])

        # concatenate s and u to pass through ds_dt
        s_aug = np.append(s, u)

        # solve the differientable equation to compute next state
        s_next = solve_ivp(self.ds_dt, (0., self.time_interval), s_aug).y[0:4, -1] # last index is the action applied

        # project state to the valid space.
        s_next[StateIndex.THETA] = wrap(s_next[StateIndex.THETA], self.angle_range[0],
                     self.angle_range[1])
        s_next[StateIndex.THETA_DOT] = np.clip(s_next[StateIndex.THETA_DOT],
                                               self.angular_velocity_range[0],
                                               self.angular_velocity_range[1])
        s_next[StateIndex.X_DOT] = np.clip(s_next[StateIndex.X_DOT],
                                            self.velocity_range[0],
                                            self.velocity_range[1])

        return s_next 

def integrate(self, t):
        # Docstring of superclass
        result = solve_ivp(
            self._system_equation, [self._t, t], self._y, t_eval=[t], args=self._f_params, **self._solver_kwargs
        )
        self._t = t
        self._y = result.y.T[-1]
        return self._y 

def __init__(self, *args, **kwargs):
        raise NotImplementedError(
            """It is not yet clear how to use scipy's `solve_ivp` with topology changes

Previous attempts where made but turned out to be clumsy..."""
        ) 

def test_no_integration():
    for method in ['RK23', 'RK45', 'Radau', 'BDF', 'LSODA']:
        sol = solve_ivp(lambda t, y: -y, [4, 4], [2, 3],
                        method=method, dense_output=True)
        assert_equal(sol.sol(4), [2, 3])
        assert_equal(sol.sol([4, 5, 6]), [[2, 2, 2], [3, 3, 3]]) 

def test_integration_const_jac():
    rtol = 1e-3
    atol = 1e-6
    y0 = [0, 2]
    t_span = [0, 2]
    J = jac_linear()
    J_sparse = csc_matrix(J)

    for method, jac in product(['Radau', 'BDF'], [J, J_sparse]):
        res = solve_ivp(fun_linear, t_span, y0, rtol=rtol, atol=atol,
                        method=method, dense_output=True, jac=jac)
        assert_equal(res.t[0], t_span[0])
        assert_(res.t_events is None)
        assert_(res.success)
        assert_equal(res.status, 0)

        assert_(res.nfev < 100)
        assert_equal(res.njev, 0)
        assert_(0 < res.nlu < 15)

        y_true = sol_linear(res.t)
        e = compute_error(res.y, y_true, rtol, atol)
        assert_(np.all(e < 10))

        tc = np.linspace(*t_span)
        yc_true = sol_linear(tc)
        yc = res.sol(tc)

        e = compute_error(yc, yc_true, rtol, atol)
        assert_(np.all(e < 15))

        assert_allclose(res.sol(res.t), res.y, rtol=1e-14, atol=1e-14) 

def test_integration_sparse_difference():
    n = 200
    t_span = [0, 20]
    y0 = np.zeros(2 * n)
    y0[1::2] = 1
    sparsity = medazko_sparsity(n)

    for method in ['BDF', 'Radau']:
        res = solve_ivp(fun_medazko, t_span, y0, method=method,
                        jac_sparsity=sparsity)

        assert_equal(res.t[0], t_span[0])
        assert_(res.t_events is None)
        assert_(res.success)
        assert_equal(res.status, 0)

        assert_allclose(res.y[78, -1], 0.233994e-3, rtol=1e-2)
        assert_allclose(res.y[79, -1], 0, atol=1e-3)
        assert_allclose(res.y[148, -1], 0.359561e-3, rtol=1e-2)
        assert_allclose(res.y[149, -1], 0, atol=1e-3)
        assert_allclose(res.y[198, -1], 0.117374129e-3, rtol=1e-2)
        assert_allclose(res.y[199, -1], 0.6190807e-5, atol=1e-3)
        assert_allclose(res.y[238, -1], 0, atol=1e-3)
        assert_allclose(res.y[239, -1], 0.9999997, rtol=1e-2) 

def test_integration_complex():
    rtol = 1e-3
    atol = 1e-6
    y0 = [0.5 + 1j]
    t_span = [0, 1]
    tc = np.linspace(t_span[0], t_span[1])
    for method, jac in product(['RK23', 'RK45', 'BDF'],
                               [None, jac_complex, jac_complex_sparse]):
        with suppress_warnings() as sup:
            sup.filter(UserWarning,
                       "The following arguments have no effect for a chosen solver: `jac`")
            res = solve_ivp(fun_complex, t_span, y0, method=method,
                            dense_output=True, rtol=rtol, atol=atol, jac=jac)

        assert_equal(res.t[0], t_span[0])
        assert_(res.t_events is None)
        assert_(res.success)
        assert_equal(res.status, 0)

        assert_(res.nfev < 25)
        if method == 'BDF':
            assert_equal(res.njev, 1)
            assert_(res.nlu < 6)
        else:
            assert_equal(res.njev, 0)
            assert_equal(res.nlu, 0)

        y_true = sol_complex(res.t)
        e = compute_error(res.y, y_true, rtol, atol)
        assert_(np.all(e < 5))

        yc_true = sol_complex(tc)
        yc = res.sol(tc)
        e = compute_error(yc, yc_true, rtol, atol)

        assert_(np.all(e < 5)) 

def integrate(self, tf, **kwargs):
        """Integrate up to a time tf.
        """
        t0 = self.t
        y0 = self._y
        rhs = self.rhs.get('rhs')

        # solve problem and silence warnings for options that don't apply to a given method
        kept_warnings = []
        with warnings.catch_warnings(record=True) as ws:
            results = solve_ivp(rhs, (t0, tf), y0,
                                method=self.options.method,
                                atol=self.options.atol,
                                rtol=self.options.rtol,
                                max_step=self.options.max_step,
                                min_step=self.options.min_step,
                                first_step=self.options.first_step,
                                **kwargs)

            # update the internal state
            self._y = results.y[:, -1]
            self._t = results.t[-1]

            # discard warnings for arguments with no effect
            for w in ws:
                if 'The following arguments have no effect' not in str(w.message):
                    kept_warnings.append(w)

        # display warnings we don't want to silence
        for w in kept_warnings:
            warnings.warn(w.message, type(w)) 

def integrate_streamline(
    X, Y, U, V, integration_direction, init_states, interpolation_num=250, average=True
):
    """use streamline's integrator to alleviate stacking of the solve_ivp. Need to update with the correct time."""
    import matplotlib.pyplot as plt

    n_cell = init_states.shape[0]

    res = np.zeros((n_cell * interpolation_num, 2))
    j = -1  # this index will become 0 when the first trajectory found

    for i in tqdm(range(n_cell), "integration with streamline"):
        strm = plt.streamplot(
            X,
            Y,
            U,
            V,
            start_points=init_states[i, None],
            integration_direction=integration_direction,
            density=100,
        )
        strm_res = np.array(strm.lines.get_segments()).reshape((-1, 2))

        if len(strm_res) == 0:
            continue
        else:
            j += 1
        t = np.arange(strm_res.shape[0])
        t_linspace = np.linspace(t[0], t[-1], interpolation_num)
        f = interpolate.interp1d(t, strm_res.T)

        cur_rng = np.arange(j * interpolation_num, (j + 1) * interpolation_num)
        res[cur_rng, :] = f(t_linspace).T

    res = res[: cur_rng[-1], :]  # remove all empty trajectories
    n_cell = int(res.shape[0] / interpolation_num)

    if n_cell > 1 and average:
        t_len = len(t_linspace)
        avg = np.zeros((t_len, 2))

        for i in range(t_len):
            cur_rng = np.arange(n_cell) * t_len + i
            avg[i, :] = np.mean(res[cur_rng, :], 0)

        res = avg

    plt.close()

    return t_linspace, res


# ---------------------------------------------------------------------------------------------------
# fate related 

def kuramoto(theta_init, frequencies, time, number_of_snapshots):
    """Kuramoto model

    Generate data for the Kuramoto model represented by the differential equation

        d/dt x_i = w_i + (2/d) * (sin(x_1 - x_i) + ... + sin(x_d - x_i)) + 0.2 * sin(x_i).

    See [1]_ and [2]_ for details.

    Parameters
    ----------
    theta_init: ndarray
        initial distribution of the oscillators
    frequencies: ndarray
        natural frequencies of the oscillators
    time: float
        integration time for BDF method
    number_of_snapshots: int
        number of snapshots

    Returns
    -------
    snapshots: ndarray(number_of_oscillators, number_of_snapshots)
        snapshot matrix containing random displacements of the oscillators in [-0.1,0.1]
    derivatives: ndarray(number_of_oscillators, number_of_snapshots)
        matrix containing the corresponding derivatives

    References
    ----------
    .. [1] P. Gelß, S. Klus, J. Eisert, C. Schütte, "Multidimensional Approximation of Nonlinear Dynamical Systems",
           arXiv:1809.02448, 2018
    .. [2] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, R. Spigler, "The Kuramoto model: A simple
           paradigm for synchronization phenomena", Rev. Mod. Phys. 77, pp. 137-185 , 2005
    """

    number_of_oscillators = len(theta_init)

    def kuramoto_ode(_, theta):
        [theta_i, theta_j] = np.meshgrid(theta, theta)
        return frequencies + 2 / number_of_oscillators * np.sin(theta_j - theta_i).sum(0) + 0.2 * np.sin(theta)

    sol = spint.solve_ivp(kuramoto_ode, [0, time], theta_init, method='BDF',
                          t_eval=np.linspace(0, time, number_of_snapshots))
    snapshots = sol.y
    derivatives = np.zeros([number_of_oscillators, number_of_snapshots])
    for i in range(number_of_snapshots):
        derivatives[:, i] = kuramoto_ode(0, snapshots[:, i])
    return snapshots, derivatives 

def test_max_step():
    rtol = 1e-3
    atol = 1e-6
    y0 = [1/3, 2/9]
    for method in [RK23, RK45, Radau, BDF, LSODA]:
        for t_span in ([5, 9], [5, 1]):
            res = solve_ivp(fun_rational, t_span, y0, rtol=rtol,
                            max_step=0.5, atol=atol, method=method,
                            dense_output=True)
            assert_equal(res.t[0], t_span[0])
            assert_equal(res.t[-1], t_span[-1])
            assert_(np.all(np.abs(np.diff(res.t)) <= 0.5))
            assert_(res.t_events is None)
            assert_(res.success)
            assert_equal(res.status, 0)

            y_true = sol_rational(res.t)
            e = compute_error(res.y, y_true, rtol, atol)
            assert_(np.all(e < 5))

            tc = np.linspace(*t_span)
            yc_true = sol_rational(tc)
            yc = res.sol(tc)

            e = compute_error(yc, yc_true, rtol, atol)
            assert_(np.all(e < 5))

            # See comment in test_integration.
            if method is not LSODA:
                assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)

            assert_raises(ValueError, method, fun_rational, t_span[0], y0,
                          t_span[1], max_step=-1)

            if method is not LSODA:
                solver = method(fun_rational, t_span[0], y0, t_span[1],
                                rtol=rtol, atol=atol, max_step=1e-20)
                message = solver.step()

                assert_equal(solver.status, 'failed')
                assert_("step size is less" in message)
                assert_raises(RuntimeError, solver.step) 

def setUp(self):
        """Consider the Fermi-Pasta-Ulam problem and Kuramoto model for testing the routines in sle.py"""

        # set tolerance
        self.tol = 1e-5

        # number of oscillators
        self.fpu_d = 4
        self.kuramoto_d = 10

        # parameters for the Fermi-Pasta_ulam problem
        self.fpu_m = 2000
        self.fpu_psi = [lambda t: 1, lambda t: t, lambda t: t ** 2, lambda t: t ** 3]

        # parameters for the Kuramoto model
        self.kuramoto_x_0 = 2 * np.pi * np.random.rand(self.kuramoto_d) - np.pi
        self.kuramoto_w = np.linspace(-5, 5, self.kuramoto_d)
        self.kuramoto_t = 100
        self.kuramoto_m = 1000
        self.kuramoto_psi = [lambda t: np.sin(t), lambda t: np.cos(t)]
        self.kuramoto_basis = [[tdt.constant_function()] + [tdt.sin(i,1) for i in range(self.kuramoto_d)], [tdt.constant_function()] + [tdt.cos(i,1) for i in range(self.kuramoto_d)]]
        self.kuramoto_initial = tt.ones([11, 11], [1, 1], 11)

        # exact coefficient tensors
        self.fpu_xi_exact = mdl.fpu_coefficients(self.fpu_d)
        self.kuramoto_xi_exact = mdl.kuramoto_coefficients(self.kuramoto_d, self.kuramoto_w)

        # generate test data for FPU
        self.fpu_x = 0.2 * np.random.rand(self.fpu_d, self.fpu_m) - 0.1
        self.fpu_y = np.zeros((self.fpu_d, self.fpu_m))
        for j in range(self.fpu_m):
            self.fpu_y[0, j] = self.fpu_x[1, j] - 2 * self.fpu_x[0, j] + 0.7 * (
                    (self.fpu_x[1, j] - self.fpu_x[0, j]) ** 3 - self.fpu_x[0, j] ** 3)
            for i in range(1, self.fpu_d - 1):
                self.fpu_y[i, j] = self.fpu_x[i + 1, j] - 2 * self.fpu_x[i, j] + self.fpu_x[i - 1, j] + 0.7 * (
                        (self.fpu_x[i + 1, j] - self.fpu_x[i, j]) ** 3 - (self.fpu_x[i, j] - self.fpu_x[i - 1, j]) ** 3)
                self.fpu_y[-1, j] = - 2 * self.fpu_x[-1, j] + self.fpu_x[-2, j] + 0.7 * (
                        -self.fpu_x[-1, j] ** 3 - (self.fpu_x[-1, j] - self.fpu_x[-2, j]) ** 3)

        # generate test data for Kuramoto
        number_of_oscillators = len(self.kuramoto_x_0)

        def kuramoto_ode(_, theta):
            [theta_i, theta_j] = np.meshgrid(theta, theta)
            return self.kuramoto_w + 2 / number_of_oscillators * np.sin(theta_j - theta_i).sum(0) + 0.2 * np.sin(theta)

        sol = spint.solve_ivp(kuramoto_ode, [0, self.kuramoto_t], self.kuramoto_x_0, method='BDF',
                              t_eval=np.linspace(0, self.kuramoto_t, self.kuramoto_m))
        self.kuramoto_x = sol.y
        self.kuramoto_y = np.zeros([number_of_oscillators, self.kuramoto_m])
        for i in range(self.kuramoto_m):
            self.kuramoto_y[:, i] = kuramoto_ode(0, self.kuramoto_x[:, i]) 

def initialize(self):
        self.options.declare('index', desc='the index of this segment in the parent phase.')

        self.options.declare('grid_data', desc='the grid data of the corresponding phase.')

        self.options.declare('method', default='RK45', values=('RK45', 'RK23', 'BDF', 'Radau'),
                             desc='The integrator used within scipy.integrate.solve_ivp. Currently '
                                  'supports \'RK45\', \'RK23\', and \'BDF\'.')

        self.options.declare('atol', default=1.0E-6, types=(float,),
                             desc='Absolute tolerance passed to scipy.integrate.solve_ivp.')

        self.options.declare('rtol', default=1.0E-6, types=(float,),
                             desc='Relative tolerance passed to scipy.integrate.solve_ivp.')

        self.options.declare('ode_class',
                             desc='System defining the ODE')

        self.options.declare('ode_init_kwargs', types=dict, default={},
                             desc='Keyword arguments provided when initializing the ODE System')

        self.options.declare('time_options', types=TimeOptionsDictionary,
                             desc='Time options for the phase')

        self.options.declare('state_options', types=dict,
                             desc='Dictionary of state names/options for the segments parent Phase')

        self.options.declare('control_options', default=None, types=dict, allow_none=True,
                             desc='Dictionary of control names/options for the segments parent Phase.')

        self.options.declare('polynomial_control_options', default=None, types=dict, allow_none=True,
                             desc='Dictionary of polynomial control names/options for the segments '
                                  'parent Phase.')

        self.options.declare('design_parameter_options', default=None, types=dict, allow_none=True,
                             desc='Dictionary of design parameter names/options for the segments '
                                  'parent Phase.')

        self.options.declare('input_parameter_options', default=None, types=dict, allow_none=True,
                             desc='Dictionary of input parameter names/options for the segments '
                                  'parent Phase.')

        self.options.declare('ode_integration_interface', default=None, allow_none=True,
                             types=ODEIntegrationInterface,
                             desc='The instance of the ODE integration interface used to provide '
                                  'the ODE to scipy.integrate.solve_ivp in the segment.  If None,'
                                  ' a new one will be instantiated for this segment.')

        self.options.declare('output_nodes_per_seg', default=None, types=(int,), allow_none=True,
                             desc='If None, results are provided at the all nodes within each'
                                  'segment.  If an int (n) then results are provided at n '
                                  'equally distributed points in time within each segment.')

        self.recording_options['options_excludes'] = ['ode_integration_interface'] 

def take_step(self, s, a):
        """Computes the next state from the current state and the action."""
        torque_1_action = a[0]
        # Clip the action to valid values.
        torque_1_action = np.clip(torque_1_action, self.action_range[0],
                                  self.action_range[1])

        torque_2_action = a[1]
        # Clip the action to valid values.
        torque_2_action = np.clip(torque_2_action, self.action_range[0],
                                  self.action_range[1])

        torque_3_action = a[2]
        # Clip the action to valid values.
        torque_3_action = np.clip(torque_3_action, self.action_range[0],
                                  self.action_range[1])

        # Add the action to the state so it can be passed to _dsdt.
        s_aug = np.append(
            s, np.array([torque_1_action, torque_2_action, torque_3_action]))

        ## Compute next state.
        # Type of integration and integration step.
        dt_in = self.time_interval
        if self.goal_range[0] < s[StateIndex.THETA_1] < self.goal_range[
            1] and self.goal_range[0] < s[StateIndex.THETA_1] + s[
            StateIndex.THETA_2] < self.goal_range[1] and self.goal_range[
            0] < s[StateIndex.THETA_1] + s[StateIndex.THETA_2] + s[
            StateIndex.THETA_3] < self.goal_range[1]:
            dt_in = self.time_interval
        else:
            dt_in = self.time_interval * 2.5

        ns = solve_ivp(self.ds_dt, (0., dt_in), s_aug).y[0:6, -1]

        # Project variables to valid space.
        theta_1 = wrap(ns[StateIndex.THETA_1], self.angle_1_range[0],
                       self.angle_1_range[1])
        ns[StateIndex.THETA_1] = np.clip(theta_1, self.angle_1_range[0],
                                       self.angle_1_range[1])
        ns[StateIndex.THETA_1_DOT] = np.clip(ns[StateIndex.THETA_1_DOT],
                                           self.angular_velocity_range[0],
                                           self.angular_velocity_range[1])

        theta_2 = wrap(ns[StateIndex.THETA_2], self.angle_2_range[0],
                       self.angle_2_range[1])
        ns[StateIndex.THETA_2] = np.clip(theta_2, self.angle_2_range[0],
                                       self.angle_2_range[1])
        ns[StateIndex.THETA_2_DOT] = np.clip(ns[StateIndex.THETA_2_DOT],
                                           self.angular_velocity_range[0],
                                           self.angular_velocity_range[1])

        theta_3 = wrap(ns[StateIndex.THETA_3], self.angle_3_range[0],
                       self.angle_3_range[1])
        ns[StateIndex.THETA_3] = np.clip(theta_3, self.angle_3_range[0],
                                       self.angle_3_range[1])
        ns[StateIndex.THETA_3_DOT] = np.clip(ns[StateIndex.THETA_3_DOT],
                                           self.angular_velocity_range[0],
                                           self.angular_velocity_range[1])

        return ns 

def test_integration():
    rtol = 1e-3
    atol = 1e-6
    y0 = [1/3, 2/9]

    for vectorized, method, t_span, jac in product(
            [False, True],
            ['RK23', 'RK45', 'Radau', 'BDF', 'LSODA'],
            [[5, 9], [5, 1]],
            [None, jac_rational, jac_rational_sparse]):

        if vectorized:
            fun = fun_rational_vectorized
        else:
            fun = fun_rational

        with suppress_warnings() as sup:
            sup.filter(UserWarning,
                       "The following arguments have no effect for a chosen solver: `jac`")
            res = solve_ivp(fun, t_span, y0, rtol=rtol,
                            atol=atol, method=method, dense_output=True,
                            jac=jac, vectorized=vectorized)
        assert_equal(res.t[0], t_span[0])
        assert_(res.t_events is None)
        assert_(res.success)
        assert_equal(res.status, 0)

        assert_(res.nfev < 40)

        if method in ['RK23', 'RK45', 'LSODA']:
            assert_equal(res.njev, 0)
            assert_equal(res.nlu, 0)
        else:
            assert_(0 < res.njev < 3)
            assert_(0 < res.nlu < 10)

        y_true = sol_rational(res.t)
        e = compute_error(res.y, y_true, rtol, atol)
        assert_(np.all(e < 5))

        tc = np.linspace(*t_span)
        yc_true = sol_rational(tc)
        yc = res.sol(tc)

        e = compute_error(yc, yc_true, rtol, atol)
        assert_(np.all(e < 5))

        tc = (t_span[0] + t_span[-1]) / 2
        yc_true = sol_rational(tc)
        yc = res.sol(tc)

        e = compute_error(yc, yc_true, rtol, atol)
        assert_(np.all(e < 5))

        # LSODA for some reasons doesn't pass the polynomial through the
        # previous points exactly after the order change. It might be some
        # bug in LSOSA implementation or maybe we missing something.
        if method != 'LSODA':
            assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)