[sage-devel] exact numerical integration
I tried doing some integrals today and the output doesn't make much sense to me: sage: f = e^(-x^2) sage: f.integrate(x, 0, 0.1) 2066*sqrt(pi)/36741 sage: f.integrate(x, 0, 1/10) sqrt(pi)*erf(1/10)/2 H. Does this mean erf(1/10) is a rational number? That's a little surprising to me. In fact: sage: RR(f.integrate(x, 0, 0.1)) 0.0996676643523801 sage: RR(f.integrate(x, 0, 1/10)) 0.0996676642903363 What's going on here? david --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@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-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Fwd: sage-devel exact numerical integration
Begin forwarded message: From: Andrzej Chrzęszczyk [EMAIL PROTECTED] Date: March 5, 2008 6:23:53 PM EST To: [EMAIL PROTECTED] Subject: sage-devel exact numerical integration Dear David Try sage: maxima_console() (%i1) integrate(%e^(-x^2),x,0,0.1); ... `rat' replaced .05623145800914245 by 2066/36741 = .05623145804414686 2066 sqrt(%pi) (%o1) -- 36741 then you will see that (behind the scene) maxima replaces more accurate result .05623145804414686 sqrt(%pi) by the less accurate one: 2066 sqrt(%pi)/36741 (default maxima behaviour) Your exact calculations are OK but why do you mix the exact-inexact. Pure numerical version using GSL: sage: numerical_integral(lambda x:e^(-x^2),0,-.1) (-0.099667664290336327, 1.1065333570574191e-15) is in a good accordance with your exact calculations Andrzej Chrzeszczyk (I'm not in sage-devel so I'm using your e-mail adress, I hope You will excuse me) Okay, I can see this makes sense from within Maxima, since you get to see the replacement message. But in Sage, it's really terrible! When I do sage: f = e^(-x^2) sage: f.integral(x, 0, 0.1) 2066*sqrt(pi)/36741 I have absolutely no idea what is going on in the background. It could be maxima, or sympy, or some other library that someone plugged in, or who knows what. How am I supposed to tell that 2066*sqrt(pi)/36741 is not an exact answer? Since it contains sqrt(pi), it's very suggestive that it's supposed to be exact. Another example: if I do sage: f = x*e^(2*x) sage: f.integral(x, 0, 1) e^2/4 + 1/4 Is that an exact answer? Or it just close enough to e^2/4? What use is the answer if I can't tell? david --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@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-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: Fwd: sage-devel exact numerical integration
David Harvey wrote: Begin forwarded message: *From: *Andrzej Chrzęszczyk [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] *Date: *March 5, 2008 6:23:53 PM EST *To: [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] *Subject: **sage-devel exact numerical integration* Dear David Try sage: maxima_console() (%i1) integrate(%e^(-x^2),x,0,0.1); ... `rat' replaced .05623145800914245 by 2066/36741 = .05623145804414686 2066 sqrt(%pi) (%o1) -- 36741 then you will see that (behind the scene) maxima replaces more accurate result .05623145804414686 sqrt(%pi) by the less accurate one: 2066 sqrt(%pi)/36741 (default maxima behaviour) According to http://www.ma.utexas.edu/maxima/maxima_11.html Function: RAT (exp, v1, ..., vn) converts exp to CRE form by expanding and combining all terms over a common denominator and cancelling out the greatest common divisor of the numerator and denominator as well as converting floating point numbers to rational numbers within a tolerance of RATEPSILON[2.0E-8]. and KEEPFLOAT[FALSE] if TRUE prevents floating point numbers from being converted to rational numbers. Maybe we should set the keepfloat thing or set ratepsilon to a much smaller value? Jason --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@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-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: Fwd: sage-devel exact numerical integration
Jason Grout wrote: David Harvey wrote: Begin forwarded message: *From: *Andrzej Chrzęszczyk [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] *Date: *March 5, 2008 6:23:53 PM EST *To: [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] *Subject: **sage-devel exact numerical integration* Dear David Try sage: maxima_console() (%i1) integrate(%e^(-x^2),x,0,0.1); ... `rat' replaced .05623145800914245 by 2066/36741 = .05623145804414686 2066 sqrt(%pi) (%o1) -- 36741 then you will see that (behind the scene) maxima replaces more accurate result .05623145804414686 sqrt(%pi) by the less accurate one: 2066 sqrt(%pi)/36741 (default maxima behaviour) According to http://www.ma.utexas.edu/maxima/maxima_11.html Function: RAT (exp, v1, ..., vn) converts exp to CRE form by expanding and combining all terms over a common denominator and cancelling out the greatest common divisor of the numerator and denominator as well as converting floating point numbers to rational numbers within a tolerance of RATEPSILON[2.0E-8]. and KEEPFLOAT[FALSE] if TRUE prevents floating point numbers from being converted to rational numbers. Maybe we should set the keepfloat thing or set ratepsilon to a much smaller value? Examples: sage: maxima_console() Maxima 5.13.0 http://maxima.sourceforge.net Using Lisp CLISP 2.41 (2006-10-13) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. This is a development version of Maxima. The function bug_report() provides bug reporting information. (%i1) keepfloat: true; (%o1)true (%i2) integrate(%e^(-x^2),x,0,0.1); (%o2).05623145800914245 sqrt(%pi) (%i3) keepfloat: false; (%o3)false (%i4) ratepsilon: 1e-20; (%o4)9.999E-21 (%i5) integrate(%e^(-x^2),x,0,0.1); `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced 0.1 by 1/10 = 0.1 `rat' replaced .05623145800914245 by 244715115/4351925482 = .05623145800914245 244715115 sqrt(%pi) (%o5) --- 4351925482 Jason --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@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-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: Fwd: sage-devel exact numerical integration
In the present example all data are exact and the results is OK In the previous example you have mixed symbolic and numerical data (0.1) and that was risky Andrzej On 5 Mar, 23:59, David Harvey [EMAIL PROTECTED] wrote: Begin forwarded message: From: Andrzej Chrzęszczyk [EMAIL PROTECTED] Date: March 5, 2008 6:23:53 PM EST To: [EMAIL PROTECTED] Subject: sage-devel exact numerical integration Dear David Try sage: maxima_console() (%i1) integrate(%e^(-x^2),x,0,0.1); ... `rat' replaced .05623145800914245 by 2066/36741 = .05623145804414686 2066 sqrt(%pi) (%o1) -- 36741 then you will see that (behind the scene) maxima replaces more accurate result .05623145804414686 sqrt(%pi) by the less accurate one: 2066 sqrt(%pi)/36741 (default maxima behaviour) Your exact calculations are OK but why do you mix the exact-inexact. Pure numerical version using GSL: sage: numerical_integral(lambda x:e^(-x^2),0,-.1) (-0.099667664290336327, 1.1065333570574191e-15) is in a good accordance with your exact calculations Andrzej Chrzeszczyk (I'm not in sage-devel so I'm using your e-mail adress, I hope You will excuse me) Okay, I can see this makes sense from within Maxima, since you get to see the replacement message. But in Sage, it's really terrible! When I do sage: f = e^(-x^2) sage: f.integral(x, 0, 0.1) 2066*sqrt(pi)/36741 I have absolutely no idea what is going on in the background. It could be maxima, or sympy, or some other library that someone plugged in, or who knows what. How am I supposed to tell that 2066*sqrt(pi)/36741 is not an exact answer? Since it contains sqrt(pi), it's very suggestive that it's supposed to be exact. Another example: if I do sage: f = x*e^(2*x) sage: f.integral(x, 0, 1) e^2/4 + 1/4 Is that an exact answer? Or it just close enough to e^2/4? What use is the answer if I can't tell? david --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@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-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
