[sage-support] Re: numerical evaluation of integral?
Dear all,
Thanks a lot for your help. It looks like there is yet another useful
function going to be implemented in sage. Just to give an example
where quad is faster than GSL (sorry for the length):
sage: def insol1(epsilon,ecc,varpi,phi,lambd):
... S0=1368
... phirad=phi*pi/180.
... lambdarad=lambd*pi/180.
...
... sindelta = sin(epsilon)*sin(lambdarad)
... cosdelta = sqrt(1-sindelta^2)
... nu = lambdarad-varpi # true anomaly
... sin_phi_sin_delta = sin(phirad)*sindelta
... cos_phi_cos_delta = cos(phirad)*cosdelta
... tan_phi_tan_delta = min(max(sin_phi_sin_delta /
cos_phi_cos_delta,-1),1)
... H0 = acos(-tan_phi_tan_delta)# solar time of sunset
... sin_H0 = sqrt(1-tan_phi_tan_delta*tan_phi_tan_delta)
#2*H0: day length
... rho = (1-ecc^2)/(1+ecc*cos(nu)) # earth-sun distance in
astronomical units
... W = S0/(pi*rho^2)*(H0*sin_phi_sin_delta
+cos_phi_cos_delta*sin_H0) # daily mean shortwave down
... dtdl = rho^2 / sqrt(1-ecc^2) # 1/(dlambda/dt)
... return W*dtdl
sage: timeit('quad(lambda x: insol(0.3,0.05,0,80,x),0,360)')
5 loops, best of 3: 711 ms per loop
sage: timeit('numerical_integral(lambda x: insol(0.3,0.05,0,80,x),
0,360)')
5 loops, best of 3: 1.06 s per loop
The difference does not seem to be very large, though.
Cheers
Stan
On Dec 5, 6:17 am, Jason Grout <[EMAIL PROTECTED]> wrote:
> William Stein wrote:
> >> Should we phase GSL out of numerical_integral too? Should we replace it
> >> with the equivalent scipy call (which would make it massively shorter
> >> and simpler)?
>
> > Yes, it is very tempting to do so. One thing is that each function
> > evaluation
> > could in theory be much faster with the GSL version, since GSL takes
> > as input a C-level callback function, whereas I think scipy's quadpack
> > wrapper doesn't.
>
> Using the example from this thread:
>
> sage: timeit('quad(ff, 0, 18)')
> 625 loops, best of 3: 157 µs per loop
> sage: timeit('numerical_integral(ff, 0, 18)')
> 625 loops, best of 3: 78.4 µs per loop
>
> So GSL is twice as fast.
>
> A few more examples pulled from thin air:
>
> sage: f=250*cos(pi*x/180)^1.8 + 170.35*sin(x*pi)+log(1+x)
> sage: ff = fast_float(f, 'x')
> sage: timeit('quad(ff, 0, 18)')
> 625 loops, best of 3: 136 µs per loop
> sage: timeit('numerical_integral(ff, 0, 18)')
> 625 loops, best of 3: 104 µs per loop
> sage: f=250*cos(pi*x/180)^1.8 +
> 170.35*sin(x*pi)+log(1+x)+1/sqrt(1+x)+x^(0.3)
> sage: ff = fast_float(f, 'x')
> sage: timeit('quad(ff, 0, 18)')
> 625 loops, best of 3: 699 µs per loop
> sage: timeit('numerical_integral(ff, 0, 18)')
> 625 loops, best of 3: 142 µs per loop
>
> In each case, GSL is better or way better.
>
>
>
> > Jason said:
>
> >> Both GSL and scipy call quadpack.
>
> > I'm not sure exactly what this means, since probably scipy's quadpack
> > is a fortran library, but GSL is definitely built 100% fortran free
> > (there's no fortran code in gsl and no fortran dependencies). Maybe
> > GSL has a C port of quadpack, or some other sort of translation of the
> > code. So I suspect GSL and scipy are calling into separate separate
> > code that is compiled using a different compiler, so there could be
> > differences in performance and capabilities.
>
> Sorry; from the GSL docs:
>
> The library reimplements the algorithms used in quadpack, a numerical
> integration package written by Piessens, Doncker-Kapenga, Uberhuber and
> Kahaner. Fortran code for quadpack is available on Netlib.
>
> So you're right, there could be performance differences in the library
> alone.
>
> Jason
--~--~-~--~~~---~--~~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
William Stein wrote:
>> Should we phase GSL out of numerical_integral too? Should we replace it
>> with the equivalent scipy call (which would make it massively shorter
>> and simpler)?
>
> Yes, it is very tempting to do so. One thing is that each function
> evaluation
> could in theory be much faster with the GSL version, since GSL takes
> as input a C-level callback function, whereas I think scipy's quadpack
> wrapper doesn't.
Using the example from this thread:
sage: timeit('quad(ff, 0, 18)')
625 loops, best of 3: 157 µs per loop
sage: timeit('numerical_integral(ff, 0, 18)')
625 loops, best of 3: 78.4 µs per loop
So GSL is twice as fast.
A few more examples pulled from thin air:
sage: f=250*cos(pi*x/180)^1.8 + 170.35*sin(x*pi)+log(1+x)
sage: ff = fast_float(f, 'x')
sage: timeit('quad(ff, 0, 18)')
625 loops, best of 3: 136 µs per loop
sage: timeit('numerical_integral(ff, 0, 18)')
625 loops, best of 3: 104 µs per loop
sage: f=250*cos(pi*x/180)^1.8 +
170.35*sin(x*pi)+log(1+x)+1/sqrt(1+x)+x^(0.3)
sage: ff = fast_float(f, 'x')
sage: timeit('quad(ff, 0, 18)')
625 loops, best of 3: 699 µs per loop
sage: timeit('numerical_integral(ff, 0, 18)')
625 loops, best of 3: 142 µs per loop
In each case, GSL is better or way better.
>
> Jason said:
>
>> Both GSL and scipy call quadpack.
>
> I'm not sure exactly what this means, since probably scipy's quadpack
> is a fortran library, but GSL is definitely built 100% fortran free
> (there's no fortran code in gsl and no fortran dependencies). Maybe
> GSL has a C port of quadpack, or some other sort of translation of the
> code. So I suspect GSL and scipy are calling into separate separate
> code that is compiled using a different compiler, so there could be
> differences in performance and capabilities.
Sorry; from the GSL docs:
The library reimplements the algorithms used in quadpack, a numerical
integration package written by Piessens, Doncker-Kapenga, Uberhuber and
Kahaner. Fortran code for quadpack is available on Netlib.
So you're right, there could be performance differences in the library
alone.
Jason
--~--~-~--~~~---~--~~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Thu, Dec 4, 2008 at 7:25 PM, Jason Grout <[EMAIL PROTECTED]> wrote: > > William Stein wrote: >> On Thu, Dec 4, 2008 at 7:11 PM, Tim Lahey <[EMAIL PROTECTED]> wrote: >>> On Dec 4, 2008, at 10:05 PM, Jason Grout wrote: >>> Tim Lahey wrote: >> Jason > Is there an easy way to get the integrand, variable and bounds out of the > integral? That way, if one has tried to analytically evaluate it, they > can pull it out and try numerically evaluating it easily. In fact, it > probably could be done automatically. > sage: a=integrate(250*cos(pi*x/180)^1.8 + 170.35 , x, 0, 18) sage: a integrate(250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) sage: a.arguments() (250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) >>> In that case, it should be simple to feed the integral to scipy. One >>> could write a simple wrapper to numerically integrate integrals which >>> can't be done analytically. >> >> It would be better to call the numerical_integral function >> that is already in Sage, which Josh Kantor wrote, which >> is pretty sophisticated. It uses GSL and a C callback function. >> Then improve the implementation of that function to also use >> scipy. To easy steps instead of one hard one. > > Should we phase GSL out of numerical_integral too? Should we replace it > with the equivalent scipy call (which would make it massively shorter > and simpler)? Yes, it is very tempting to do so. One thing is that each function evaluation could in theory be much faster with the GSL version, since GSL takes as input a C-level callback function, whereas I think scipy's quadpack wrapper doesn't. Jason said: > Both GSL and scipy call quadpack. I'm not sure exactly what this means, since probably scipy's quadpack is a fortran library, but GSL is definitely built 100% fortran free (there's no fortran code in gsl and no fortran dependencies). Maybe GSL has a C port of quadpack, or some other sort of translation of the code. So I suspect GSL and scipy are calling into separate separate code that is compiled using a different compiler, so there could be differences in performance and capabilities. -- William --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
Tim Lahey wrote: > > On Dec 4, 2008, at 10:14 PM, William Stein wrote: >> >> It would be better to call the numerical_integral function >> that is already in Sage, which Josh Kantor wrote, which >> is pretty sophisticated. It uses GSL and a C callback function. >> Then improve the implementation of that function to also use >> scipy. To easy steps instead of one hard one. > > Scipy probably just calls quadpack so it's likely best to have > the numerical_integral function call that instead of scipy. I > was thinking from a UI perspective, though. If one already has > an integral defined and they attempted to integrate it and failed, > one could pass that to a function to evaluate it numerically. I don't > know if that's what the numerical_integral function does. Both GSL and scipy call quadpack. I don't see any reason for us to interface directly with quadpack if it's already been wrapped in scipy. I doubt there will be much of a speed difference, and I'd rather not maintain a wrapper when the scipy guys are perfectly willing to maintain a wrapper. The numerical_integral function calls various routines in GSL, depending on the integral. Jason --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
William Stein wrote: > On Thu, Dec 4, 2008 at 7:11 PM, Tim Lahey <[EMAIL PROTECTED]> wrote: >> On Dec 4, 2008, at 10:05 PM, Jason Grout wrote: >> >>> Tim Lahey wrote: > Jason Is there an easy way to get the integrand, variable and bounds out of the integral? That way, if one has tried to analytically evaluate it, they can pull it out and try numerically evaluating it easily. In fact, it probably could be done automatically. >>> sage: a=integrate(250*cos(pi*x/180)^1.8 + 170.35 , x, 0, 18) >>> sage: a >>> integrate(250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) >>> sage: a.arguments() >>> (250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) >>> >> In that case, it should be simple to feed the integral to scipy. One >> could write a simple wrapper to numerically integrate integrals which >> can't be done analytically. > > It would be better to call the numerical_integral function > that is already in Sage, which Josh Kantor wrote, which > is pretty sophisticated. It uses GSL and a C callback function. > Then improve the implementation of that function to also use > scipy. To easy steps instead of one hard one. Should we phase GSL out of numerical_integral too? Should we replace it with the equivalent scipy call (which would make it massively shorter and simpler)? I'm asking because that seemed to be the preference for numpy vs. gsl for RDF/CDF matrices and vectors. Jason --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Dec 4, 2008, at 10:14 PM, William Stein wrote: It would be better to call the numerical_integral function that is already in Sage, which Josh Kantor wrote, which is pretty sophisticated. It uses GSL and a C callback function. Then improve the implementation of that function to also use scipy. To easy steps instead of one hard one. Scipy probably just calls quadpack so it's likely best to have the numerical_integral function call that instead of scipy. I was thinking from a UI perspective, though. If one already has an integral defined and they attempted to integrate it and failed, one could pass that to a function to evaluate it numerically. I don't know if that's what the numerical_integral function does. Cheers, Tim. smime.p7s Description: S/MIME cryptographic signature
[sage-support] Re: numerical evaluation of integral?
On Thu, Dec 4, 2008 at 7:11 PM, Tim Lahey <[EMAIL PROTECTED]> wrote: > > On Dec 4, 2008, at 10:05 PM, Jason Grout wrote: > >> >> Tim Lahey wrote: >>> Jason >>> >>> Is there an easy way to get the integrand, variable and bounds out of the >>> integral? That way, if one has tried to analytically evaluate it, they >>> can pull it out and try numerically evaluating it easily. In fact, it >>> probably could be done automatically. >>> >> >> sage: a=integrate(250*cos(pi*x/180)^1.8 + 170.35 , x, 0, 18) >> sage: a >> integrate(250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) >> sage: a.arguments() >> (250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) >> > > In that case, it should be simple to feed the integral to scipy. One > could write a simple wrapper to numerically integrate integrals which > can't be done analytically. It would be better to call the numerical_integral function that is already in Sage, which Josh Kantor wrote, which is pretty sophisticated. It uses GSL and a C callback function. Then improve the implementation of that function to also use scipy. To easy steps instead of one hard one. William --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Dec 4, 2008, at 10:05 PM, Jason Grout wrote: Tim Lahey wrote: Jason Is there an easy way to get the integrand, variable and bounds out of the integral? That way, if one has tried to analytically evaluate it, they can pull it out and try numerically evaluating it easily. In fact, it probably could be done automatically. sage: a=integrate(250*cos(pi*x/180)^1.8 + 170.35 , x, 0, 18) sage: a integrate(250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) sage: a.arguments() (250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) In that case, it should be simple to feed the integral to scipy. One could write a simple wrapper to numerically integrate integrals which can't be done analytically. Cheers, Tim. --- Tim Lahey PhD Candidate, Systems Design Engineering University of Waterloo http://www.linkedin.com/in/timlahey smime.p7s Description: S/MIME cryptographic signature
[sage-support] Re: numerical evaluation of integral?
Tim Lahey wrote: > >> Jason > > Is there an easy way to get the integrand, variable and bounds out of the > integral? That way, if one has tried to analytically evaluate it, they > can pull it out and try numerically evaluating it easily. In fact, it > probably could be done automatically. > sage: a=integrate(250*cos(pi*x/180)^1.8 + 170.35 , x, 0, 18) sage: a integrate(250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) sage: a.arguments() (250*cos(pi*x/180)^1.8 + 170.35, x, 0, 18) Jason --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Dec 4, 2008, at 9:38 PM, Jason Grout wrote:
Robert Dodier wrote:
On Dec 4, 2:04 pm, "William Stein" <[EMAIL PROTECTED]> wrote:
sage: f.n()
and get back a floating point number. This is surprisingly not
implemented in Sage, but it isn't.
(That's basically because Maxima itself doesn't seem to have such
functionality.)
I'm guessing that f.n() just turns on the numer flag for Maxima.
numer causes any literal numbers or symbolic constants
to be replaced by floating point values. However the integrate
function is called as without numer. If you want a numerical
integration, call quad_qags or some other Quadpack function.
FYI, scipy has numerical integration, based on quadpack:
http://docs.scipy.org/doc/scipy/reference/integrate.html
sage: from scipy.integrate import quad
sage: f(x) = 250*cos(pi*x/180)^1.8 + 170.35
sage: from sage.ext.fast_eval import fast_float
sage: ff = fast_float(f, 'x')
sage: quad(ff,0,18)
(7435.2795815640284, 8.2548185859776835e-11)
sage: timeit('quad(ff,0,18)')
625 loops, best of 3: 118 µs per loop
There are lots of options you can pass. If you want an infinite
limit,
then, use scipy.integrate.Inf.
It sounds like it would be good to use this if we wanted a numerical
approximation of an integral.
Jason
Is there an easy way to get the integrand, variable and bounds out of
the
integral? That way, if one has tried to analytically evaluate it, they
can pull it out and try numerically evaluating it easily. In fact, it
probably could be done automatically.
Cheers,
Tim.
---
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://www.linkedin.com/in/timlahey
smime.p7s
Description: S/MIME cryptographic signature
[sage-support] Re: numerical evaluation of integral?
Robert Dodier wrote:
> On Dec 4, 2:04 pm, "William Stein" <[EMAIL PROTECTED]> wrote:
>
>> sage: f.n()
>>
>> and get back a floating point number. This is surprisingly not
>> implemented in Sage, but it isn't.
>> (That's basically because Maxima itself doesn't seem to have such
>> functionality.)
>
> I'm guessing that f.n() just turns on the numer flag for Maxima.
> numer causes any literal numbers or symbolic constants
> to be replaced by floating point values. However the integrate
> function is called as without numer. If you want a numerical
> integration, call quad_qags or some other Quadpack function.
FYI, scipy has numerical integration, based on quadpack:
http://docs.scipy.org/doc/scipy/reference/integrate.html
sage: from scipy.integrate import quad
sage: f(x) = 250*cos(pi*x/180)^1.8 + 170.35
sage: from sage.ext.fast_eval import fast_float
sage: ff = fast_float(f, 'x')
sage: quad(ff,0,18)
(7435.2795815640284, 8.2548185859776835e-11)
sage: timeit('quad(ff,0,18)')
625 loops, best of 3: 118 µs per loop
There are lots of options you can pass. If you want an infinite limit,
then, use scipy.integrate.Inf.
It sounds like it would be good to use this if we wanted a numerical
approximation of an integral.
Jason
--~--~-~--~~~---~--~~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Dec 4, 2:04 pm, "William Stein" <[EMAIL PROTECTED]> wrote: > sage: f.n() > > and get back a floating point number. This is surprisingly not > implemented in Sage, but it isn't. > (That's basically because Maxima itself doesn't seem to have such > functionality.) I'm guessing that f.n() just turns on the numer flag for Maxima. numer causes any literal numbers or symbolic constants to be replaced by floating point values. However the integrate function is called as without numer. If you want a numerical integration, call quad_qags or some other Quadpack function. FWIW Robert Dodier --~--~-~--~~~---~--~~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-support] Re: numerical evaluation of integral?
On Wed, Dec 3, 2008 at 12:44 AM, Stan Schymanski <[EMAIL PROTECTED]> wrote:
>
> Dear all,
>
> I would like to evaluate a symbolic equation containing an integral
> numerically:
> ((integrate(250*cos(pi*x/180)^1.8 + 170.35,x,0,18)/a_v)(a_v=1)).n()
> does not work. Is there a way of doing this? The real equation is a
> lot longer than the above, so I am looking for a simple automatic way.
> Usually, if I substitute all symbolic variables with numerical values,
> the function just evaluates to a number, but if I have a symbolic
> integral in there, it does not work.
>
> Thanks for your help!
I'm just going to try to clarify your question. You want to be able
to take, e.g., f defined below:
sage: f = integrate(cos(x)^1.3,x,0,18) + 2/3
sage: f
integrate(cos(x)^1.30, x, 0, 18) + 2/3
and do
sage: f.n()
and get back a floating point number. This is surprisingly not
implemented in Sage, but it isn't.
(That's basically because Maxima itself doesn't seem to have such
functionality.)
Of course you can compute integrals numerically as illustrated below:
sage: var('a_v')
sage: ((numerical_integral(250*cos(pi*x/180)^1.8 +
170.35,0,18)[0]/a_v)(a_v=1)).n()
7435.27958156403
This should be enough "power" for you to hopefully be able to solve
the problem you
have. But obviously it would be desirable that float(foo) would work
when foo is a
symbolic expression involving formal integrals. It probably gets old,
but that's the
sort of thing we'll almost certainly fix in the pynac based symbolics...
-- William
--~--~-~--~~~---~--~~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
