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 = 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 feedback. 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.

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 code now reads,

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