[Top][All Lists]
[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