+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: HP Prime Numerical Integration (/thread-18631.html)
HP Prime Numerical Integration - Ruda975 - 08-03-2022 06:15 PM
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.
RE: HP Prime Numerical Integration - KeithB - 08-04-2022 12:18 PM
Look at the thread just below yours:
https://www.hpmuseum.org/forum/thread-18566.html
RE: HP Prime Numerical Integration - Albert Chan - 08-04-2022 04:50 PM
(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