HP Prime Numerical Integration

[Ruda975]

Hello,
I'd like to know the method of numerical integration the hp prime uses, because it's great.
For example, I tried to compute the incomplete elliptic integral of the first kind:

EllipticF(x,m) = integrate 1/sqrt(1 - msin^2(t)) dt from 0 to x
For parameters x=2, m=3, so:
EllipticF(2,3) = integrate 1/sqrt(1 - 3sin^2(t)) dt from 0 to 2

Wolfram Alpha tells me that the result is:
EllipticF(2,3) = 1.001077-1.490278*i

On HP Prime the command would be:

int(1/SQRT(1-3*(SIN(T)^2)),T,0,2)

When approximated in HP Prime this results in 1.00107718634 - 1.49027809758*i

When I tried the same thing in Wolfram Alpha:

N[integrate (1/sqrt(1-3sin^2(t)) dt from 0 to 2)]
the result is 0.9996911316703074 - 1.4848975811276912*i

Then I tried other methods, such as programming my own methods in python. I tried the Romberg integration, Gauss Legandre Quadrature, Simpson Rule, Gauss Kronrod.
None of which worked as good. Even romberg in HP Prime can't calculate it well enough.

romberg(1/sqrt(1-3*(sin(t)^2)),t,0,2) returns the last approximation 0.991198869818-1.47660568998*i

The scipy.integrate.quad in python works well enough, but it's not as good as HP Prime int.
Since 1/sqrt(1-3sin^2(x)) has poles from 0 to 2, I suspect that some pole finding algorithm might be used.

[KeithB]

Look at the thread just below yours:
https://www.hpmuseum.org/forum/thread-18566.html

[Albert Chan]

(08-03-2022 06:15 PM)Ruda975 Wrote:  int(1/SQRT(1-3*(SIN(T)^2)),T,0,2)

When approximated in HP Prime this results in 1.00107718634 - 1.49027809758*i

XCAS int(...) is using the same algorithm, with full 53 bits precision (HP Prime use 48 bits)

XCAS> int(1/sqrt(1-3*sin(x)^2), x, 0., 2.)

Low accuracy, error estimate 3.0628046608e-06.
Error might be underestimated if initial boundary was +/-infinity

1.00107391075 - 1.49027871636*i

XCAS only get 6 digits accuracy, HP Prime getting 7 is just lucky.

With a pole within integral limits, even XCAS getting 6 correct digits is luck.
Error estimate seems to be in the right ballpark ... but, don't bet on it.

---

To accurately calculate EllipticF by integration, we locate the pole, and not touch it.

pole = ± asin(sqrt(1/3)) + n*pi ≈ ± 0.6155 + n*pi
For pi/2 .. 2, there is no pole.

N[integrate (1/sqrt(1-3sin^2(t)) dt from pi/2 to 2)] = -0.31885796059775706 i

Using AGM, complete elliptic integral of the first kind does not care where pole is.
EllipticK(3) = pi/2 / AGM(1, sqrt(1-3)) = 1.0010773804561062 - 1.1714200841467699 i

Sum 2 terms, we have EllipticF(2,3) = 1.0010773804561062 - 1.490278044744527 i