Consider the initial value problem given below where the behaviour of
the system is studied with respect to given set of inputs and initial
condition. Lets make it parameteric by adding a scalar parameter to the
ode and setting its value to zero.

ODE 1
-----

.. math::
   :nowrap:

   \begin{aligned} 
   (t_0,t_f,N)&=(0,10,25)\\ 
   x(t_0)&=\begin{bmatrix} 
   0\\
   0 
   \end{bmatrix} \\ 
   u(t)&=0.2t-1\\ 
   p&=0\\ 
   \begin{bmatrix}
   \dot{x}_0\\
   \dot{x}_1 
   \end{bmatrix}&= 
   \begin{bmatrix} 
   (1-x_1^2)x_0-x_1+u\\ 
   x_0
   \end{bmatrix}+p 
   \end{aligned}


The problem can be solved in 2 ways using Simulate.

1. eliminate controls from the ode
2. supply controls evaluated on the grid points

Elimination
+++++++++++

* Firstly, the ode needs to modelled for which we call the :code:`get_symbols()`
  function. The arguments are the number of states, controls, parameters and intervals
  for the grid. This also has an additional option to model the ode using symbolics 
  of the SX or MX class. MX is usually preferred due to its generality over SX.
  The :code:`get_symbols()` return a dictionary with required symbolics to construct
  the ode. This dictionary will keep updating as we pass it through each function.

* Using the symbolics stored in the dictionary, the ode is constructed. The
  ode will be passed to the :code:`scale_ode()` method where, a time domain transformation
  takes place. In the scaled ode, the initial and final time are simply parameters.

* Post domain transformation, the numerical solver for integration and its options must be set.
  Currently, CVODES, rk4, IRK methods are supported. For detailed options list, refer to integrator section of
  CasADi API. Calling the :code:`integrator()` creates an integrator object for the scaled ode.

* Lastly, :code:`simulate()` is called where the initial state, initial and final time, control and
  parameter if any is passed. The system is integrated and results are stored. This can be accessed using keys
  'x_res','u_res' and 't_grid' if required. 

.. code:: python

   import simulatetraj as a
   import casadi as cs
   sym_dict=a.get_symbols(n_x=2,n_u=0,n_p=1,N=25)
   t=sym_dict['t']
   x=sym_dict['x']
   u=sym_dict['u']
   p=sym_dict['p']
   f=cs.vertcat((1-x[1]**2)*x[0]-x[1]+0.2*t-1,x[0])+p
   _=a.scale_ode(sym_dict=sym_dict,f=f)
   _=a.integrator(sym_dict=sym_dict,options=None)
   a.simulate(t0=0,tf=10,x0=cs.DM([0,0]),p=cs.DM(0),sym_dict=sym_dict)
   a.plot_sol(sym_dict=sym_dict)

.. image:: ../../tests/Figure_1.svg
   :alt: Alternate text

.. image:: ../../tests/Figure_2.svg
   :alt: Alternate text

without control elimination
+++++++++++++++++++++++++++

The control input is approximated as a piecewise constant function.

.. code:: python

   import simulatetraj as a
   import casadi as cs
   sym_dict=a.get_symbols(n_x=2,n_u=1,n_p=1,N=25)
   t=sym_dict['t']
   x=sym_dict['x']
   u=sym_dict['u']
   p=sym_dict['p']
   f=cs.vertcat((1-x[1]**2)*x[0]-x[1]+u,x[0])+p
   _=a.scale_ode(sym_dict=sym_dict,f=f)
   _=a.integrator(sym_dict=sym_dict,options=None)
   a.simulate(t0=0,tf=10,x0=cs.DM([0,0]),u=cs.linspace(-1,1,25).T,p=cs.DM(0),sym_dict=sym_dict)
   a.plot_sol(sym_dict=sym_dict)   

.. image:: ../../tests/Figure_3.svg
   :alt: Alternate text

.. image:: ../../tests/Figure_4.svg
   :alt: Alternate text

.. image:: ../../tests/Figure_5.svg
   :alt: Alternate text


ODE 2
-----

In this example, the derivative of the square function is integrated.

.. math::
   :nowrap:   

   \begin{aligned} 
   (t_0,t_f,N)&=(0,10,100000)\\ 
   x(t_0)&=0\\
   \dot{x}&=2t
   \end{aligned}

.. code:: python

   sym_dict=get_symbols(n_x=1,n_u=0,n_p=0,N=100000)
   t=sym_dict['t']
   f=2*t
   _=scale_ode(sym_dict=sym_dict,f=f)
   _=integrator(sym_dict=sym_dict,options=None)
   x0=cs.DM([0])
   t0,tf=0,10
   simulate(t0=t0,tf=tf,x0=x0,sym_dict=sym_dict)
   plot_sol(sym_dict=sym_dict)
   print('Global error x(N+1) for xdot=2t:',sym_dict['x_res'][:,-1]-100)

.. code:: bash

   Global error x(N+1) for xdot=2t: DM(1.49132e-08)

.. image:: ../../tests/Figure_6.svg
   :alt: Alternate text

Lotka voltera/prey-predator model
---------------------------------

.. math::
   :nowrap:

   \begin{aligned} 
   (t_0,t_f,N)&=(0,15,1000)\\
   (\alpha,\beta)&=(0.01,0.02)\\
   x(t_0)&=
   \begin{bmatrix} 
   20\\ 
   20
   \end{bmatrix}\\ 
   \begin{bmatrix} 
   \dot{x}_0\\ 
   \dot{x}_1
   \end{bmatrix}&= 
   \begin{bmatrix} 
   x_0-\alpha x_0 x_1\\ 
   -x_1+\beta x_0 x_1 
   \end{bmatrix} 
   \end{aligned}

.. code:: python

   sym_dict=get_symbols(n_x=2,n_u=0,n_p=2,N=1000)
   x=sym_dict['x']
   p=sym_dict['p']
   f=cs.vertcat(x[0]-p[0]*x[0]*x[1],-x[1]+p[1]*x[0]*x[1])
   _=scale_ode(sym_dict=sym_dict,f=f)
   _=integrator(sym_dict=sym_dict,options=None)
   x0=cs.DM([20,20])
   t0,tf=0,15
   p=cs.DM([0.01,0.02])
   simulate(t0=t0,tf=tf,x0=x0,p=p,sym_dict=sym_dict)
   plot_sol(sym_dict=sym_dict)
   plt.plot(sym_dict['x_res'][0,:].full().flatten(),sym_dict['x_res'][1,:].full().flatten(),'o')
   plt.show()

.. image:: ../../tests/Figure_8.svg
   :alt: Alternate text