HP42s first major program (Double Integral) Best way to approach?

[Albert Chan]

(05-27-2020 06:15 AM)Werner Wrote:  Then the inner integral becomes

int(a,b,r/(sqrt(r^2-x^2),dr)

which is simply sqrt(b^2-x^2) - sqrt(a^2-x^2) ...

Continued on, we have

\(I = \int_a^b\int_0^c\frac{r\sqrt{c^2-x^2}}{\sqrt{r^2-x^2}}\;dx\;dr
= \int _0 ^c \sqrt{c^2-x^2}\left( \sqrt{b^2-x^2} - \sqrt{a^2-x^2} \right)\;dx\)

The shape look like an ellipse, with semi-axis, c, c*(b-a)

a, b, c = 24, 29, 12
I ≈ pi/4 * (b-a) * c^2 ≈ 565.4866776461628

If integrand can be transformed to bell-shape, integral is easy to integrate.
see https://www.hpmuseum.org/forum/thread-13...#pid127590

Let \(x = c \cos(t)\; → dx = -c \sin(t)\;dt\)

\(I = \int _0 ^{\pi /2} (c^2-x^2) \left( \sqrt{b^2-x^2} - \sqrt{a^2-x^2}\right)\;dt \)

Shape of transformed curve: https://www.wolframalpha.com/input/?i=pl...0+to+pi%29

Code:

function I(a,b,c, n)    -- n = number of trapezoids
    local A, B = (a+c)*(a-c), (b+c)*(b-c)
    local function f(y) return y/(sqrt(B+y)+sqrt(A+y)) end
    local function g(t) t=c*sin(t); return f(t*t) end
    n = n or 8          -- default = 8 trapezoids
    local h, t = pi/(n+n), 0
    for i=1,n-1 do t = t + g(i*h) end
    return h * (b-a) * ((b+a)*t + c*c/2)
end

Convergence with trapezoid rule is amazingly fast:

lua> a, b, c = 24, 29, 12
lua> for n=1,8 do print(n, I(a,b,c, n)) end
1       565.4866776461628
2       581.3714770410786
3       581.40570463529
4       581.405804245352
5       581.4058045759407
6       581.4058045771336
7       581.4058045771384
8       581.4058045771384