\documentclass{beamer}
%\documentclass[handout]{beamer} 


\usetheme{Boadilla}

\usepackage{hyperref}
\usepackage{ulem}
\usepackage{fontawesome}
\usepackage{comment}
%\usepackage{listings}

\newcommand{\scalefont}[1]{%
  \edef\scale@fontsize{\fpeval{#1*\f@size}}%
  \edef\scale@fontbaselineskip{\fpeval{1.2*\scale@fontsize}}%
  \fontsize{\scale@fontsize}{\scale@fontbaselineskip}\selectfont
}

\begin{document}

\title[The power of Series]{Stability via computed Power Series}
\subtitle{JMM 2022 in Seattle}
\author[thelwerj@jmu.edu]{Roger Thelwell}
\institute[JMU]{James Madison University}
\date{\today}
%\date{January 7, 2022}

\frame{\titlepage}

\title[IVODEs: Error, Convergence, and Control]{The power of Series}
\subtitle{JMM 2022 in \sout{Seattle} online}
\author[thelwerj@jmu.edu]{Steve Lucas, Joseph Rudmin Jr.,  Jim Sochacki, Roger Thelwell*,  Anthony Tongen, Paul Warne }
\institute[JMU]{James Madison University}
\date{\today}
%\date{October 10, 2015}

\frame{\titlepage}
%%%%%%
\begin{frame}
\frametitle{Abstract}

\only<1>{\color{gray} We exploit the power of an iterative symbolic computation to efficiently compute the Lyapunov spectrum of highly nonlinear systems of ODEs. The technique constructs a truncated Taylor series representation (a jet of degree $k$) at  a generic point of the separation of the flow. An approximate spectrum is then recovered and analyzed. Several examples will be discussed. (circa Sept 2, 2021)}
\only<2>{We exploit the power of an iterative symbolic computation to explore ODEs. We will talk about solutions, stability, error bounds, convergence, and control. The methods and theory are accessible to undergraduate students at the end of a first course in ODE, and the techniques are useful extensions that can be applied to power series solutions, and applicable to a wide range complex questions. (\today)}
 
\end{frame}

% SLIDE

%%%%%%

\begin{frame}
\frametitle{Outline}
    \tableofcontents
    \setcounter{tocdepth}{0}

\end{frame}

% SLIDE




%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The (toy) problem}
\begin{frame}
\frametitle{$y' = \alpha y^2$}

We'll begin with a toy problem...
\[ 
y' = \alpha y^2 \qquad y(0) = y_0
\]
\pause
An analytic solution to this one is easy:
\[
y(t) = - \frac{y_0}{\alpha  y_0 t - 1}
\]

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A (power series) solution}
\begin{frame}
\frametitle{$y' = \alpha y^2$ : series solution?}

What happens if we try (formal, for now) series?\\
\pause
Let
\[
y(t) = \sum_{k=0}^\infty y_k t^k
\]
then
\[
\sum_{k=0}^\infty (k+1) y_{k+1} t^k = \alpha \sum_{k=0}^\infty \left( \sum_{i,j\ge0}^{i+j = k} y_i y_j\right) t^k
\]
so we equate coefficients to get
\[
y_{k+1} = \frac{\alpha}{k+1} \sum_{i+j = k} y_i y_j
\]
\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{$y' = \alpha y^2$: series solution?}

Or, with \textsc{MAPLE}:

\begin{verbatim}
> ODE1 := diff(y(t),t) = alpha*y(t)^2;
> IC := y(0) = y0;
> y1 := dsolve({ODE1,IC},y(t))
\end{verbatim}
and
\begin{verbatim}
> Y1 := dsolve({ODE1,IC},y(t),series);
       Y1 := 1+2*alpha*t*y0+3*alpha^2*y0^2*t^2+...
\end{verbatim}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The power of Power Series}
\subsection{Sensitivity}
\begin{comment}
\begin{frame}
\frametitle{$y' = \alpha y^2$ : sensitivity?}

Given either solution, we see
\[
\partial_{\alpha} y1 = \frac{y_0^2 t}{\left( \alpha y_0 t - 1 \right)^2}
\]
or 
\[
\partial_{y_0} Y1 =1+2 \alpha t y_0+3 \alpha^2 y_0^2 t^2+ ...
\]

Computing sensitivity to perturbation in initial condition and/or parameter(s) seems easy!\\

%\pause
%Applications: stability and inverse problems. 
\end{frame}
\end{comment}


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{$y' = \alpha y^2$: sensitivity?}
Consider the ODE system
\[
\mathbf{\dot x} = \mathbf{F(x)} \qquad \mbox{ with } \mathbf{x}(0) = \mathbf{x}_0.
\]
and let
\[
\phi(t;\mathbf{x}_0)
\]
represent the flow of this system through the initial point $\mathbf{x}_0$.\\
\pause
Suppose we want to understand the divergence of two nearby trajectories:
\[
\phi(t;\mathbf{y}_0) - \phi(t;\mathbf{x}_0) \approx D_\mathbf{x} \phi(t;\mathbf{x}_0) (\mathbf{y}_0 - \mathbf{x}_0), 
\]
but $D_\mathbf{x} \phi(t;\mathbf{x}_0)$ is usually hard to compute! 

\pause

\url{https://en.wikipedia.org/wiki/Lyapunov_exponent}

\end{frame}



%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{Aside: The Lyapunov exponent}

For any curve of initial conditions $\mathbf{x}_s$, define
\[
\mathbf{v}(t) = \partial_{\mathbf s} \phi(t;\mathbf{x}_s) \big|_{s=0},
\]
then $\mathbf{v}(t)$ satisfies the first variation equation
\[
\dot{\mathbf{v}}= D_\mathbf{x} \mathbf{F}(\phi(t;\mathbf{x}) \big|_{\mathbf{x}_0}  \mathbf{v}_0 \quad \mbox{ with } \quad \mathbf{v}_0 = \partial_{\mathbf s} \mathbf{x}_s. 
\]
\pause
The Lyapunov exponent captures the (exponential) growth rate of \alert{ $\mathbf{v}(t)$ }. 
 
 \pause
The largest \textit{Lyapunov} exponent, $\lambda$, can be computed as: 
%\[
%|| \delta \mathbf v(t) || \approx exp(\lambda(t) \cdot t)  \mathbf ||\delta \mathbf{v_0}||,
%\]
\[
|| \mathbf v(t) || \approx exp(\lambda \cdot t)  \mathbf || \mathbf{v_0}||,
\]
%where $\delta \mathbf{v}$ is the (first) variation of trajectory $\mathbf{v}_0$.\\
\pause More generally, the full Lyapunov spectrum might be computed, but this is more difficult.

\pause
Can \alert{power series} help?

\end{frame}


%%%%%%%%%%%%%%%%%
%%%%%%%
\begin{frame}[fragile]
\frametitle{ $\displaystyle y' = \alpha(t) y$}

With \textsc{Maple} \\
\begin{verbatim}
>> restart:
>> Order := 4:
>> alpha := t -> sum(a[k]*t^k,k=0..Order):
>> GROWTH := diff(y(t),t) = alpha(t)*y(t):
>> Yseries := dsolve({GROWTH,y(0)=y[0]},y(t),type='series');
\end{verbatim}
\begin{align*}
y \left( t \right) = & y_0 + a_0 \, y_0t+ \left( 1/2\,{a_0}^2 \, y_0+1/2\,a_1 \,y_0 \right) t^2 + \\
 &  \left( 1/6\,{a_0}^3 \, y_0+1/2\,a_1 \, a_0 \, y_0 +1/3\,a_2 \, y_0 \right) t^3 + O \left( t^4 \right)
\end{align*}
\pause
which we can check
\begin{verbatim}
>> SOLN1 := y[0] * exp(int(alpha(tau),tau=0..t));
>> check := taylor(SOLN1,t=0) - Yseries;

\end{verbatim}

\end{frame}

%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{Why?}
From
\begin{align*}
y \left( t \right) = & y_0 + a_0 \, y_0t+ \left( 1/2\,{a_0}^2 \, y_0+1/2\,a_1 \,y_0 \right) t^2 + \\
 &  \left( 1/6\,{a_0}^3 \, y_0+1/2\,a_1 \, a_0 \, y_0 +1/3\,a_2 \, y_0 \right) t^3 + O \left( t^4 \right)
\end{align*}
we can find
\begin{align*}
\partial_{y_0} y(t) = & 1+a_0 t+ \left( 1/2\,{a_0}^{2}+1/2\,a_1 \right) {t}^{2}+\\
 & \left( 1/6\,{a_0}^{3}+1/2\,a_1 \, a_0+1/3\,a_2 \right) t^3 +O \left( t^4 \right)
\end{align*}
Sensitivity to initial conditions! 
\pause
Which we can easily compute...\\
\begin{verbatim}
>> Yy0 := taylor(diff(Yseries,y_0),t=0);
\end{verbatim}
%\begin{align*} & = 1+a_{{0}}t+ \left( 1/2\,{a_{{0}}}^{2}+1/2\,a_{{1}} \right) {t}^{2}+ \\
%& \qquad \left( 1/6\,{a_{{0}}}^{3}+1/2\,a_{{1}}a_{{0}}+1/3\,a_{{2}} \right) t^{3}+O \left( {t}^{4} \right) 
%\end{align*}
\end{frame}

%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{The Lyapunov exponent}
We have
\[
\partial_{y_0} y \left( t \right) =  \sum_{n=0}^\infty f_n(y_0,\ldots,y_{n-1}) t^n 
\]
and so $\lambda(t)$ is easy to compute: \\
\texttt{>> simplify(taylor(ln(Yy0))/t);}
\[
\lambda(t) = a_0 + 1/2 a_1 t + 1/3 a_2 t^2 + O (t^3)
\]
For our problem,  a direct calculation verifies this:
\[
\frac{1}{t} \int_0^t \alpha(\tau) d\tau).
\]
This time average is the mean coefficient on $[0,t]$\\
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Error}
\begin{frame}[fragile]
\frametitle{$y' = 2 * y^2, \quad y(0) = 1$: error?}

What about the error?
\begin{verbatim}
> alpha := 2; y0 := 1; plot(abs(Y1-y1),t=0..0.5);
\end{verbatim}
\includegraphics[scale=0.3]{JMM1_error.png}

Can we quantify this?

\pause
Let's take a little diversion...

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\subsubsection{An Aside}
\begin{frame}
\frametitle{$y' = \alpha y^m$}

The solution to the constant coefficient nonlinear IVODE
\[ y' = \alpha y^m \qquad y(0) = y_0
\]
is messy:
\[
y(t) = \left( (\alpha -\alpha m) t + y_0^{1-m} \right) ^{-\left( m-1 \right) ^{-1}}.
\]
But the ratio of $\frac{y}{y'}$ isn't!
\[
\frac{y}{y'} =  \frac{(\alpha -\alpha m ) t+y_0^{1-m}}{\alpha}
\]
and so
\[
y'(t) = \underbrace{\frac{\alpha}{(\alpha -\alpha m)t +  y_0^{1-m} }}_{K(t)} y(t),
\]
a non-constant coefficient \alert{LINEAR} ode.

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{$y' = \alpha y^m$: important aside}

So
\[
y'(t) = \underbrace{\frac{\alpha}{(\alpha -\alpha t m) t + y_0^{1-m} }}_{K(t)} y(t),
\]
has solution
\[
y(t) = y_0 exp\left( \int_0^t K(\tau) d \tau \right) 
\]
or, via series, 
\[ 
Y_{k+1} = \frac{\alpha ( 1 + (m-1)k)}{y_0^{1-m}(k+1)} Y_k
\]

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{$y' = \alpha y^m$: important aside}
From 
\[ 
Y_{k+1} = \frac{\alpha ( 1 + (m-1)k)}{y_0^{1-m}(k+1)} Y_k
\]
and for $m \ge 2$,
\[
Y_{k+1} \le  (m-1)|y_0|^{m-1} Y_k := C_\infty Y_k.
\]
This leads directly to a geometric series bounding $y(t)$:
\[
y(t) \le \frac{|y_0|}{1-C_\infty} = |y_0| \sum_{k=0} (C_\infty t)^k 
\]
Now for the bound...
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{$y' = \alpha y^2$: back to error}

From 
\[
y(t) \le \frac{|y_0|}{1-C_\infty} = |y_0| \sum_{k=0} (C_\infty t)^k 
\]
we see that the absolute error is
\[
|Y1 - y1| \le |y_0| \sum_{k=n+1}^\infty C_\infty |t|^k \le \frac{|y_0| C_\infty^{n+1}}{1-C_\infty |t|}
\]
where $C_\infty = |y_0 \alpha|$.

An \alert{ERROR} bound!
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{$y' = 2 * y^2, \quad y(0) = 1 $: error plots}

\begin{verbatim}
> ee :=  abs(Y1-y1);
> m := 2; Cinf := y0*alpha;
> EE := N -> abs(y0)*(Cinf*t)^(N+1)/( 1 - Cinf*abs(t));
> plot({ee,EE(5)},t=0..0.3);
> plot({ee-EE(5)},t=0..0.2);
\end{verbatim}

\includegraphics[scale=0.3]{JMM3_error.png}
\includegraphics[scale=0.3]{JMM4_error.png}
 
 Great!

\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Convergence}
\begin{frame}[fragile]
\frametitle{$y' = 2 y^2, \quad y(0) = 1$: radius of convergence?}

\begin{verbatim}
>plot({ee-EE(5)},t=0..0.48);
\end{verbatim}
\includegraphics[scale=0.3]{JMM5_error.png}

Hmmmm.... 
\pause
\faLightbulbO
\pause 
\[
y(t) = - \frac{y_0}{\alpha y_0 t - 1}
\]

\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{$y' = \alpha y^2$: radius of convergence}

But what if we only have this form?
 \begin{verbatim}
> Y20 := dsolve({ODE1,IC},y(t),series);
       Y20 := 1+2*alpha*t*y0+3*alpha^2*y0^2*t^2+...
\end{verbatim}
\pause
Approximate using the root or ratio test. Or calculate the Pad\'{e} approximant.

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{$y' = 2 y^2, y(0) = 1$: radius of convergence?}

 \begin{verbatim}
> Order := 20: alpha := 2; y0 := 1;
> Y20 := rhs(dsolve({ODE1,IC},y(t),series)):
> Ycoeff := [seq(coeff(convert(Y20,polynom),t,i),i=1..20)];
> RatioT := i -> abs(a[i+1]/a[i]): RootT := i -> abs(a[i])^(1/i):
> RootEst := [seq(1/RootT(i),i=1..Order)]: RatioEst := [seq(1/RatioT(i),i=1..Order)]: 
> plot ....      
\end{verbatim}

\includegraphics[scale=0.3]{JMM6_radius_of_cvg.png}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Control}
\begin{frame}
\frametitle{$y' = \alpha y^2$: for control?}

Suppose we want to control 
\[ y' = \alpha y^2 \qquad y(0) = y(0) \]
so that $y(T) = \beta$.\\
If we apply frictional damping to the system 
\[ y' - \alpha y^2 = u, \]
where $u = k t y'$,
can we drive the system to the desired state?
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{$y' = \alpha y^2$: for control?}

Let's try to drive the system so that $Y(0.2) = 3.2.$\\


\begin{verbatim}
> ODE := diff(y(t),t) = alpha*y(t)^2 + k*t*diff(y(t),t);
> IC := y(0) = 2;
> Y := convert(rhs(dsolve({ODE1,IC},y(t),series)),polynom);
> kvals := solve(subs(t=0.2,Y)=3.2);
       -2.051118481+8.300750459*I, -.6504685102 , ...
> subs({k=kvals[2],t=0.2},Y);
        3.200000000
\end{verbatim}
     
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{$y'=\alpha y^2 + kt y'$}

It looks like we can.
\begin{verbatim}
> kvals := solve(subs(t=0.2,Y)=3.2);
       -2.051118481+8.300750459*I, -.6504685102 , ...
>plot({Yk(0),Yk(kvals[2]),Yk(kvals[3])},t=0..0.2);
\end{verbatim}
\includegraphics[scale=0.35]{JMM2_control.png}

Repeated application allows trajectory control, and our error bound still applies to the forced system!

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add some numerical analysis

\begin{comment}
\begin{frame}
\frametitle{$y' = \alpha y^2$ : numerical analysis}

Can we measure accumulated numerical error? Just add it to the symbolic recursion!

\end{frame}
\end{comment}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Some (non-polynomial) problems}
\begin{frame}
\frametitle{Torsional Deformation of Blatz-Ko material?}

Torsional deformation of a compressible elastic solid cylinder can be modelled as: 
\begin{equation*}
\begin{small}
y''\left[ {1 + \frac{{\sqrt 2 }} 
{{{{y'}^2}}}{{\left( {\frac{x} 
{{yy'}}} \right)}^{\sqrt 2  - 1}}} \right] = \left[ {\sqrt 2 {{\left( {\frac{x} 
{{yy'}}} \right)}^{\sqrt 2 }} + 1} \right]\left[ {\frac{y} 
{{{x^2}}} - \frac{{y'}} 
{x}} \right] + \frac{{{\pi ^2}}} 
{{16}}y
\end{small} 
\end{equation*}
A solution to arbitrary order can also be computed via power series.
\pause
Just introduce auxiliary variables to generate a polynomial system.
All the same ideas apply!

\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Six degree-of-freedom flight mechanics}

\includegraphics[scale=0.35]{JMM7_6dof_schematic.png}
from: \url{http://plantarchy.us/katko/projects/dope/DSTO-TR-0931-PR.pdf}

A solution to arbitrary order can be computed via power series.
\pause
Just introduce auxiliary variables to generate a polynomial system.
All the same ideas apply!
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
\begin{frame}
\frametitle{Conclusions}

We considered a toy problem already cast as a polynomial ode. Extension and application of these methods will rely on the use of auxiliary variables to build a system of polynomial IVODEs. Once the system is polynomial, series methods allow remarkably direct analysis.
\begin{itemize}
\item Analytic approximation of solution
\item Stability and sensitivity
\item Easy error \alert{BOUND}
\item Radius of Convergence estimate
\item Simple control?
\end{itemize}
These techniques apply to a broad range of highly nonlinear ODE.

\end{frame} 

%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{\faQuestion}

\begin{center}
\Large{ Thanks! }\\[10mm] 

Questions? \\

\large{ thelwerj@jmu.edu }\\
\end{center}

\url{http://educ.jmu.edu/~sochacjs/PSM.html}

 \end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{document}

\begin{frame}
\frametitle{For more information}, 

Thelwell et al.
\newblock Cauchy Kowalevski and Polynomial ODE 
\newblock {\em EJDE}, 11, 1--8, 2012.


\end{frame}

\end{document}

\begin{frame}
\frametitle{References}
\begin{thebibliography}{10}

\bibitem{PaperXXFehlberg1964}
E.~Fehlberg
\newblock Numerical integration of differential equations by power series expansions, illustrated by physical examples.
\newblock Technical Report NASA-TN-D-2356, NASA, 1964.


%\bibitem{PaperXXNeidinger2010}
%R.D. Neidinger
%\newblock Introduction to automatic differentiation and Matlab object-oriented programming.
%\newblock {\em SIAM Review}, 52(3):545--563, 2010.

\bibitem{PaperXXCarothers2005}
David Carothers et. al.
\newblock Some properties of solutions to polynomial systems of differential equations.
\newblock {\em Electronic Journal of Differential Equations}, 2005:1--18, 2005.

\bibitem{PaperXXWarne2006}
P.~G. Warne et. al.
\newblock Explicit a-priori error bounds and adaptive error control for approximation of nonlinear initial value differential systems.
\newblock {\em Comput. Math. Appl.}, 52(12):1695--1710, 2006.

\end{thebibliography}

\end{frame}

\end{document}

plot({subs(k=0,Y),subs(k=kvals[2],Y),subs(k=kvals[3],Y)},t=0..0.2);
plot({subs(k=0,Y),subs(k=kvals[2],Y),subs(k=kvals[3],Y)},t=0..0.2);






\end{document}





\begin{frame}
\frametitle{The Butterfly}

\url{file://Users/faculty/Downloads/Lorenz.ogv.480p.webm}
\vfill
\verb+https://commons.wikimedia.org/wiki/File:Lorenz.ogv+

\pause
\vfill

How do we measure the sensitivity to perturbation in initial condition?

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Preliminaries}
%%%%%%%%%%
\subsection{Some notation}

\begin{frame}
\frametitle{A general ODE system}

Consider the ODE system
\[
\mathbf{\dot x} = \mathbf{F(x)} \qquad \mbox{ with } \mathbf{x}(0) = \mathbf{x}_0.
\]
and let
\[
\phi(t;\mathbf{x}_0)
\]
represent the flow of this system through the initial point $\mathbf{x}_0$.\\
\pause
Then the ratio
\[
\frac{|| \phi(t;\mathbf{y}_0) - \phi(t;\mathbf{x}_0) || }{ || \mathbf{y}_0 - \mathbf{x}_0 || } 
\] quantifies the rate at which trajectories spread. 

\end{frame}


%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Calculating sensitivity...}
isn't so easy. 
\pause
We want to understand the divergence of two nearby trajectories:


\[
\phi(t;\mathbf{y}_0) - \phi(t;\mathbf{x}_0) \approx D_x \phi(t;\mathbf{x}_0) (\mathbf{y}_0 - \mathbf{x}_0), 
\]

but $D_x \phi(t;\mathbf{x}_0)$ is usually hard to compute! 

\bigskip
\pause
For any curve of initial conditions $\mathbf{x}_s$, define

\[
\mathbf{v}(t) = \partial_{\mathbf s} \phi(t;\mathbf{x}_s) \big|_{s=0},
\]
then $\mathbf{v}(t)$ satisfies the first variation equation
\[
\dot{\mathbf{v}}= D_x \mathbf{F}(\phi(t;\mathbf{x}) \big|_{\mathbf{x}_0}  \mathbf{v}_0 \quad \mbox{ with } \quad \mathbf{v}_0 = \partial_{\mathbf s} \mathbf{x}_s. 
\]
\pause
 \alert{It is $\mathbf{v}(t)$ that will ultimately give us our growth rate.}

\end{frame}


%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{The Lyapunov exponent}

The \textit{Lyapunov} exponent, $\lambda(t)$, is the exponential growth rate measuring sensitivity to initial conditions. It is classically computed as: 
%\[
%|| \delta \mathbf v(t) || \approx exp(\lambda(t) \cdot t)  \mathbf ||\delta \mathbf{v_0}||,
%\]
\[
|| \mathbf v(t) || \approx exp(\lambda \cdot t)  \mathbf || \mathbf{v_0}||,
\]
%where $\delta \mathbf{v}$ is the (first) variation of trajectory $\mathbf{v}_0$.\\

But what happens if we use \alert{power series} approach?

\end{frame}


%%%%%%%%%%%%%%%%%
\section{Examples}

\subsection{Scalar: $y'=a(t) y$}

%%%%%%%
\begin{frame}[fragile]
\frametitle{Ex 1: $\displaystyle y' = \alpha(t) y$}

With \textsc{Maple} \\
\begin{verbatim}
>> restart:
>> Order := 4:
>> alpha := t -> sum(a[k]*t^k,k=0..Order):
>> GROWTH := diff(y(t),t) = alpha(t)*y(t):
>> Yseries := dsolve({GROWTH,y(0)=y[0]},y(t),type='series');
\end{verbatim}
\begin{align*}
y \left( t \right) = & y_0 + a_0 \, y_0t+ \left( 1/2\,{a_0}^2 \, y_0+1/2\,a_1 \,y_0 \right) t^2 + \\
 &  \left( 1/6\,{a_0}^3 \, y_0+1/2\,a_1 \, a_0 \, y_0 +1/3\,a_2 \, y_0 \right) t^3 + O \left( t^4 \right)
\end{align*}
\pause
which we can check
\begin{verbatim}
>> SOLN1 := y[0] * exp(int(alpha(tau),tau=0..t));
>> check := taylor(SOLN1,t=0) - Yseries;

\end{verbatim}

\end{frame}

%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{Why?}
From
\begin{align*}
y \left( t \right) = & y_0 + a_0 \, y_0t+ \left( 1/2\,{a_0}^2 \, y_0+1/2\,a_1 \,y_0 \right) t^2 + \\
 &  \left( 1/6\,{a_0}^3 \, y_0+1/2\,a_1 \, a_0 \, y_0 +1/3\,a_2 \, y_0 \right) t^3 + O \left( t^4 \right)
\end{align*}
we can find
\begin{align*}
\partial_{y_0} y(t) = & 1+a_0 t+ \left( 1/2\,{a_0}^{2}+1/2\,a_1 \right) {t}^{2}+\\
 & \left( 1/6\,{a_0}^{3}+1/2\,a_1 \, a_0+1/3\,a_2 \right) t^3 +O \left( t^4 \right)
\end{align*}
Sensitivity to initial conditions! 
\pause
Which we can easily compute...\\
\begin{verbatim}
>> Yy0 := taylor(diff(Yseries,y_0),t=0);
\end{verbatim}
%\begin{align*} & = 1+a_{{0}}t+ \left( 1/2\,{a_{{0}}}^{2}+1/2\,a_{{1}} \right) {t}^{2}+ \\
%& \qquad \left( 1/6\,{a_{{0}}}^{3}+1/2\,a_{{1}}a_{{0}}+1/3\,a_{{2}} \right) t^{3}+O \left( {t}^{4} \right) 
%\end{align*}
\end{frame}

%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{The Lyapunov exponent}
We have
\[
\partial_{y_0} y \left( t \right) =  \sum_{n=0}^\infty f_n(y_0,\ldots,y_{n-1}) t^n 
\]
and so $\lambda(t)$ is easy to compute: \\
\texttt{>> simplify(taylor(ln(Yy0))/t);}
\[
\lambda(t) = a_0 + 1/2 a_1 t + 1/3 a_2 t^2 + O (t^3)
\]
For our problem,  a direct calculation verifies this:
\[
\frac{1}{t} \int_0^t \alpha(\tau) d\tau).
\]
This time average is the mean coefficient on $[0,t]$\\
\end{frame}


%%%%%%%%%%%%%%%%%%%%

\subsection{2D: Planer system}

\begin{frame}
\frametitle{Ex: A 2D planar system}

Consider the system: 
\begin{align*}
x' & = -y + x(1-x^2-y^2) \\
y' & = \quad x + y(1 - x^2 -y^2)
\end{align*}
which can be recast (and decoupled!) as 
\begin{align*}
r' & = r(1-r^2) \\
\theta' & = 1
\end{align*}
with a stable orbit at $r=1.$

\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Ex: A 2D planar system: exponents?}

With $ \dot{r} = r(1-r^2)$ and $\dot{\theta} = 1$, the variational problem is simple...
\[ \dot{\mathbf v}= 
\begin{bmatrix}
 1-3 r^2  & 0 \\
0 & 0
\end{bmatrix}_{r=1} \mathbf{v_0}
\]
Remember that it's $\mathbf v$ that we need to find, and it provides an (approximate) basis of $D_x \phi$. \bigskip
\pause

So $v_r = \exp(-2t) v_0$, and the lyapunov exponent is -2 in the radial direction (and 0 in the tangential direction).  

\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[fragile]
\frametitle{Ex: A 2D planar system: series}

We can use power series to explore the original cartesian system.

\begin{verbatim}
>> restart:
>> Order := 2:
>> F1 := (x,y) -> -y+x*(1-x^2-y^2);
>> F2 := (x,y) -> x+y*(1-x^2-y^2);
>> Planar := [diff(x(t),t) = F1(x(t),y(t)) , ...
              diff(y(t),t) = F2(x(t),y(t))];
>> ART := dsolve(Planar,[x(t),y(t)],type = 'series);
>> BOB := convert(subs({x(0)=x,y(0)=y},ART),polynom);
>> JacP := {x,y} -> Jacobian(BOB,[x,y]);
>> SingularValues(JacP(1,0));
\end{verbatim}

\end{frame}

%%%%%%%%%%%%%%%%%%%



\begin{frame}[fragile]
\frametitle{Ex: A 2D planar system: series}
We were able to compute  singular values as a time series, 

\begin{verbatim}
>> check := SingularValues(JacP(1,0)); 
>> plot( ln(check[1](t))/t , t=0..0.0001);
\end{verbatim}
\begin{align*}
exp(\lambda_1 t) \approx \sqrt{3*t^2+1-2*t+2*t*\sqrt{2*t^2-2*t+1}  } \\
exp(\lambda_2 t) \approx \sqrt{3*t^2+1-2*t-2*t*\sqrt{2*t^2-2*t+1}  }
\end{align*}

\pause
or evaluate the Jacobian at  a specific time and then compute the singular values.  
Either way:\\
at \texttt{t = 1e-7}, $\lambda_1 \approx$ \texttt{ 0.00000005} and $\lambda_2 \approx $ \texttt{ -2.000000015}
\end{frame}



%%%%%%%%%%%%%%%%%%%
\subsection{3D: Lorenz}

\begin{frame}
\frametitle{3D: Lorenz}

And now for the real work...

\end{frame}

%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile]
\frametitle{3D: Lorenz}

With \textsc{Maple}  and \textsc{Matlab}\\


\end{frame}
%%%%%%%%%%%%%%%%%%%




\begin{frame}
\frametitle{Ex: Lorenz}

Uhoh....


\includegraphics[width=5in,height=2.5in]{spec_lyap}

Sprott gets: $\left{ 0.966, 0, -14.6\right} $.  At least our trace, 13.6667, is right.

\end{frame}

\begin{frame}
\frametitle{Ex: Lorenz}

\includegraphics[width=5in,height=2.5in]{growth_1}

\end{frame}




%%%%%%%%%%%%%%%%%

\begin{frame}<1->
\frametitle{The PROCESS}

\begin{itemize}
\item[] Generate a polynomial (quadratic) system.(Reverse Polish)\\
\item[] Compute analytic series approximation $Y(t;ic)$
\item[] Differentiate to find  sensitivity to IC $S(t;ic) := \partial_{ic} Y(t;ic)$
\item[] Build IC series $\left{ y_0(T) \right}$ 'on' attractor.
\item[] Choose tau for local continuation length
\item[20] WHILE T$<$ Tmax
\item[\phantom{20}] \quad Evaluate $D = svd(S(tau;y_0(T))$
\item[\phantom{20}] \quad Compute $\lambda(T)$ (or mean) for local time interval
\item[\phantom{20}]\quad Advance: $T = T + r$
\item[] GOTO 20
\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%

\section{Conclusion}

\begin{frame}
\frametitle{Questions?}


The community typically evolves trajectory and  linearized variational problem for many time steps to accumulate growth, then re-orthogonalizes. 
\bigskip


Instead
 of answering questions, we now have many more:  \\
\pause

\medskip 
A direct TNB frame calculation seems clunky, and the resulting power series computation messy.

\medskip

Is there a trick we can play to rotate our local SVD frame onto the physical frame? It seems so, but ...
\medskip

\pause

Another series approach might be to compute  
\[
|| \phi(t;\mathbf{x}_0) - \phi(t;\mathbf{y}_0)|| 
\]
directly using the analytic (approximate) integrator.

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Conclusions}

\begin{itemize}
\item Easy to find approximate solution operator as a function of IV.
\item Seems to provide information about local growth of variations.
\item No need to evolve tangent space (a la Wolf)?
\item Non-autonomous? No fear!
\item Non-linear? No problem! 
\end{itemize}
\end{frame}



% SLIDE

\begin{frame}
\begin{center}
\Large{ Thanks! }\\[10mm] 

Questions? \\
\pause
(I sure have a lot!)\\[5mm]

\large{ thelwerj@jmu.edu }\\
\end{center}
 \end{frame}


\end{document}

\begin{frame}
\frametitle{For more information}
\begin{thebibliography}{10}

James Sochacki
\newblock Polynomial ordinary differential equations - examples, solutions, properties.
\newblock {\em Neural Parallel \& Scientific Computations}, 18(3-4):441--450, 2010.


\end{frame}

\end{document}

\begin{frame}
\frametitle{References}
\begin{thebibliography}{10}

\bibitem{PaperXXFehlberg1964}
E.~Fehlberg
\newblock Numerical integration of differential equations by power series expansions, illustrated by physical examples.
\newblock Technical Report NASA-TN-D-2356, NASA, 1964.

%\bibitem{PaperXXNeidinger2010}
%R.D. Neidinger
%\newblock Introduction to automatic differentiation and Matlab object-oriented programming.
%\newblock {\em SIAM Review}, 52(3):545--563, 2010.

\bibitem{PaperXXCarothers2005}
David Carothers et. al.
\newblock Some properties of solutions to polynomial systems of differential equations.
\newblock {\em Electronic Journal of Differential Equations}, 2005:1--18, 2005.

\bibitem{PaperXXWarne2006}
P.~G. Warne et. al.
\newblock Explicit a-priori error bounds and adaptive error control for approximation of nonlinear initial value differential systems.
\newblock {\em Comput. Math. Appl.}, 52(12):1695--1710, 2006.

\end{thebibliography}

\end{frame}

%  SLIDE

\begin{frame}
\frametitle{For more information}
\begin{thebibliography}{10}

\bibitem{PaperXXSochacki2010}
James Sochacki
\newblock Polynomial ordinary differential equations - examples, solutions, properties.
\newblock {\em Neural Parallel \& Scientific Computations}, 18(3-4):441--450, 2010.


\end{thebibliography}

\end{frame}

\end{document}



%%%%%%%%%%


\subsection{$\displaystyle \frac{d}{dt} y = y^3$}

\begin{frame}[fragile]
\frametitle{Ex 2: $\displaystyle \frac{d}{dt} y = y^3$}

Consider 
\[ \frac{d}{dt} y(t) = y^3(t). \]
Now let $v(t) = y^2$ and notice that
\begin{align*}
y' & = y^3  = v \cdot y\\
v' & = 2y \cdot y'  = 2 v^2
\end{align*}
The variational problem for this would be 
%(\textbf{w} is the 2-tensor of 2nd variations)
\[
\delta\mathbf{w}' = \begin{bmatrix}v&y\\0&4v\end{bmatrix} \delta\mathbf{w}
\]
from which we find the exponents to be $\lambda_1 = \frac{1}{t} \int_0^t v(\tau), \lambda_2 = \frac{1}{t} \int_0^t 4v(\tau)$.
\end{frame}

\begin{frame}
\frametitle{Stability}
If $y(0) = y_0,$ then $v(0) = y_0^2,$
and our series solution is (from \textsc{Maple}) 
\begin{align*}
y(t)=&y_0+{y_0}^{3}t+3/2\,{y_0}^{5}{t}^{2}+5/2\,{y_0}^{7}t^3+O\left(t^4\right)\\
v(t)=&y_0^2+2\,{y_0}^4t+4\,{y_0}^6 t^2+8\,{y_0}^{8} t^3+O\left(t^4\right),
\end{align*}
and we have complete stability information.
\[
\partial_{y_0} y(t) = 1+3 y_0^2 t+ (15/2)y_0^4 t^2+(35/2) y_0^6 t^3+O(t^4)
\]
with the corresponding exponent ($\lambda(t) = \frac{1}{t} \ln (\partial_{y_0} y(t) )$),
\[
\lambda(t) = y_0^2 + y_0^4 t + (4/3) y_0^6 t^2 + O (t^3)
\]
at \alert{any point} in the flow!
\end{frame}
