Aaron,
allow me to reply to your Numerical.cc question instead of Chi as
I am to blame for that part of the code. Yes, the second piece is
the "erorr term" which takes the averages of second derivatives along
the whole integration interval. True, formally one should look for
the maximum of y2 at any given sub-interval, and by finding them,
adjust the mesh spacing such that the whole error becomes smaller
than eps. In practice, we are making the assumption that the integrands
are reasonably behaved funcs for which second derivatives exist in all
points of the integration interval, and are smooth as well. If you sum up
the error contributions from all subintervals in the full range [a,b],
you end up with just adding up
h^3 / 12 * ( f''(c1) + f''(c2) + ... + f''(cn) ) / n == h^3 / 12 * f''(c)
(i.e. with averaging over all second derivatives) where
x0 < c1 < x1, x1 < c2 < x2, etc. and x0 < c < xn-1. The code does just that.
Again, in principle, you should seek for c such that the above formula is
exactly fulfilled. In practice, if you work with non-pathological functions,
it is almost "true by construction".
At any rate, this is a h^3 contribution. I would be surprised to see
any effects other than, say, 6th digit stuff, from this term even
if you cut it completely.
Simon
PS: no absolute signs here, minus at the front, the sign of f''() does
the job. The construction with the absolute sign you quote could be
used in the hypothetical true evaluation of f''(c) and corresponding
adjustment of the mesh spacings.
On Fri, 6 Aug 2004, Aaron Joseph Maschinot wrote:
>
> hi, chi:
>
> in the splint function inside of Numerical.h, there is the following line
> for implementing trapezoidal numerical integration:
>
> s += (hi/2.0) * (y[i]+y[i+1] - (hi*hi/12.)*(y2[i]+y2[i+1]));
>
> the first contribution is the trapezoidal area; the second is the error
> term, i assume. am i interpreting this right?
>
> if so, is this formula correct? my reference says that the error term for
> the trapezoidal rule is:
>
> err = - h1 * h1 * h1 * y2[c] / 12.0
>
> where c is the value such that y2[c] is maximum on the interval. i wonder
> about the relevance of taking the end-point average second derivative and
> assuming that this is the maximum.
>
> also, can you even add the error term like this? doesn't the error
> derivation require an absolute value sign:
>
> abs( int(f(x),x=a,b) - trap(f(x),x=a,b) ) < err
>
>
> aaron
>
-- Simon Sirca Dept of Physics, University of Ljubljana Tel: +386 1 4766-574 Jadranska 19 Fax: +386 1 2517-281 1000 Ljubljana, Slovenia
This archive was generated by hypermail 2.1.2 : Mon Feb 24 2014 - 14:07:31 EST