summary.tex

\documentclass[10pt]{article}
\usepackage{amsmath}
\begin{document}

\section{System}
\begin{align}
u_t &= u_{xx} + (1-v)F(w) \\
v_t &= \epsilon v_{xx} + (1-v)F(w)\\
w &= hu+(1-h)v
\end{align}
with $h \in [0,1], \epsilon \in (0,1]$.

We take a discontinuous nonlinearity
\begin{equation}
F(w) = \left\{\begin{array}{ll}
\mathrm{exp}\left(Z\left\{\frac{w-h}{\sigma+(1-\sigma)w} \right\} \right) & w \geq w^*\\
0 & w < w^*
\end{array} \right.
\end{equation}


\section{Waves}
\begin{align}
&U_{\xi\xi} + cU_\xi + (1-V)F(W) = 0\\
&\epsilon V_{\xi\xi} + cV_\xi + (1-V)F(W) = 0\\
&W = hU+(1-h)V
\end{align}

At $-\infty$, $(U,V) \to (1,1)$, and at $\infty$, $(U,V) \to (0,0)$.

In the region $W < w^*$ where $F(W) = 0$, these equations do not involve $U$ and $V$ at all---they are differential equations for $U_\xi$ and $V_\xi$.
As a result, integration may not hit $(0,0)$. To fix this, we integrate the system.

Let
\begin{equation}
z = \frac{1}{c}\int_\xi^\infty (1-V)F(W)
\end{equation}
Clearly $z=0$ at $\infty$.  After integration we will get $z=1$ at $-\infty$.

This gives
\begin{equation}
z_\xi = -\frac{1}{c}(1-V)F(W)
\end{equation}
so
\begin{align}
&U_{\xi\xi}+cU_\xi-cz_\xi = 0\\
&\epsilon V_{\xi\xi}+cV_\xi-cz_\xi = 0
\end{align}
These may now be integrated.
\begin{align}
&U_\xi + cU - cz = q \\
&\epsilon V_\xi + cV - cz = r
\end{align}
where $q$ and $r$ are constants.
At $\infty$, each of $U_\xi$, $U$, and $z$ go to 0, so $q=0$.
Likewise $r=0$.
At $-\infty$, $U_\xi = 0$ and $U = 1$, so we get $z=1$.
\begin{align}
&U_\xi + cU - cz = 0 \\
&\epsilon V_\xi + cV - cz = 0
\end{align}

Thus, our end behavior is $(U,V,z) \to (1,1,1)$ at $-\infty$, and $(U,V,z) \to (0,0,0)$ at $\infty$.

\subsection{$\epsilon = 0$}
In this case, we get $V = z$, and the system is two-dimensional.
\begin{align}
U_\xi &= c(V-U) \label{front_eps0} \\
V_\xi &= - \frac{1}{c}(1-V)F(W)
\end{align}

For shooting, we want to find the positive eigendirection at $(1,1)$.
The system linearizes to
\[
\begin{bmatrix}
P\\Q
\end{bmatrix}'
=
\begin{bmatrix}
-c & c \\
-K_u/c & -K_v/c
\end{bmatrix}
\begin{bmatrix}
P\\Q
\end{bmatrix}
\]
where $K = (1-v)F(w)$.

At $(U,V) = (1,1)$ we have $K_u = 0$ and $K_v = - e^{(1-h)Z}$, so the matrix is
\[
\begin{bmatrix}
-c & c \\
0 & e^{(1-h)Z}/c
\end{bmatrix}
\]
The positive eigenvalue is $e^{(1-h)Z}/c$, with eigenvector $[1,1+e^{(1-h)Z}/c^2]$.

\subsection{$\epsilon \neq 0$}
Here the system is
\begin{align}
U_\xi &= -cU + cz \\
V_\xi &= \frac{1}{\epsilon}(-cV + cz) \\
z_\xi &= -\frac{1}{c}(1-V)F(W)
\end{align}
This linearizes to
\[
\begin{bmatrix}
P\\Q\\R
\end{bmatrix}'
=
\begin{bmatrix}
-c & 0 & c \\
0 & -c/\epsilon & c/\epsilon \\
-K_u/c & -K_v/c & 0
\end{bmatrix}
\begin{bmatrix}
P\\Q\\R
\end{bmatrix}
\]

At $(1,1)$ the matrix becomes
\[
\begin{bmatrix}
-c & 0 & c \\
0 & -c/\epsilon & c/\epsilon \\
0 & e^{(1-h)Z}/c & 0
\end{bmatrix}
\]
This has one positive eigenvalue, which is quite a mess.
We found these with computer algebra.
\[
    -1/2 (c e^{Z h - Z} - \sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon e^{Z h -
    Z}}) e^{-Z h + Z}/\epsilon
\]
Its eigenvector is:
\begin{multline*}
    [1,\\
        -\frac{1}{2} (c^2 \epsilon e^{3 Z h - Z} - c^2 e^{3 Z h -
        Z} - ((\sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon e^{Z h - Z}} c + 2) \epsilon \\
        - \sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon e^{Z h -
        Z}} c) e^{2 Z h}) e^{-2 Z h}/\epsilon^2, \\
        \frac{1}{2} ((2 c \epsilon - c) e^{Z h} +
\sqrt{c^2 e^{Z h} + 4 \epsilon e^Z} e^{1/2 Z h}) e^{-Z h}/(c \epsilon)]
\end{multline*}


\section{Finding $c$}
We find $c$ by integrating the ODE from a point close to $(1,1,1)$, in the positive eigendirection.
At different values of $c$, we end up at different points, and there's only one correct value where we hit $(0,0,0)$.
We can tell whether our value is too high or low by where the terminus is: if $c$ is too low, we finish with $U,V < 0$, and vice versa.  This lets us perform a simple binary search to approach the correct value of $c$.

Starting in an eigendirection lets us approximate the ODE backward in time with simple exponential decay.
Since the $(1,1,1)$ fixed point is a saddle, backward integration eventually leaves along another eigendirection, so we do this instead.


\section{Linearization}
First, moving frame version:
\begin{align}
u_t &= u_{\xi\xi} + c u_\xi + (1-v)F(w) \\
v_t &= \epsilon v_{\xi\xi} + c v_\xi+ (1-v)F(w)
\end{align}
We linearize about the front with perturbation $(p,q)$.  This gives:
\begin{align}
p_t &= p_{\xi\xi} + c p_\xi + K_u p + K_v q \\
q_t &= \epsilon q_{\xi\xi} + c q_\xi+ K_u p + K_u q
\end{align}
where $K = (1-v)F(w)$.
$K_u$ and $K_v$ are evaluated at the front.

For $w \geq w^*$, $K_u$ and $K_v$ are each messy functions of $u$ and $v$; we find them with computer algebra.
For $w < w^*$, since $F = 0$ uniformly we have $K_u = K_v = 0$.

This gives us the operator
\[
L = \begin{bmatrix}
\partial_{\xi\xi} + c \partial_\xi + K_u  &  K_v  \\
K_u  &  \epsilon \partial_{\xi\xi} + c \partial_\xi + K_v
\end{bmatrix}
\]
We wish to solve $LX = \lambda X$, i.e.
\begin{align}
\lambda p &= p_{\xi\xi} + c p_\xi + K_u p + K_v q \\
\lambda q &= \epsilon q_{\xi\xi} + c q_\xi+ K_u p + K_v q
\end{align}

\subsection{$\epsilon = 0$}
We will study the $\epsilon=0$ case for now, to simplify things.
We can rewrite this as a first-order system $X' = AX$ by introducing $p_x = s$.
\begin{align}
p_\xi &= s \\
s_\xi &= \lambda p - c s - K_u p - K_v q \\
q_\xi &= (\lambda q - K_u p - K_v q)/c
\end{align}
Writing as a matrix, where $K_u$ and $K_v$ depend on the front,
\[
X = \begin{bmatrix}p\\s\\q\end{bmatrix},
A = \begin{bmatrix}
0 & 1 & 0 \\
\lambda - K_u  &  -c  &  -K_v \\
- K_u/c  &  0  &  (\lambda - K_v)/c
\end{bmatrix}
\]

\subsection{$\epsilon \neq 0$}
In this case, the system becomes four-dimensional, with $p_x = s$ and $q_x = r$.
(these are in non-alphabetical order\ldots)
\begin{align}
p_\xi &= s \\
s_\xi &= \lambda p - c s - K_u p - K_v q \\
q_\xi &= r \\
r_\xi &= (\lambda q - c r - K_u p - K_v q)/\epsilon
\end{align}
Writing the system as $X' = AX$, where $K_u$ and $K_v$ depend on the front,
\[
X = \begin{bmatrix}p\\s\\q\\r\end{bmatrix},
A = \begin{bmatrix}
0 & 1 & 0 & 0 \\
\lambda - K_u  &  -c  &  -K_v  & 0 \\
0 & 0 & 0 & 1 \\
- K_u/\epsilon &  0  &  (\lambda - K_v)/\epsilon  &  - c / \epsilon
\end{bmatrix}
\]



\section{Eigenvalues}
\subsection{$\epsilon = 0$}
\subsubsection{$+\infty$}
At $(u,v) = 0$, we have $K_u = K_v = 0$, so this simplifies to
\[
A = \begin{bmatrix}
0 & 1 & 0 \\
\lambda &  -c  &  0 \\
0 &  0  &  \lambda/c
\end{bmatrix}
\]
This has eigenvalues $(-c \pm \sqrt{c^2 + 4 \lambda})/2$ and $\lambda/c$.  For $\Re \lambda > 0$, the one with negative real part is
\[\lambda_1 = \frac{-c - \sqrt{c^2 + 4 \lambda}}{2}\]
One eigenvector for it is
\[v_1 = [1,(-c - \sqrt{c^2 + 4 \lambda})/2,0] = [1,\lambda_1,0] \]

\subsubsection{$-\infty$}
At $(u,v) = (1,1)$ we get
\[
A = \begin{bmatrix}
0 & 1 & 0 \\
\lambda &  -c  &  e^{(1-h)Z} \\
0 &  0  &  (\lambda + e^{(1-h)Z})/c
\end{bmatrix}
\]
with eigenvalues $(-c \pm \sqrt{c^2 + 4 \lambda})/2$ and $(\lambda + e^{(1-h)Z})/c$.

We're actually interested in the adjoint $-A^T$, whose eigenvalues are the negation of those.
Thus, the unique positive eigenvalue is
\[\lambda_2 = \frac{c + \sqrt{c^2 + 4 \lambda}}{2} = -\lambda_1 \]
One eigenvector (of $-A^T$) for that is
\[
v_2 = \left[1, -\frac{c + \sqrt{c^2 + 4 \lambda}}{2\lambda}, \frac{(c + \sqrt{c^2+4\lambda})}{\lambda e^{Z(h-1)}(c + \sqrt{c^2+4\lambda} + 2 \lambda/c)+ 2 \lambda / c}\right]
\]

\subsection{$\epsilon \neq 0$}
\subsubsection{$+\infty$}
At $(u,v) = 0$, we have $K_u = K_v = 0$, so this simplifies to
\[
A = \begin{bmatrix}
0 & 1 & 0 & 0 \\
\lambda &  -c  &  0  & 0 \\
0 & 0 & 0 & 1 \\
0 &  0  &  \lambda /\epsilon  &  - c / \epsilon
\end{bmatrix}
\]
The eigenvalues with negative real part are:
\begin{align*}
    \lambda_1 &= -\frac{1}{2} (c + \sqrt{c^2 + 4 \epsilon \lambda})/\epsilon\\
    \lambda_2 &= -\frac{1}{2} (c + \sqrt{c^2 + 4 \lambda})
\end{align*}
Their respective eigenvectors are:
\begin{align*}
    v_1 &= [0,0,1,-\frac{1}{2} (c + \sqrt{c^2 + 4 \epsilon \lambda})/\epsilon] = [0,0,1,\lambda_1] \\
    v_2 &= [1, -\frac{1}{2} (c + \sqrt{c^2 + 4 \lambda}), 0, 0] = [1,\lambda_2,0,0]
\end{align*}
We actually need the two with positive real part as well, in order to use the Wronskian to test for mode mixing.
\begin{align*}
    \lambda_3 &= -\frac{1}{2} (c - \sqrt{c^2 + 4 \epsilon \lambda})/\epsilon\\
    \lambda_4 &= -\frac{1}{2} (c - \sqrt{c^2 + 4 \lambda})
\end{align*}
Their associated eigenvectors are:
\begin{align*}
    v_3 &= [0,0,1,\lambda_3] \\
    v_4 &= [1,\lambda_4,0,0]
\end{align*}


\subsubsection{$-\infty$}
At $(u,v) = (1,1)$, we have $K_u = 0$ and $K_v = -e^{(1-h)Z}$.
The matrix becomes:
\[
A = \begin{bmatrix}
0 & 1 & 0 & 0 \\
\lambda &  -c  &  e^{(1-h)Z}  & 0 \\
0 & 0 & 0 & 1 \\
0 &  0  &  (\lambda + e^{(1-h)Z})/\epsilon  &  - c / \epsilon
\end{bmatrix}
\]
%Again, we look for eigenvalues and eigenvectors of $-A^T$, and choose positive real part.
We look for eigenvalues with positive real part.  In this case, we don't work with the adjoint.
We follow the method from Lafortune, Lega, and Madrid.

The first eigenvector is particularly horrific.
It might be worth seeing if another CAS can produce a simpler version, maybe by writing it in terms of the eigenvalue.
\begin{align*}
    \nu_1 &=-1/2 \left(c e^{Z h - Z} - \sqrt{(c^2 + 4 \epsilon \lambda) e^{2 Z h - 2 Z} +
4 \epsilon e^{Z h - Z}}\right) e^{-Z h + Z}/\epsilon
\\
    \nu_2 &= -1/2 c + 1/2 \sqrt{c^2 + 4 \lambda}
\end{align*}
\begin{align*}
    w_1 =& [1, ((\sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} c \epsilon - \sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} c) e^{3/2 Z h} \\
& \quad - (c^2 \epsilon - c^2 -
2 \epsilon \lambda) e^{2 Z h} + 2 \epsilon e^{Z h + Z})/(\sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} \epsilon e^{3/2 Z h} \\
& \quad + (2 c \epsilon^2 -
    c \epsilon) e^{2 Z h}), \\
    &\!-1/2 (c^2 \epsilon e^{3 Z h - Z} - c^2 e^{3 Z h - Z} -
    ((\sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon \lambda e^{2 Z h - 2 Z} + 4 \epsilon e^{Z h
    - Z}} c + 2) \epsilon \\
    & \quad - \sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon \lambda e^{2 Z h -
    2 Z} + 4 \epsilon e^{Z h - Z}} c) e^{2 Z h} + 2 (\epsilon^2 e^{3 Z h - Z} -
    \epsilon e^{3 Z h - Z}) \lambda) e^{-2 Z h}/\epsilon^2, \\
    & 1/2 (2 c \epsilon^2 +
    (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} e^{-1/2 Z h} -
    3 c) \epsilon \\
    & \quad - (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h} +
    4 \epsilon e^Z} c^2 \epsilon e^{-Z h - Z} - \sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} c^2 e^{-Z h - Z} \\ 
    & \quad + (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} \epsilon^2 e^{-Z h - Z} - \sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} \epsilon e^{-Z h - Z}) \lambda) e^{3/2 Z h} \\
    & \quad + (c^3 \epsilon e^{-Z h -
    Z} - c^3 e^{-Z h - Z} + 3 (c \epsilon^2 e^{-Z h - Z} - c \epsilon e^{-Z h -
Z}) \lambda) e^{2 Z h})/\epsilon^3]
\\
w_2 =& [1, \nu_2, 0, 0]
\end{align*}
We also need the eigenvalues with negative real part, for the Wronskian mode-mixing test.
\begin{align*}
    \nu_3 &=-1/2 \left(c e^{Z h - Z} + \sqrt{(c^2 + 4 \epsilon \lambda) e^{2 Z h - 2 Z} +
4 \epsilon e^{Z h - Z}}\right) e^{-Z h + Z}/\epsilon
\\
    \nu_4 &= -1/2 c - 1/2 \sqrt{c^2 + 4 \lambda}
\end{align*}
\begin{align*}
    w_3 =& [1, ((\sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} c \epsilon - \sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} c) e^{3/2 Z h} \\
& \quad + (c^2 \epsilon - c^2 -
2 \epsilon \lambda) e^{2 Z h} - 2 \epsilon e^{Z h + Z})/(\sqrt{c^2 e^{Z h} +
4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} \epsilon e^{3/2 Z h} \\
& \quad - (2 c \epsilon^2 -
    c \epsilon) e^{2 Z h}), \\
    &\!-1/2 (c^2 \epsilon e^{3 Z h - Z} - c^2 e^{3 Z h - Z} +
    ((\sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon \lambda e^{2 Z h - 2 Z} + 4 \epsilon e^{Z h
    - Z}} c - 2) \epsilon \\
    & \quad - \sqrt{c^2 e^{2 Z h - 2 Z} + 4 \epsilon \lambda e^{2 Z h -
    2 Z} + 4 \epsilon e^{Z h - Z}} c) e^{2 Z h} + 2 (\epsilon^2 e^{3 Z h - Z} -
    \epsilon e^{3 Z h - Z}) \lambda) e^{-2 Z h}/\epsilon^2, \\
    & 1/2 (2 c \epsilon^2 -
    (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h} + 4 \epsilon e^Z} e^{-1/2 Z h} +
    3 c) \epsilon \\
    & \quad + (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h} +
    4 \epsilon e^Z} c^2 \epsilon e^{-Z h - Z} - \sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} c^2 e^{-Z h - Z} \\ 
    & \quad + (\sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} \epsilon^2 e^{-Z h - Z} - \sqrt{c^2 e^{Z h} + 4 \epsilon \lambda e^{Z h}
    + 4 \epsilon e^Z} \epsilon e^{-Z h - Z}) \lambda) e^{3/2 Z h} \\
    & \quad + (c^3 \epsilon e^{-Z h -
    Z} - c^3 e^{-Z h - Z} + 3 (c \epsilon^2 e^{-Z h - Z} - c \epsilon e^{-Z h -
Z}) \lambda) e^{2 Z h})/\epsilon^3]
\\
w_4 =& [1, \nu_4, 0, 0]
\end{align*}

\section{Evans function}
\subsection{$\epsilon = 0$}
We perform two integrations.
These depend on a fixed front, a solution of the ODE \eqref{front_eps0}.
We pick left, middle, and right points $L, M, R$, along this front.
The left endpoint is just the point we started shooting from (in our implementation, $\xi=0$).
Currently, we pick the middle as the first point where the average of $u$ and $v$ is below 0.3,
and the right is the first point where $u < 0.001$.
These are somewhat arbitrary.

As $\xi$ changes during the integration, the front gives us different values of $u,v$ in $A$.

We integrate $X' = AX$ from $\xi = R$ to $M$, with initial condition $e^{\lambda_1 (R-M)}v_1$.

We integrate $X' = -A^T X$ from $L$ to $M$, with initial condition $e^{\lambda_1 (M-L)}v_2' = e^{-\lambda_2 (M-L)}v_2'$. The vector $v_2'$ is another eigenvector collinear with $v_2$, but normalized so that $v_1 \cdot v_2' = 1$.
\[v_2' = \frac{v_2}{v_2 \cdot v_1}\]

The Evans function is the dot product of these results at $\xi = M$.
(being careful to do this the correct way: \verb|a.' * b| is what we want,
not \verb|dot(a,b)|, which does \verb|a' * b|)


\end{document}