In [1]:

    

%matplotlib inline
%load_ext autoreload
%autoreload 2
import matplotlib
matplotlib.rcParams['text.latex.unicode'] = True


    

    




/Users/jslater/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py:582: UserWarning: matplotlibrc text.usetex can not be used unless ghostscript-8.60 or later is installed on your system
  % gs_req)

In [2]:

    

# -*- coding: utf-8 -*-
"""
Created on Wed Mar  7 14:10:25 2018

@author: kirby
"""


import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import mousai as ms
from scipy import pi, sin
from scipy.fftpack import *

plt.close('all')
def duff_osc_ss(x, params):
    omega = params['omega']
    t = params['cur_time']
    xdot = np.array([[x[1]],[-x[0]-.1*x[0]**3-.1*x[1]+1*sin(omega*t)]])
    return xdot

def duff_osc2(x, v, params):
    omega = params['omega']
    t = params['cur_time']
    return np.array([[-x-0.1*x**3-.1*v+1*sin(omega*t)]])

print('second order')
tfso, xfso, e, amps, phases = ms.hb_freq(duff_osc2, sp.array([[.1,-.1,0]]), omega = 0.1, \
                            eqform='second_order', num_harmonics=1, num_time_steps = 51, f_tol = 1e-10)

plt.figure()
#plt.plot(time, xc.T, t, x.T, '*')
time, xc = ms.time_history(tfso,xfso)
plt.plot(time, xc.T, label = 'second order 1 harmonic')


print('on to state space')
t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), omega = 0.1, \
                            eqform='first_order', num_harmonics=1, num_time_steps = 51, f_tol = 1e-7)

#print(x,e)
print(e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
time, xc = ms.time_history(t,x)


#plt.plot(time, xc.T, t, x.T, '*')
plt.plot(time, xc.T, label = 'state space 1 harmonic')


print('on to state space')
t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), omega = 0.1, \
                            eqform='first_order', num_harmonics=3, num_time_steps = 51, f_tol = 1e-7)

print(e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
time, xc = ms.time_history(t,x)

plt.plot(time, xc.T, label = '3 harmonics')



t, x, e, amps, phases = ms.hb_time(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), omega = 0.1, eqform='first_order', 
                                   num_harmonics=5)

print(x,e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
time, xc = ms.time_history(t,x)

#plt.figure()
#plt.plot(time, xc.T, t, x.T, '*')
plt.plot(time, xc.T,'-.', label = '5 harmonics')
#plt.plot(time, xc.T, '3 harmonics')

time, xc = ms.time_history(tfso,xfso)
plt.plot(time, xc.T, '--', label = 'second order 1 harmonic')
#plt.grid(True)
#plt.legend(['hb_freq Displacement','hb_freq Velocity','hb_time Displacement','hb_time Velocity'])
plt.legend()
plt.grid()


    

    




second order
on to state space
[[ 6.11014807e-12 -4.09835182e-12]
 [ 1.42057024e-12 -7.02622938e-12]]
Constant term of FFT of signal should be zero:  (1.9081958235744878e-17+0j)
on to state space
[[-2.14516624e-12 -4.96904744e-12  3.85635721e-14 -6.51966255e-15
  -2.41996735e-13 -5.27259020e-13]
 [-1.56015062e-12 -1.02396267e-11 -9.73592067e-15 -2.05416694e-13
  -7.08425755e-13 -1.89628532e-12]]
Constant term of FFT of signal should be zero:  (-7.632783294297951e-17+0j)
[[-0.01039999  0.52526596  0.85219917  0.92080925  0.73014714  0.29325872
  -0.27439117 -0.72047079 -0.91897883 -0.85775937 -0.5396823 ]
 [ 0.10151528  0.07897289  0.03477871 -0.0108248  -0.05572575 -0.09438379
  -0.09561399 -0.05728763 -0.0121241   0.03345823  0.07723494]] [[ 1.28577704e-13 -4.66016115e-14 -2.08721929e-14 -9.70386965e-14
  -4.72538675e-15  4.86277685e-14 -7.23032745e-15 -6.25471896e-14
  -9.42822209e-14  5.72805692e-14 -3.35009798e-14]
 [-3.41100845e-15  1.38338126e-13 -5.91036421e-11 -2.92572452e-10
  -8.80753803e-13  1.00032829e-14  1.70332498e-14  1.09752658e-12
   2.84168106e-10  5.94388774e-11  3.53379652e-13]]
Constant term of FFT of signal should be zero:  (-2.214069885344422e-06+0j)






    

In [3]:

    

t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), .1, \
                            eqform='first_order', num_harmonics=1, num_time_steps = 51, f_tol = 1e-7)
time, xc = ms.time_history(t,x)
plt.plot(time, xc.T)
plt.grid(True)


    

In [4]:

    

e


    

    
Out[4]:





array([[ 6.11014807e-12, -4.09835182e-12],
       [ 1.42057024e-12, -7.02622938e-12]])

In [5]:

    

tfso, xfso, e, amps, phases = ms.hb_freq(duff_osc2, sp.array([[.1,-.1,0]]), .1, \
                            eqform='second_order', num_harmonics=1, num_time_steps = 51)
print('Equation errors (should be zero): ', e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])


    

    




Equation errors (should be zero):  [[-4.51508908e-10 -3.68637431e-09]]
Constant term of FFT of signal should be zero:  (1.9081958235744878e-17+0j)

In [6]:

    

tfso, xfso, e, amps, phases = ms.hb_freq(duff_osc2, num_variables=1, omega = 0.2, \
                            eqform='second_order', num_harmonics=5, mask_constant=False, num_time_steps = 501, f_tol = 1e-10)
print('Equation errors (should be zero): ', e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])


    

    




Equation errors (should be zero):  [[-4.21601663e-11  2.06073875e-11  3.91625760e-11  1.49842769e-11
  -2.75224080e-11 -2.00321577e-11 -5.70806888e-11 -7.32032897e-12
   1.69858586e-11  2.56796522e-13  1.09424955e-11]]
Constant term of FFT of signal should be zero:  (1.9081958235744878e-17+0j)

In [7]:

    

e


    

    
Out[7]:





array([[-4.21601663e-11,  2.06073875e-11,  3.91625760e-11,
         1.49842769e-11, -2.75224080e-11, -2.00321577e-11,
        -5.70806888e-11, -7.32032897e-12,  1.69858586e-11,
         2.56796522e-13,  1.09424955e-11]])

In [8]:

    

xfso


    

    
Out[8]:





array([[-0.03242912,  0.54981788,  0.86040844,  0.94737891,  0.74788442,
         0.32695025, -0.2816328 , -0.73791066, -0.93460228, -0.88170109,
        -0.56416394]])

In [9]:

    

time, xc = ms.time_history(tfso,xfso)

plt.plot(time,xc.T)
plt.grid(True)
plt.title('Second order with 3 harmonics')


    

    
Out[9]:





Text(0.5,1,'Second order with 3 harmonics')






    

In [10]:

    

t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0,1,-1],[.1,-.1,0]]), .1, \
                            eqform='first_order', num_harmonics=1, num_time_steps = 51)#, f_tol = 1e-12)
print('x = ', x)


    

    




x =  [[-0.0089479   0.82363988 -0.81469198]
 [ 0.09458913 -0.04651966 -0.04806948]]

In [11]:

    

e


    

    
Out[11]:





array([[ 6.11014807e-12, -4.09835182e-12],
       [ 1.42057024e-12, -7.02622938e-12]])

In [12]:

    

time, xc = ms.time_history(t,x)

plt.plot(time,xc.T)
plt.grid(True)
plt.title('State space with 3 harmonics')


    

    
Out[12]:





Text(0.5,1,'State space with 3 harmonics')






    

In [13]:

    

t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, num_variables = 2, omega = 0.1,
                                   eqform='first_order', num_harmonics=5, num_time_steps=51, mask_constant=True, f_rtol = 1e-8)#, f_rtol=1e-8)
time, xc = ms.time_history(t, x)
plt.plot(time, xc.T, label='5 harmonics')


    

    
Out[13]:





[<matplotlib.lines.Line2D at 0x121dd26d8>,
 <matplotlib.lines.Line2D at 0x121dd2898>]






    

In [14]:

    

t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, x0 = x, omega = 0.1,
                                   eqform='first_order', num_harmonics=5, num_time_steps=501, mask_constant=True)#, f_rtol = 1e-8)#, f_rtol=1e-8)
time, xc = ms.time_history(t, x)
plt.plot(time, xc.T, label='5 harmonics')


    

    




/Users/jslater/anaconda3/lib/python3.6/site-packages/scipy/optimize/nonlin.py:476: RuntimeWarning: invalid value encountered in double_scalars
  and dx_norm/self.x_rtol <= x_norm))






    
Out[14]:





[<matplotlib.lines.Line2D at 0x121f176d8>,
 <matplotlib.lines.Line2D at 0x121f17898>]






    

In [15]:

    

e[0,:]


    

    
Out[15]:





array([ 2.97210699e-16,  2.99660039e-18,  3.00766780e-16, -5.46025033e-17,
        2.01847324e-16,  6.24921018e-17, -3.24441405e-16,  1.35786207e-16,
       -3.23564866e-16, -2.18315442e-16])

In [16]:

    

print('ss 1')
t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0, 1, -1], [.1, -.1, 0]]), .1,
                                   eqform='first_order', num_harmonics=1, num_time_steps=501)
time, xc = ms.time_history(t, x)
plt.plot(time, xc.T, label='1 harmonics')

print('ss 3')
t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, sp.array([[0, 1, -1], [.1, -.1, 0]]), .1,
                                   eqform='first_order', num_harmonics=3, num_time_steps=501)  # , f_tol = 1e-8)
time, xc = ms.time_history(t, x)
print(e)
plt.plot(time, xc.T, label='3 harmonics')

print('ss 5')
t, x, e, amps, phases = ms.hb_freq(duff_osc_ss, x, .1,
                                   eqform='first_order', num_harmonics=5, num_time_steps=501, mask_constant=True, f_tol=1e-10)

print(e)
time, xc = ms.time_history(t, x)
plt.plot(time, xc.T, label='5 harmonics')

plt.legend()


    

    




ss 1
ss 3
[[-2.14311653e-12 -4.96362386e-12  3.00848948e-14 -4.12724722e-15
  -2.95001565e-13 -5.55438119e-13]
 [-1.47733220e-12 -1.06241971e-11 -4.67820026e-14 -1.24719700e-13
  -7.98423837e-13 -1.84443363e-12]]
ss 5
[[-1.34859807e-14 -7.25206451e-16  7.08049606e-15  2.67639781e-15
   4.90583682e-15  3.45813679e-15  4.89511001e-15 -2.87054656e-15
   7.11233377e-15 -4.88032750e-15]
 [-1.07502715e-14 -6.67308769e-13  1.33500129e-14  2.37640135e-14
   1.52139239e-14 -1.17099220e-13 -1.31167025e-14 -1.04777743e-14
  -1.14137115e-14 -1.00101420e-14]]






    
Out[16]:





<matplotlib.legend.Legend at 0x121feacf8>






    

In [17]:

    

plt.plot(np.abs(fft(x).T))


    

    
Out[17]:





[<matplotlib.lines.Line2D at 0x122076748>,
 <matplotlib.lines.Line2D at 0x122076898>]






    

In [18]:

    

e


    

    
Out[18]:





array([[-1.34859807e-14, -7.25206451e-16,  7.08049606e-15,
         2.67639781e-15,  4.90583682e-15,  3.45813679e-15,
         4.89511001e-15, -2.87054656e-15,  7.11233377e-15,
        -4.88032750e-15],
       [-1.07502715e-14, -6.67308769e-13,  1.33500129e-14,
         2.37640135e-14,  1.52139239e-14, -1.17099220e-13,
        -1.31167025e-14, -1.04777743e-14, -1.14137115e-14,
        -1.00101420e-14]])

In [19]:

    

X = fft(x)
X[:,2:9]= 0*X[:,2:9]
x = ifft(X)


    

In [20]:

    

x


    

    
Out[20]:





array([[-0.0090069 +9.63875444e-16j,  0.5044076 +9.00765094e-16j,
         0.85767626-2.15492758e-16j,  0.93863876-1.07980295e-15j,
         0.72159009-6.81639953e-16j,  0.27544167+5.13476008e-16j,
        -0.25815754+1.10825124e-15j, -0.70979356+4.07292397e-16j,
        -0.93607513-7.69860485e-16j, -0.86515946-1.04691560e-15j,
        -0.51956178-9.99484334e-17j],
       [ 0.0946996 -2.35961818e-17j,  0.08015332-2.30142801e-16j,
         0.04015893-1.67613368e-16j, -0.01258564+9.08845826e-17j,
        -0.06133435+2.43123008e-16j, -0.09060985+1.11109312e-16j,
        -0.09111735-1.50810055e-16j, -0.06269575-2.36406834e-16j,
        -0.01436868-4.56038413e-17j,  0.03852034+1.98517794e-16j,
         0.07917942+2.10538385e-16j]])

In [21]:

    

time, xc = ms.time_history(t, x)
plt.plot(time, xc.T, label='5 harmonics')

plt.legend()


    

    
Out[21]:





<matplotlib.legend.Legend at 0x1220ee8d0>






    

In [22]:

    

e


    

    
Out[22]:





array([[-1.34859807e-14, -7.25206451e-16,  7.08049606e-15,
         2.67639781e-15,  4.90583682e-15,  3.45813679e-15,
         4.89511001e-15, -2.87054656e-15,  7.11233377e-15,
        -4.88032750e-15],
       [-1.07502715e-14, -6.67308769e-13,  1.33500129e-14,
         2.37640135e-14,  1.52139239e-14, -1.17099220e-13,
        -1.31167025e-14, -1.04777743e-14, -1.14137115e-14,
        -1.00101420e-14]])

In [23]:

    

np.max(e)


    

    
Out[23]:





2.3764013454246053e-14

In [24]:

    

import scipy.integrate as odeint


    

In [25]:

    

help(odeint.solve_ivp)


    

    




Help on function solve_ivp in module scipy.integrate._ivp.ivp:

solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)
    Solve an initial value problem for a system of ODEs.
    
    This function numerically integrates a system of ordinary differential
    equations given an initial value::
    
        dy / dt = f(t, y)
        y(t0) = y0
    
    Here t is a one-dimensional independent variable (time), y(t) is an
    n-dimensional vector-valued function (state), and an n-dimensional
    vector-valued function f(t, y) determines the differential equations.
    The goal is to find y(t) approximately satisfying the differential
    equations, given an initial value y(t0)=y0.
    
    Some of the solvers support integration in the complex domain, but note that
    for stiff ODE solvers, the right-hand side must be complex-differentiable
    (satisfy Cauchy-Riemann equations [11]_). To solve a problem in the complex
    domain, pass y0 with a complex data type. Another option is always to
    rewrite your problem for real and imaginary parts separately.
    
    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). The
        vectorized implementation allows a faster approximation of the Jacobian
        by finite differences (required for stiff solvers).
    t_span : 2-tuple of floats
        Interval of integration (t0, tf). The solver starts with t=t0 and
        integrates until it reaches t=tf.
    y0 : array_like, shape (n,)
        Initial state. For problems in the complex domain, pass `y0` with a
        complex data type (even if the initial guess is purely real).
    method : string or `OdeSolver`, optional
        Integration method to use:
    
            * 'RK45' (default): Explicit Runge-Kutta method of order 5(4) [1]_.
              The error is controlled assuming accuracy of the fourth-order
              method, 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.
            * 'RK23': Explicit Runge-Kutta method of order 3(2) [3]_. 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.
            * 'Radau': Implicit Runge-Kutta method of the Radau IIA family of
              order 5 [4]_. The error is controlled with a third-order accurate
              embedded formula. A cubic polynomial which satisfies the
              collocation conditions is used for the dense output.
            * 'BDF': Implicit multi-step variable-order (1 to 5) method based
              on a backward differentiation formula for the derivative
              approximation [5]_. The implementation follows the one described
              in [6]_. A quasi-constant step scheme is used and accuracy is
              enhanced using the NDF modification. Can be applied in the complex
              domain.
            * 'LSODA': Adams/BDF method with automatic stiffness detection and
              switching [7]_, [8]_. This is a wrapper of the Fortran solver
              from ODEPACK.
    
        You should use the 'RK45' or 'RK23' method for non-stiff problems and
        'Radau' or 'BDF' for stiff problems [9]_. If not sure, first try to run
        'RK45'. If needs unusually many iterations, diverges, or fails, your
        problem is likely to be stiff and you should use 'Radau' or 'BDF'.
        'LSODA' can also be a good universal choice, but it might be somewhat
        less convenient to work with as it wraps old Fortran code.
    
        You can also pass an arbitrary class derived from `OdeSolver` which
        implements the solver.
    dense_output : bool, optional
        Whether to compute a continuous solution. Default is False.
    t_eval : array_like or None, optional
        Times at which to store the computed solution, must be sorted and lie
        within `t_span`. If None (default), use points selected by the solver.
    events : callable, list of callables or None, optional
        Types of events to track. Each is defined by a continuous function of
        time and state that becomes zero value in case of an event. Each function
        must have the signature ``event(t, y)`` and return a float. The solver will
        find an accurate value of ``t`` at which ``event(t, y(t)) = 0`` using a
        root-finding algorithm. Additionally each ``event`` function might have
        the following attributes:
    
            * terminal: bool, whether to terminate integration if this
              event occurs. Implicitly False if not assigned.
            * direction: float, direction of a zero crossing. If `direction`
              is positive, `event` must go from negative to positive, and
              vice versa if `direction` is negative. If 0, then either direction
              will count. Implicitly 0 if not assigned.
    
        You can assign attributes like ``event.terminal = True`` to any
        function in Python. If None (default), events won't be tracked.
    vectorized : bool, optional
        Whether `fun` is implemented in a vectorized fashion. Default is False.
    options
        Options passed to a chosen solver. All options available for already
        implemented solvers are listed below.
    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). But if a component of `y`
        is approximately below `atol`, the error only needs to fall within
        the same `atol` threshold, and the number of correct digits is not
        guaranteed. 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`.
    jac : {None, array_like, sparse_matrix, callable}, optional
        Jacobian matrix of the right-hand side of the system with respect to
        y, required by the 'Radau', 'BDF' and 'LSODA' method. The Jacobian matrix
        has shape (n, n) and its element (i, j) is equal to ``d f_i / d y_j``.
        There are three ways to define the Jacobian:
    
            * If array_like or sparse_matrix, the Jacobian is assumed to
              be constant. Not supported by 'LSODA'.
            * If callable, the Jacobian is assumed to depend on both
              t and y; it will be called as ``jac(t, y)`` as necessary.
              For the 'Radau' and 'BDF' methods, the return value might be a
              sparse matrix.
            * If None (default), the Jacobian will be approximated by
              finite differences.
    
        It is generally recommended to provide the Jacobian rather than
        relying on a finite-difference approximation.
    jac_sparsity : {None, array_like, sparse matrix}, optional
        Defines a sparsity structure of the Jacobian matrix for a
        finite-difference approximation. Its shape must be (n, n). This argument
        is ignored if `jac` is not `None`. If the Jacobian has only few non-zero
        elements in *each* row, providing the sparsity structure will greatly
        speed up the computations [10]_. A zero entry means that a corresponding
        element in the Jacobian is always zero. If None (default), the Jacobian
        is assumed to be dense.
        Not supported by 'LSODA', see `lband` and `uband` instead.
    lband, uband : int or None
        Parameters defining the bandwidth of the Jacobian for the 'LSODA' method,
        i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting
        these requires your jac routine to return the Jacobian in the packed format:
        the returned array must have ``n`` columns and ``uband + lband + 1``
        rows in which Jacobian diagonals are written. Specifically
        ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used
        in `scipy.linalg.solve_banded` (check for an illustration).
        These parameters can be also used with ``jac=None`` to reduce the
        number of Jacobian elements estimated by finite differences.
    min_step, first_step : float, optional
        The minimum allowed step size and the initial step size respectively
        for 'LSODA' method. By default `min_step` is zero and `first_step` is
        selected automatically.
    
    Returns
    -------
    Bunch object with the following fields defined:
    t : ndarray, shape (n_points,)
        Time points.
    y : ndarray, shape (n, n_points)
        Values of the solution at `t`.
    sol : `OdeSolution` or None
        Found solution as `OdeSolution` instance; None if `dense_output` was
        set to False.
    t_events : list of ndarray or None
        Contains for each event type a list of arrays at which an event of
        that type event was detected. None if `events` was None.
    nfev : int
        Number of evaluations of the right-hand side.
    njev : int
        Number of evaluations of the Jacobian.
    nlu : int
        Number of LU decompositions.
    status : int
        Reason for algorithm termination:
    
            * -1: Integration step failed.
            *  0: The solver successfully reached the end of `tspan`.
            *  1: A termination event occurred.
    
    message : string
        Human-readable description of the termination reason.
    success : bool
        True if the solver reached the interval end or a termination event
        occurred (``status >= 0``).
    
    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.
    .. [3] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
           Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
    .. [4] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
           Stiff and Differential-Algebraic Problems", Sec. IV.8.
    .. [5] `Backward Differentiation Formula
            <https://en.wikipedia.org/wiki/Backward_differentiation_formula>`_
            on Wikipedia.
    .. [6] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
           COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
    .. [7] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
           Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
           pp. 55-64, 1983.
    .. [8] L. Petzold, "Automatic selection of methods for solving stiff and
           nonstiff systems of ordinary differential equations", SIAM Journal
           on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
           1983.
    .. [9] `Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>`_ on
           Wikipedia.
    .. [10] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
            sparse Jacobian matrices", Journal of the Institute of Mathematics
            and its Applications, 13, pp. 117-120, 1974.
    .. [11] `Cauchy-Riemann equations
             <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
             Wikipedia.
    
    Examples
    --------
    Basic exponential decay showing automatically chosen time points.
    
    >>> from scipy.integrate import solve_ivp
    >>> def exponential_decay(t, y): return -0.5 * y
    >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
    >>> print(sol.t)
    [  0.           0.11487653   1.26364188   3.06061781   4.85759374
       6.65456967   8.4515456   10.        ]
    >>> print(sol.y)
    [[2.         1.88836035 1.06327177 0.43319312 0.17648948 0.0719045
      0.02929499 0.01350938]
     [4.         3.7767207  2.12654355 0.86638624 0.35297895 0.143809
      0.05858998 0.02701876]
     [8.         7.5534414  4.25308709 1.73277247 0.7059579  0.287618
      0.11717996 0.05403753]]
    
    Specifying points where the solution is desired.
    
    >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
    ...                 t_eval=[0, 1, 2, 4, 10])
    >>> print(sol.t)
    [ 0  1  2  4 10]
    >>> print(sol.y)
    [[2.         1.21305369 0.73534021 0.27066736 0.01350938]
     [4.         2.42610739 1.47068043 0.54133472 0.02701876]
     [8.         4.85221478 2.94136085 1.08266944 0.05403753]]
    
    Cannon fired upward with terminal event upon impact. The ``terminal`` and
    ``direction`` fields of an event are applied by monkey patching a function.
    Here ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts at
    position 0 with velocity +10. Note that the integration never reaches t=100
    because the event is terminal.
    
    >>> def upward_cannon(t, y): return [y[1], -0.5]
    >>> def hit_ground(t, y): return y[1]
    >>> hit_ground.terminal = True
    >>> hit_ground.direction = -1
    >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
    >>> print(sol.t_events)
    [array([ 20.])]
    >>> print(sol.t)
    [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
     1.11088891e-01 1.11098890e+00 1.11099890e+01 2.00000000e+01]

In [26]:

    

def duff_osc_ss(t, x):
    omega = 0.1    
    xdot = np.array([x[1],-x[0]-.1*x[0]**3-.1*x[1]+1*sin(omega*t)])
    return xdot


    

In [27]:

    

sol = odeint.solve_ivp(fun = duff_osc_ss, t_span=(0,600), y0 = np.array([0,0]))


    

In [28]:

    

plt.plot(sol.t,sol.y.T)
plt.axis([500,600,-1,1])
plt.grid(True)


    

In [29]:

    

time
y = sin(time)
plt.plot(time, sin(time))


    

    
Out[29]:





[<matplotlib.lines.Line2D at 0x121b26eb8>]






    

In [30]:

    

import scipy.fftpack as fftp


    

In [31]:

    

ffty = fftp.fft(y)


    

In [32]:

    

plt.plot(np.abs(ffty))


    

    
Out[32]:





[<matplotlib.lines.Line2D at 0x11ea4fa58>]






    

In [33]:

    

plt.plot(ms.fft_to_rfft(ffty))


    

    
Out[33]:





[<matplotlib.lines.Line2D at 0x121ba0dd8>]






    

In [34]:

    

plt.plot(rfft(y))


    

    
Out[34]:





[<matplotlib.lines.Line2D at 0x1224c8240>]






    

In [35]:

    

plt.plot(rfft(y+1)[:2],'*')


    

    
Out[35]:





[<matplotlib.lines.Line2D at 0x1225266d8>]






    

In [36]:

    

y.size


    

    
Out[36]:





200

In [37]:

    

plt.plot(rfft(np.cos(time))[:25],'*',rfft(y)[:25],'x')
plt.grid()


    

In [38]:

    

rfft = ms.fft_to_rfft(ffty)
rfft


    

    
Out[38]:





array([-6.23085483e-15,  3.65510101e-16,  7.69951025e-16, -6.36491753e-15,
        1.05094900e-15,  4.34409923e-16,  4.57159117e-15, -3.67719809e-15,
        3.78997887e-15, -2.10270340e-15,  8.11578815e-15, -7.34600178e-15,
        7.81693320e-15, -2.00710084e-15,  1.32811324e-14, -5.37005458e-15,
        1.90136542e-14,  6.64379137e-16,  3.63745285e-14, -1.40211486e-13,
       -1.00000000e+02, -1.03886241e-15, -3.96452581e-14, -1.03191619e-14,
       -2.13367569e-14,  1.03808448e-15, -1.43332521e-14, -3.07580320e-14,
       -1.23109705e-14, -9.82541347e-15, -6.94832314e-14,  2.72928007e-14,
       -8.39982935e-15, -4.66247286e-16, -8.39996108e-15,  6.41440650e-15,
       -6.92582565e-15,  5.74587678e-16, -5.39378006e-15,  8.45678455e-15,
       -4.57075312e-15,  7.01400090e-16, -6.11213130e-15,  1.25393487e-15,
       -5.16210277e-15, -2.55473457e-17, -4.49127876e-15,  2.08623362e-15,
       -4.36686410e-15,  1.44873246e-15, -7.29861897e-15, -1.90295910e-15,
       -4.58957064e-15,  8.53052472e-16, -1.23202687e-15,  5.20797970e-16,
       -3.34363360e-15,  4.93226508e-16, -6.60196758e-15, -7.70402240e-15,
       -2.26430247e-15, -2.01025521e-15, -1.74821920e-16, -5.04252669e-15,
       -1.57873299e-15, -4.17965806e-17, -4.46884755e-15, -2.84899673e-14,
       -1.84113198e-15, -4.02403466e-16, -5.22972573e-14,  3.19183251e-14,
       -2.21056151e-15, -4.26327749e-15, -6.69217246e-15, -2.28252594e-15,
       -1.37045226e-15, -1.35709917e-14,  1.40471501e-14,  1.31413422e-14,
       -7.07848918e-14,  1.34722612e-14,  1.03624630e-14,  1.98514301e-14,
       -1.13251513e-15,  5.08516169e-15, -4.93769208e-15,  6.03559429e-15,
       -1.00200641e-15,  2.60252866e-15, -1.85775943e-15,  7.88971930e-15,
       -1.18039189e-15,  2.16284690e-15, -1.70887445e-15,  5.27649253e-15,
       -1.02806372e-15,  1.01465866e-15,  5.84593330e-16,  3.09501858e-15,
       -4.44089210e-16,  2.89757875e-15, -1.34188531e-15,  4.67539226e-15,
       -1.68665448e-15, -4.88111384e-16, -2.64352643e-15,  1.15708228e-15,
       -2.90447029e-15, -4.37921095e-16,  1.35094283e-15,  5.23618038e-15,
       -2.94871148e-16, -1.48176598e-15, -5.31255689e-15, -5.89388173e-15,
       -1.06404707e-15, -1.26888916e-14,  1.35223415e-14,  1.19690618e-14,
       -7.07848918e-14,  1.56377990e-14,  1.49081933e-14,  2.72653138e-14,
       -5.15742671e-16,  5.83303210e-15, -4.99051077e-15,  3.17663670e-14,
       -1.43780955e-15,  2.42016849e-15,  5.78368667e-14, -1.81206033e-14,
       -1.04062703e-15,  3.61586612e-15, -3.76263196e-15, -6.52467693e-15,
       -2.81550268e-15,  1.67715785e-15,  1.61606510e-15,  6.94646322e-15,
        1.22889755e-15,  1.15715706e-15, -1.53044917e-15, -9.70688286e-16,
        6.68981067e-16,  2.76170539e-15,  2.04071272e-15,  3.37625374e-15,
        3.77506735e-16,  2.07676944e-15, -4.63408371e-15,  1.26085390e-15,
       -3.74655385e-16,  1.73628448e-15,  3.48747067e-16,  5.59636964e-15,
       -1.39958834e-15,  3.34725385e-15, -3.07379604e-15,  6.36617312e-15,
        8.87410592e-17,  2.50920051e-15,  8.10765009e-16,  1.07519133e-14,
       -5.27376251e-16,  2.63604741e-15, -2.60536557e-15,  2.80697736e-14,
       -9.24399445e-16,  1.18431785e-14,  4.06508927e-14, -2.67908848e-14,
        3.27742021e-16,  2.99556377e-15, -1.48826145e-15, -9.89850593e-15,
        4.77976461e-16,  5.29729435e-15, -3.81505258e-18, -5.20836643e-15,
       -7.10542736e-15,  2.09413330e-15, -2.65933699e-15, -4.01455659e-15,
       -1.45315897e-15,  1.35635766e-15,  4.59204351e-16, -1.73636868e-15,
        1.15733455e-15,  6.76581965e-16, -1.72007016e-15, -2.77423292e-15,
       -1.55876305e-15,  2.75837049e-15,  1.36280024e-15, -4.31618944e-16,
        4.94006840e-16,  1.89833881e-15, -3.02324444e-16, -1.36058882e-17])

In [39]:

    

plt.plot(rfft/np.max(np.abs(rfft)))


    

    
Out[39]:





[<matplotlib.lines.Line2D at 0x1226b7668>]






    

In [40]:

    

from scipy.fftpack import *


    

In [41]:

    

a = np.array([[9, -9, 1, -1],[9, -9, 1, -1]])
fft(a)


    

    
Out[41]:





array([[ 0.+0.j,  8.+8.j, 20.+0.j,  8.-8.j],
       [ 0.+0.j,  8.+8.j, 20.+0.j,  8.-8.j]])

In [42]:

    

rfft(a)


    

    
Out[42]:





array([[ 0.,  8.,  8., 20.],
       [ 0.,  8.,  8., 20.]])

In [44]:

    

fft(a)


    

    
Out[44]:





array([ 0. +0.j        , 13.5+7.79422863j, 13.5-7.79422863j])

In [45]:

    

rfft(a)


    

    
Out[45]:





array([ 0.        , 13.5       ,  7.79422863])

In [46]:

    

time = np.array([0,1,2,3])


    

In [ ]:

    

plt.plot(time_e, x.T*3.2/5, time, a.T, '*')


    

In [ ]:

    

x[0,0]


    

In [ ]:

    

a[0,0]


    

In [ ]:

    

rfft(a)


    

In [ ]:

    

rfft(x)


    

In [ ]:

    

plt.plot(irfft(x).T)


    

In [ ]:

    

x.shape


    

In [ ]:

    

time


    

In [ ]:

    

np.sin(time)


    

In [ ]:

    

tt = np.linspace(0,np.pi*2,101)


    

In [ ]:

    

np.zeros((2,2)).shape


    

In [ ]:

    

plt.plot(tt,sin(tt))


    

In [ ]:

    

rft = rfft(np.cos(tt))
rft


    

In [ ]:

    

rft[3:]= rft[3:]*0
rft


    

In [ ]:

    

plt.plot(irfft(rft))


    

In [ ]:

    

a


    

In [ ]:

    

time


    

In [ ]:

    

plt.plot(time, a.T, '*')


    

In [ ]:

    

b = irfft(np.hstack((rfft(a), a*0)))


    

In [ ]:

    

plt.plot(b.T)


    

In [ ]: