Advanced Mathematics for Engineers and Scientists/Method of Lines

From Wikibooks, open books for an open world
Jump to: navigation, search

Method of Lines[edit]





Example Problem[edit]

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\,


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}\,


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)}\,

 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)}

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.



A Third Order TVD RK Scheme in C[edit]



Use n = 0 to free memory.
 	static double *x1 = NULL, *x2 = NULL;
	static unsigned int N = 0;
	if(n == 0){
		x1 = NULL;
		x2 = NULL;
	if(n > N){
		N = n;
	unsigned int i;

Fortran 90 equivalent:

   double precision, external:: f
   integer         :: n
   double precision:: x(n)
   double precision:: delta_t
   ! Locals

   if (n .eq. 0) return
end subroutine RK3

f() should be defined as

integer:: n
double precision:: f(n)
double precision:: pass(*)
double precision:: r(n)