Numerical integration methods

[Albert Chan]

(07-18-2022 06:51 PM)Tonig00 Wrote:  Int((x^2*((e^(1/x)^3)-1),x,1,100)

This, performed with a G2 with ROM 14603, takes around a minute giving 4.805 (in home mode).

I think it is a long time for this hardware.

Does anybody knows which numerical method the Prime uses?

Turns out the issue is less related to integration routine.
It is more related to inability to evaluate integrand accurately.
Without this fuzziness, Gaussian Quadrature should not have trouble.

Rewrite expression using expm1() does not help. It get "simplified" away.
We might as well split up integrals, with both ∫ accurately calculated.

CAS> ∫(x^2*exp(1/x^3), x, 1., 100.), ∫(-x^2, x, 1., 100.)

[333337.805043, -333333.]

CAS> sum(Ans)

4.80504345335

log2(333333) ≈ 18.3 bits ≈ 5.5 decimal digits.
Decimal accuracy ≈ 14 - 5.5 = 8 to 9 digits

---

If expm1() work as advertised, apply x = exp(y) substitution is better.
This is because integrand looks very similar to 1/x = ln(x)'

∫(x^2*expm1(1/x^3),x,1.,100.) = ∫(exp(3y)*expm1(1/exp(3y)),y,0.,ln(100.))

lua> Q = require'quad'
lua> f = function(x) return x^2 * expm1(1/x^3) end
lua> g = function(x) x=exp(3*x); return x*expm1(1/x) end

It is difficult to fit 1/x well as polynomial.
For example, Q.gauss() use n=20 Gaussian Quadrature Weights and Abscissae

lua> Q.gauss(f, 1, 100)
4.778930417963698
lua> Q.gauss(g, 0, log(100))
4.805043460316733

lua> Q.gauss(f, 1, 10) + Q.gauss(f, 10, 100)
4.805043229893461
lua> Q.gauss(g, 0, log(10)) + Q.gauss(g, log(10), 2*log(10))
4.805043460319849

The transformation also helped tanh-sinh quadrature method
Note: Q.quad() have default eps = 1e-9

lua> Q.quad(f, 1, 100)
4.805043460309022         4.0108968271012166e-011     111
lua> Q.quad(g, 0, log(100))
4.8050434603198475       1.1409030351596591e-009     58