(35S) Quick integration

[Roberto Volpi]

Hi all,

the native integrator of our HP35S works fine, but it is not as quick as we would like.

Inputing functions in RPN mode instead of EQN makes us save approx 30% of time, but waiting 5 minutes for a result is still possible, and here a short programme to avoid it.

The algorithm is based upon Simpson's Rule, with n=2. So just 3 sampling points are taken, but the value we obtain can be an acceptable approximation most of the times, and it is really quick, as the f(x) to be integrated needs to be evaluated just 3 times.

After all, in case f(x) is a polynomial of degree 3 or lower, we obtain an exact value.

NOTE: it assumes that there is a LBL F to be used to store the f(x) to be analysed.

LBL I
STO B
X<>Y
STO A
STO X
XEQ F001
RCL B
STO X
XEQ F
X<>Y
R down
+
X<>Y
RCL B
2
/
STO X
XEQ F001
4
X
R up
+
RCL B
RCL -A
6
/
*
RTN


INSTRUCTIONS:

- Input the integrand f(x) in LBL F

- As with the native integrator, input "a" value, press ENTER, and input "b" value

- press XEQ I ENTER

we quickly obtain:

stack y: midpoint between a and b
stack x: numerical integration


Enjoy!

[PedroLeiva]

An example please, just to check if I am doing things well
Pedro

[Roberto Volpi]

Dear Pedro,

here a short and quick example.

f(x)=∫√(1+4X^2)dx in interval [0,1]

which can be interpreted as the arc length of the parabola y=x^2 in interval [0,1]

With the native HP35S integrator, after approx 1 minute and 37 seconds:
=1.47894285754

With this program, after 1 second:
=
stack y: 0.5
stack x: 1.48215537045

Please take into account that also the value of the native HP35S integrator is approximate.

[PedroLeiva]

I will try. Thanks a lot, Pedro

[Albert Chan]

With only 3 sample points, I would use Gaussian quadrature.
It is more accurate (exact for quintic polynomial or less), and does not touch end-points.

Code is slightly more complicated, but I think it is worth it.

Code:

def gauss3(f,a,b, w=0.6**0.5):
    c, d = (b+a)/2, (b-a)/2
    g = d*w
    return (f(c-g) + f(c+g) + 1.6*f(c))/1.8 * d

Your example, ∫(√(1+4*x^2), x = 0..1) ≈ 1.47894285754

Simpson  3-pts: 1.48215537045, relative error ≈ -0.217%
Gaussian 3-pts: 1.47910755366, relative error ≈ -0.011%

Update: integral formula corrected.

[Thomas Klemm]

This program for the HP-15C uses a 3-point Gaussian quadrature:

Code:

001 {    42 21 11 } f LBL A
002 {       44  1 } STO 1
003 {          30 } −
004 {           2 } 2
005 {          10 } ÷
006 {       44  0 } STO 0
007 {    44 40  1 } STO + 1
008 {          48 } .
009 {           6 } 6
010 {          11 } √x̅
011 {          20 } ×
012 {       44  2 } STO 2
013 {    45 40  1 } RCL + 1
014 {       32 15 } GSB E
015 {           5 } 5
016 {          20 } ×
017 {       44  3 } STO 3
018 {       45  1 } RCL 1
019 {    45 30  2 } RCL − 2
020 {       32 15 } GSB E
021 {           5 } 5
022 {          20 } ×
023 {    44 40  3 } STO + 3
024 {       45  1 } RCL 1
025 {       32 15 } GSB E
026 {           8 } 8
027 {          20 } ×
028 {    45 40  3 } RCL + 3
029 {    45 20  0 } RCL × 0
030 {           9 } 9
031 {          16 } CHS
032 {          10 } ÷
033 {       43 32 } g RTN
034 {    42 21 15 } f LBL E
035 {          23 } SIN
036 {       43 36 } g LSTx
037 {          10 } ÷
038 {       43 32 } g RTN

It yields an exact result for polynomials of degree \(2n − 1\) (in this case \(5\)) or less.

Example

\(
\begin{align}
\int_0^2 \frac{\sin(x)}{x} \, dx
\end{align}
\)

0 ENTER 2
GSB A

1.605418622

The correct value is:
1.605412977

It should be easy to translate the program for the HP-35S.

References
Inspired by: (11C) Gaussian integration

[Roberto Volpi]

I am pretty sure there are plenty of algorithms out there for numerical integration, just waiting to be mentioned.

And, if other calculators, or more powerful machines, come into discussion, this list could be endless.

Any proposal for the HP35S?

It would really be a nice challenge make a shorter and hence quicker programme.

[rawi]

Roberto Volpi wrote:

Quote:Inputing functions in RPN mode instead of EQN makes us save approx 30% of time, but waiting 5 minutes for a result is still possible, and here a short programme to avoid it.

Could you give an example or a hint to a thread how to use a function in RPN mode? In the manual I did not find any information about that.
Thank you very much!

[Roberto Volpi]

(11-13-2022 05:27 AM)rawi Wrote:  Roberto Volpi wrote:

Quote:Inputing functions in RPN mode instead of EQN makes us save approx 30% of time, but waiting 5 minutes for a result is still possible, and here a short programme to avoid it.

Could you give an example or a hint to a thread how to use a function in RPN mode? In the manual I did not find any information about that.
Thank you very much!

Well, it is essencially inputing the same operations in a programme label (ex. F) in RPN mode.

As I understand not every participant in this forum is an engineer with 30 years RPN practice, let's give a quick example.

Given f(x)=∫√(1+4X^2), as in the given example, we could input in RPN mode as following.

1) using X variable, we could do:

LBL F
RCL X
squared
4
*
1
+
squared root
RTN

every time this subroutine is executed, the f(x) is evaluated according to the content if the independent variable X, corresponding to X register.

2) just using x stack value:

LBL F
squared
4
*
1
+
squared root
RTN

[rawi]

Hi Roberto,

thank you very much for your fast answer.
I know how to code functions in RPN but I do not know how to tell the intergration function of the HP 35s to use the program F for integration instead of an EQN-equation. If I put in lower limit Enter upper limit and then integration it asks for the equation to integrate. Is there a way to tell the calculator to use a program instead of an equation put in in EQN? This would be a little more like with the HP 15 or the HP 42s.
Thanks once more and best regards
Raimund

[Roberto Volpi]

(11-13-2022 12:49 PM)rawi Wrote:  Hi Roberto,

thank you very much for your fast answer.
I know how to code functions in RPN but I do not know how to tell the intergration function of the HP 35s to use the program F for integration instead of an EQN-equation. If I put in lower limit Enter upper limit and then integration it asks for the equation to integrate. Is there a way to tell the calculator to use a program instead of an equation put in in EQN? This would be a little more like with the HP 15 or the HP 42s.
Thanks once more and best regards
Raimund

In the HP35S when you want to solve or integrate a function other than an equation the "FN=" instruction must be used (up left, in yellow).

So, in case you have a function in, say LBL F, before solving or integrating just press FN=, it will prompt a letter, press F, so the HP35S knows that it will use LBL F when asked to solve or integrate.

[Roberto Volpi]

Taking advantage of the suggestion of Thomas Klemm, here a version for the HP35S of a programme for quick integral calculation using as algorithm Gaussian integration.

A few keystrokes are necessary but overall it is an advancement. It is in fact little more accurate and permits to evaluate integrals like that given in the example (∫(sinX/X) dx in interval [0,2] that may not be evaluated by the native HP35S integrator).

It assumes that LBL F contains the f(x) to be integrated in RPN mode without using X register)

NOTE: IT ALREADY INCLUDES ERRATA POINTED OUT BY THOMAS KLEMM.

LBL I
STO B
-
2
/
STO A
STO+B
.6
sqr root
*
STO C
RCL +B
XEQ F001
5
*
STO D
RCL B
RCL -C
XEQ F001
5
*
STO +D
RCL B
XEQ F001
8
*
RCL +D
RCL *A
-9
/
RTN

The integral in the examble should be

LBL F
SIN
last x
/
RTN

INSTRUCTIONS:

- Input lower and upper integral value as usual (in this case, press "0" "ENTER" "2"

- Press "XEQ I ENTER"

- result: 1.6054186216


Enjoy!

[Thomas Klemm]

It appears that the command on line 2 got lost in translation:
STO B

Also RCL +A should rather be:
RCL *A

(11-14-2022 05:05 AM)Roberto Volpi Wrote:  ∫(sinX/X) dx in interval [0,2] that may not be evaluated by the native HP35S integrator

From the HP 35s scientific calculator user's guide (p. 8-5):

Quote:If the calculator attempted to evaluate this function at x = 0, the lower limit of integration, an error (DIVIDE BY 0) would result. However, the integration algorithm normally does not evaluate functions at either limit of integration, unless the endpoints of the interval of integration are extremely close together or the number of sample points is extremely large.

The given result is:
1.6054

---

(11-12-2022 05:13 PM)Albert Chan Wrote:  Your example, ∫(√(1+x^2), x = 0..1) ≈ 1.47894285754

That should be:
∫(√(1+4x^2), x = 0..1) ≈ 1.47894285754

[Roberto Volpi]

(11-14-2022 07:39 AM)Thomas Klemm Wrote:  It appears that the command on line 2 got lost in translation:
STO B

Also RCL +A should rather be:
RCL *A

(11-14-2022 05:05 AM)Roberto Volpi Wrote:  ∫(sinX/X) dx in interval [0,2] that may not be evaluated by the native HP35S integrator

From the HP 35s scientific calculator user's guide (p. 8-5):

Quote:If the calculator attempted to evaluate this function at x = 0, the lower limit of integration, an error (DIVIDE BY 0) would result. However, the integration algorithm normally does not evaluate functions at either limit of integration, unless the endpoints of the interval of integration are extremely close together or the number of sample points is extremely large.

The given result is:
1.6054

---

(11-12-2022 05:13 PM)Albert Chan Wrote:  Your example, ∫(√(1+x^2), x = 0..1) ≈ 1.47894285754

That should be:
∫(√(1+4x^2), x = 0..1) ≈ 1.47894285754

Thanks Thomas,

I have already amended the programme as per your instruction.

I love this forum...

[Thomas Klemm]

(11-12-2022 05:16 PM)Thomas Klemm Wrote:  It yields an exact result for polynomials of degree \(2n − 1\) (in this case \(5\)) or less.

This made me wonder if we can trick the Gaussian quadrature to make the integral 0 when we make the function 0 at the selected values.
But then the 5th degree polynomial must be of the form:

\(
\begin{align}
p(x) = a (5x^3 - 3x)(x^2 + px + q)
\end{align}
\)

In fact, it follows that:

\(
\begin{align}
\int_{-1}^{1} p(x) \, dx
&= a \int_{-1}^{1} (5x^3 - 3x)(x^2 + px + q) \, dx \\
&= a \int_{-1}^{1} 5 x^5 + 5 p x^4 + (5 q - 3) x^3 - 3 p x^2 - 3 q x \, dx \\
&= a \left[ \left.\tfrac{5}{6}x^6 + p x^5 + \frac{5 q - 3}{4} x^4 - p x^3 - \tfrac{3}{2} q x^2 \right|_{-1}^{1} \; \right] \\
&= a [2p - 2p] =0
\end{align}
\)

In retrospect, this is obvious since \(p_3(x) = \tfrac{1}{2}\left(5x^{3} - 3x \right) \) is orthogonal to \(x^k\) for all \(k \in \{0, 1, 2\}\). (Cf. Fundamental theorem)
But at first it wasn't to me.

[Albert Chan]

We can take advantage of odd/even function:

∫((5*x^3-3*x) * (x^2+p*x+q), x=-1 .. 1)      // odd function * even function = odd function
= ∫((5*x^3-3*x) * (p*x), x=-1 .. 1)               // odd function * odd function = even function
= 2p * ∫(5*x^4 - 3*x^2, x=0 .. 1)       
= 2p * (5/5 - 3/3)
= 0

This explained coefficients (5, -3) for Gaussian quadrature 3-points abscissa

(5*x^3-3*x) = 5*x*(x^2-0.6) = 0      → x = 0, ±√(0.6)

[J-F Garnier]

(11-12-2022 03:56 PM)Roberto Volpi Wrote:  here a short and quick example.

f(x)=∫√(1+4X^2)dx in interval [0,1]

With the native HP35S integrator, after approx 1 minute and 37 seconds:
=1.47894285754

With this program, after 1 second:
=
stack y: 0.5
stack x: 1.48215537045

Please take into account that also the value of the native HP35S integrator is approximate.

Your comparison with the HP-35S integration full 12-digit result is not fair.
On my HP-35S using the built-in integration and the default FIX 4 display, I get 1.4789 in about 3s.
If you need even less accuracy, you can use FIX 1 or FIX 2, or better SCI 1 or SCI 2 to specify a relative target accuracy.

(11-12-2022 05:29 PM)Roberto Volpi Wrote:  Any proposal for the HP35S?

It would really be a nice challenge make a shorter and hence quicker programme.

So here is mine (enter the integration limits on the stack):
LBL I
SCI 2 ; set target accuracy
FN= F ; select function
∫FN d X ; compute the integral
FIX 4 ; restore your preferred display mode
RTN

J-F

[Roberto Volpi]

(11-14-2022 02:50 PM)J-F Garnier Wrote:  

(11-12-2022 03:56 PM)Roberto Volpi Wrote:  here a short and quick example.

f(x)=∫√(1+4X^2)dx in interval [0,1]

With the native HP35S integrator, after approx 1 minute and 37 seconds:
=1.47894285754

With this program, after 1 second:
=
stack y: 0.5
stack x: 1.48215537045

Please take into account that also the value of the native HP35S integrator is approximate.

Your comparison with the HP-35S integration full 12-digit result is not fair.
On my HP-35S using the built-in integration and the default FIX 4 display, I get 1.4789 in about 3s.
If you need even less accuracy, you can use FIX 1 or FIX 2, or better SCI 1 or SCI 2 to specify a relative target accuracy.

(11-12-2022 05:29 PM)Roberto Volpi Wrote:  Any proposal for the HP35S?

It would really be a nice challenge make a shorter and hence quicker programme.

So here is mine (enter the integration limits on the stack):
LBL I
SCI 2 ; set target accuracy
FN= F ; select function
∫FN d X ; compute the integral
FIX 4 ; restore your preferred display mode
RTN

J-F

Your proposal is pretty interesting.

In this way we don't even need to upload a programme: we just fix 3 or 4 decimal figures accuracy and the native HP35S integrator become much much faster.

As for me, the suggestion of Thomas Klemm is still valuable because it permits numerical calculation as those given in the example, with f(x)= sinX/X interval [0,2] that native integrator may not process.

We are putting together lot of interesting details...

[Thomas Klemm]

(11-12-2022 05:29 PM)Roberto Volpi Wrote:  It would really be a nice challenge make a shorter and hence quicker programme.

We could also just calculate the integral:

\(
\begin{align}
\int_0^1 \sqrt{1 + 4 x^2} \, dx
&= \frac{2 \sqrt{5} + \sinh^{-1}(2)}{4} \\
\\
&\approx 1.47894285754460… \\
\end{align}
\)

Code:

001 {           2 } 2
002 {           0 } 0
003 {          11 } √x̅
004 {           2 } 2
005 {    43 22 23 } g HYP⁻¹ SIN
006 {          40 } +
007 {           4 } 4
008 {          10 } ÷

R/S

1.478942858

[Roberto Volpi]

(11-15-2022 09:49 AM)Thomas Klemm Wrote:  

(11-12-2022 05:29 PM)Roberto Volpi Wrote:  It would really be a nice challenge make a shorter and hence quicker programme.

We could also just calculate the integral:

\(
\begin{align}
\int_0^1 \sqrt{1 + 4 x^2} \, dx
&= \frac{2 \sqrt{5} + \sinh^{-1}(2)}{4} \\
\\
&\approx 1.47894285754460… \\
\end{align}
\)

Code:

001 {           2 } 2
002 {           0 } 0
003 {          11 } √x̅
004 {           2 } 2
005 {    43 22 23 } g HYP⁻¹ SIN
006 {          40 } +
007 {           4 } 4
008 {          10 } ÷

R/S

1.478942858

Well, one step forward and no calculator is needed at all ...