Detailed discussion of Gosse's numerical integration question below.
A few points are worth pulling out:
-- Problem reports should include a complete procedure/configuration
for reproducing the problem and an explanation of the derivation
of the putatively correct answer. (I could not reproduce some of
Gosse's issues.)
-- We need usage-oriented documentation in addition to pure reference
documentation; in this case, simply setting rombergmin to a higher
value resolves the problem. (Easy to do for this case.)
-- Error messages should be more helpful. (Easy to do for this case.)
-- Handling NaNs etc. properly would be a good thing. (But this is
a big project.)
-- Defaults should be more user-friendly, and proportional to current
machine speeds; in particular, rombergit should probably be 16 or
more, and rombergmin should be 5. (This is easy to do for this
case.)
-- It would be nice to have a more intelligent, adaptive numerical
integration routine, which perhaps does some symbolic
preprocessing.
This wouldn't be the routine used by power users for production
runs,
but would be useful to casual and new users. (A project.)
Best,
-s
---------------------------
> Is it possible to calculate with maxima :
> integrate(t**(10.84)*2.71828**(-t*t),t,0,1000);
> maxima 5.9 returns nothing.
It returns the unevaluated integral because the 'integrate' command is
symbolic, not numeric, integration, and Maxima can't find the result in
closed form.
If you integrate numerically using romberg(...), it returns 0.0, because
by default the minimum number of iterations Romberg performs is 0; so it
doesn't even 'see' the non-zero part of the function. Try setting
rombergmin to some larger value, and adjusting rombergit so that it
converges:
romberg(t**(10.84)*2.71828**(-t*t),t,0,1000),
rombergmin:4;
=> ROMBERG failed to converge
Let's increase the maximum number of iterations to see if it will
converge:
romberg(t**(10.84)*2.71828**(-t*t),t,0,1000),
rombergmin:4,rombergit:15;
=> 52.3756 (3.85 sec)
Rombergmin is documented, but it is not clear from the main body of
Romberg documentation how relevant it might be in a case like this. As
usual, the documentation is OK for reference (looking up specific
facts), but not very useful for helping users actually use the system
effectively. It would be nice, too, if the "ROMBERG failed to converge"
error message were more helpful, e.g. "Romberg failed to converge:
adjusting Rombergit and Rombergtol may be helpful".
It seems to me that the default for rombergmin ought to be larger than
0. Perhaps something like 5 would be more reasonable; that still takes
negligeable time in the linear case where iteration isn't needed.
Also, rombergit could be increased. Given that machines are several
hundred times faster now than when this default was originally set, I'd
increase it by at least 5 (factor of 32 slower) to 16.
As a sanity check, you can look at some close-by closed-form integrals:
i10: integrate(t^10/exp(t^2),t,0,1000);
float(i10) => 0.0 oops! (underflow*overflow: should really be NaN)
bfloat(i10) => ... erf(1.0b3) ... oops! bfloat doesn't handle erf
float(expand(i10)) => 26.17 OK
i11: integrate(t^11/exp(t^2),t,0,1000);
float(i10) => 60.0 OK
float(expand(i10)) => 60.0 OK
As for your other example:
> integrate(f**(2/9.84),f,0,1000);
> returns 299. However, it is 1218.
When I try this, it returns 30750*2^(25/41)*5^(25/41)/37, the exact
result after rationalizing 9.84; float(%) => 3383.8. Integrating
numerically, I get:
romberg(f**(2/9.84),f,0,1000) => 3383.685794569355 (.27 sec)
for a relative error of 2.2e-5, less than rombergtol=1.0e-4. With
rombergtol=1.0e-6, rombergit:22, I get a result with a relative error of
-3.4e-7, though it does take 9.4 sec.
What was the exact procedure/configuration which gave 299, and what was
the calculation which gave 1218 as the correct answer?
I tested under Maxima 5.9.0/GCL 2.5.0/mingw32/W2k/Athlon.
-s