[VA] SRC #009 - Pi Day 2021 Special

[robve]

Quote:a. Find a real root x in [3, 4] of the following equation:

$$ \int_0^x \left( \frac{\sin t}{t} \mathrm{e}^{t/\tan t} \right)^x\,dt = \frac{x\,x^x}{x!} $$

The LHS is an improper integral with the integrand defined on (0,π). The improper integral can be evaluated for x in [3,pi], for example with the excellent HP PRIME or by using a vintage HP or SHARP BASIC calculator using Romberg with midpoint quadrature (for open intervals, see e.g. NR2 Ch4.4). For example, Romberg open interval integration gives 15.93599534 for x=3.141592654 on my SHARP PC-1350. Edit: I should add that this value is the same as x^x/Gamma(x)=15.93599533 for x=3.141592654 on my SHARP PC-1350 with X=PI: GOSUB "GAMMA": PRINT X^X/Y.

Rewrite the equation to

$$ \int_0^x \left( \frac{\sin t}{t} \mathrm{e}^{t/\tan t} \right)^x\,dt - \frac{x^x}{\Gamma x} = 0 $$

After some hunting on the interval [3,π] we find the root x=π.

That makes this a remarkable equation, which I am not yet sure where it came from.

SHARP BASIC:


' ROMBERG QUADRATURE FOR IMPROPER INTEGRALS WITH OPEN INTERVALS
' Functions to integrate are defined with label "F1", "F2",... should return Y given X
' VARIABLES
' A,B range
' F$,F function label (or number) to integrate
' Y result
' E relative error: integral = Y with precision E (attempts E = 1E-10)
' H step size
' N max number of Romberg steps (=10)
' I iteration step
' U current row
' O previous row
' J,S,X scratch
' A(27..26+2*N) scratch auto-array

100 "QROMO" INPUT "f=F";F
110 INPUT "a=";A
120 INPUT "b=";B
' init and first midpoint step
130 E=1E-9,N=7,F$="F"+STR$ F,H=B-A,X=A+H/2: GOSUB F$: O=27,U=O+N,A(O)=H*Y,I=1
' next midpoint step
140 H=H/3,S=0
150 FOR J=1 TO 3^I STEP 3: X=A+(J-.5)*H: GOSUB F$: S=S+Y,X=A+(J+1.5)*H: GOSUB F$: S=S+Y: NEXT J
' integrate
160 A(U)=H*S+A(O)/3,S=1
170 FOR J=1 TO I: S=9*S,A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1): NEXT J
' loop until convergence
180 IF I>5 LET Y=A(U+I): IF ABS(Y-A(O+I-1))<=E*Y+E PRINT Y: END
190 J=O,O=U,U=J,I=I+1: IF I<N GOTO 140
' no convergence, output result with error estimate
200 E=ABS((Y-A(U+N-2))/(Y+E)): PRINT Y,E: END

300 "F1" Y=(SIN X/X*EXP(X/TAN X))^B: RETURN

400 "GAMMA" IF X<=0 LET Y=9E99: RETURN
410 Y=1+76.18009173/(X+1)-86.50532033/(X+2)+24.01409824/(X+3)
420 Y=Y-1.231739572/(X+4)+1.208650974E-3/(X+5)-5.395239385E-6/(X+6)
430 Y=EXP(LN(Y*2.506628275/X)+(X+.5)*LN(X+5.5)-X-5.5): RETURN


b. Find a real root x in [3, 4] of the following equation:

Again we have an improper integral. We could follow the same strategy as a. to numerically solve it.

However, due to Poisson we already know that

$$ W(x) = \frac{1}{\pi}\int_0^\pi \log\left(1+x\frac{\sin t}{t}\mathrm{e}^{t\cot t} \right) dt $$

Hence,

$$ \Omega = W(1) = \frac{1}{\pi}\int_0^\pi \log\left(1+\frac{\sin t}{t}\mathrm{e}^{t\cot t} \right) dt $$

Rewrite the given equation and take x=π. By Poisson we have

$$ \frac{1}{\Omega} \int_0^\pi \log\left( 1+\frac{\sin t}{t} \mathrm{e}^{t\cot t} \right) dt = \pi $$

Hence the root of the given equation is at x=π. QED