Normal forms
============

Background
~~~~~~~~~~

A **normal form** of :math:`\dot{x}=f(x)` (:math:`\dot{x}` means :math:`dx/dt`) centered at an equilibrium :math:`x_0` is an equation :math:`\dot{x}=h(x)` such that :math:`h(x)` is homeomorphic to :math:`f(x)` in a neighborhood of :math:`x_0` and :math:`h'(x)f'(x_0)^Tx-f'(x_0)^Th(x)=0`. The normal form :math:`\dot{x}=h(x)` has the same qualitative dynamics as the starting equation :math:`\dot{x}=f(x)` near :math:`x_0` since the vector fields :math:`f,h` are topologically equivalent in a neighborhood of :math:`x_0`. Normal forms enjoy neat geometric properties such as: the flow :math:`\varphi_h^t(x)` generated by :math:`\dot{x}=h(x)` preserves the foliation :math:`F` induced by the flow :math:`e^{f'(x)t}` near :math:`x_0`; that is, if :math:`\ell` is a leaf of :math:`F`, then so is :math:`\varphi_h^t(\ell)` for appropriate :math:`\ell` and :math:`t`.

Normal forms can be computed by a sequence of substitutions in the context of asymptotic series. A substitution :math:`x=\varphi(y)` transforms :math:`\dot{x}=f(x)` to :math:`\dot{y}=(\varphi')^{-1}(y)f(\varphi(y))`. Differentiate the substitution with respect to :math:`t` and solve for :math:`\dot{y}` to check the last statement. This transformation can be computed using

The fundamental theorem of Lie series for vector fields
    If :math:`\varphi^t_g(x)` is the flow generated by the vector field :math:`g(x)`, then the *similarity transformation* :math:`S_{g}f=(\varphi_g')^{-1}f\circ\varphi_g` is formally expanded as

    .. math::

       S_{g}f\sim e^{L_g}f = \left(I + L_g + \frac1{2!}L^2_g + \cdots \right)f

    where :math:`L_g(f(x))=[f(x),g(x)]=f'(x)g(x)-g'(x)f(x)`.

Notice that if :math:`f_i,g_j` are degree :math:`i,j` homogenous polynomial vector fields, respectively, then :math:`L_{j}f_i(x)`, where :math:`L_j` denotes :math:`L_{g_j}`, is degree :math:`i+j-1`. So, the subsitution :math:`x=\varphi_{g_j}(y)` leaves the Taylor series of :math:`f` unaltered up to degree :math:`j-1` by

.. math::

   e^{L_j}f(x) &= \left(I + L_j + \cdots \right)\left(f_1(x) + \cdots + f_{j-1}(x) + f_j(x) + O(x^{j+1}) \right) \\
   &= \underbrace{f_1(x) + \cdots + f_{j-1}(x)}_{\text{deg }<j} + \underbrace{L_jf_1(x) + f_{j}(x)}_{h_j(x)} + O(x^{j+1})

where :math:`f_i(x)` denotes the degree :math:`i` term in the Taylor series of :math:`f` centered at :math:`x_0`.

A sequence of polynomial vector fields :math:`g_j(x)` (called *generators*) for :math:`j>1` can be chosen to ensure that the equation :math:`\dot{y}=f_1(y)+h_2(y)+h_3(y)+\cdots` resulting from the substitution :math:`x=\cdots\circ\varphi_3\circ\varphi_2(y)` is a normal form of :math:`\dot{x}=f(x)` by the following.

Rearranging the formula for the modified degree :math:`j` term :math:`h_j(x)` labeled above and using :math:`[f,g]=-[g,f]`, we have

.. math::

  L_1g_j(x) = f_j(x) - h_j(x).

where :math:`L_1` denotes :math:`L_{f_1}`. If :math:`h_j(x)` is the part of :math:`f_j(x)` in the nullspace of :math:`L_1^T`, so that :math:`L_1^T(h_j(x))=h'_j(x)f'(x_0)^Tx-f'(x_0)^Th_j(x)=0`, then the last equation can be solved for :math:`g_j(x)`.

The main class of this package, :class:`normal_forms.normal_form.normal_form`, implements the following algorithm for computing a normal form :math:`\dot{x}=f_1(x)+h_2(x)+\cdots+h_k(x)+O(x^{k+1})` of :math:`\dot{x}=f(x)` up to a specified degree :math:`k`.

0. Select an orthonormal basis :math:`N` for the nullspace of :math:`L_1^T`

For :math:`j=2,\ldots,k`

1. Set :math:`h_j(x)=N^Tf_j(x)` to the projection of :math:`f_j(x)` onto :math:`N`
2. Solve :math:`L_1g_j(x)=f_j(x)-h_j(x)` for :math:`g_j(x)`
3. Compute the transformation :math:`e^{L_j}f(x)` up to degree :math:`k` and relabel it as :math:`f(x)`

.. note::

   In some theoretically-oriented presentations of normal forms, a preliminary linear substitution :math:`x=Ty` is performed to bring :math:`f_1(y)` into a canonical form such as Jordan form, and a basis :math:`N` for the nullspace of :math:`L_1^T` is chosen to simplify an ensuing analysis. This package deals directly with :math:`f_1(x)=f'(x_0)x` derived from the supplied righthand side :math:`f` since computing Jordan forms is ill-conditioned. And, this package chooses :math:`N` agnostically with respect to any ensuing analysis corresponding to the singular value decomposition modes of :math:`L_1^T` with zero singular values.

The terms :math:`h_j(x)` that appear in the normal form of :math:`\dot{x}=f(x)` are fully determined by :math:`f_1(x)=f'(x_0)x` and :math:`N`. If :math:`f_1'(x_0)` is diagonalizable over :math:`\mathbb{C}` with eigenvalues :math:`\lambda_i`, then so is :math:`L_1:\mathcal{V}_j^n\rightarrow\mathcal{V}_j^n` but with eigenvalues :math:`m_1\lambda_1+\cdots+m_n\lambda_n-\lambda_i` for all nonnegative :math:`m_1,\ldots,m_n` such that :math:`|m_1|+\cdots+|m_n|=j`, where :math:`\mathcal{V}_j^n` is the space of :math:`n`-dim homogenous degree :math:`j` polynomial vector fields. If :math:`L_1` has all nonzero eigenvalues, that is has no integer relations :math:`m_1\lambda_1+\cdots+m_n\lambda_n=\lambda_i` with :math:`m_1,\ldots,m_n` as before, for :math:`j>1`, then :math:`h_j(x)=0` since :math:`L_1^T:\mathcal{V}_j^n\rightarrow\mathcal{V}_j^n` has a trivial nullspace. Otherwise, :math:`h_j(x)` may appear as a nontrivial vector field in the nullspace of :math:`L_1^T`.


A 1-D example: saddle-node bifurcation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The equation

.. math::

   \dot{x} = r+1-x-e^{-x}

after expanding the righthand side is

.. math::

   \dot{x} &= r+1-x-\left(1-x+\frac12x^2+O(x^3)\right) \\
   &= r-\frac12x^2+O(x^3).

Since :math:`f_1=0`, every polynomial is in the nullspace of :math:`L_1^T=0` and therefore the equation is already in normal form. More pedantically, a basis for the nullspace of :math:`L_1^T:\mathcal{V}_j^n\rightarrow\mathcal{V}_j^n` is :math:`\{x^j\}` and clearly every degree :math:`j` term for :math:`j>1` in the Taylor series of :math:`f` is a multiple of one of these basis monomials.

The above equation undergoes a saddle-node bifurcation at :math:`r=0`: a saddle and a node that exist for :math:`r>0` coalesce at :math:`x=0` when :math:`r=0` and disappear as :math:`r` is decreased [1]_. The above equation is then topologically equivalent in a neighborhood of :math:`x=0` to any equation :math:`\dot{x}=f(x)` with :math:`f:\mathbb{R}\rightarrow\mathbb{R}` that has a saddle node equilibrium at :math:`x=0`.

To compute the normal form with this package, supply the equation righthand side :math:`f`, the center of the expansion :math:`x_0`, and the maximum degree :math:`k` of the expansion. The ``sympy`` variants of non-algebraic functions such as ``sin``, ``cos``, and ``exp`` should be used when defining ``f``.


.. ipython::

   In [0]: import sympy

   In [0]: f = lambda x, r=0: r + 1 - x - sympy.exp(-x)

   In [0]: from normal_forms import normal_form

   In [0]: h = normal_form(f, x=0, k=2)

You can see a symbolic representation of the normal form righthand side using the normal form object's ``fun`` attribute. You can also call the normal form object to evaluate :math:`f_1(x)+h_2(x)+\cdots+h_k(x)`.

.. ipython::

   In [0]: h.fun

   In [0]: h(2)

A 2-D example: Duffing's equation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. plot:: normal_forms/examples/normal_form/05.py

.. [1] Strogatz, Nonlinear Dynamics and Chaos, Example 3.1.2