Python integrate.solve_ivp方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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))) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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)) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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)) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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]]) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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..."""
        ) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import solve_ivp [as 别名]
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)