commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <Luc.Maison...@free.fr>
Subject Re: [math] Getting second derivative for discrete dataset
Date Sat, 31 Mar 2012 17:13:41 GMT
Le 31/03/2012 01:31, SUJIT PAL a écrit :
> Thanks Luc, thanks for getting back so quickly, and will wait for your information tomorrow.

Hi Sujit,

Sorry for the delay.

If you know your sample points (xi, yi) are all separated by the same
interval h (i.e x(i+1) = x(i) + h for all i), then you can compute any
derivatives using the following method:

First write down the series expansion of the function y = f(x):

  y(x + k * h) = y(x)
               + k * h * y'(x)
               + (k * h)^2 / 2 y''(x)
               + ...
               + (k * h)^n / n! yn(x)
               + ...

Then you truncate this expression to at least include the derivative you
want, and apply it to a number of points equal to your truncation order.
Here, as an example, we will truncate just to order 2 in order to have
y'', so we have three components (y, y' and y''), so we need three
points. Let's take k = -1, k = 0 and k = 2. We have the following
expressions:

 y(x - h) = y(x) - h y'(x) + h^2/2 y''(x)
 y(x)     = y(x)
 y(x+h)   = y(x) + h y'(x) + h^2/2 y''(x)

Note that the first and last expression differ only in the sign of the
odd terms h y'(x), the even terms (h^0 and h^2) have the same sign.

You can rewrite these expressions as a simple square matrix vector
multiplication:

 [ y(x-h) ]   [ 1  -h  h^2/2 ]   [  y(x)  ]
 [ y(x)   ] = [ 1   0  0     ] * [  y'(x) ]
 [ y(x+h) ]   [ 1   h  h^2/2 ]   [ y''(x) ]

Then you invert this system and get:

 [  y(x)  ]    [    0,        1,        0    ]   [ y(x-h) ]
 [  y'(x) ]  = [ -1 / 2h,     0,      1 / 2h ] * [ y(x)   ]
 [ y''(x) ]    [  1 / h^2, -2 / h^2,  1 / h^2]   [ y(x+h) ]

The third row is the one you want:

 y''(x) = (y(x-h) - 2 * y(x) + y(x+h)) / h^2

You can check this formula using again the initial expansion, and you
will see this combination basically canels out the y(x) and y'(x) part
and only y''(x) remains.

You could also have used more points, which could improve accuracy as it
fits your function to an higher degree. Beware that you should *not*
increase the number of points too much, otherwise you would get Gibbs
oscillation problems (i.e a polynomial that has very large changes
between the sampling points).

Hope this helps
Luc


> 
> -sujit
> 
> On Mar 30, 2012, at 2:36 PM, Luc Maisonobe wrote:
> 
>>
>>
>>
>> SUJIT PAL <sujit.pal@comcast.net> a écrit :
>>
>>> Hi,
>>>
>>> I have a (newbie) question about how to go about solving the problem
>>> below with commons-math.
>>>
>>> 1) I have histogram data (equal x intervals) based off some
>>> distribution.
>>> 2) I need a way to calculate the second differential between two given
>>> x points in the distribution.
>>>
>>> Here is what I think I should do:
>>> 1) Feed the (x,y) values represented by the top of each histogram entry
>>> within my x-range into a UnivariateRealInterpolator between ranges
>>> (x1,y1) and (x2,y2).
>>> 2) The interpolate method gives me back a PolynomialSplineFunction
>>> (call if psf).
>>> 3) I get the second differential UnivariateRealFunction (call it urf)
>>> by calling psf.getPolynomialSplineDerivative().derivative()
>>> 4) I compute the second derivative value as the difference of
>>> urf.value(x2) - urf.value(x1)
>>>
>>> Does this sound like a reasonable/accurate approach? Any
>>> suggestions/gotchas or pointers to better approaches?
>>
>> No, this is not an approach I would recommend. A much more straightforward
>> way to compute this derivative is to apply a dedicated finite differences formula
>> on three consecutive points. I'll give ou such a formula tomorrow (I am on a mobile
>> device right now and it is difficult to write a complète mail...)
>>
>> Luc
>>
>>>
>>> Thanks very much,
>>>
>>> Sujit
>>>
>>>
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: user-help@commons.apache.org
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>> For additional commands, e-mail: user-help@commons.apache.org
>>
> 
> 
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
> For additional commands, e-mail: user-help@commons.apache.org
> 


---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message