.. _ch05-python-sympy_scypy:


===============================================================
SymPy
===============================================================

So far we have mostly discussed scientific computing from a numerical
perspective. Symbolic computing can also done readily in Python using
`SymPy <http://www.sympy.org/>`_. The SymPy website has many great 
tutorials as well as the documentation for SymPy.

-------------------------------------------------------
An example: Finding wavespeeds of a hyperbolic system
-------------------------------------------------------

In physics, wavelike phenomena often arise from systems of hyperbolic
conservation laws. Generally, these take the form

.. math::

   \frac{\partial\textbf{U}}{\partial t} + \nabla\cdot\textbf{F}(\textbf{U}) = 0

where :math:`\textbf{U}` is a vector of conserved quantities and
:math:`\textbf{F}` is a flux tensor describing the flow of these quantities.
In two space dimensions the above becomes

.. math::

   \frac{\partial\textbf{U}}{\partial t} + \frac{\partial\textbf{F}}{\partial x} + \frac{\partial\textbf{G}}{\partial y} = 0

where :math:`\textbf{F} = \left(\textbf{F}(\textbf{U}),\textbf{G}(\textbf{U})\right)`
consists of the :math:`x` and :math:`y` direction fluxes. We can use the
chain rule to decompose the flux derivatives, which gives the so called
quasi-linear form of the equation:

.. math::

   \frac{\partial\textbf{U}}{\partial t} + \frac{\partial\textbf{F}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial x} + \frac{\partial\textbf{G}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial y} = 0

where :math:`\frac{\partial\textbf{F}}{\partial\textbf{U}}` and
:math:`\frac{\partial\textbf{G}}{\partial\textbf{U}}` are Jacobian
matrices of each flux. These matrices are interesting both mathematically
and physically. In particular, the eigenvalues of these matrices yield
the wavespeeds in each direction.

Let's consider the shallow water equations (in conservation form)

.. math::

   \frac{\partial h}{\partial t} &+& \frac{\partial}{\partial x}&(hu) &+& \frac{\partial}{\partial y}(hv) = 0 \\
   \frac{\partial (hu)}{\partial t} &+& \frac{\partial}{\partial x}&\left(hu^2 + \frac{1}{2}gh^2\right) &+& \frac{\partial}{\partial y}(huv) = 0 \\
   \frac{\partial (hv)}{\partial t} &+& \frac{\partial}{\partial x}&(huv) &+& \frac{\partial}{\partial y}\left(hv^2 + \frac{1}{2}gh^2\right) = 0

where the conserved variables are the height of the fluid :math:`h`, and
the momenta :math:`hu` and :math:`hv`. Additionally, :math:`g` is the
gravitational acceleration. The following SymPy code generates the flux
Jacobian in the :math:`x` direction, and finds the corresponding
wavespeeds.

.. literalinclude:: ./codes/sweEigs.py

:download:`Download this code <./codes/sweEigs.py>`

**Exercises:**

* Generate :math:`y` direction wavespeeds as well

* Generate eigenvectors (waves) of the flux Jacobian

===============================================================
SciPy
===============================================================

SciPy is a Python package that contains numerous scientific toolboxes
that are commonly used in the scientific computing community. SciPy is
built on top of NumPy, and two work together quite well.

Included in SciPy are implementations of numerical interpolation and 
quadrature rules, optimization, image processing, statistics, special
functions, etc. The routines available in SciPy are highly optimized,
robust, and well documented. To see examples of how to use SciPy for
scientific computing check out the `SciPy users guide <https://docs.scipy.org/doc/scipy/tutorial/index.html#user-guide>`_

.. _ch05-scipy-lorenz:
   
-----------------------------------------------
An example: Solving a nonlinear, chaotic, ODE
-----------------------------------------------

To demonstrate one small piece of SciPy we'll solve the
`Lorenz system <https://en.wikipedia.org/wiki/Lorenz_system>`_
of equations:

.. math::

   \frac{dx}{dt} &= \sigma (y - x) \\
   \frac{dy}{dt} &= x (\rho - z) - y \\
   \frac{dz}{dt} &= xy - \beta z

where :math:`\rho,\sigma,\beta` are parameters which we'll fix at
:math:`28`, :math:`10`, and :math:`8/3` respectively. These equations
arise as a highly simplified, reduced order model, for convection as
applicable to weather prediction. While quite simple to write down,
these equations are nonlinear and can not be solved analytically.
Furthermore, they give rise to complex and chaotic behavior, and are
arguably the most well known example of a (continuous) chaotic system.

As such, these are a perfect candidate for numerical treatment. SciPy
makes solving equations like this nearly trivial, as the following code
demonstrates.

.. literalinclude:: ./codes/lorenz_demo.py
   :linenos:

:download:`Download this code <./codes/lorenz_demo.py>`

This should have produced the following two figures.
	  
.. _figLorenzTime:

.. figure:: ./_figures/lorenzTime.png
   :align: center

   Time history of each component in the Lorenz system
	  
.. _figLorenzPhase:

.. figure:: ./_figures/lorenzPhase.png
   :align: center

   Phase space trajectory of the Lorenz system from the given
   initial condition.