appendixA.tex

\chapter{Further Mathematical Details}\label{section:appA}
\section{Finite Difference Scheme}
With the $m=500$ spatially discretised points given as $\textbf{x}=[x_1,\cdots,x_m]^T$ where $x_i=(i-1)\Delta x$, with $\Delta x=\frac{1}{m-1}$. This ensures $x_1=0$ and $x_m=1$. Letting $U_i^t$ denote the numerical approximation to $u(x_i,t)$, we use a second-order central difference approximation \cite{finitediff} to evaluate the second-order derivatve $\frac{\partial^2 u}{\partial x^2}(x_i,t)$. Namely,
\begin{equation}
\frac{\partial^2 u}{\partial x^2}(x_i,t)\approx\frac{1}{\Delta x^2}\left(U_{i+1}^t-2U_i^t+U_{i-1}^t\right).
\end{equation}
By denoting $\textbf{U}^t$ as the vector of numerical approximations to $u(x,t)$ at some time $t$ across the whole spatial domain $\textbf{x}$, so that $\textbf{U}^t=\left[U_1^t,\cdots,U_m^t\right]^T$, the numerical approximation of the second-order derivative across the whole spatial domain at some time $t$ can be computed in matrix form as
\begin{equation}
    \frac{\partial^2 u}{\partial \textbf{x}^2}\approx A\textbf{U}^t.
\end{equation}
$A$ is the discrete second-order differential operator, and is given by
\begin{equation}\label{A}
A=\frac{1}{\Delta x^2}\begin{bmatrix}
   -2&  1&  &  & 0\\
   1&  -2&  1&  & \\
   &  \ddots&  \ddots&  \ddots& \\
   &  &  1&  -2& 1\\
   0&  &  &  1& -2
  \end{bmatrix}.
\end{equation}
The sparse nature of $A$ allows for computational advantages when implementing the finite-difference scheme. In order to implement homogeneous Neumann boundary conditions, a first-order central difference approximation is used with `ghost' nodes appended at $x_{-1}$ and $x_{m+1}$. This results in altering entries $A_{1,2}$ and $A_{m,m-1}$ from a $1$ to a $2$, as seen in \eqref{Aneumann}. To implement homogeneous Dirichlet conditions, the first and last rows of $A$ are set to $0$, as seen in \eqref{Adirichlet}. Since homogeneous Dirichlet conditions require $u(x)=0$ at $x=0,1$ for all $t>0$, we also set the initial conditions and kinetic functions to equal $0$ at the end nodes.
\begin{multicols}{2}
\begin{equation}\label{Aneumann}
    \begin{split}
A&=\frac{1}{\Delta x^2}\begin{bmatrix}
   -2&  2&  &  & 0\\
   1&  -2&  1&  & \\
   &  \ddots&  \ddots&  \ddots& \\
   &  &  1&  -2& 1\\
   0&  &  &  2& -2
  \end{bmatrix}.\\
  A &\textit{ with homogeneous Neumann conditions}
    \end{split}
\end{equation}
\break
\begin{equation}\label{Adirichlet}
    \begin{split}
A&=\frac{1}{\Delta x^2}\begin{bmatrix}
   -0&  0& \cdots &  & 0\\
   1&  -2&  1&  & \\
   &  \ddots&  \ddots&  \ddots& \\
   &  &  1&  -2& 1\\
   0&\cdots  &  &  0& 0
  \end{bmatrix}.\\
  A & \textit{ with homogeneous Dirichlet conditions}
    \end{split}
\end{equation}
\end{multicols}


\section{Functional Form of $\text{IC}_1$.}
The form of $\text{IC}_1$ is taken from those used in \cite{gaffmonk}, and given by
\begin{equation}\label{ic1}
    \begin{split}
    u_0(x)&=u_\star+E_ux^7(1-x^2)[A_ux^3+B_ux^2+C_ux+D_u]\\
    v_0(x)&=v_\star+E_vx^5(1-x^2)[A_vx^3+B_vx^2+C_vx+D_v],
    \end{split}
\end{equation}
for $x\in[0,1]$. The parameter values $(a,b)=(0.1,0.9)$ are used, yielding $(u_\star,v_\star)=(1,0.9)$. The coefficients in \eqref{ic1} are given as
\begin{align*}
A_u&=-130.8444445,\ \ B_u=337.0666669,\ \ C_u=-281.6000002,\ \ D_u=75.3777778,\\
A_v&=-170.6666682,\ \ B_v=412.4444479,\ \ C_v=-312.8888910,\ \ D_v=71.1111113,\\
\end{align*}
and $E_u=0.00600$, $E_v=0.00125$.


\section{Derivation of Mean of Skewed Truncated Gaussian Distribution}

Using a result from \cite{skewed}, we have that the $m$-th moment of a random variable $X$ following a truncated skewed Gaussian distribution, on the domain $[a,b]$ is given by
\begin{equation}\label{mth}
\mathbb{E}[X^m]=\sum_{r=0}^mC_m^r\mu^{m-r}\omega^rs_{\rho,r}(u,v).
\end{equation}
Here, $u=\frac{a-\mu}{\omega}$, $v=\frac{b-\mu}{\omega}$, and $C_m^r=\begin{pmatrix}m\\r\end{pmatrix}$ is a binomial coefficient. The function $s_{\rho,r}$ is given by
\begin{equation}\label{s}
    s_{\rho,r}(u,v)=(r-1)s_{\rho,r-2}(u,v)+q_{\rho,r}(u,v),
\end{equation}
with $s_{\rho,0}(u,v)=1$, and the function $q_{\rho,r}$ defined by
\begin{equation}\label{q}
    q_{\rho,r}(u,v)=-\frac{\left[x^{r-1}f_\rho(x)\right]|_u^v}{\left[F_\rho(x)\right]|_u^v}+\frac{2}{\sqrt{2\pi}}\frac{\rho}{\hat{\rho}^r}\frac{\left[\phi(\hat{\rho x})\right]|_u^v}{\left[F_\rho(x)\right]|_u^v}m_{r-1}(\hat{\rho}u, \hat{\rho}v).
\end{equation}
The functions $f_\rho$ and $F_\rho$ denote the pdf and cdf of the standard skew Gaussian distribution,

 To compute the expection, $\mathbb{E}[X]$, we take the first moment, namely $m=1$. Equation \eqref{mth} can therefore be considerably simplified to
\begin{equation}
    \begin{split}
\mathbb{E}[X]&=C_1^0\mu s_{\rho,0}(u,v)+C_1^1\omega s_{\rho,1}(u,v)\\
&=\mu+\omega q_{\rho,1}(u,v),
\end{split}
\end{equation}
since $s_{\rho,0}(u,v)=1$ and $s_{\rho,1}(u,v)=q_{\rho,1}(u,v)$. We also use the fact that $m_1(u,v)=1$, from \cite{skewed}. Evaluating $q_{\rho,1}(u,v)$ from \eqref{q} leads to the desired result presented in \eqref{computetau}.