Irregularly gridded Discrete Laplacian Operator

[Daniel J Sebald]

> The definition for the interior points along an axis of 1D problem for
> the octave-forge version of del2 can be written as
>
> D (2:end-1) =  (M(3:end) - 2 * M(2:end  - 1) + M(1:end -2)) ./
> (dx(1:end-1) .* dx(2:end)) ./ 2

I think the above definition is a bit suspect.  Think of deriving the 
approximation to the Laplacian operator, which I believe the second formula:

> whereas the equivalent in Matlab appears to be
>
> D(2:end-1) = ((M(3:end) - M(2:end-1)) ./  dx(2:end) + (M(1:end-2) -
> M(2:end-1)) ./ dx(1:end-1)) ./ (dx(1:end-1) + dx(2:end))

comes closer to.  Imagine first taking the first derivative, and then taking 
another derivative on the first derivative, i.e., two steps.  Taking the first 
order derivative:

(M(3:end) - M(2:end-1)) ./  dx(2:end)

is the derivative approximation at one location, and

(M(1:end-2) - M(2:end-1)) ./ dx(1:end-1)

is *minus* the derivative approximation at the second location.  Taking the 
difference (i.e., adding because the terms are switched around in the second 
case) and then dividing by the interval (dx(1:end-1) + dx(2:end)) gives the 
second derivative.

Note, if we carry out the math for the second formula we get

D(2:end-1) =
( (M(3:end)-M(2:end-1))*dx(1:end -1) + (M(1:end-2)-M(2:end-1))*dx(2:end) )
./ ( dx(1:end-1)*dx(2:end)*(dx(1:end-1) + dx(2:end)) )

Now, compare the above and the very first formula and you will see that the 
first formula assumes equispacing, i.e., dx(1:end -1) = dx(2:end) in the 
numerator so that the terms factor out and cancel with (dx(1:end-1) + 
dx(2:end)) of the denominator, but that is mistaken algebra.

Unless the spacing is severerly skewed--i.e., dx(1:end -1) != dx(2:end) and in 
a major way--the difference between the appoximations shouldn't appear too bad; 
but yes, convergence might show the flaw.  So, if you want to compare the two, 
replace the very first formula with what is written above and I suspect the 
results will pretty much match.

Dan