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

\usetheme{Boadilla}
%\usepackage{comment}
%\usepackage{listings}

\begin{document}

\title[Power series and Nonlinear D.E.]{Differential Equations}
\subtitle{through the lens of power series}
\author[thelwerj@jmu.edu]{Roger Thelwell}
\institute[JMU]{James Madison University}
\date{\today}
\date{October 27, 2014}

\frame{\titlepage}

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

% SLIDE

\section{Introduction}


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

\begin{frame}
\frametitle{What are series?}

If you have a sequence $ \{ a_n \}$, then add the terms: 
\[ a_0 + a_1 + a_2 + a_3 + \ldots = \sum_{n=0}^\infty a_n\]
 Sometimes it converges, sometimes it doesn't.

We are usually interested in the series that coverge.

\end{frame}

\begin{frame}
\frametitle{Examples}


\begin{itemize}
\item[]
$\displaystyle \sum_{n=1}^\infty \frac{1}{n}  = 1 + \frac{1}{2} + \frac{1}{3} + \ldots $ %harmonic series
\item[]
$\displaystyle \sum_{n=1}^\infty \frac{1}{n^2}  = 1 + \frac{1}{4} + \frac{1}{9} + \ldots $% Basel problem
\item[]
$\displaystyle \sum_{n=0}^\infty x^n  = 1 + x^1 + x^2 + x^ 3 + \ldots $% geometric |x| <1
\item[]
$\displaystyle \sum_{n=0} a_n x^n = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \ldots$ % power series
\item[]
$\displaystyle \sum_{n=0}^\infty \frac{f^{(n)}(0)}{n!} x^n = f(0) + \frac{f'(0)}{1!} x + \frac{f''(0)}{2!} x^2 + \frac{f^{(3)}(0)}{3!} x^3 + \ldots$
%\item[]
%$\displaystyle \sum_{n=0}^\infty \frac{f^{(n)}(c)}{n!} (x-c)^n = f(c) + \frac{f'(c)}{1!} (x-c) + \frac{f''(c)}{2!} (x-c)^2 + \frac{f(3)(c)}{3!} (x-c)^3 + \ldots$
\end{itemize}


\end{frame}

%%%

\begin{frame}
\frametitle{And more examples}

And series for some common functions:
\[\exp(x) = \sum_{n=0}^\infty \frac{1}{n!} x^n\]
\[\sin(x) = \sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)!} x^{(2n+1)}\]
\[\frac{1}{1-x} = \sum_{n=0}^\infty x^n\]

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Example 1: $\displaystyle \frac{dy}{dx}= y$}
%%%%%%%%%%
%\subsection{Intuition}

%%
\begin{frame}
\frametitle{Example 1}
Consider \[\frac{d}{dx} y(x) = y(x) \quad y(0) = y_0 \]
\end{frame}
 
\begin{frame}
\frametitle{Phase Portrait}
What do you think the solutions of 
\[\frac{d}{dx} y(x) = y(x) \quad y(0) = y_0 \]
to look like?\\

Let's use a \alert{phase portrait} to get some intuition about typical solutions.
\end{frame}

%%

\begin{frame}<1->
\frametitle{An exact solution}

We've got
\only<1>{ \begin{align*}
\frac{dy}{dx} &= y  \\
\end{align*}
Let's try to integrate it.
}
\only<2>{
\begin{align*}
\frac{dy}{dx} &= y  \\
\frac{1}{y} \, \frac{dy}{dx} &= 1 \\
\int \frac{1}{y} \, \frac{dy}{dx} \, dx &= \int 1 \, dx \\
\int \frac{1}{y} \, dy &= \int 1 \, dx \\
\ln(y) &=  x + C  \quad \alert{careful!} \\
\mbox{So  } y(x) & = \exp(x+C) = K \exp(x). \\
\end{align*}
Using the initial condition $y(x=0) = y_0$,\\ we have 
\[ \alert{y(x) = y_0 \exp(x) = y_0  e^x} \]
}
\end{frame}


%%

%\subsection{Power Series}

\begin{frame}
\frametitle{A power series solution?}
How do we solve using power series?

Find coefficients $a_n$ so that
\[y(x) = \sum_{n=0}^\infty a_n x^n = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \ldots\] 

to satisfy 
\[\frac{d}{dx} y(x) = y(x)\]

\end{frame}

\begin{frame}<1->
\frametitle{Power series}

If \[y(x) =  a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \ldots =  \sum_{n=0}^\infty a_n x^n, \qquad \]
 then 
\only<1>{\[ \frac{d}{dx} y = \quad a_1 + 2 a_2 x + 3 a_3 x^2 +  \ldots  = \sum_{n=0}^\infty \alert{???} x^n  \]} 
\only<2->{\[ \frac{d}{dx} y = \quad a_1 + 2 a_2 x + 3 a_3 x^2 +  \ldots = \sum_{n=0}^\infty (n+1) a_{n+1} x^n\]}
\only<3>
{And substitute into \[\frac{d}{dx} y(x) = y(x)\]}

\end{frame}

%%

\begin{frame}
\frametitle{Power series}
We get
\[ \sum_{n=0}^\infty (n+1) a_{n+1} x^n = \sum_{n=0}^\infty a_n x^n \]
\pause
If we compare coefficients for each power of $x,$ we see
 \[ (n+1) a_{n+1} = a_n \implies \alert{a_{n+1} = \frac{1}{n+1} a_n}\]
So,
\begin{align*}
a_1 &= \frac{1}{1} a_0  = a_0 \\
a_2 &= \frac{1}{2} a_1 = \frac{1}{2} ( a_0) = \frac{a_0}{2!} \\
a_3 &= \frac{1}{3} a_2 = \frac{1}{3} \left( \frac{a_0}{2!} \right) = \frac{a_0}{3!} \\
\end{align*}
and \ldots
\pause 
\[ a_n = \frac{1}{n!} a_0\]

\end{frame}

\begin{frame}<1->
\frametitle{Power series}
Since \[ a_n = \frac{1}{n!} a_0,\]
then 
\only<1>{\[y(x) = \sum_{n=0}^\infty a_n x^n = \sum_{n=0}^\infty a_0 \frac{1}{n!} x^n = \alert{???}\]}
 \only<2>{\[y(x) = \sum_{n=0}^\infty a_n x^n = \sum_{n=0}^\infty a_0 \frac{1}{n!} x^n =a_0 \exp(x)\]}
\pause
Again using $y(x=0) = y_0,$ we see that $a_0 = y_0$, and recover the exact solution
\[
y(x) = y_0 \exp(x)
\]
Usually, though, the function is harder identify.
\end{frame}


\section{Example 2: $\displaystyle \frac{d}{dx} y = y (1-y)$}

\begin{frame}
\frametitle{Example 2:}
Now consider 
\[\frac{d}{dx} y = y (1-y),\]
the \alert{logistic equation}.
\end{frame}

\begin{frame}
\frametitle{Phase Portrait}

What do you think solutions of 
\[\frac{d}{dx} y = y ( 1-y ) \quad y(0) = y_0 \]
to look like?\\

Let's use a \alert{phase portrait} again to get some intuition about typical solutions.
\end{frame}


% \subsection{Exact solution}

\begin{frame}[fragile]
\frametitle{Exact solution}
\only<1>{
Let's solve $\frac{d}{dx} y = y ( 1-y )$ by integrating directly.
\begin{align*}
\frac{d}{dx} y &= y ( 1-y ) \\
\frac{1}{y( 1-y )} \frac{d}{dx} y   &= 1 \\
\int \frac{1}{y ( 1-y )} \, dy  &= \int 1 \, dy \\
%\int \frac{A}{y} + \frac{B}{1-y} \, dy &= \int 1 \, dy \\
%\int \frac{1}{y} + \frac{1}{1-y} \, dy &= \int 1 \, dy \\
%\ln(y) - \ln(1-y) & = x + C, \\
%\mbox{so  }\frac{y}{1-y} & = \exp(x+C) = K \exp(x) \\
%\mbox{and  }y(x=0) = y_0 \implies  K = \frac{y_0}{1-y_0}
\end{align*}
}
\only<2>{
Lets solve $\frac{d}{dx} y = y ( 1-y )$ by integrating directly.
\begin{align*}
\frac{d}{dx} y &= y ( 1-y ) \\
\frac{1}{y( 1-y )} \frac{d}{dx} y   &= 1 \\
\int \frac{1}{y ( 1-y )} \, dy  &= \int 1 \, dy \\
\int \frac{\alert{A}}{y} + \frac{\alert{B}}{1-y} \, dy &= \int 1 \, dy \\
%\int \frac{1}{y} + \frac{1}{1-y} \, dy &= \int 1 \, dy \\
%\ln(y) - \ln(1-y) & = x + C, \\
%\mbox{so  }\frac{y}{1-y} & = \exp(x+C) = K \exp(x) \\
%y(x) = \frac{1}{1+ K \exp(-x)} \mbox{ with } K = \frac{y_0}{1-y_0}
\end{align*}
}
\only<3>{
Lets solve $\frac{d}{dx} y = y ( 1-y )$ by integrating directly.
\begin{align*}
\frac{d}{dx} y &= y ( 1-y ) \\
\frac{1}{y( 1-y )} \frac{d}{dx} y   &= 1 \\
\int \frac{1}{y ( 1-y )} \, dy  &= \int 1 \, dy \\
%\int \frac{A}{y} + \frac{B}{1-y} \, dy &= \int 1 \, dy \\
\int \frac{1}{y} + \frac{1}{1-y} \, dy &= \int 1 \, dy \\
\ln(y) - \ln(1-y) & = x + C, \\
\mbox{so  }\frac{y}{1-y} & = \exp(x+C) = K \exp(x) \\
y(x) = \frac{1}{1+ K \exp(-x)} & \mbox{ with } K = \frac{y_0}{1-y_0}
\end{align*}
}
\end{frame}

%\subsection{Power Series solution}

\begin{frame}
\frametitle{Power series}
Let's now solve $\frac{d}{dx} y = y ( 1-y )$ using power series. Again, assume 
\[
y(x) = \sum_{n=0}^\infty a_n x^n.
\]
Substituting,
\[ \sum_{n=0}^\infty (n+1) a_{n+1} x^n = \left( \sum_{n=0}^\infty a_n x^n \right) \cdot \left( 1 - \sum_{n=0}^\infty a_n x^n \right)\]
\vfill

\end{frame}

\begin{frame}
\frametitle{aside: Products}
Power series are easy to add, subtract, differentiate and integrate  - do it term by term.\\
\bigskip
\pause
If $\displaystyle A = \sum_{n=0} a_n x^n $ and $\displaystyle B = \sum_{n=0} b_n  x^n $, what is $A \cdot B$?
\pause

\begin{align*}
(a_0 + & a_1 x + a_2 x^2+  a_3 x^3 + \ldots) \cdot (b_0 + b_1 x+ b2x^2 + b3 x^3+ \ldots)  \\
& = (a_0+b_0) +(a_0 b_1 + a_1 b_0) x +  (a_0 b_2 + a_1 b_1 + a_2 b_0) x^2  \ldots  \\
& \qquad + (a_0 b_3+ a_1 b_2 + a_2 b_1 + a_3 b_0) x^3 + \ldots \\
&= \sum_{n=0}  \left[ \sum_{i=0}^n a_i b_{n-i} \right] x^n 
\end{align*}
We'll call this the \alert{Cauchy Product}
 
\end{frame}

\begin{frame}[fragile]
\frametitle{aside: \textsc{Matlab}}

In \textsc{Matlab}, \[\sum_{n=0}^{degree} \sum_{i=0}^n a_i b_{n-i} \]
is
\begin{verbatim}
function cvec = cauchy_product(avec,bvec,degree)

for n = 1:degree + 1
    for i = 1:n
        j = n-i+1;
        cvec(n) = avec(i)*bvec(j)+cvec(n); 
    end
end
\end{verbatim}
\end{frame}

\begin{frame}
\frametitle{Power series}
Back to it:
\begin{align*}
\sum_{n=0}^\infty (n+1) a_{n+1} x^n &= \left( \sum_{n=0}^\infty a_n x^n \right) \cdot \left( 1 -\sum_{n=0}^\infty a_n x^n \right) \\
&= \sum_{n=0}^\infty a_n x^n - \sum_{n=0}^\infty \sum_{i=0}^n a_i a_{n-i} x^n)  \\
%& = \sum_{n=0}^\infty  \left( a_n - \sum_{i=0}^n a_i a_{n-i}  \right) x^n  \\
& =  \sum_{n=0}^\infty  \left( a_n - c_n \right) x^n 
\end{align*}
where $c_n$ is the result of the Cauchy product.
\pause
We can again write  a recursion relation:
\[
(n+1) a_{n+1} = a_n - \sum_{i=0}^n a_i a_{n-i} 
\]
or 
\[
\alert{a_{n+1} =\frac{1}{n+1} \left( a_n - \sum_{i=0}^n a_i a_{n-i} \right) }
\]


\end{frame}

\begin{frame}[fragile]
\frametitle{\textsc{Matlab} solution}

\textsc{Matlab} makes it easy.  Since $a_0 = y_0$, we need
\[
a_{n+1} =\frac{1}{n+1} \left( a_n - \sum_{i=0}^n a_i a_{n-i} \right)
\]
\begin{verbatim}
   function a = solve_logistic(a0,degree)

   for n =1:degree+1
      c(n) = cauchy_product(a,a,n);
      a(n+1) = 1/(n+1) * ( a(n) - c(n) );
    end

\end{verbatim}
\end{frame}



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

\section{Example 3: $\displaystyle \frac{d}{dx} y = \sin(y)$}
\begin{frame}
\frametitle{Example 3}

Let's explore
\[
\frac{d}{dx} y = \sin(y)
\]
\end{frame}

% \subsection{Intuition}
\begin{frame}
\frametitle{Intuition}
Consider \[\frac{d}{dx} y = \sin(y) \quad y(0) = y_0 \]
What do we expect?\\

\alert{Phase Portrait} 

\end{frame}


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

%\subsection{Exact solution}

\begin{frame}<1->
\frametitle{Exact solution}
\only<1->{
\begin{align*}
\frac{d}{dx}y = sin(y)\\
\int \frac{1}{\sin(y)} \, dy &= \int  1 \, dx
\end{align*}
\vfill
}
\only<2>{
These are getting a little tedious:\\
\textsc{Wolfram Alpha} says:
\texttt{ y(x) = 2 arccot(exp(C-x))}
\vfill
}
\end{frame}

\begin{frame}
\frametitle{ok...}
\begin{align*}
\frac{d}{dx} y &= \sin(y)\\
\int \frac{1}{\sin(y)} \, dy &= \int  1 \, dx\\
% \int \frac{\sin(y)}{\sin^2(y)} \, dy &= \int  1 \, dx\\
\int \frac{\sin(y)}{1-\cos^2(y)} \, dy &= \int  1 \, dx \\
\int \frac{-1}{1-u^2} \, du &= \int  1 \, dx \\
\int \frac{A}{1-u} + \frac{B}{1+u} \, du &= \int  1 \, dx \\
1/2 \ln((1-u)/(1+u)) &= x + C \\
\frac{1-\cos(y)}{1+\cos(y)} &= K \exp(2x) \\
\end{align*}
\end{frame}

%\subsection{Power Series solution}

\begin{frame}
\frametitle{Power series}

We have \[ \frac{d}{dx} y = \sin(y) \quad y(0) = y_0\]

We do have a series for $\sin( \bullet ),$ but ...

Let's turn this into a polynomial system,

Let \[
v_1(x) = y(x), \quad v_2(x) = \sin(y(x)),  \mbox{\quad and \quad}
v_3(x) = \cos(y(x))
\]
\end{frame}

%SLIDE

\begin{frame}
\frametitle{The polynomial system}
Taking $v_1(x)= y, \quad v_2(x) = \sin(y)$ \quad and  \quad $v_3(x) = \cos(y)$, then

\begin{align*}
v_1'  &=  1 \cdot  y'  \\
        &=    v_2   \qquad \qquad v_1(0) =  y_0 \\
v_2'  &= v_3 \cdot  y' \\  
        &=  v_2 \; v_3   \qquad v_2(0) = \sin(y_0)  \\
v_3'  &=  -v_2 \cdot  y'  \\
        &=  -v_2^2   \qquad v_3(0) = \cos(y_0)
\end{align*}
We can solve this system with series recursion, just as before.\\  
\pause
\medskip
We can also consider the geometry....

\end{frame}


%SLIDE



\begin{frame}
\frametitle{Geometry}
\textsc{Matlab} allows us to see the structure.
%\begin{figure}[b]
%\vspace*{-4mm}
%\includegraphics[width=4in,height=2.8in]{helix.png}
%\end{figure}
\end{frame}

\begin{frame}
\frametitle{Decoupling}
Since $v_1(x)= y, \quad v_2(x) = \sin(y)$ \, and  \, $v_3(x) = \cos(y)$,
\begin{align*}
v_1'' = y' &= v_2 \; v_3 \\
v_1''' = (v_2 \; v_3)' &= v_2 v_3' + v_2' v_3  \\
 &= v_2(-v_2)^2+ (v_2  v_3) v_3  \\
&=  -v_2^3 + v_2 v_3^2 \\
&=  - v_2^3 + v_2 (1-v_2^2) \\
&= -2(v_2)^3  + v_2 
\end{align*}
but $v_2 = v_1'=y',$ and we've recast  
\[ y' = \sin(y) \quad  \mbox{as} \quad  y''' - y' + 2(y')^3 = 0\]

Gr\"{o}bner basis theory says this can always be done.
\end{frame}
% SLIDE

\section{Example 4: $\displaystyle \frac{d}{dx} y = y \, \sin(y)$}

\begin{frame}
\frametitle{Example 4}
And one more to think about:
\[
\frac{d}{dx} y = y \sin(y)
\]
\end{frame}

\begin{frame}
\frametitle{Intuition?}
\[
\frac{d}{dx} y = y \sin(y)
\]
\end{frame}

\begin{frame}[fragile]
\frametitle{Analytic soln?}

\textsc{Maple}

\begin{verbatim}
> dsolve(diff(y(x),x) = y(x)*sin(y(x)));
             y(x)
           /
           |          1
     x -   |      ---------- d_a + _C1 = 0
           |      _a sin(_a)
           /
\end{verbatim}
\end{frame}


\begin{frame}
\frametitle{Analytic soln?}
\begin{center}
\textsc{Wolfram Alpha}\\
\includegraphics[width=11cm,height=6cm]{wolfram_yp_siny.png}
\end{center}
\end{frame}

%%%%%

\begin{frame}
\frametitle{Series?}
\only<1>{Taking $v_1= y  , v_2 = y \sin(y), \quad v_3 = \alert{????} , v_4 = \alert{????}$, \\ }
\only<2->{Taking $v_1= y  , v_2 = y \sin(y), \quad v_3 = \sin(y), v_4 = \cos(y)$, } then
\only<3->{
\pause
\begin{align*}
v_1'  &=   1  \cdot  y'  &\\
        &= v_2  & v_1(0) &=  y_0 \\
v_2'  &=  (y \cos(y) + \sin(y) )\cdot  y' &\\ 
         &=  (v_1 \; v_4 + v_3) \; v_2 & v_2(0) &= y_0 \sin(y_0)  \\
v_3' & = \cos(y)   \cdot  y'  &\\
       &=  v_2 v_3  & v_3(0) &= \sin(y_0) \\
v_4' &= -\sin(y) \cdot y' & \\
       &= -v_2^2  & v_4(0) &= \cos(y_0) 
\end{align*}
Let $v_5 = v_1 v_4$ to reduce to a quadratic system.
}
\end{frame}


\section{Theory}

\begin{frame}<1->
\frametitle{THEORY}

We've been able to recast EVERY ODE that we've considered as a polynomial system.\\
\bigskip
\textbf{U\$D 50 CASH PRIZE} to the first person to send me an ODE with analytic solution that \textbf{CAN'T} be recast as a polynomial system.

\end{frame}

% SLIDE 15

%\subsection{Decoupling}
\begin{frame}
\frametitle{THEORY}
%\framesubtitle{Decoupling}

Carothers et. al. 2005 \cite{PaperXXCarothers2005}
\begin{theorem} 

A function $v$ is the
solution to an arbitrary component of  a polynomial system
of differential equations  if and only if for some $n$ there is a polynomial $Q$
in $n+1$ variables so that $Q(v, v^\prime, \cdots, v^{(n)})=0.$ 
\end{theorem}

\bigskip

\only<2->{\begin{center}
\alert{This implies that the motion of our moon may be described \\
\textbf{without}\\
reference to the earth, sun, or any other planets!} \end{center} }

\end{frame}



% SLIDE 15

%\subsection{Error Bound}
\begin{frame}
\frametitle{THEORY}
\framesubtitle{Error Bound}

Warne et. al. 2006  \cite{PaperXXWarne2006}  \\
\bigskip
If we have  (at $a = 0$) a system
${\bf x} ' (t) = {\bf h}({\bf x} (t)),$  ${\bf x} (0) = {\bf b},$
then
\begin{equation}
\left\| {\bf x}(t) - \sum\limits_{k = 0}^n {\bf x}_k t^k  \right\|_{\infty} \leq \frac{\left\| {\bf b} \right\|_{\infty} \left| {Kt} \right|^{n + 1} }{1 - \left| {Mt} \right|} ~~~{\text{for}}\,\,\,\,\,\,\,m \geq 2\
\end{equation}

Where the parameters K and M  depend on immediately observable quantities of the original system; \\[3mm]

$M$ is the largest row sum of coefficients, and
$K=(m-1) c^{m-1}$, where $c = \max \{ 1,||\mathbf{b}||_\infty \}$ and $m = deg(\bf{h}).$ 

\end{frame}

% SLIDE 15

\section{Conclusion}

\begin{frame}
\frametitle{Conclusions}
\begin{itemize}
\item<1-> Easily compute \textit{arbitrarily high order Taylor coefficients}
\item<2-> The tools can \textit{solve highly nonlinear and stiff problems}
\item<3-> Semi-analytic methods and
\item<4-> \textit{interpolation free} to machine capability (error and calculation)
\end{itemize}
\end{frame}

% SLIDE

\begin{frame}
\begin{center}
\Large{Thanks to all the faculty and students of}\\
\includegraphics[width=5cm,height=5cm]{VMI_logo.png}\\
 \Large{and in particular to LTC Dr. Chalishajar.}\\
\end{center}
 \end{frame}

\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}
