paper.tex

\documentclass{article}
\usepackage{iclr2021_conference,times}

\input{math_commands.tex}
\usepackage[colorlinks=True]{hyperref}
\usepackage{url}
\usepackage{graphicx}
\graphicspath{{./figures/}}
\DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\usepackage{lipsum}
\usepackage[caption=false,font=footnotesize]{subfig}
\usepackage{optidef}
\usepackage{tabularx}
%\usepackage{IEEEtrantools}

% Front matter
\title{A Finite Difference Approach to Solving the Transmission Line Telegraph Equation}

\author{Christian Y. Cahig
	\thanks{Under the supervision of Engr. Michael S. Villame.} \\
	Department of Electrical Engineering and Technology\\
	Mindanao State University - Iligan Institute of Technology\\
	Iligan City, Philippines \\
	\texttt{\{christian.cahig\}@g.msuiit.edu.ph} \\
}

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\iclrfinalcopy

\begin{document}

\maketitle

\begin{abstract}

In power system dynamics, the propagation of a voltage signal through a transmission line
is described by a hyperbolic partial differential equation known as a telegraph equation.
Although closed form analytic solutions are tractable for special cases,
numerical approaches are preferred in more general scenarios.
This work introduces the use of finite difference approximation of derivatives
such that solving the telegraph equation reduces to iteratively solving a linear algebraic equations.
To aid discussion, we also provide illustrative examples that investigate three aspects of this numerical approach.
First, we empirically show the computational advantage,
in terms of reduced execution time,
in vectorizing the iterative update scheme.
Second, we demonstrate how the discretization of the spatial and the temporal domains
affect the accuracy of the approximation as well as the stability of the numerical method.
Lastly, we simulate the voltage variation of a scenario where a fault occurs at the receiving-end bus.
Codes and files are available in \url{https://github.com/christian-cahig/FDATE}.

\end{abstract}

\section{Introduction}
\label{sec: Introduction}

Consider a transmission line of length $X$ characterized by per-unit-length
series resistance $R$,
series inductance $L$,
shunt conductance $G$,
and
shunt capacitance $C$.
Let $u \left(x,t\right)$ be the instantaneous voltage signal (referred to ground)
at point $x$ along the length of the line at time $t$,
where $0 \leq x \leq X$ and $0 \leq t \leq T$.
We refer to $x = 0$ as the \textit{sending end} and $x=X$ as the \textit{receiving end} of the line.
From elementary transmission line theory,
the propagation of a voltage signal through the line is described by the \textit{telegraph equation}:
\begin{equation}
   \label{eqn: Telegraph eqn full}
   \frac{1}{LC} \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}}
   =
   \frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}}
   +
   \left(\frac{G}{C} + \frac{R}{L}\right) \frac{\partial u \left(x,t\right)}{\partial t}
   +
   \left(\frac{RG}{LC}\right) u \left(x,t\right)
\end{equation}
which is a hyperbolic partial differential equation (PDE) \citep{Chapra2015}.
Letting
\begin{equation*}
   c^{2} = \frac{1}{LC}, \quad \alpha = \frac{G}{C}, \quad \beta = \frac{R}{L}
\end{equation*}
we can rewrite \Eqref{eqn: Telegraph eqn full} more succinctly as
\begin{equation}
   \label{eqn: Telegraph eqn short}
   c^{2} \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}}
   =
   \frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}}
   +
   \left(\alpha + \beta\right) \frac{\partial u \left(x,t\right)}{\partial t}
   +
   \alpha \beta u \left(x,t\right)
\end{equation}
The first-order term on the right-hand side of \Eqref{eqn: Telegraph eqn short} is the \textit{dissipation term},
while the zeroth-order term is the \textit{dispersion term}.
In the absence of losses, \textit{i.e.}, $R=G=0$, the telegraph equation reduces into one describing a wave
that propagates at a velocity $c$ and an angular frequency $\omega$ given as
\begin{equation}
   \label{eqn: Omega}
   \omega = \frac{1}{X \sqrt{LC}}
\end{equation}

Solving the telegraph equation is valuable in the analysis of power system dynamics.
However, deriving the expression for the exact analytic solution may not always be tractable nor the most efficient course;
in which case numerical approaches that approximate the PDE as a combination of algebraic operations are used.
This work presents a basic finite difference method for numerically solving the transmission line telegraph equation.
The remainder of the paper proceeds as follows.
Section \ref{sec: Finite Difference Approximation} details how the telegraph equation is approximated as a linear equation via discretization and finite differences.
Section \ref{sec: Illustrative Examples} presents and discusses results from select worked examples.
Section \ref{sec: Conclusion} concludes the work.

\section{Finite Difference Approximation}
\label{sec: Finite Difference Approximation}

\subsection{Discretization Scheme}
\label{subsec: Discretization Scheme}

We transform the continuous spatial domain into a set of equally separated discrete points, \textit{i.e.},
\begin{equation*}
   0 \leq x \leq X \quad
   \longrightarrow \quad
   x_{k} = k \Delta x,\ \ 0 \leq k \leq K \in \mathbb{Z}
\end{equation*}
In other words, we approximate the spatial domain by sampling $K+1$ points spaced $\Delta x$ apart.
Note that $x_{0}$ corresponds to $x=0$ just as $x_{K}$ to $x=X$.
Similarly, for the temporal domain:
\begin{equation*}
   0 \leq t \leq T \quad
   \longrightarrow \quad
   t_{n} = n \Delta t,\ \ 0 \leq n \leq N \in \mathbb{Z}
\end{equation*}
where $t_0$ corresponds to $t=0$ as $t_{N}$ to $t=T$.
The voltage defined on the continuous domain is likewise discretized,
and is parametrized by $k$ and $n$:
\begin{equation*}
   u \left(x,t\right) \quad
   \longrightarrow \quad
   u \left(x_{k},t_{n}\right)
\end{equation*}
For notational convenience, $u_{k}^{n} = u \left(x_{k},t_{n}\right)$.

\subsection{Difference Equation}
\label{subsec: Difference Equation}

We can approximate the continuous derivatives as central divided differences:
\begin{alignat*}{4}
   \frac{\partial u \left(x,t\right)}{\partial t}
   &\quad\longrightarrow\quad&
   \frac{\partial u_{k}^{n}}{\partial t}
   &=\ 
   \frac{u_{k}^{n+1} - u_{k}^{n-1}}{2 \Delta t} \\
   \frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}}
   &\quad\longrightarrow\quad&
   \frac{\partial{2} u_{k}^{n}}{\partial t^{2}}
   &=\ 
   \frac{u_{k}^{n+1} - 2 u_{k}^{n} + u_{k}^{n-1}}{\left(\Delta t\right)^{2}} \\
   \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}}
   &\quad\longrightarrow\quad&
   \frac{\partial^{2} u_{k}^{n}}{\partial x^{2}}
   &=\ 
   \frac{u_{k+1}^{n} - 2 u_{k}^{n} + u_{k-1}^{n}}{\left(\Delta x\right)^{2}}
\end{alignat*}
Substituting these into their continuous counterparts, we approximate the telegraph as a difference equation:
\begin{equation}
   \label{eqn: Difference eqn full}
   c^{2} \frac{u_{k+1}^{n} - 2 u_{k}^{n} + u_{k-1}^{n}}{\left(\Delta x\right)^{2}}
   =
   \frac{u_{k}^{n+1} - 2 u_{k}^{n} + u_{k}^{n-1}}{\left(\Delta t\right)^{2}}
   +
   \left(\alpha + \beta\right) \frac{u_{k}^{n+1} - u_{k}^{n-1}}{2 \Delta t}
   +
   \alpha \beta u_{k}^{n}
\end{equation}

This numerical approximation of \Eqref{eqn: Telegraph eqn short} suggests that we can estimate the voltage
at point $x_{k}$ at the next time instant $t_{n+1}$ given the voltages at $x_{k}$ and
at the neighbouring points at the current time instant
(that is, $u_{k}^{n}$, $u_{k-1}^{n}$, and $u_{k+1}^{n}$) and
the voltage at $x_{k}$ at the preceding time instant
(that is, $u_{k}^{n-1}$).

\subsection{Update Scheme}
\label{subsec: Update Scheme}

From \Eqref{eqn: Difference eqn full}, we can obtain the ``update'' $u_{k}^{n+1}$
given $u_{k}^{n}$, $u_{k-1}^{n}$, $u_{k+1}^{n}$, and $u_{k}^{n-1}$.
To express this more explicity, we can rewrite \Eqref{eqn: Difference eqn full} as
\begin{equation}
   \label{eqn: Difference eqn short}
   A u_{k}^{n+1} = E u_{k-1}^{n} + F u_{k}^{n} + E u_{k+1}^{n} - B u_{k}^{n-1}
\end{equation}
where
\begin{alignat}{2}
   \label{eqn: A}
   A &=\ 1 + \frac{\Delta \left(\alpha + \beta\right)}{2} \\
   \label{eqn: B}
   B &=\ 1 - \frac{\Delta \left(\alpha + \beta\right)}{2} \\
   \label{eqn: E}
   E &= \left(c \frac{\Delta t}{\Delta x}\right)^{2} \\
   \label{eqn: F}
   F &= 2 - 2 \left(c \frac{\Delta t}{\Delta x}\right)^{2} - \alpha \beta \left(\Delta t\right)^{2}
\end{alignat}

\subsection{Some Remarks}
\label{subsec: Some Remarks}

\subsubsection{Encoding Initial and Boundary Conditions}
\label{subsubsec: Encoding Initial and Boundary Conditions}

Notice that the difference equation approximation applies for $k=1,2,\ldots,K-1$ and $n=1,2,\ldots,N-1$.
It requires initial (\textit{i.e.}, at $t_{0}$) and boundary (\textit{i.e.}, at $x_{0}$ and $x_{K}$) values to be specified separately.

In general, initial voltage values are exressed as a function of $x$:
\begin{equation*}
   u \left(x,0\right) = \mu \left(x\right)
   \quad\longrightarrow\quad
   u_{k}^{0} = \mu \left(x_{k}\right),\ \forall k.
\end{equation*}
It is also common to have predetermined initial time rate of change of voltage,
which can then be approximated by a forward finite divided difference:
\begin{equation*}
   \frac{\partial u \left(x,0\right)}{\partial t} = \xi^{0}
   \quad\longrightarrow\quad
   \frac{\partial u_{k}^{0}}{\partial t} =
   \frac{u_{k}^{1} - u_{k}^{0}}{\Delta t} = \xi^{0}
   \quad\longrightarrow\quad
   u_{k}^{1} = u_{k}^{0} + \xi^{0} \Delta t,\ \forall k.
\end{equation*}

The sending- and receiving-end voltages can be expressed as functions of $t$:
\begin{alignat*}{4}
   & u \left(0,t\right) &&= \nu_{0} \left(t\right)
   &&\quad\longrightarrow\quad
   u_{0}^{n} &&= \nu_{0} \left(t_n\right),\ \forall n \\
   & u \left(X,t\right) &&= \nu_{X} \left(t\right)
   &&\quad\longrightarrow\quad
   u_{K}^{n} &&= \nu_{X} \left(t_n\right),\ \forall n.
\end{alignat*}
Information at the boundaries may also be expressed in terms of space-derivatives,
which can be approximated usign forward and backward finite divided differences:
\begin{alignat*}{8}
   & \frac{\partial u \left(0,t\right)}{\partial x} &&= \gamma_{0}
   &&\quad\longrightarrow\quad
   \frac{\partial u_{0}^{n}}{\partial x} &&= \frac{u_{1}^{n} - u_{0}^{n}}{\Delta x} = \gamma_{0}
   &&\quad\longrightarrow\quad
   u_{0}^{n} &&= u_{1}^{n} - \gamma_{0} \Delta x,\ \forall n \\
   & \frac{\partial u \left(X,t\right)}{\partial x} &&= \gamma_{X}
   &&\quad\longrightarrow\quad
   \frac{\partial u_{K}^{n}}{\partial x} &&= \frac{u_{K}^{n} - u_{K-1}^{n}}{\Delta x} = \gamma_{X}
   &&\quad\longrightarrow\quad
   u_{K}^{n} &&= u_{K-1}^{n} + \gamma_{X} \Delta x,\ \forall n.
\end{alignat*}

\subsubsection{Vectorizing the Update Scheme}
\label{subsubsec: Vectorizing the Update Scheme}

The update scheme \Eqref{eqn: Difference eqn short} is essentially a system of $K-1$ linear equations:
\begin{equation*}
   A
   \begin{bmatrix}
      u_{1}^{n+1} \\
      u_{2}^{n+1} \\
      \vdots \\
      u_{K-1}^{n+1}
   \end{bmatrix}
   =
   \begin{bmatrix}
      E & F & E & 0 & \cdots & 0 & 0 & 0 \\
      0 & E & F & E & \cdots & 0 & 0 & 0 \\
      \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
      0 & 0 & 0 & 0 & \cdots & E & F & E
   \end{bmatrix}
   \begin{bmatrix}
      u_{0}^{n} \\
      u_{1}^{n} \\
      u_{2}^{n} \\
      \vdots \\
      u_{K-1}^{n} \\
      u_{K}^{n}
   \end{bmatrix}
   - B
   \begin{bmatrix}
      0 & 1 & 0 & \cdots & 0 & 0 \\
      0 & 0 & 1 & \cdots & 0 & 0 \\
      \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
      0 & 0 & 0 & \cdots & 1 & 0
   \end{bmatrix}
   \begin{bmatrix}
      u_{0}^{n-1} \\
      u_{1}^{n-1} \\
      u_{2}^{n-1} \\
      \vdots \\
      u_{K-1}^{n-1} \\
      u_{K}^{n-1}
   \end{bmatrix}
\end{equation*}
Letting
\begin{alignat}{3}
   \tilde{\mathbf{u}}^{n} &=\
   \begin{bmatrix}
      u_{1}^{n}, u_{2}^{n}, \ldots, u_{K-1}^{n}
   \end{bmatrix}^{\intercal}\ \in \mathbb{R}^{K-1} \nonumber \\
   \mathbf{u}^{n} &=\
   \begin{bmatrix}
      u_{0}^{n}, u_{1}^{n}, \ldots, u_{K-1}^{n}, u_{K}^{n}
   \end{bmatrix}^{\intercal}\ \in \mathbb{R}^{K+1} \nonumber \\
   \label{eqn: E matrix}
   \mathbf{E} &=
   \begin{bmatrix}
      E & F & E & 0 & \cdots & 0 & 0 & 0 \\
      0 & E & F & E & \cdots & 0 & 0 & 0 \\
      \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
      0 & 0 & 0 & 0 & \cdots & E & F & E
   \end{bmatrix}\ \in \mathbb{R}^{\left(K-1\right) \times \left(K+1\right)} \\
   \label{eqn: B matrix}
   \mathbf{B} &=\
   B
   \begin{bmatrix}
      \mathbf{0} & \mathbf{I} & \mathbf{0}
   \end{bmatrix}\ \in \mathbb{R}^{\left(K-1\right) \times \left(K+1\right)}
\end{alignat}
where $\mathbf{0}$ is an $\left(K-1\right)$-vector of zeros and
$\mathbf{I}$ is the identity matrix of size $\left(K-1\right)$,
the update equations can be expressed compactly as
\begin{alignat}{2}
   \label{eqn: Difference eqn vectorized Pt. 1}
   A \tilde{\mathbf{u}}^{n+1} &=\ \mathbf{E} \mathbf{u}^{n} - \mathbf{B} \mathbf{u}^{n-1} \\
   \label{eqn: Difference eqn vectorized Pt. 2}
   \mathbf{u}^{n+1} &=\ 
   \begin{bmatrix}
      u_{0}^{n+1} \\
      \tilde{\mathbf{u}}^{n+1} \\
      u_{K}^{n+1}
   \end{bmatrix}
\end{alignat}
where $u_{0}^{n+1}$ and $u_{K}^{n+1}$ are determined from the boundary conditions.

\subsubsection{Courant-Friedrichs-Levy (CFL) Condition}
\label{subsubsec: Courant-Friedrichs-Levy (CFL) Condition}

In general, the finer the spatial and temporal domains are discretized, the better $u_{k}^{n}$ approximates $u \left(x,t\right)$.
However, as $\Delta x$ and $\Delta t$ get smaller, the number of gridpoints at which $u\left(x,t\right)$ is to be approximated increases,
and so does the computational burden.
Moreover, a judicious choice of the spatial and temporal steps helps avoid instability
(\textit{i.e.}, when the error drastically accumulates).
The CFL condition is a commonly used guide for selecting $\Delta x$ or $\Delta t$ (given the other):
\begin{equation}
   \label{eqn: CFL condition}
   \epsilon = c \frac{\Delta t}{\Delta x} \leq 1
\end{equation}
where $\epsilon$ is called the \textit{CFL number}.
In other words, for a particular $\Delta x$, the time step should be
\begin{equation*}
   \Delta t \leq \frac{\Delta x}{c}
\end{equation*}
to avoid an unstable approximation.
Intuitively, this upper limit on $\Delta t$ says that the simulation cannot be incremented any more than the time required
for a wave to travel one grid step in space.

\section{Illustrative Examples}
\label{sec: Illustrative Examples}

To supplement the discussion in the preceding sections, we present a few examples that focus on different aspects of
the finite difference numerical approximation of the transmission line telegraph equation.
We first investigate the benefits of using the vectorized update scheme
(\Eqref{eqn: E matrix}--\ref{eqn: Difference eqn vectorized Pt. 2})
in Section \ref{subsec: On Vectorization}.
Then, in Sections \ref{subsec: On Spatial Domain Discretization} and \ref{subsec: On Temporal Domain Discretization},
we look into the effects of varying how ``fine'' the spatial or the temporal domains are discretized.
%Lastly, we simulate faulted bus and faulted lines in Sections \ref{subsec: Simulating a Faulted Bus} and
%\ref{subsec: Simulating a Faulted Line}, respectively.
Lastly, we simulate a bus fault scenario in Section \ref{subsec: Simulating a Faulted Bus}.
All procedures are accomplished using NumPy \citep{NumPy2020} on an
Intel\textsuperscript{\textregistered} Core\textsuperscript{\texttrademark} i7-10750H with 16GB of RAM.
Source codes and related files are available in
{\tt illustrative examples/}
within the project repository.

\subsection{On Vectorization}
\label{subsec: On Vectorization}

Consider the scenario modelled as
\begin{equation}
   \label{eqn: Telegraph eqn for On Vectorization}
   \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}}
   =
   \frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}}
   +
   3.00001 \frac{\partial u \left(x,t\right)}{\partial t}
   +
   3 \times 10^{-5} u \left(x,t\right)
\end{equation}
where $0 \leq x \leq 1$, $0 \leq t \leq 1$,
and subject to
\begin{alignat}{3}
   \label{eqn: mu for On Vectorization}
   u \left(x,0\right) &=\ && \sin \left(5 \pi x\right) + 2 \sin \left(7 \pi x\right),\ \forall x \\
   \label{eqn: xi for On Vectorization}
   \frac{\partial u \left(x,0\right)}{\partial t} &=\ && 0,\ \forall x \\
   \label{eqn: nu_0 for On Vectorization}
   u \left(0,t\right) &=\ && 0,\ \forall t \\
   \label{eqn: nu_X for On Vectorization}
   u \left(1,t\right) &=\ && 0,\ \forall t
\end{alignat}
Ten $x$-domain discretizations are examined;
we start from $K+1=100$ grid points and increment by 50.
To avoid numerical instability, we use $\epsilon = 0.10$ to set the corresponding time steps.
For each domain discretization, the telegraph equation is solved using both
the basic (\Eqref{eqn: Difference eqn short}--\ref{eqn: F})
and the vectorized (\Eqref{eqn: E matrix}--\ref{eqn: Difference eqn vectorized Pt. 2}) update schemes.
Twenty independent runs are performed, where execution times are the metrics of interest.
All procedures are documented in
{\tt illustrative examples/on vectorization/on vectorization.ipynb}.

\begin{figure}[t!]
   \centering
   \includegraphics[scale=0.7]{on vectorization.png}
   \caption{Execution times for basic and vectorized update schemes at different domain discretizations.}
	\label{fig: On Vectorization results}
\end{figure}

Runtimes for the different discretizations and update schemes are summarized in \Figref{fig: On Vectorization results}.
The vectorized update scheme is significantly faster than the naive loop approach,
and the speedup gained increases as the domain discretization gets finer.
The disparity in the execution times is attributed to NumPy's support for array-based computing;
in fact, one can expect to arrive at similar findings when using other numerical computing tools such as
MATLAB\textsuperscript{\textregistered} \citep{MATLABR2020b}.
This suggests that, for the same computational resources and time duration,
vectorization enables one to iterate and model more scenarios
than using the explicit loop-based approach.
Henceforth, we will use the vectorized update scheme.

\subsection{On Spatial Domain Discretization}
\label{subsec: On Spatial Domain Discretization}

Consider the scenario
\begin{alignat}{3}
   \label{eqn: Telegraph eqn for On Spatial Domain Discretization}
   \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}} &=\ && \frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}} + 3 \frac{\partial u \left(x,t\right)}{\partial t} \\
   \label{eqn: mu for On Spatial Domain Discretization}
   u \left(x,0\right) &=\ && \sin \left(5 \pi x\right) + 2 \sin \left(7 \pi x\right),\ \forall x \\
   \label{eqn: xi for On Spatial Domain Discretization}
   \frac{\partial u \left(x,0\right)}{\partial t} &=\ && 0,\ \forall x \\
   \label{eqn: nu_0 for On Spatial Domain Discretization}
   u \left(0,t\right) &=\ && 0,\ \forall t \\
   \label{eqn: nu_X for On Spatial Domain Discretization}
   u \left(1,t\right) &=\ && 0,\ \forall t
\end{alignat}
where $0 \leq x \leq 1$ and $0 \leq t \leq 1$.
From \cite{Zhang2019}, the corresponding analytic solution is
\begin{alignat}{3}
   \label{eqn: Analytic soln for On Spatial Domain Discretization Part 1}
   u \left(x,t\right) &=\ && e^{-1.5t}
   \sin\left(5 \pi x\right) \left[\cos\left(\theta\left(5\right)t\right) + \frac{1.5}{\theta\left(5\right)} \sin \left(\theta\left(5\right)t\right)\right] \nonumber \\
   & &&+ 2e^{-1.5t} \sin\left(7 \pi x\right) \left[\cos\left(\theta\left(7\right)t\right) + \frac{1.5}{\theta\left(7\right)} \sin \left(\theta\left(7\right)t\right)\right] \\
   \label{eqn: Analytic soln for On Spatial Domain Discretization Part 2}
   \theta \left(z\right) &=\ && \sqrt{\left(z \pi\right)^{2} - 1.5^{2}}
\end{alignat}
To illustrate how discretization of the spatial domain affects the quality of the numerical approximation,
we consider varying $K$ (and hence, the number of grid points)
while fixing $N=999$ (that is, 1000 grid points along the temporal domain).
We compare the results obtained using the finite difference approximations
and those obtained by the analyic solution
(\Eqref{eqn: Analytic soln for On Spatial Domain Discretization Part 1}--\ref{eqn: Analytic soln for On Spatial Domain Discretization Part 2});
in particular, we look at the mean squared error (MSE) calculated as
\begin{equation*}
   \sqrt{\sum_{n=0}^{N} \sum_{k=0}^{K} \left\{u \left(x_{k}, t_{n}\right) - \hat{u} \left(x_{k}, t_{n}\right)\right\}^{2} }
\end{equation*}
where $u \left(x_{k}, t_{n}\right)$ and $\hat{u} \left(x_{k}, t_{n}\right)$ are
the respective analytic and approximate instantaneous voltages evaluated at the grid points along $x-$ and $t-$domains.
All procedures are documented in
{\tt illustrative examples/on spatial domain discretization/on spatial domain discretization.ipynb}.

\begin{figure}[t!]
   \centering
   \includegraphics[scale=0.7]{on spatial domain discretization - mse and cfl.png}
   \caption{MSE and $\epsilon$ for various discretization schemes of the $x$-domain.}
	\label{fig: On Spatial Domain Discretization results: MSE and CFL}
\end{figure}

\begin{figure}[t!]
   \centering
   \includegraphics[scale=0.7]{on spatial domain discretization - finite time blowup.png}
   \caption{Numerical instability when the $x-$domain is approximated by 1025 grid points and the $t-$domain by 1000 grid points.}
	\label{fig: On Spatial Domain Discretization results: Finite time blowup}
\end{figure}

\begin{figure}[t!]
   \centering
   \includegraphics[scale=0.7]{on spatial domain discretization - analytic vs numerical.png}
   \caption{Comparing analytic solution with finite difference approximations at various $K$.}
	\label{fig: On Spatial Domain Discretization results: Analytic vs numerical solutions}
\end{figure}

Referring to \Figref{fig: On Spatial Domain Discretization results: MSE and CFL},
we find a general trend of the approximation accuracy improving
as more grid points are sampled from the spatial domain.
One can explain this by considering the limit definition of the partial derivative with respect to $x$:
\begin{equation}
   \label{eqn: Limit definition of the partial derivative wrt to x}
   \frac{\partial u \left(x,t\right)}{\partial x}
   =
   \lim_{\Delta x \longrightarrow 0}
   \frac{u \left(x + \Delta x, t\right)}{\Delta x}
\end{equation}
Notice also the diminishing returns in accuracy with increasing $K$:
there is a significant drop in MSE as one increases the number of grid points from 50 to 150,
but the improvement tapers off and is insignificant from about 450 to 990 grid points.
From a practical viewpoint, one can specify a ``good enough'' accuracy level so that the number of grid points--and hence,
the computational expenses--can be kept to a manageable value.
Another salient observation from \Figref{fig: On Spatial Domain Discretization results: MSE and CFL}
is that the MSE increases without bounds when 1025 grid points are sampled from the $x-$domain.
This ``finite time blowup'' (see \Figref{fig: On Spatial Domain Discretization results: Finite time blowup})
is expected to occur since
\begin{equation*}
   \epsilon = c \frac{\Delta t}{\Delta x} = \frac{\frac{1 - 0}{1000 - 1}}{\frac{1 - 0}{1025 - 1}} = \frac{1024}{999} > 1
\end{equation*}
Nevertheless, the finite difference approximation performs well when the number of spatial grid points yields a numerically stable solution
(see \Figref{fig: On Spatial Domain Discretization results: Analytic vs numerical solutions}).

\subsection{On Temporal Domain Discretization}
\label{subsec: On Temporal Domain Discretization}

Consider again the scenario in
\Eqref{eqn: Telegraph eqn for On Spatial Domain Discretization}--\ref{eqn: nu_X for On Spatial Domain Discretization}.
This time we illustrate how the discretization of the temporal domain affects the quality of the numerical approximation.
Specifically, we vary $N$ while fixing $K=99$ (that is, 100 grid points along the spatial domain).
We compare the results obtained using the finite difference approximations and those obtained by the analytic solution
(\Eqref{eqn: Analytic soln for On Spatial Domain Discretization Part 1}--\ref{eqn: Analytic soln for On Spatial Domain Discretization Part 2})
in terms of the MSE.
All procedures are documented in
{\tt illustrative examples/on temporal domain discretization/on temporal domain discretization.ipynb}.

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{on temporal domain discretization - mse and cfl.png}
	\caption{MSE and $\epsilon$ for various discretization schemes of the $t$-domain.}
	\label{fig: On Temporal Domain Discretization results: MSE and CFL}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{on temporal domain discretization - finite time blowup.png}
	\caption{Numerical instability when the $x-$domain is approximated by 100 grid points and the $t-$domain by 95 grid points.}
	\label{fig: On Temporal Domain Discretization results: Finite time blowup}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{on temporal domain discretization - analytic vs numerical.png}
	\caption{Comparing analytic solution with finite difference approximations at various $N$.}
	\label{fig: On Temporal Domain Discretization results: Analytic vs numerical solutions}
\end{figure}

We see in \Figref{fig: On Temporal Domain Discretization results: MSE and CFL} that,
for a fixed discretization of the spatial domain,
the finite difference approximation generally improves
(\textit{i.e.}, the MSE decreases)
as more grid points are sampled from the temporal domain.
This can be explained by considering the limit definition of the partial time derivative:
\begin{equation}
	\label{eqn: Limit definition of the partial derivative wrt to t}
	\frac{\partial u \left(x,t\right)}{\partial t}
	=
	\lim_{\Delta t \longrightarrow 0}
	\frac{u \left(x, t + \Delta t\right)}{\Delta t}
\end{equation}
However, we note that even if we decrease $\Delta t$ to zero
(\textit{i.e.}, the continuous limit)
the interpolation error is still nonzero;
hence, the MSE never quite reaching null.
Similar to the case with $x-$domain discretization,
the accuracy gain tapers off with increasing $N$.
From a practical perspective, one can attain approximately the same MSE when using 750 and 1000 grid points,
with the former being computationally lighter.
We also observe a sudden increase in MSE when 95 grid points are sampled from the $t-$domain.
This error accumulation,
better illustrated in \Figref{fig: On Temporal Domain Discretization results: Finite time blowup},
is an expected outcome since
\begin{equation*}
	\epsilon = c \frac{\Delta t}{\Delta x}
	= \frac{\frac{1 - 0}{95 - 0}}{\frac{1 - 0}{100-0}}
	= \frac{100}{95}
	> 1
\end{equation*}
When the number of temporal grid points yields $\epsilon <= 1$,
the finite difference approximation performs well (see \Figref{fig: On Temporal Domain Discretization results: Analytic vs numerical solutions}).

\subsection{Simulating a Faulted Bus}
\label{subsec: Simulating a Faulted Bus}

Consider a 150-km transmission line with
$R = 1.5 \times 10^{-2}$ $\Omega$/m,
$L = 1 \times 10^{-3}$ H/m,
$G = 1 \times 10^{-10}$ S/m, and
$C = 1 \times 10^{-8}$ F/m.
From \Eqref{eqn: Telegraph eqn short},
\begin{equation*}
	c^{2} = 1 \times 10^{11}, \quad
	\alpha = 1 \times 10^{-2}, \quad
	\beta = 1.5
\end{equation*}
which gives the telegraph equation as
\begin{equation}
	\label{eqn: Telegraph eqn for Simulating a Faulted Bus}
	1 \times 10^{11} \frac{\partial^{2} u \left(x,t\right)}{\partial x^{2}}
	=
	\frac{\partial^{2} u \left(x,t\right)}{\partial t^{2}}
	+
	1.51 \frac{\partial u \left(x,t\right)}{\partial t}
	+
	0.015 u \left(x,t\right) 
\end{equation}
For all $x$ at $t = 0$,
\begin{alignat}{3}
	\label{eqn: mu for Simulating a Faulted Bus}
	u \left(x,0\right) &=\ && \sin \left(120 \pi x\right) \\
	\label{eqn: xi for Simulating a Faulted Bus}
	\frac{\partial u \left(x,0\right)}{\partial t} &=\ && 0
\end{alignat}
The sending-end bus injects power at a time-varying voltage described by $1.5 \sin \left(120 \pi t\right)$.
We assume a simulation period of $0 \leq t \leq 1.5$ s.
At $t = 0.3$ s, the receiving-end bus is faulted, \textit{i.e.},
it is shorted to ground.
As part of the network protection scheme,
the receiving-end bus is disconnected from the line at $t = 0.7$ s
and the sending-end bus 0.10 s later.
The fault is said to be ``cleared'' after both buses are disconnected from the line.
We describe the sending-end boundary condition as
\begin{alignat}{4}
	\label{eqn: Dirichlet condition for sending-end bus for Simulating a Faulted Bus}
	u \left(0,t\right) &=\ && 1.5 \sin \left(120 \pi t\right),\ \quad &0 \leq t < 0.8 \\
	\label{eqn: Neumann condition for sending-end bus for Simulating a Faulted Bus}
	\frac{\partial u \left(0,t\right)}{\partial x} &=\ && 0,\ \quad &0.8 \leq t \leq 1.5
\end{alignat}
The boundary condition for the receiving-end bus is described as
\begin{alignat}{4}
	\label{eqn: Dirichlet condition for receiving-end bus for Simulating a Faulted Bus}
	u \left(150 \times 10^{3}, t\right) &=\ &&0,\ \quad &0.3 \leq t \leq 0.7 \\
	\label{eqn: Neumann condiiton for receiving-end bus for Simulating a Faulted Bus}
	\frac{\partial u \left(150 \times 10^{3},t\right)}{\partial x} &=\ && 0,\ \quad &0 \leq t < 0.3,\ 0.8 < t \leq 1.5
\end{alignat}
We discretize the $x-$domain by sampling 100 grid points from it (that is, $K = 99$).
The $t-$domain is discretized in such a way that the time step $\Delta t$ is a tenth of the maximum time step derived from \Eqref{eqn: CFL condition}:
\begin{equation*}
	\Delta t = 0.10 \frac{\Delta x}{c} \longrightarrow N = \frac{T}{\Delta t}
\end{equation*}
Plugging in known values, the number of grid points along the $t-$domain is $N + 1$, which--rounded to the nearest thousand--is 15000.
The discretization corresponds to $\epsilon = 0.105205$.
All procedures for this simulation are documented in
{\tt illustrative examples/simulating a faulted bus/simulating a faulted bus.ipynb}.

The simulation results are summarized in Figures
\ref{fig: Simulating a Faulted Bus results: Sending- and receiving-end voltages}--\ref{fig: Simulating a Faulted Bus results: Voltage profiles at post-fault}.
The time-evolution of the sending- and receiving-end voltages are shown in \Figref{fig: Simulating a Faulted Bus results: Sending- and receiving-end voltages}.
The voltage profiles across the length of the line at different time instants before, during, and after the fault event
are shown in Figures
\ref{fig: Simulating a Faulted Bus results: Voltage profiles at pre-fault},
\ref{fig: Simulating a Faulted Bus results: Voltage profiles at during fault}, and
\ref{fig: Simulating a Faulted Bus results: Voltage profiles at post-fault},
respectively.
We emphasize two main observations.
First, in the absence of fault, there tends to be a decreasing trend in voltage from $x=0$ to $x=150 \times 10^{3}$
(see \Figref{fig: Simulating a Faulted Bus results: Voltage profiles at pre-fault}).
This is a direct result of the presence of nonzero $R$ and $G$, which means that the line is not lossless.
Second, when the fault is cleared, the voltage along the line acts like a wave that bounces back and forth between the two ends,
where in each reflection the amplitude decreases, and will eventually reach zero.
The gradual decrease in voltage is also attributable to the losses of the line.

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.7]{simulating a faulted bus - sending- and receiving-end voltages.png}
	\caption{Variation of the sending- and receiving-end voltages through time.}
	\label{fig: Simulating a Faulted Bus results: Sending- and receiving-end voltages}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{simulating a faulted bus - voltage profiles at pre-fault.png}
	\caption{Voltage profile across the length of the line before the fault event.}
	\label{fig: Simulating a Faulted Bus results: Voltage profiles at pre-fault}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{simulating a faulted bus - voltage profiles during fault.png}
	\caption{Voltage profile across the length of the line during the fault event.}
	\label{fig: Simulating a Faulted Bus results: Voltage profiles at during fault}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.7]{simulating a faulted bus - voltage profiles at post-fault.png}
	\caption{Voltage profile across the length of the line after the fault is cleared.}
	\label{fig: Simulating a Faulted Bus results: Voltage profiles at post-fault}
\end{figure}

\section{Conclusion}
\label{sec: Conclusion}

Finite difference approximations of differential equations gained patronage due to its simplicity.
We have shown how this approach enables us to approximate the solution to the transmission line telegraph equation
by iteratively solving systems of linear algebraic equations.
We also demonstrated through various experiments that
(i) vectorizing the update scheme offers significantly less execution times than that of the explicit loop-based scheme;
(ii) the accuracy of the approximation generally increases as the discretization of the spatial and of the temporal domains get finer;
and
(iii) the CFL condition can be used as a guide for choosing domain discretizations that yield numerical stability.

\bibliography{references}
\bibliographystyle{iclr2021_conference}

\end{document}