forum.alglib.net

ALGLIB forum
It is currently Sat Apr 27, 2024 7:14 am

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 5 posts ] 
Author Message
 Post subject: autogkintegrate error estimate
PostPosted: Mon Jan 16, 2012 6:04 pm 
Offline

Joined: Mon Jan 16, 2012 5:58 pm
Posts: 3
Hello!

From the documentation I take it that alglib's numerical integrators can give estimates on the integration error. But I can't find where you get this; is it contained in autgoresults() for example?

Here is the code I want errors from, adapted from the example in the manual:

Code:
// function taking double variable to mpfr variable
void dobtomp(mpfr_t mp, double d) {
   int i;
   ostringstream strs;
   strs << scientific << setprecision(15) << d;
   mpfr_set_str(mp, strs.str().c_str(), 10, GMP_RNDN);
}

using namespace alglib;

// function to be integrated; bessel function of the first kind J_3(x)
void int_function_1_func(double x, double xminusa, double bminusx, double &y, void *ptr)
{
   mpfr_t my, mx;
   mpfr_inits2(200, mx, my, (mpfr_ptr)0);
   dobtomp(mx, x);
   mpfr_jn(my, 3, mx, GMP_RNDN);
   y=mpfr_get_d(my, GMP_RNDN);
}

// uses autogk to calculate integral
int main(int argc, char **argv)
{

    double a = 0;
    double b = 3;
    autogkstate s;
    double v;
    autogkreport rep;

    autogksmooth(a, b, s);
    alglib::autogkintegrate(s, int_function_1_func);
    autogkresults(s, v, rep);
   // prints result to screen
   printf("%.5f\n", double(v));
   
   // calculates analytical result and prints to screen
   mpfr_t m1, m2, m4;
   mpfr_inits2(200, m1, m2, m4, (mpfr_ptr)0);
   mpfr_set_str(m1, "3e0", 10, GMP_RNDN);
   mpfr_set_str(m4, "2e0", 10, GMP_RNDN);
   mpfr_jn(m2, 2, m1, GMP_RNDN);
   mpfr_mul(m2, m2, m4, GMP_RNDN);
   mpfr_j0(m1, m1, GMP_RNDN);
   mpfr_set_str(m4, "1e0", 10, GMP_RNDN);
   mpfr_sub(m4, m4, m1, GMP_RNDN);
   mpfr_sub(m4, m4, m2, GMP_RNDN);
   mpfr_out_str(stdout, 10, 10, m4, GMP_RNDN);
   cout << endl;
   mpfr_clears(m1, m2, m4, (mpfr_ptr)0);

    return 0;
}


Best,
Jorgen


Top
 Profile  
 
 Post subject: Re: autogkintegrate error estimate
PostPosted: Tue Jan 17, 2012 12:14 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 906
Documentation is a bit misleading here. It tells that it is possible to estimate integration error with Gauss-Kronrod formula, and that such estimates are used by algorithm internally, but it does not mean that this estimate is returned.

If you want, you can peek into State.SumErr field, which contains estimate of the total absolute error (it is better to do so after algorithm completion). However, this feature is considered undocumented and may change in future without warning, so you should not rely on it if you want to write production code which will work with future versions of ALGLIB.


Top
 Profile  
 
 Post subject: Re: autogkintegrate error estimate
PostPosted: Thu Jan 19, 2012 12:32 pm 
Offline

Joined: Mon Jan 16, 2012 5:58 pm
Posts: 3
Thank you for the reply!

I'm only using the code for myself so future versions shouldn't be a problem. I'm wondering though, how exact is this error estimate? Is it a ballpark figure or the actual error? And if not the latter, do you know of any library with a numerical integrator (capable of integrating expressions containing Bessel functions, should be smooth enough) which does give you the error?

Help is appreciated!

(I know I could write one myself but my project is starting to drag on)


Top
 Profile  
 
 Post subject: Re: autogkintegrate error estimate
PostPosted: Fri Jan 20, 2012 6:08 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 906
This error estimate is calculated using Gauss-Kronrod formula. Gauss-Kronrod formula is a very accurate quadrature formula with two sets of nodes - one for rough estimation of the results, another one for more accurate value. Best nodes are used as integral value, difference between values obtained on best and worst nodes is used as estimate of the error.

This formula gives accurate error estimates for moderate precisions (relative error from 1.0 to 1.0E-8). However, ALGLIB tries to reduce numerical error down to the numerical noise level. Thus, estimates become only order of magnitude correct, because of noise in the input data and numerical errors. So I think that if total error is 1.0E-14, you should not think that it is actually 1E-14. It is as small as possible with given amount of numerical noise, but may be larger than that.


Top
 Profile  
 
 Post subject: Re: autogkintegrate error estimate
PostPosted: Mon Jan 23, 2012 7:05 pm 
Offline

Joined: Mon Jan 16, 2012 5:58 pm
Posts: 3
Ok, thank you! That was helpful.

J?rgen


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 5 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 231 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group