a
    BCCfX                     @   s   d dl ZddlmZmZ ddlmZmZmZm	Z	m
Z
mZ ddlmZ dZdZdZd	d
 ZG dd deZG dd deZG dd deZG dd deZG dd deZG dd deZdS )    N   )	OdeSolverDenseOutput)validate_max_stepvalidate_tolselect_initial_stepnormwarn_extraneousvalidate_first_step)dop853_coefficientsg?皙?
   c	                 C   s   ||d< t t|dd |dd ddD ]H\}	\}
}t|d|	 j|
d|	 | }| |||  || ||	< q*||t|dd j|  }| || |}||d< ||fS )a8  Perform a single Runge-Kutta step.

    This function computes a prediction of an explicit Runge-Kutta method and
    also estimates the error of a less accurate method.

    Notation for Butcher tableau is as in [1]_.

    Parameters
    ----------
    fun : callable
        Right-hand side of the system.
    t : float
        Current time.
    y : ndarray, shape (n,)
        Current state.
    f : ndarray, shape (n,)
        Current value of the derivative, i.e., ``fun(x, y)``.
    h : float
        Step to use.
    A : ndarray, shape (n_stages, n_stages)
        Coefficients for combining previous RK stages to compute the next
        stage. For explicit methods the coefficients at and above the main
        diagonal are zeros.
    B : ndarray, shape (n_stages,)
        Coefficients for combining RK stages for computing the final
        prediction.
    C : ndarray, shape (n_stages,)
        Coefficients for incrementing time for consecutive RK stages.
        The value for the first stage is always zero.
    K : ndarray, shape (n_stages + 1, n)
        Storage array for putting RK stages here. Stages are stored in rows.
        The last row is a linear combination of the previous rows with
        coefficients

    Returns
    -------
    y_new : ndarray, shape (n,)
        Solution at t + h computed with a higher accuracy.
    f_new : ndarray, shape (n,)
        Derivative ``fun(t + h, y_new)``.

    References
    ----------
    .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
           Equations I: Nonstiff Problems", Sec. II.4.
    r   r   Nstart)	enumeratezipnpdotT)funtyfhABCKsacdyy_newf_new r%   S/var/www/html/django/DPS/env/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.pyrk_step   s    /."r'   c                       s   e Zd ZU dZeZejed< eZ	ejed< eZ
ejed< eZejed< eZejed< eZeed< eZeed< eZeed	< ejd
dddf fdd	Zdd Zdd Zdd Zdd Z  ZS )
RungeKuttaz,Base class for explicit Runge-Kutta methods.r   r   r   EPordererror_estimator_ordern_stagesMbP?ư>FNc
              	      s   t |
 t j|||||dd d | _t|| _t||| j\| _| _	| 
| j| j| _|	d u rt| j
| j| j| j| j| j| j| j	| _nt|	||| _tj| jd | jf| jjd| _d| jd  | _d | _d S )NT)Zsupport_complexr   dtyper   )r	   super__init__y_oldr   max_stepr   nrtolatolr   r   r   r   r   	directionr,   h_absr
   r   emptyr-   r1   r   error_exponent
h_previousselfr   t0Zy0t_boundr5   r7   r8   Z
vectorizedZ
first_stepZ
extraneous	__class__r%   r&   r3   U   s"    
 zRungeKutta.__init__c                 C   s   t |j| j| S N)r   r   r   r)   )r?   r   r   r%   r%   r&   _estimate_errori   s    zRungeKutta._estimate_errorc                 C   s   t | ||| S rD   )r   rE   )r?   r   r   scaler%   r%   r&   _estimate_error_norml   s    zRungeKutta._estimate_error_normc              
   C   s  | j }| j}| j}| j}| j}dtt|| jtj	 |  }| j
|krP|}n| j
|k r`|}n| j
}d}d}	|s||k rd| jfS || j }
||
 }| j|| j  dkr| j}|| }
t|
}t| j||| j|
| j| j| j| j	\}}|tt|t||  }| | j|
|}|dk rh|dkr6t}nttt|| j  }|	rZtd|}||9 }d}qn|ttt|| j  9 }d}	qn|
| _|| _|| _ || _|| _
|| _dS )Nr   Fr   r   T)TN)r   r   r5   r7   r8   r   abs	nextafterr9   infr:   ZTOO_SMALL_STEPrA   r'   r   r   r   r   r   r   maximumrG   
MAX_FACTORminSAFETYr<   max
MIN_FACTORr=   r4   )r?   r   r   r5   r7   r8   Zmin_stepr:   Zstep_acceptedZstep_rejectedr   Zt_newr#   r$   rF   Z
error_normfactorr%   r%   r&   
_step_implo   s`    "




 


zRungeKutta._step_implc                 C   s$   | j j| j}t| j| j| j|S rD   )r   r   r   r*   RkDenseOutputt_oldr   r4   )r?   Qr%   r%   r&   _dense_output_impl   s    zRungeKutta._dense_output_impl)__name__
__module____qualname____doc__NotImplementedr   r   Zndarray__annotations__r   r   r)   r*   r+   intr,   r-   rJ   r3   rE   rG   rR   rV   __classcell__r%   r%   rB   r&   r(   J   s"   
Cr(   c                   @   s   e Zd ZdZdZdZdZeg dZ	eg dg dg dgZ
eg dZeg d	Zeg d
g dg dg dgZdS )RK23a  Explicit Runge-Kutta method of order 3(2).

    This uses the Bogacki-Shampine pair of formulas [1]_. The error is controlled
    assuming accuracy of the second-order method, but steps are taken using the
    third-order accurate formula (local extrapolation is done). A cubic Hermite
    polynomial is used for the dense output.

    Can be applied in the complex domain.

    Parameters
    ----------
    fun : callable
        Right-hand side of the system: the time derivative of the state ``y``
        at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
        scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must
        return an array of the same shape as ``y``. See `vectorized` for more
        information.
    t0 : float
        Initial time.
    y0 : array_like, shape (n,)
        Initial state.
    t_bound : float
        Boundary time - the integration won't continue beyond it. It also
        determines the direction of the integration.
    first_step : float or None, optional
        Initial step size. Default is ``None`` which means that the algorithm
        should choose.
    max_step : float, optional
        Maximum allowed step size. Default is np.inf, i.e., the step size is not
        bounded and determined solely by the solver.
    rtol, atol : float and array_like, optional
        Relative and absolute tolerances. The solver keeps the local error
        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
        relative accuracy (number of correct digits), while `atol` controls
        absolute accuracy (number of correct decimal places). To achieve the
        desired `rtol`, set `atol` to be smaller than the smallest value that
        can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
        allowable error. If `atol` is larger than ``rtol * abs(y)`` the
        number of correct digits is not guaranteed. Conversely, to achieve the
        desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
        than `atol`. If components of y have different scales, it might be
        beneficial to set different `atol` values for different components by
        passing array_like with shape (n,) for `atol`. Default values are
        1e-3 for `rtol` and 1e-6 for `atol`.
    vectorized : bool, optional
        Whether `fun` may be called in a vectorized fashion. False (default)
        is recommended for this solver.

        If ``vectorized`` is False, `fun` will always be called with ``y`` of
        shape ``(n,)``, where ``n = len(y0)``.

        If ``vectorized`` is True, `fun` may be called with ``y`` of shape
        ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave
        such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
        the returned array is the time derivative of the state corresponding
        with a column of ``y``).

        Setting ``vectorized=True`` allows for faster finite difference
        approximation of the Jacobian by methods 'Radau' and 'BDF', but
        will result in slower execution for this solver.

    Attributes
    ----------
    n : int
        Number of equations.
    status : string
        Current status of the solver: 'running', 'finished' or 'failed'.
    t_bound : float
        Boundary time.
    direction : float
        Integration direction: +1 or -1.
    t : float
        Current time.
    y : ndarray
        Current state.
    t_old : float
        Previous time. None if no steps were made yet.
    step_size : float
        Size of the last successful step. None if no steps were made yet.
    nfev : int
        Number evaluations of the system's right-hand side.
    njev : int
        Number of evaluations of the Jacobian.
        Is always 0 for this solver as it does not use the Jacobian.
    nlu : int
        Number of LU decompositions. Is always 0 for this solver.

    References
    ----------
    .. [1] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
           Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
          )r         ?      ?)r   r   r   )rb   r   r   )r   rc   r   )gqq?gUUUUUU?gqq?)grqǱ?gUUUUUUgqqg      ?)r   gUUUUUUgrq?)r   r   gUUUUUU)r   gUUUUUU?gqq)r   r   r   NrW   rX   rY   rZ   r+   r,   r-   r   arrayr   r   r   r)   r*   r%   r%   r%   r&   r_      s"   \
r_   c                
   @   s   e Zd ZdZdZdZdZeg dZ	eg dg dg dg d	g d
g dgZ
eg dZeg dZeg dg dg dg dg dg dg dgZdS )RK45a  Explicit Runge-Kutta method of order 5(4).

    This uses the Dormand-Prince pair of formulas [1]_. The error is controlled
    assuming accuracy of the fourth-order method accuracy, but steps are taken
    using the fifth-order accurate formula (local extrapolation is done).
    A quartic interpolation polynomial is used for the dense output [2]_.

    Can be applied in the complex domain.

    Parameters
    ----------
    fun : callable
        Right-hand side of the system. The calling signature is ``fun(t, y)``.
        Here ``t`` is a scalar, and there are two options for the ndarray ``y``:
        It can either have shape (n,); then ``fun`` must return array_like with
        shape (n,). Alternatively it can have shape (n, k); then ``fun``
        must return an array_like with shape (n, k), i.e., each column
        corresponds to a single column in ``y``. The choice between the two
        options is determined by `vectorized` argument (see below).
    t0 : float
        Initial time.
    y0 : array_like, shape (n,)
        Initial state.
    t_bound : float
        Boundary time - the integration won't continue beyond it. It also
        determines the direction of the integration.
    first_step : float or None, optional
        Initial step size. Default is ``None`` which means that the algorithm
        should choose.
    max_step : float, optional
        Maximum allowed step size. Default is np.inf, i.e., the step size is not
        bounded and determined solely by the solver.
    rtol, atol : float and array_like, optional
        Relative and absolute tolerances. The solver keeps the local error
        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
        relative accuracy (number of correct digits), while `atol` controls
        absolute accuracy (number of correct decimal places). To achieve the
        desired `rtol`, set `atol` to be smaller than the smallest value that
        can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
        allowable error. If `atol` is larger than ``rtol * abs(y)`` the
        number of correct digits is not guaranteed. Conversely, to achieve the
        desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
        than `atol`. If components of y have different scales, it might be
        beneficial to set different `atol` values for different components by
        passing array_like with shape (n,) for `atol`. Default values are
        1e-3 for `rtol` and 1e-6 for `atol`.
    vectorized : bool, optional
        Whether `fun` is implemented in a vectorized fashion. Default is False.

    Attributes
    ----------
    n : int
        Number of equations.
    status : string
        Current status of the solver: 'running', 'finished' or 'failed'.
    t_bound : float
        Boundary time.
    direction : float
        Integration direction: +1 or -1.
    t : float
        Current time.
    y : ndarray
        Current state.
    t_old : float
        Previous time. None if no steps were made yet.
    step_size : float
        Size of the last successful step. None if no steps were made yet.
    nfev : int
        Number evaluations of the system's right-hand side.
    njev : int
        Number of evaluations of the Jacobian.
        Is always 0 for this solver as it does not use the Jacobian.
    nlu : int
        Number of LU decompositions. Is always 0 for this solver.

    References
    ----------
    .. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
           formulae", Journal of Computational and Applied Mathematics, Vol. 6,
           No. 1, pp. 19-26, 1980.
    .. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
           of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
             )r   r   g333333?g?gqq?r   )r   r   r   r   r   )r   r   r   r   r   )g333333?g?r   r   r   )gII?ggqq@r   r   )gq@g 1'gR<6R#@gE3ҿr   )g+@g>%gr!@gE]t?g/pѿ)gUUUUUU?r   gVI?gUUUUU?gϡԿg10?)g2Tr   gĿ
UZkq?ggX
?g{tg?)r   g#
!gJ<@gFC)r   r   r   r   )r   gF@gFj'NgDg@)r   gdDgaP#$@g2)r   g<p@g@갘g,@)r   gRq#g_40g.
@gF)r   g'?g'gK@Nrd   r%   r%   r%   r&   rf   %  s0   Srf   c                       s   e Zd ZdZejZdZdZej	dedef Z	ej
Z
ejde ZejZejZejZej	ed d Zejed d Zejddddf fd	d
	Zdd Zdd Zdd Z  ZS )DOP853a"  Explicit Runge-Kutta method of order 8.

    This is a Python implementation of "DOP853" algorithm originally written
    in Fortran [1]_, [2]_. Note that this is not a literal translation, but
    the algorithmic core and coefficients are the same.

    Can be applied in the complex domain.

    Parameters
    ----------
    fun : callable
        Right-hand side of the system. The calling signature is ``fun(t, y)``.
        Here, ``t`` is a scalar, and there are two options for the ndarray ``y``:
        It can either have shape (n,); then ``fun`` must return array_like with
        shape (n,). Alternatively it can have shape (n, k); then ``fun``
        must return an array_like with shape (n, k), i.e. each column
        corresponds to a single column in ``y``. The choice between the two
        options is determined by `vectorized` argument (see below).
    t0 : float
        Initial time.
    y0 : array_like, shape (n,)
        Initial state.
    t_bound : float
        Boundary time - the integration won't continue beyond it. It also
        determines the direction of the integration.
    first_step : float or None, optional
        Initial step size. Default is ``None`` which means that the algorithm
        should choose.
    max_step : float, optional
        Maximum allowed step size. Default is np.inf, i.e. the step size is not
        bounded and determined solely by the solver.
    rtol, atol : float and array_like, optional
        Relative and absolute tolerances. The solver keeps the local error
        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
        relative accuracy (number of correct digits), while `atol` controls
        absolute accuracy (number of correct decimal places). To achieve the
        desired `rtol`, set `atol` to be smaller than the smallest value that
        can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
        allowable error. If `atol` is larger than ``rtol * abs(y)`` the
        number of correct digits is not guaranteed. Conversely, to achieve the
        desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
        than `atol`. If components of y have different scales, it might be
        beneficial to set different `atol` values for different components by
        passing array_like with shape (n,) for `atol`. Default values are
        1e-3 for `rtol` and 1e-6 for `atol`.
    vectorized : bool, optional
        Whether `fun` is implemented in a vectorized fashion. Default is False.

    Attributes
    ----------
    n : int
        Number of equations.
    status : string
        Current status of the solver: 'running', 'finished' or 'failed'.
    t_bound : float
        Boundary time.
    direction : float
        Integration direction: +1 or -1.
    t : float
        Current time.
    y : ndarray
        Current state.
    t_old : float
        Previous time. None if no steps were made yet.
    step_size : float
        Size of the last successful step. None if no steps were made yet.
    nfev : int
        Number evaluations of the system's right-hand side.
    njev : int
        Number of evaluations of the Jacobian. Is always 0 for this solver
        as it does not use the Jacobian.
    nlu : int
        Number of LU decompositions. Is always 0 for this solver.

    References
    ----------
    .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
           Equations I: Nonstiff Problems", Sec. II.
    .. [2] `Page with original Fortran code of DOP853
            <http://www.unige.ch/~hairer/software.html>`_.
          Nr   r.   r/   Fc
              
      sZ   t  j|||||||||	f	i |
 tjtj| jf| jjd| _	| j	d | j
d  | _d S )Nr0   r   )r2   r3   r   r;   r   ZN_STAGES_EXTENDEDr6   r   r1   
K_extendedr-   r   r>   rB   r%   r&   r3     s    zDOP853.__init__c                 C   st   t |j| j}t |j| j}t t |dt | }t |}|dk}t || ||  ||< || | S )Ng?r   )r   r   r   E5E3hypotrH   Z	ones_like)r?   r   r   err5err3denomZcorrection_factormaskr%   r%   r&   rE     s    
zDOP853._estimate_errorc           	      C   s   t |j| j| }t |j| j| }t j|d }t j|d }|dkr\|dkr\dS |d|  }t || t |t	|  S )Nra   r   g        g{Gz?)
r   r   r   rn   ro   Zlinalgr   rH   sqrtlen)	r?   r   r   rF   rq   rr   Zerr5_norm_2Zerr3_norm_2rs   r%   r%   r&   rG     s    zDOP853._estimate_error_normc           
      C   s  | j }| j}tt| j| j| jd dD ]N\}\}}t|d | j	|d | | }| 
| j||  | j| ||< q(tjtj| jf| jjd}|d }| j| j }	|	|d< || |	 |d< d|	 || j|   |d< |t| j| |dd < t| j| j| j|S )Nr   r   r0   r   ra   r`   )rm   r=   r   r   A_EXTRAC_EXTRAr-   r   r   r   r   rT   r4   r;   r   ZINTERPOLATOR_POWERr6   r1   r   r   DDop853DenseOutputr   )
r?   r   r   r   r    r!   r"   FZf_oldZdelta_yr%   r%   r&   rV     s"    ""zDOP853._dense_output_impl)rW   rX   rY   rZ   r   ZN_STAGESr-   r+   r,   r   r   r   ro   rn   ry   rw   rx   r   rJ   r3   rE   rG   rV   r^   r%   r%   rB   r&   rj     s&   Q		
rj   c                       s$   e Zd Z fddZdd Z  ZS )rS   c                    s8   t  || || | _|| _|jd d | _|| _d S )Nr   )r2   r3   r   rU   shaper+   r4   )r?   rT   r   r4   rU   rB   r%   r&   r3   )  s
    
zRkDenseOutput.__init__c                 C   s   || j  | j }|jdkr8t|| jd }t|}n$t|| jd df}tj|dd}| jt| j| }|jdkr|| j	d d d f 7 }n
|| j	7 }|S )Nr   r   )Zaxisra   )
rT   r   ndimr   Ztiler+   Zcumprodr   rU   r4   )r?   r   xpr   r%   r%   r&   
_call_impl0  s    


zRkDenseOutput._call_implrW   rX   rY   r3   r   r^   r%   r%   rB   r&   rS   (  s   rS   c                       s$   e Zd Z fddZdd Z  ZS )rz   c                    s(   t  || || | _|| _|| _d S rD   )r2   r3   r   r{   r4   )r?   rT   r   r4   r{   rB   r%   r&   r3   B  s    
zDop853DenseOutput.__init__c                 C   s   || j  | j }|jdkr(t| j}n0|d d d f }tjt|t| jf| jjd}t	t
| jD ]2\}}||7 }|d dkr||9 }qf|d| 9 }qf|| j7 }|jS )Nr   r0   ra   r   )rT   r   r}   r   Z
zeros_liker4   Zzerosrv   r1   r   reversedr{   r   )r?   r   r~   r   ir   r%   r%   r&   r   H  s    
 

zDop853DenseOutput._call_implr   r%   r%   rB   r&   rz   A  s   rz   )numpyr   baser   r   commonr   r   r   r   r	   r
    r   rN   rP   rL   r'   r(   r_   rf   rj   rS   rz   r%   r%   r%   r&   <module>   s    <mnr 