Partial Differential Equations/Method of Lines

Method of Lines

Introduction

The method of lines is an interesting numeric method for solving partial differential equations. The idea is to semi-discretize the PDE into a huge system of (continuous and interdependent) ODEs, which may in turn be solved using the methods you know and love, such as forward stepping Runge-Kutta. This is where the name "method of lines" comes from: the solution is composed of a juxtaposition of lines (curves, more properly).

The method of lines is applicable to IBVPs only. The variables (or variable) that have boundary constraints are discretized, and the (one) variable that is associated with the initial value(s) becomes the independent variable for the ODE solution. Put more formally, the method of lines is obtained when all of the variables are discretized except for one of them, resulting in a coupled system of ODEs.

Pure BVPs, such as Laplace's equation on some terrible boundary, can be solved in a manner similar to Jacobi iteration known as the method of false transients.

Example Problem

Consider the following nonlinear IBVP, where $h = h(z, t)$:

$\frac{\partial h}{\partial t} = 2 \left(\frac{\partial h}{\partial z}\right)^2 + h\frac{\partial^2 h}{\partial z^2}\,$

$h(z, 0) = 0\,$

$h(0, t) = 1 \qquad \mbox{and} \qquad h(1, t) = 0\,$

This intimidating nondimensional system describes the flow of liquid in an interior corner ("groove"), where $h$ is the height of the liquid and $z$ is the distance along the axis of the corner. The IBVP states that the corner is initially empty, and that at one end ($z = 0$) the height is fixed to 1. Naturally, fluid will flow along $z$, increasing $h$ until the other boundary ($z = 1$), where $h$ is always zero, is reached. It's worth noting that though this IBVP can't be solved analytically, similarity solutions to the nonlinear PDE exist for some other situations, such as fixed height at z = 0 and no constraint to the right of that.

The boundary values are associated with the spatial variable, $z$. Suppose that $h$ is discretized in space but not time to $n$ points on a uniform grid. Then $h(z, t)$ is approximated as follows:

$h(i \Delta z, t) \approx h_i(t)\,$

So $h_i(t)$ is a sequence (or vector) that has an index $i$ instead of a continuous dependence on $z$; note however that the whole sequence still has a continuous dependence on time. $i$ is an index that runs from 0 to $n - 1$, note that zero based indexing is being used, so that $i = 0$ corresponds to $z = 0$ and $i = n - 1$ corresponds to $z = 1$. Looking at the boundaries, we know that:

$h_0(t) = 1 \qquad \mbox{and} \qquad h_{n - 1}(t) = 0\,$

What about the points of the sequence inside the boundaries? Initially (at $t = 0$),

$h_i(0) = 0\,$

Suppose now that the right side of the PDE is discretized however you like, so that:

$f_i = \mbox{The spatial discretization of} \quad 2 \left(\frac{\partial h}{\partial z}\right)^2 + h\frac{\partial^2 h}{\partial z^2}\,$

If, for example, central differences were used, one would eventually arrive at:

$f_i(h_{i - 1}, h_{i}, h_{i + 1}) = \frac{1}{2 (\Delta z)^2} (h_{i + 1} h_{i + 1} - 2 h_{i + 1} h_{i - 1} + h_{i - 1} h_{i - 1} + 2 h_i (h_{i + 1} - 2 h_i + h_{i - 1}))\,$

To construct the method of lines approximation for this problem, first $h(z, t)$ is replaced with $h_i(t)$, and then the right side of the equation is replaced with its discretization (the exact equality is a slight abuse of notation):

$\frac{\partial}{\partial t}(h_i(t)) = f_i\,$

Since $h_i(t)$ depends continuously on time and nothing else, this differential equation becomes:

$\frac{d h_i}{d t} = f_i\,$

Putting everything together in one place, the solution to $h_i(t)$ is the solution to the following IVP:

$h_i(0) = 0 \qquad \mbox{(this is the initial value)}\,$

$\begin{cases} h_0 = 1 \qquad \mbox{(this is the left boundary value)} \\ \cfrac{d h_1}{d t} = f_1 \\ \quad \vdots \\ \cfrac{d h_i}{d t} = f_i \\ \quad \vdots \\ \cfrac{d h_{n - 2}}{d t} = f_{n - 2} \\ h_{n - 1} = 0 \qquad \mbox{(this is the right boundary value)} \end{cases}$

Solving this problem gives the approximation for $h(x, t)$. Note that if a second order forward stepping method is used, such as second order Runge Kutta, this solution method will have approximately the same accuracy as the Crank-Nicolson method, but without simultaneous solution of a mess of algebraic equations, which would make a nonlinear problem prohibitively difficult since a tidy matrix solution wouldn't work.

The method of lines is especially popular in electromagnetics, for example, using the Helmholtz equation to simulate the passage of light through a lens for better lens design; the reason for the popularity is that the system of ODEs can (in this case) be solved analytically, so that the accuracy is limited only by the spatial discretization.

A word of caution: an explicit forward stepping method of lines solution bears similarity to a forward stepping finite difference method; so there is no reason to believe that the method of lines doesn't suffer the same stability issues.

A Third Order TVD RK Scheme in C

What follows is an efficient TVD RK scheme implemented in C. Memory is automatically allocated per need (but it will just assume that there is no problem in memory allocation) and freed when n = 0. Note that the solution will not be TVD (total variation diminishing) unless the discretization of $f(x)$ is TVD also.

This isn't strictly a method of lines solver, of course; it may find use whenever some large, interdependent vector ODE must be stepped forward.

/*
Third order TVD autonomous RK solution. Steps x'(t) = f(x) forward, where x is some vector of size n.
 
x is the address of the first element. x must be an array of size n. f is a function as shown above.
It's fed the address of the first element of x and the index of the element being worked on, n.
Use n = 0 to free memory.
*/
void RK3(double (*f)(double *x, unsigned int i, unsigned int n), double *x, unsigned int n, double delta_t){
        static double *x1 = NULL, *x2 = NULL;
        static unsigned int N = 0;
        if(n == 0){
                free(x1);
                free(x2);
                x1 = NULL;
                x2 = NULL;
                return;
        }
        if(n > N){
                N = n;
                x1 = realloc(x1, n*sizeof(double));
                x2 = realloc(x2, n*sizeof(double));
        }
 
        unsigned int i;
 
        for(i = 0; i != n; i++)
                x1[i] = x[i] + delta_t*f(x, i, n);
 
        for(i = 0; i != n; i++)
                x2[i] = 3.0/4.0*x[i] + 1.0/4.0*(x1[i] + delta_t*f(x1, i, n));
 
        for(i = 0; i != n; i++)
                x[i] = 1.0/3.0*x[i] + 2.0/3.0*(x2[i] + delta_t*f(x2, i, n));
 
}

Fortran 90 equivalent:

subroutine RK3(f, x, n, delta_t, pass)
   double precision, external:: f
   integer         :: n
   double precision:: x(n)
   double precision:: delta_t
   double precision:: pass(*) !-- Parameters for the routine f()
   ! Locals
   double precision:: x1(n), x2(n)

   if (n .eq. 0) return
   x1 = x + delta_t *  f(x, n, pass)
   x2 = 3d0/4d0 * x + 1d0/4d0 * (x1 + delta_t * f(x1, n, pass)) 
   x  = 1d0/3d0 * x + 2d0/3d0 * (x2 + delta_t * f(x2, n, pass)) 
end subroutine RK3

f() should be defined as

function f(x, n, pass) result (r)
integer:: n
double precision:: f(n)
double precision:: pass(*)
double precision:: r(n)
! START OF CODE