Crash when Solve for Integral with 15C

[Gerson W. Barbosa]

(09-03-2017 06:47 PM)Gerson W. Barbosa Wrote:  Only the first 22 digits are correct. 600+ terms would be necessary for 100 digits, but a clean 50g will run out of memory after around 380 terms.

EXPANDing inside the loop is slower at first, but becomes faster as the number of terms increases, not to mention a lesser memory requirement. Still, not a recommended method for the task (727 seconds on the emulator).

600

« 0 DUP ROT
FOR n n 3 * 1 + DUP n ! * 6 ROT ^ SWAP / + EXPAND
NEXT →STR "/" " " SREPL DROP "'" " " SREPL DROP OBJ→ SWAP 12 ALOG * SWAP IQUOT →STR DUPDUP SIZE DUP 11 - SWAP SUB "." OVER + SREPL DROP
»

EVAL

=>

59639380918976775196293204202056658471019738405733913
533679285187391995176580099542319971838.280912881192



PS:

A stack-only version hasn't helped, at least on the emulator, mostly due to my clumsy stack handling:

6 600

« DUP2 3 * 1 + DUP UNROT OVER 3 ^ UNROT ^ 4 PICK ! 1 6 PICK
START DUP2 6 PICK * / 7 ROLL + EXPAND 6 ROLLD 5 PICK / UNROT OVER / ROT 5 ROLL 1 - 5 ROLLD 4 ROLL 3 - 4 ROLLD
NEXT 5 DROPN →STR "/" " " SREPL DROP "'" " " SREPL DROP OBJ→ SWAP 12 ALOG * SWAP IQUOT →STR DUPDUP SIZE DUP 11 - SWAP SUB "." OVER + SREPL DROP
»

EVAL

=>

59639380918976775196293204202056658471019738405733913
533679285187391995176580099542319971838.280912881192

(1293 seconds)

PPS:

It's better to sum the terms in ascending order, otherwise the time-consuming expansions of long expressions are done since the beginning. Execution time down to 769 seconds despite more command inside the loop. Stack handling optimization might make it faster than the first version.

6 601

« OVER 3 ^ ROT 1 0 1 1 6 ROLLD 6 ROLLD 6 ROLLD DUP 8 ROLL
START DUP2 6 PICK * / 7 ROLL + EXPAND 6 ROLLD 5 PICK * UNROT OVER * ROT 5 ROLL 1 + 5 ROLLD 4 ROLL 3 + 4 ROLLD
NEXT 5 DROPN →STR "/" " " SREPL DROP "'" " " SREPL DROP OBJ→ SWAP 12 ALOG * SWAP IQUOT →STR DUPDUP SIZE DUP 11 - SWAP SUB "." OVER + SREPL DROP
»

EVAL

=>

59639380918976775196293204202056658471019738405733913
533679285187391995176580099542319971838.280912881192

The arguments are

Y: upper limit of integration
X: number of terms (n-1, in the previous version)

PPPS:

A version using local variables isn't any faster:

6 601

« 0 UNROT OVER 3 ^ ROT 1 DUPDUP → k b d f n
« 1 SWAP
START b d f * / + EXPAND k 'b' STO* 3 'd' STO+ n 'f' STO* 1 'n' STO+
NEXT
» DUP TYPE 28 - NOT SWAP →STR SWAP { "'" SWAP + "/1'" + } IFT "/" " " SREPL DROP "'" " " SREPL DROP OBJ→ SWAP 12 ALOG * SWAP IQUOT →STR "#" + DUPDUP SIZE DUP 12 - SWAP SUB "." OVER + SREPL DROP "#" "" SREPL DROP
»

EVAL

=>

59639380918976775196293204202056658471019738405733913
533679285187391995176580099542319971838.280912881192

786 seconds on the emulator.

[emece67]

(09-01-2017 11:51 AM)emece67 Wrote:  \(5.963938091897677519629320420205\underline{665}×10^{91}\) in a wp34s after 2'20''.
...

I used for this an wp34s that runs the new double exponential integration method. Can anybody with a wp34s still using the Romberg method repeat this calculation abd report the results?

[Pekis]

Hello,

After u=x^3 substitution, is this the correct form ?

Integ(e^x^3 dx, 0, 6) = Integ(e^u/(3*u^(2/3)), 0, 216)

u=x^3 => dx=du/(3*x^2)
x=u^(1/3) => 3*x^2=3*u^(2/3)
=> dx=du/(3*u^(2/3))

Thanks

[Gerson W. Barbosa]

(09-04-2017 09:30 AM)Pekis Wrote:  Hello,

After u=x^3 substitution, is this the correct form ?

Integ(e^x^3 dx, 0, 6) = Integ(e^u/(3*u^(2/3)), 0, 216)

u=x^3 => dx=du/(3*x^2)
x=u^(1/3) => 3*x^2=3*u^(2/3)
=> dx=du/(3*u^(2/3))

Thanks

Yes, it is.

Elsewhere in this thread I'd wrongly written

Integ(e^x^3 dx, 0, 6) = Integ(e^u/u^(2/3), 0, 216)

neglecting the 1/3 constant in my own hand-written note:

Integ(e^x^3 dx, 0, 6) = 1/3* Integ(e^u/u^(2/3), 0, 216)

Sorry for the lack of attention and thanks for the correction.

Gerson.

[Gamo]

Found another problem to solve for integral this time HP 15C gave error message as error 0 have to do with the 'Complex Results from Real Numbers'

HP Prime gave answer as 3.64899739786 E -2
The answer should be 0.03662

The integral equation attached as jpg.

Here is the program:

LBL A
ENTER
SQR
LSTx
1
-
+
X<>Y
LN
1/x
-
RTN
P/R


Gamo

[Joe Horn]

(09-06-2017 11:15 AM)Gamo Wrote:  Found another problem to solve for integral this time HP 15C gave error message as error 0 have to do with the 'Complex Results from Real Numbers'

HP Prime gave answer as 3.64899739786 E -2
The answer should be 0.03662

Wolfram Alpha's result agrees with Prime: the answer is exactly 2-ɣ-log(4) which is ≈ 0.03648997397857652....

ɣ is the Euler-Mascheroni constant.

[Dieter]

(09-06-2017 11:15 AM)Gamo Wrote:  Found another problem to solve for integral this time HP 15C gave error message as error 0 have to do with the 'Complex Results from Real Numbers'

Where should a complex result appear here?

Hmhmhm... you do not happen to have user mode on, have you?
Check the keycode for the SQRT and 1/x commands. What do you see? 11 and 15 (correct) or 32 11 and 32 15? The latter is GSB A resp. GSB E due to user mode.

(09-06-2017 11:15 AM)Gamo Wrote:  The answer should be 0.03662

And that's what the 15C returns.

(09-06-2017 11:15 AM)Gamo Wrote:  Here is the program:

LBL A
ENTER
SQR
LSTx
1
-
+

Here is another error. This plus should be a divide.

With a correct program you will get the correct result:

FIX 5
0 [ENTER] 1 [∫yx] [A] => 0,03662

Edit: this is not the correct result, for more details see the posts below and the mentioned HP-Journal article.

Dieter

[Didier Lachieze]

(09-06-2017 11:15 AM)Gamo Wrote:  Here is the program:

LBL A
ENTER
SQR
LSTx
1
-
+
X<>Y
LN
1/x
-
RTN
P/R

The ENTER on line 2 is not necessary as the stack is filled with the x value.
On line 7 it should be a division and not an addition.
With the following program:

Code:

001-42.21.11 LBL A
002-      11 SQR
003-      34 X<>Y
004-       1 1
005-      30 -
006-      10 /
007-      34 X<>Y
008-   43 12 LN
009-      15 1/x
010-      30 -
011-   43 32 RTN

FIX 9 0 ENTER 1 f x A returns: 0,036489974 on the HP 15C Windows Emulator.

[Gerson W. Barbosa]

(09-06-2017 11:15 AM)Gamo Wrote:  Found another problem to solve for integral this time HP 15C gave error message as error 0 have to do with the 'Complex Results from Real Numbers'

HP Prime gave answer as 3.64899739786 E -2
The answer should be 0.03662

The integral equation attached as jpg.

Here is the program:

LBL A
ENTER
SQR
LSTx
1
-
+
X<>Y
LN
1/x
-
RTN
P/R


Gamo

The '+' sign in your listing should be a '/'.

More about this integral at

http://www.hpl.hp.com/hpjournal/pdfs/Iss...980-08.pdf

starting at page 23.

Gerson.

[Dieter]

(09-06-2017 12:10 PM)Joe Horn Wrote:  

(09-06-2017 11:15 AM)Gamo Wrote:  ...
The answer should be 0.03662

Wolfram Alpha's result agrees with Prime: the answer is exactly 2-ɣ-log(4) which is ≈ 0.03648997397857652....

FTR: The WP34s integration routine seems to perform better here:

FIX 5:  ∫ ≈ 0,03649
ALL 3:  ∫ ≈ 0,03648997398363982

The latter is at least correct in all displayed digits. ;-)

(09-06-2017 12:29 PM)Gerson W. Barbosa Wrote:  More about this integral at ...

Ah, now I know why this integral looked so familiar! It was mentioned in the HP-Journal article on the then-new 34C and its integrate function (which seems to be more or less the same in the 15C). The indeed wrong answer in FIX 5 mode is explained there, just as the correct result 0,0364900 in FIX 7 mode: "The answers are wrong but the calculator is right". ;-) Again, the error is caused by problems with the numeric evaluation of the integrand which has an very sharp spike near x=0, and it is shown how this can be avoided and the correct result can be obtained.

Dieter

[Didier Lachieze]

(09-06-2017 12:29 PM)Gerson W. Barbosa Wrote:  More about this integral at

http://www.hpl.hp.com/hpjournal/pdfs/Iss...980-08.pdf

starting at page 23.

Gerson.

It's interesting that the program provided in Dr Kahan's article ignores the fact that the stack is filled with the value of the parameter when the function is called.

The first two lines are:

Code:

LBL A   Begin with u in the X register 
ENTER   Save u in the stack

instead of:

Code:

LBL A   Begin with u in the X, Y, Z and T registers

[Dieter]

(09-06-2017 01:02 PM)Didier Lachieze Wrote:  It's interesting that the program provided in Dr Kahan's article ignores the fact that the stack is filled with the value of the parameter when the function is called.

Right, the initial ENTER is not required as Solve and Integrate fill the stack with x. But the ENTER allows to calculate the function for different x (without having to manually fill the stack beforehand) in order to get an impression of its shape and the problems following from this (cf. fig. 10 a...c).

Dieter

[Gamo]

Right that was from HP Journal magazine I try that formula and miss enter the last program from - and suppose to be divide. No wonder it gave out as Error 0

Gamo

[emece67]

(09-06-2017 12:45 PM)Dieter Wrote:  

(09-06-2017 12:10 PM)Joe Horn Wrote:  Wolfram Alpha's result agrees with Prime: the answer is exactly 2-ɣ-log(4) which is ≈ 0.03648997397857652....

FTR: The WP34s integration routine seems to perform better here:

FIX 5:  ∫ ≈ 0,03649
ALL 3:  ∫ ≈ 0,03648997398363982

The latter is at least correct in all displayed digits. ;-)

The last version of the wp34s integration routine returns:
FIX 5: \(\int\approx0.03649\) in 4''
ALL 3: \(\int\approx0.03648997397857\underline{909}\) in 9''
DBLMODE: \(\int\approx0.036489973978576520559023667001244\underline{13}\) in 24''

[Dieter]

(09-06-2017 02:02 PM)emece67 Wrote:  The last version of the wp34s integration routine returns: ...

I must have missed something here. What method was implemented in this new integration routine? And what is the first WP34s version that uses this one instead of the Romberg method?

Dieter

[emece67]

(09-06-2017 08:57 PM)Dieter Wrote:  

(09-06-2017 02:02 PM)emece67 Wrote:  The last version of the wp34s integration routine returns: ...

I must have missed something here. What method was implemented in this new integration routine? And what is the first WP34s version that uses this one instead of the Romberg method?

Dieter

The new method is the double-exponential one. It was uploaded to the svn repository some months ago.

[SlideRule]

ALL

I found the book L’Hôpital's Analyse des infiniments petits …
[attachment=5161]
… to be VERY interesting & EDIFYING with respect to the underlying foundation/basis of the Calculus.

from the preface:
"The Analyse consists of ten chapters, which L’Hôpital called “sections.” We consider it to have three parts. The first part, an introduction to the differential calculus, consists of the first four chapters:
1. In which we give the Rules of this calculus.
2. Use of the differential calculus for finding the Tangents of all kinds of curved lines.
3. Use of the differential calculus for finding the greatest and the least ordinates, to which are reduced questions De maximis & minimis.
4. Use of the differential calculus for finding inflection points and cusps.
Taken together, these chapters provide a thorough introduction to the differential calculus in about 70 pages. The next five chapters are devoted to what can only be described as an advanced text on differential geometry, motivated in part by what were then cutting-edge research problems in optics and other fields. The final chapter is mildly polemical, demonstrating the superiority of Leibniz’ new calculus, when compared to the methods of René Descartes (1596–1660) and Johann van Waveren Hudde (1628–1704)."

BEST!
SlideRule

[johanw]

To revive this old thread: the recently upgraded official HP15C app, version 1.7.3, does not crash on it any more.

In FIX4 I got after 6 minutes 5.9639 * 10^{91} on a reasonably modern phone (Samsung A50).

[Albert Chan]

(09-02-2017 08:55 PM)Gerson W. Barbosa Wrote:  BTW, as an attempt to improve the execution time. I did

u=x^3 -> x=u^(1/3)

du/dx=3*x^2 -> dx=du/(3*x^2)

-> Integ(e^x^3 dx, 0, 6) = 1/3*Integ(e^u/u^(2/3), 0, 216)

A better way might be integrate from the right, with x=u^3, dx=3 u^2 du
The transformation produce a nice bell-shaped curve.

e^x^3 is increasing function, and \(\displaystyle{e^{6^3}\over e^{5^3}} = e^{91} ≈ 3.3×10^{39}\), we can do \(\int_5^6 e^{x^3}\,dx\)

Sharp EL-506W, n=100 (Simpson's Rule, 200 points), time ≈ 14 sec.

∫(e^u/u^(2/3)/3, u=125. .. 216) = 5.965307107 × 10⁹¹
∫(e^(6-u³)³×3u², u=0. .. 1)        = 5.963938092 × 10⁹¹