Re: [Help-gsl] Origin of 2nd round-off term "dy" in function central_deriv.c

[Rene Girard]

Hi Brian,

Thanks for the quick reply. 

"I think the dy contribution is trying to capture the rounding error from terms 
like x+h which is O(|x|*DBL_EPSILON),"  Is "dy" not rather trying to express 
the round-off due to cancellation caused by taking difference like "f(x+h) - 
f(x-h)" ?
Note that in the expression for "dy" we have max between absolute value of r3 
and r5 which are differences. I do not agree that the expression for "dy" 
should be divided by h^2. In fact, the round off error is directly propotional 
to 1/h because on the denominator of the expression for the derivatives of both 
O(h^2) and O(h^4) we find only "h". To substantiate this, let us note here that 
to obtain the expression for h_opt we have to start with 

   T(h) = To *h^2/(ho^2)  where To is the truncation error based on "ho" and T 
is
                                          the truncation error based on "h". 
Also the
                                          stepsize is squared because  the 
derivative 
                                          is at least of O(h^2) 
   R(h) = Ro*ho/h where Ro is the round-off error based on ho 
                                    and R is based on h.

To find h_opt we must minimize the function g(h) = T(h) + R(h). Taking the 
derivative with respect to h and equation the resulting expression to 0 we have:

        2*To*h/(ho^2) - Ro*ho/h^2 = 0

       from which we get :

          h^3 = (ho)^3 * (Ro/(2*To))

and h_opt = ho*(Ro/(2*To))^(1/3)

So it difficult for me to see why the "dy" term should be divided by "h^2" 
Note that to obtain the expression for h_opt, the "dy" term is not included in 
the term for roud-off error in the function g(h). Even if it was included it 
would not change the expression for h_opt as the "dy" term is not a function of 
"h" and its derivative with respect to "h" would be zero. 

Also note that dividing the "dy" by h^2 would increase further the round-off 
error.

I did not see in any of the textbooks I consulted (listed in my previous 
e-mail) an expression for "dy" and yet that it should be divided by h^2.

I am sorry to be somehow persistant but what is the exact origin of 2nd 
round-off term "dy" and how is it derived ?

This is important because I am trying to develop functions to compute 
cross-derivative d^2f(x,y)/dxdy and where the expression
for h_opt = h*(Ro/To)^(1/4) where the derivative will be of O(h^2) but
what will be the correct expression for "dy" ?

Thanks you in advance for your collaboration.

Regards

Rene Girard



        

    



Brian Gough <address@hidden> wrote: At Fri, 16 Feb 2007 21:00:49 -0500 (EST),
Rene Girard wrote:
> On line 24 of the function "central_deriv.c" we have the following
> line:
> 
> double dy = GSL_MAX(fabs(r3),fabs(r5))* fabs(x)*GSL_DBL_EPSILON
> 
> I understand the rest of the function very well and I am able to
> derive the equation for the optimal stepsize "h_opt" in function
> "gsl_deriv_central.c"; however, I am trying to look for a derivation
> for the expression of "dy". I have look in numerous textbooks on
> numerical analysis, in particular the book of Conte and de Boor
> given as a reference in the GSL Reference Manual. 

Thanks for the email.  I think the dy contribution is trying to
capture the rounding error from terms like x+h which is
O(|x|*DBL_EPSILON), but it looks like there is a factor of 1/h^2
missing.

-- 
Brian Gough



                
---------------------------------
Share your photos with the people who matter at Yahoo! Canada Photos