On Dec. 20, 2011, Barton Willis wrote:
----------------------
>Numerically evaluating this expression results in subtractive cancellation;
>try:
>
>(%i9) e : integrate(integrate(cos(x*y*log(2)),y,1,10),x,1,10)$
>
>(%i10) makelist(block([fpprec : 15], bfloat(block([fpprec : 10 * k],
>expand(bfloat(e))))), k,1,10);
----------------------
On Dec. 21, 2011, Raymond Toy wrote:
-------------------------
> I think this is a pretty good example where, if the goal is a numerical
>answer, it's better to use a numerical algorithm instead of taking the
>analytic solution and numerically evaluating that. :-)
-------------------------
Following Barton's suggestion, mfloat uses bfloat methods
to go for an answer with a certain number of digits of precision.
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
(%i1) load(nintmd);
(%o1) "c:/work2/nintmd.mac"
(%i2) e:integrate(integrate(cos(x*y*log(2)),y,1,10),x,1,10);
(%o2) -(%i*gamma_incomplete(0,100*%i*log(2))
-2*%i*gamma_incomplete(0,10*%i*log(2))+%i*gamma_incomplete(0,%i*log(2))
-%i*gamma_incomplete(0,-%i*log(2))
+2*%i*gamma_incomplete(0,-10*%i*log(2))
-%i*gamma_incomplete(0,-100*%i*log(2)))
/(2*log(2))
(%i3) mfloat(e,6);
(%o3) -0.959759
(%i4) e:integrate(integrate(sin(x*y*log(2)),y,1,10),x,1,10);
(%o4) (gamma_incomplete(0,100*%i*log(2))-2*gamma_incomplete(0,10*%i*log(2))
+gamma_incomplete(0,%i*log(2))
+gamma_incomplete(0,-%i*log(2))
-2*gamma_incomplete(0,-10*%i*log(2))
+gamma_incomplete(0,-100*%i*log(2)))
/(2*log(2))
(%i5) mfloat(e,6);
(%o5) 0.0613141
------------------
the relevant code for mfloat and mfloat1 is:
---------------------------
/* mfloat accepts a complex expression */
mfloat(expr,mdigits):=
block([er,ei,err],
er : realpart(expr),
err : mfloat1(er,mdigits),
if not complex_number_p(err) then
(print(" realpart not evaluated"),
return(false)),
if expr - er = 0 then err
else
(ei : imagpart(expr),
ei : mfloat1(ei,mdigits),
if not complex_number_p(ei) then
(print(" imagpart not evaluated"),
return(false)),
err + %i*ei))$
/* mfloat1 expects a real expression */
mfloat1(ex,ndigits) :=
block([fpprec:16,old,new,kk,small,goodans:false],
small:bfloat(10^(-(ndigits+1))),
old : expand(bfloat(ex)),
for kk thru 10 do
(new:realpart(block([fpprec:kk*10],expand(bfloat(ex)))),
if abs(new-old) < small then
(goodans:true,return()),
old:new),
if not goodans then false
else if not numberp(new) then false
else float(new))$
--------------------------------------
Ted