[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Irregularly gridded Discrete Laplacian Operator
From: |
Daniel J Sebald |
Subject: |
Irregularly gridded Discrete Laplacian Operator |
Date: |
Thu, 19 Jul 2007 22:39:11 -0500 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.7.3) Gecko/20041020 |
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