Re: Irregularly gridded Discrete Laplacian Operator

[David Bateman]

Daniel J Sebald wrote:
>> 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:
That is what was in the del2.m version on octave-forge/

>
>> 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.

Thanks Dan, Marco pretty much convinced me with a norm of a perturbed
mesh which basically demonstrated the slowed convergence.. I sent my
matlab compatible version of del2 as a patch against the CVS to the list
yesterday..

D.


>
> Dan
>


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary