lines 8-195 of file: include/cppad/example/atomic_four/lin_ode/reverse.hpp

{xrst_begin atomic_four_lin_ode_reverse.hpp}
{xrst_spell
   lagrangian
   lll
   nnz
}

Atomic Linear ODE Reverse Mode: Example Implementation
######################################################

Purpose
*******
The ``reverse`` routine overrides the virtual functions
used by the atomic_four base; see
:ref:`reverse<atomic_four_reverse-name>` .

First Order Theory
******************
We are given a vector :math:`w \in \B{R}^m` and need to compute

.. math::

   \partial_x w^\R{T} z(r, x)

see the definition of :ref:`atomic_four_lin_ode@z(t, x)` .
Consider the Lagrangian corresponding to
:math:`w^\R{T} z(r, x)` as the objective and the ODE as the constraint:

.. math::

   L(x, \lambda)
   =
   w^\R{T} z(r, x) +
      \int_0^r \lambda(t, x)^\R{T}
         [ A(x) z(t, x) - z_t (t, x) ] \R{d} t

where :math:`\lambda : \R{R} \times \B{R}^n \rightarrow \B{R}^m`
is a smooth function.
If :math:`z(t, x)` satisfies its ODE, then

.. math::

   \partial_x w^\R{T} z(r, x)
   =
   \partial_x L(x, \lambda)

We use the following integration by parts to replace the :math:`z_t (t, x)`
term in the integral defining :math:`L(x, \lambda)`:

.. math::

   - \int_0^r \lambda(t, x)^\R{T} z_t (t, x) \R{d} t
   =
   - \left. \lambda(t, x)^\R{T} z(t, x) \right|_0^r
   +
   \int_0^r \lambda_t (t, x)^\R{T} z(t, x) \R{d} t

Adding the condition :math:`\lambda(r, x) = w`,
and noting that :math:`z(0, x) = b(x)`, we have

.. math::

   L(x, \lambda)
   =
   \lambda(0, x)^\R{T} z(0, x)
   +
   \int_0^r \lambda_t (t, x)^\R{T} z(t, x) \R{d} t
   +
   \int_0^r \lambda(t, x)^\R{T} A(x) z(t, x) \R{d} t

.. math::

   L(x, \lambda)
   =
   \lambda(0, x)^\R{T} b (x)
   +
   \int_0^r [ \lambda_t (t, x)^\R{T} + \lambda(t, x)^\R{T} A(x) ]
      z(t, x) \R{d} t

.. math::

   L(x, \lambda)
   =
   \lambda(0, x)^\R{T} b (x)
   +
   \int_0^r z(t, x)^\R{T}
      [ \lambda_t (t, x) + A(x)^\R{T} \lambda(t, x) ] \R{d} t

The partial derivative
of :math:`L(x, \lambda)` with respect to :math:`b_j`,
(not including the dependence of :math:`\lambda(t, x)` on :math:`x`)
is :

.. math::

   \partial_{b(j)} L(x, \lambda)
   =
   \lambda_j (0, x)

The partial derivative
of :math:`L(x, \lambda)` with respect to :math:`A_{i,j}`
(not including The dependence of :math:`\lambda(t, x)` on :math:`x`)
is :

.. math::

   \partial_{A(i,j)} L(x, \lambda)
   =
   \int_0^r \partial_{A(i,j)} z(t, x)^\R{T}
      [ \lambda_t (t, x) + A(x)^\R{T} \lambda(t, x) ] \R{d} t
   +
   \int_0^r z_j (t, x) \lambda_i (t, x) \R{d} t

If :math:`\lambda(t, x)` satisfies the ODE

.. math::

   0 = \lambda_t (t, x) + A(x)^\R{T} \lambda(t, x)

The partial derivative with respect to :math:`A_{i,j}` is

.. math::

   \partial_{A(i,j)} L(x, \lambda)
   =
   \int_0^r z_j (t, x) \lambda_i (t, x) \R{d} t

In summary, we can compute
an approximate solution for the initial value ODE:

.. math::

   z_t (t, x) = A(x) z(t, x) \W{,} z(0, x) = b(x)

and approximate solution for the final value ODE:

.. math::

   \lambda_t (t, x) = - A(x)^\R{T} \lambda(t, x)
   \W{,}
   \lambda(r, x) = w

Using the notation
:ref:`atomic_four_lin_ode@pattern@nnz` ,
:ref:`atomic_four_lin_ode@pattern@row` , and
:ref:`atomic_four_lin_ode@pattern@col` ,
We can compute an approximation for

.. math::

   \partial_{x(k)} w^\R{T} z(r, x)
   =
   \left\{ \begin{array}{lll}
   \int_0^r \lambda_i (t, x) z_j (r, x) \R{d} t
   & \R{where} \; i = \R{row} [k] \W{,} j = \R{col}[k]
   & \R{if} \; k < nnz
   \\
   \lambda_i (0, x)
   & \R{where} \; i = k - nnz
   & \R{otherwise}
   %
   \end{array} \right.

{xrst_toc_hidden
   include/cppad/example/atomic_four/lin_ode/reverse_2.xrst
}
Second Order Theory
*******************
:ref:`atomic_four_lin_ode_reverse_2-name` .

Simpson's Rule
**************
This example uses Simpson's rule to approximate the integral

.. math::

   \int_0^r \lambda_i (t, x) z_j (t, x) \R{d} t

Any other approximation for this integral can be used.

Source
******
{xrst_literal
   // BEGIN C++
   // END C++
}

{xrst_end atomic_four_lin_ode_reverse.hpp}