# 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 imporant 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 if formalized by the so-called method of lines. Consider a dummy variable $$\phi (t, x_1, x_2, ..., x_n)$$ that statisfies 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 presrcribed 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 in our numerical solver. It only exists on a finite number of points in time. This also means that the proper solution $$\phi$$ maybe truly 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 assymptotic solution with increasing effort. This assymptotic solution should match the solution to the original equation.

## 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.

$\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 circlular 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 imporant 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 = x1 - dt*x2;
x2 = x2 + dt*tmp_x1;
}
...

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 feadback. As such, we print the solution:

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

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. The result stored as svg 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 = x1 - dt*x2; x2 = x2 + dt*tmp_x1; } // Errors logging double error1 = x1 - 1, error2 = x2 - 0; double L2 = sqrt(error1*error1 + error2*error2); printf ("%d %g\n", Nstps, L2); } } 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 logarithmic scale, and reveal the dependency of the $$L_2$$ error on the number of steps.

\$ gnuplot
...
> set logscale xy
...
> 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 can state that the forward Euler method is first-order accurate, which corresponds with the literature.