help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## Re: [Help-gsl] Origin of 2nd round-off term "dy" in function central_der

 From: Rene Girard Subject: Re: [Help-gsl] Origin of 2nd round-off term "dy" in function central_deriv.c Date: Mon, 19 Feb 2007 20:28:09 -0500 (EST)

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

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]