# 2. The method of lines

The Navier-Stokes equations are an example of a so-called partial differential equation: The solution exists in a multi-dimensional parameter space, and as such the equations contain various corresponding partial derivatives. In the previous chapter we can indeed see derivatives with respect to $x, y, z$ and $t$. For our solver we choose to delineate between the spatial dimensions and the time dimension for two important reasons:

The time dimension is the only one associated with causality.

The equations are isotropic for the group $x,y$ and $z$ (i.e. they maybe rotated)

The actual separation in the treatment of the variables is formalized
by the so-called *method of lines*. Consider a dummy variable
$\phi (t, x_1, x_2, ..., x_n)$
that satisfies an equation of the form,

$\frac{\partial \phi}{\partial t} = \mathcal{L}(\phi),$

where
$\mathcal{L}$
is some differential operator that does *not* include
differentiation with respect to
$t$.
Because time and physical causality are related concepts, it makes sense
to solve for
$\phi$
as an initial value problem. Here
$\phi (t = 0, x_1, x_2, ..., x_n)$
is prescribed as an initial condition, and the equation can be solved by
using a time-integration method which advances the solution at the
current time
($\phi_n$)
to the next point in time
($\phi_{n+1}$)
which differ by a timestep size `dt`

. It is imporant to note
that the continuously evolving field
$\phi$
is now *discretized* as it only exists on a finite number of
points in time. This also means that the proper solution
$\phi$
is fundamentally different from the discrete solution
$\phi_n$,
and this introduces the concepts of convergence, consistency and
stability.

### Consistency

The timestep paramter is artificially introduced. It is imporant that for increasingly smaller timesteps, the discretization error (due to truncation) of the numerical schemes vanishes.

### Stability

The inevitable errors introduced by discretization should not grow in time.

### Convergence

By reducing
$\mathrm{dt} \rightarrow 0$,
the number or timesteps (and effort) required to obtained a solution at
some given time goes to infinity. As such, we hope to see the solution
*converges* towards an asymptotic solution with increasing
effort. This asymptotic solution should match the solution to the
original equation at the discrete intervals.

## Practice with an example

For many non-linear problems (like flow problems) the above is hard to formalize in numerical algorithms. For example, it is not even known if the Navier-Stokes equations have smooth solutions for the general case. So we will have to think of test cases to gain confidence in our digitally-generated data.

For now, we start with an example problem, namely,

$\frac{\mathrm{d}x_1}{\mathrm{d} t} = -x_2,$ $\frac{\mathrm{d}x_2}{\mathrm{d} t} = x_1.$

with initial conditions $x_1(t = 0) = 1$ and $x_2(t = 0) = 0$. It can be verified that the solution reads,

$x_1 = \mathrm{cos}(t),$ $x_2 = \mathrm{sin}(t).$

Which describes a circular trajectory. Lets start coding a time
integrator for this problem in a file called `circle.c`

. We
will need a few functions defined in some system header files. We begin
with,

```
#include <stdio.h>
#include <math.h>
int main() {
}
```

Whoa! those are a lot of lines of code, with ample of room for typos. It is important to check your syntax early on, as spotting a bug becomes harder as we proceed.

```
$ gcc circle.c
$ ./a.out
```

Lets add some variables, including the final time $t_{end} = 2 \pi$, when $x_1$ and $x_2$ have returned to their original positions.

```
...
int main() {
double x1 = 1, x2 = 0;
double Nstps = 50, tend = 6.2832; // ~2*pi
double dt = tend/Nstps;
}
```

The `...`

dots are added to hint at the ommited code (the
`#include`

directives in this case). Any errors? It is time
for an important design choice, *How do we advance the solution in
time?*. For the sake of minimalism, we choose the forward Euler
method:

$\phi_{n+1} = \phi_{n} + \mathrm{dt} \mathcal{L}(\phi_n).$

If we wrap that in a time loop,

```
...
double dt = tend/Nstps;
for (int i = 0; i < Nstps; i++) {
double tmp_x1 = x1;
= x1 - dt*x2;
x1 = x2 + dt*tmp_x1;
x2 }
...
```

Note that we needed to store the value for
$x_{1,n}$
because it would have been at stage `n+1`

at the time of
updating
$x_2$,
which would not correspond to the Forward Euler method. Compiling and
running this code does not give us any feedback. As such, we print the
solution:

```
...
= x2 + dt*tmp_x1;
x2 ("%g %g\n", x1, x2);
printf }
...
```

Running this code produces

```
$ ./a.out
1 0.125664
...
...
1.4787 -0.0484434
```

We can store the output in a file for analysis

`$ ./a.out > solution`

Now we can plot the data using `Gnuplot`

.

```
$ gnuplot
...
gnuplot> plot 'solution'
```

Which should show you a bunch of markers. We can add the analytical solution

```
gnuplot> plot sample [t=0:2*pi] '+' using (cos(t)) : (sin(t)) t 'circle' w l lw 2,\
'solution'
```

It appears that the numerical solution is moving away for the analytical solution. For a more presentable plot, use the following gnuplot code

```
set size square
set xlabel 'x_1'
set ylabel 'x_2'
set grid
plot sample [t=0:2*pi] '+' using (cos(t)) : (sin(t)) t 'circle' w l lw 2,\
'solution'
q
```

where the last command closes Gnuplot.

Finally, we check the convergence of this method for increasing
`Nstps`

, by wrapping the time loop in a new loop. Also at
$t = t_{end}$,
the distance to the analytical solution is computed and printed. The
total code note is:

```
#include <stdio.h>
#include <math.h>
int main() {
double tend = 6.2832; // ~2*pi
int Nstps_max = 4000;
// Convergence loop
for (int Nstps = 50; Nstps < Nstps_max; Nstps *= 2) {
double dt = tend/Nstps;
double x1 = 1, x2 = 0;
// Time loop
for (double t = 0; t <= tend; t += dt) {
double tmp_x1 = x1;
= x1 - dt*x2;
x1 = x2 + dt*tmp_x1;
x2 }
// Errors logging
double error1 = x1 - 1, error2 = x2 - 0;
double L2 = sqrt(error1*error1 + error2*error2);
("%d %g\n", Nstps, L2);
printf }
}
```

Lets give it a go:

```
$ gcc circle.c
Some error message referencing `sqrt`
```

The compiler (or rather the linker) may trow an error, because it
does not know how to compute the square root (`sqrt()`

). We
must link the math library during compilation, like so;

```
$ gcc circle.c -lm
$ ./a.out > errors
```

Now we can plot the data on a doubly-logarithmic scale, and reveal the dependency of the $L_2$ error on the number of steps.

```
$ gnuplot
...
gnuplot> set logscale xy
gnuplot> plot 'errors', 20*x**-1 lw 2 t 'x^{-1}'
```

Which may produce something like:

Because the so-called *scaling* of the error with
$L_2 \propto N^{-1}$,
fits the data well, we can conclude that the solution indeed converges.
More specifically, the forward Euler method is *first-order
accurate*, which corresponds with the literature.

Continue with Navier-Stokes related programming in chapter 3.