\documentclass [11pt] {article}
\usepackage{verbatim}
\usepackage{epsfig}
%
%
\topmargin=-0.75in
\hoffset=-0.75in
\textheight=9.0in \vsize\textheight
\textwidth=6.5in \hsize\textwidth
\def\lb{\left(}
\def\rb{\right)}
\def\lbb{\left[}
\def\rbb{\right]}
%
\begin {document}
\newpage
\setcounter{page}{1}
\setcounter{equation}{0}
%
\begin{center}
{\bf PHYS 410/555 Computational Physics}\\
{\em The Method of Lines for the Wave Equation}
\end  {center}
%
One approach to the numerical solution of time-dependent {\em partial 
differential equations} (PDEs) is to use a discretization technique, such 
as finite-differencing, but only apply it explicitly to the 
{\em spatial} part(s) of the PDE operator(s) under consideration.  
Following the spatial discretization, one is left with a set of 
coupled {\em ordinary differential equations} in $t$, which can 
then often be solved by a ``standard'' ODE integrator such as 
{\tt LSODA}.

As an example of this technique, consider the {\em wave equation} in 
one space dimension (often called the ``1D wave equation):
\begin{equation}
\label{waveA}
	\frac{\partial^2}{\partial t^2} u(x,t) = 
   c^2 \frac{\partial^2}{\partial x^2} u(x,t)
\end{equation}
Introducing the notation that a subscript denotes partial differentiation,
and suppressing the explicit $x$ and $t$ dependence, (\ref{waveA}) can also 
be written as
\begin{equation}
\label{waveb}
	u_{tt} = c^2 u_{xx}
\end{equation}
As you probably know, the wave equation describes propagation of 
disturbances, or waves, at a speed $c$: waves can either travel 
to the right (velocity $+c$), or to the left (velocity $-c$).
Without loss of generality, we can always choose units such that 
$c=1$, and, for convenience, we will do so.  Our wave equation then
becomes:
\begin{equation}
\label{waveC}
	u_{tt} = u_{xx}
\end{equation}

As with any differential equation, boundary conditions play a crucial 
role in fixing a solution of~(\ref{waveC}).  Here, we
will solve the wave equation on the domain
\begin{equation}
\label{domain}
	0 \le x \le 1 \quad\quad t \ge 0
\end{equation}
and will thus have to provide boundary conditions at $x=0$ and $x=1$, as 
well as initial conditions at $t=0$.

For concreteness, we will prescribe {\em Dirichlet boundary conditions}:
\begin{equation}
\label{bcs}
	u(0,t) = u(1,t) = 0
\end{equation}
as well as the following {\em initial conditions}:
\begin{eqnarray}
\label{ics}
	u(x,0) & = & u_0(x) = 
      \exp\left(-\left(\frac{x-x_0}{\Delta}\right)^2\right)\\
\label{icsv}
	u_t(x,0) & = & 0 
\end{eqnarray}
where $x_0$ ($0 < x_0 < 1$) and $\Delta$ are specified constants.

If we think in terms of small-amplitude waves propagating on a string, 
then the Dirichlet conditions correspond to keeping the ends of the 
string fixed.  The interpretation of the initial conditions is as 
follows: In solving~(\ref{waveC}) we have the freedom to specify 
the amplitude of the disturbance for all values of $x$, as well as 
the time-rate of change of that amplitude, again for all values of 
$x$.  

We thus have two functions worth of freedom in specifying 
our initial conditions.  We set the initial amplitude to some 
functional form given by $u_0(x)$; here we use a ``gaussian 
pulse'' that is centred at $x_0$, and that has an overall 
effective width of a few $\times \,\, \Delta$. 
We also set the initial time rate of change of the amplitude to be 
0 for all $x$.  

Such data is known as {\em time symmetric}, since it 
defines an instant in the evolution of the wave equation where there 
is a $t \to -t$ symmetry.  In other words, with time symmmetric 
initial data, if we integrate {\em backward} in time, we will see 
exactly the same solution as a function of $-t$ as we see integrating 
forward in time.

Since the wave equation describes propagating disturbances, and given 
that the initial conditions {\em are} time symmetric, 
a little reflection might convince you that the 
initial conditions~(\ref{ics}) and~(\ref{icsv}) must represent
a superposition of equal amplitude right-moving and left-moving pulses.
Thus, we should expect the solution of~(\ref{waveC}), subject to~(\ref{bcs}),
(\ref{ics}) and~(\ref{icsv}) to describe the propagation of two equal-amplitude 
pulses that are initially coincident, but that subsequently move
apart reflect off $x=0$ and $x=1$ respectively, move together, through
each other, then apart, etc. etc.  Indeed, this is precisely the 
behaviour we will observe in our subsequent numerical solution.

As mentioned above, the {\em method of lines}, involves an explicit 
discretization {\em only of the spatial part of the PDE operator}.
Here we will use the familiar $O(h^2)$ finite-difference approach
to the treatment of $u_{xx} \equiv \partial^2 u / \partial x^2$.

Howver, before proceeding to the spatial discretization, we first 
note that~(\ref{waveC}) is a second-order-in-time equation.
In order that our approach eventually produce a set of {\em first order}
ODEs in $t$, we introduce an auxiliary variable, $v(x,t)$,
\begin{equation}
	v(x,t) \equiv u_t(x,t) \equiv \frac{\partial u}{\partial t}(x,t)
\end{equation}
and then rewrite~(\ref{waveC}) as the {\em system:}
\begin{eqnarray}
\label{ut}
	u_t & = & v \\
\label{vt}
   v_t & = & u_{xx}
\end{eqnarray}
The boundary conditions become 
\begin{equation}
\label{fobcs}
	u(0,t) = u(1,t) = v(0,t) = v(1,t) = 0
\end{equation}
while the initial conditions are now
\begin{eqnarray}
\label{foics}
	u(x,0) & = & u_0(x) = 
      \exp\left(-\left(\frac{x-x_0}{\Delta}\right)^2\right)\\
\label{foicsv}
	v(x,0) & = & 0 
\end{eqnarray}

We can now proceed with the spatial discretization.
To that end, we replace the continuum spatial domain $0 \le x \le 1$
by a uniform finite difference mesh, $x_j$:
\begin{equation}
 x_j \equiv (j - 1) h \quad\quad j = 1, 2, \cdots N \quad\quad\quad
    h \equiv (N - 1)^{-1}
\end{equation}
and introduce the discrete unknowns, $u_j$ and $v_j$:
\begin{eqnarray}
	u_j & \equiv & u_j(t) \equiv u(x_j,t) \\
	v_j & \equiv & v_j(t) \equiv v(x_j,t) 
\end{eqnarray}
Using the usual centred, $O(h^2)$ approximation for the second spatial
derivative,
\begin{equation}
	u_{xx}(x_j) = \frac{u_{j+1} - 2 u_j + u_{j-1}}{h^2} + O(h^2)
\end{equation}
eqs.~(\ref{ut}) and (\ref{vt}) become a set of $2(N-2)$ coupled ODEs for the 
$2(N-2)$ unknowns $u_j(t)$ and $v_j(t)$,~~$j = 2,\cdots N-1$:
\begin{eqnarray}
\label{interior}
	\frac{du_j}{dt} & = & v_j
   \quad\quad\quad\quad\quad\quad\quad\quad\quad j = 2,\cdots N-1 \\
	\frac{dv_j}{dt} & = & \frac{u_{j+1} - 2 u_j + u_{j-1}}{h^2} \quad\quad
   j = 2,\cdots N-1
\end{eqnarray}
We can implement the Dirichlet boundary conditions as follows: if 
the boundary conditions are satisfied at the initial time, $t=0$, 
then they will be satisfied at all future times provided that 
the time derivatives of $u$ and $v$ vanish at the boundaries.  Using
this observation, we
can now write down a complete set of $2 N$ coupled ODEs in the $2 N$
unknowns $u_j(t)$ and $v_j(t)$ which can then be solved using {\tt LSODA}:
\begin{eqnarray}
	\frac{du_1}{dt} &=& 0 \\
	\frac{du_j}{dt} &=& v_j \quad\quad\quad\quad\quad\quad\quad\quad\quad 
   j = 2,\cdots N-1 \\
	\frac{du_N}{dt} &=& 0 \\[0.1in]
	\frac{dv_1}{dt} &=& 0 \\
	\frac{dv_j}{dt} &=& \frac{u_{j+1} - 2 u_j + u_{j-1}}{h^2} \quad\quad
   j = 2,\cdots N-1 \\
	\frac{dv_N}{dt} &=& 0 
\end{eqnarray}

This solution procedure is implemented by the the program 
{\tt wave} (See {\tt $\sim$phys410/ode/wave}).  You will follow 
an analagous approach to solve the {\em diffusion equation} in 
the final homework.
\end {document}
