Fortran - Integration advice

Author

Integration advice

Madhusudan Singh

2005-05-16, 3:58 am

Hi

I have a function that remains zero for a large fraction of the domain, but
then goes to order ~ 1e14 - 1e15 in one region.

I am using the routine dqags from quadpack, with limit=4000 and
epsabs=1e-16 and epsrel=1e-08. I am getting integrals of the order of 1e8
with abserr set to 1e-02 - 1e+00.

Now, if I look at the function that is being integrated, I can see that the
answer for the integration is not quite correct in some cases. Is it
possible that dividing the function by 1e15 before returning the value to
dqags, and then multiplying the result by that, could improve the accuracy
of the integration (the change from 0 to 1e13-1e15 is quite sudden) ?

Thanks.

David Flower

2005-05-16, 8:57 am

Its obviously worth trying.

Also worth trying is to compare results in single and double precision.

What I suspect is that your error is related to exactly where in the
'one region' the function changes value. Perhaps a separate call to
your integration routine for that one region would help.

kaci_tizi_ouzou2000@yahoo.ca

2005-05-17, 3:59 am

Hi,

I do not believe that dividing by a constante factor would change your
result had it been wrong or right. This is typical of integramds
exhibiting shapr peaks. What I would suggest you to do is the
following:
[1] Sol. 1: Transform so to obtain a function in which the peak is more
"extended" and then integrate.

[2] Sol. 2: Locate your peak and concentrate and choose a high order
quadrature in the area contributing the most.

[3] Sol. 3: Use adaptive quadrature (I think it is called Konrod
quadrature). It is time consuming since it adapts the number of points
to be used based on the variation of the function (i.e. value of the
derivative).

Hope this helps,

Cheers

Madhusudan Singh wrote:
> Hi
>
> I have a function that remains zero for a large fraction of the
domain, but
> then goes to order ~ 1e14 - 1e15 in one region.
>
> I am using the routine dqags from quadpack, with limit=4000 and
> epsabs=1e-16 and epsrel=1e-08. I am getting integrals of the order of
1e8
> with abserr set to 1e-02 - 1e+00.
>
> Now, if I look at the function that is being integrated, I can see
that the
> answer for the integration is not quite correct in some cases. Is it
> possible that dividing the function by 1e15 before returning the
value to
> dqags, and then multiplying the result by that, could improve the
accuracy
> of the integration (the change from 0 to 1e13-1e15 is quite sudden) ?
>
> Thanks.

Peter Spellucci

2005-05-17, 3:59 pm

In article <42881900$0$79465$14726298@news.sunsite.dk>,
Madhusudan Singh <spammers-go-here@spam.invalid> writes:
>Hi
>
> I have a function that remains zero for a large fraction of the domain, but
>then goes to order ~ 1e14 - 1e15 in one region.
>
> I am using the routine dqags from quadpack, with limit=4000 and
>epsabs=1e-16 and epsrel=1e-08. I am getting integrals of the order of 1e8
>with abserr set to 1e-02 - 1e+00.
>
> Now, if I look at the function that is being integrated, I can see that the
>answer for the integration is not quite correct in some cases. Is it
>possible that dividing the function by 1e15 before returning the value to
>dqags, and then multiplying the result by that, could improve the accuracy
>of the integration (the change from 0 to 1e13-1e15 is quite sudden) ?
>
>Thanks.

a pure scaling should have no effect on the precision (if there are not some
fixed built in termination criteria, which I do not assume) better use
a subdivision of the interval and call dqags only on the subintervals (if you
have a rrude knowledge of the position of the peak.) it may be that the automatic
mesh refinement of dqags sometimes doesn't get the peak precise enough if there
are longer regions with a rather flat integrand
hth
peter

Sponsored Links