Advanced Mathematics for Engineers and Scientists/Method of Lines

Method of LinesEdit


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 ProblemEdit

Consider the following nonlinear IBVP, where  :


This intimidating nondimensional system describes the flow of liquid in an interior corner ("groove"), where   is the height of the liquid and   is the distance along the axis of the corner. The IBVP states that the corner is initially empty, and that at one end ( ) the height is fixed to 1. Naturally, fluid will flow along  , increasing   until the other boundary ( ), where   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,  . Suppose that   is discretized in space but not time to   points on a uniform grid. Then   is approximated as follows:


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


What about the points of the sequence inside the boundaries? Initially (at  ),


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


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


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


Since   depends continuously on time and nothing else, this differential equation becomes:


Putting everything together in one place, the solution to   is the solution to the following IVP:


Solving this problem gives the approximation for  . 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 CEdit

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   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){
		x1 = NULL;
		x2 = NULL;
	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)