Discussion:
[SciPy-User] Problems with scipy.integrate.quad and sinusoidal weights
Per Nielsen
2012-09-28 08:51:58 UTC
Permalink
Hi all,

I am getting some strange (wrong or highly unaccurate) results when I try
to use scipy.integrate.quad to integrate highly oscillatory functions.
Please consider the following code:

In [40]: from scipy.integrate import quad

In [41]: from math import exp, sin, cos

In [42]: fsin = lambda x: exp(-x) * sin(100*x)

In [43]: quad(fsin, 0., 10., args=(), weight='sin', wvar=100.)
Out[43]: (0.4974066723952844, 0.0005303917238325333)

In [44]: quad(fsin, 0., 10., args=())
Out[44]: (0.0099990111781435, 0.0006814046027542903)

where the last line, without the weight function gives the correct result.

I have looked at the documentain for quad at:

https://github.com/scipy/scipy/blob/master/scipy/integrate/quadpack.py#L134

and to me this should be the way to integrate fsin. Am I misunderstanding
the arguments or whats going on? :)

Cheers,
Per
j***@gmail.com
2012-09-28 12:13:07 UTC
Permalink
Post by Per Nielsen
Hi all,
I am getting some strange (wrong or highly unaccurate) results when I try to
use scipy.integrate.quad to integrate highly oscillatory functions. Please
In [40]: from scipy.integrate import quad
In [41]: from math import exp, sin, cos
In [42]: fsin = lambda x: exp(-x) * sin(100*x)
In [43]: quad(fsin, 0., 10., args=(), weight='sin', wvar=100.)
Out[43]: (0.4974066723952844, 0.0005303917238325333)
In [44]: quad(fsin, 0., 10., args=())
Out[44]: (0.0099990111781435, 0.0006814046027542903)
where the last line, without the weight function gives the correct result.
https://github.com/scipy/scipy/blob/master/scipy/integrate/quadpack.py#L134
and to me this should be the way to integrate fsin. Am I misunderstanding
the arguments or whats going on? :)
the weight function is automatically added
Post by Per Nielsen
integrate.quad(lambda x: np.exp(-x)*np.sin(100*x), 0., 10., args=())
Warning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
(0.0099990111781434986, 0.00068140460275433325)
Post by Per Nielsen
integrate.quad(lambda x: np.exp(-x), 0., 10., args=(), weight='sin', wvar=100.)
(0.0099987410521618428, 7.0983734843630013e-10)

(I never used wvar before, just guessing

Josef
Post by Per Nielsen
Cheers,
Per
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user
Per Nielsen
2012-09-28 14:02:36 UTC
Permalink
Post by j***@gmail.com
the weight function is automatically added
Post by Per Nielsen
integrate.quad(lambda x: np.exp(-x)*np.sin(100*x), 0., 10., args=())
Warning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
(0.0099990111781434986, 0.00068140460275433325)
Post by Per Nielsen
integrate.quad(lambda x: np.exp(-x), 0., 10., args=(), weight='sin',
wvar=100.)
(0.0099987410521618428, 7.0983734843630013e-10)
(I never used wvar before, just guessing
Josef
Ahh, you are completely right. I was thinking that weight function meant
what basis function it used in the numerical integration, but of course it
simply multiply it to the integrand.

Thanks for making that clear to me :)

Per
Post by j***@gmail.com
Post by Per Nielsen
Cheers,
Per
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user
Loading...