HP-35s SOLVE with integration

[Rick314]

I want to find the first null in Bessel function J0(x) using an HP-35s. The code below finds Jn(x). But trying to SOLVE it for J0(x) = 0 (at x = 2.4048) results in a "SOLVE ACTIVE" error at the "FN= I" instruction. I finally found in the manual under "Restrictions on Solving and Integrating" that "SOLVE and <integral>FN cannot call a routine that contains an FN=label instruction." This surprised me because older calculators such as the HP-15C can SOLVE an expression involving an integral. Is there any way to SOLVE J0(x) = 0 (at x = 2.4048) using an HP-35s?

Code:

LBL J  Jn(x) = (integral 0 to pi of LBL I)/pi
0      low lim
pi     hi lim
RAD
FN= I
<integrate>FN d T
pi
/
RTN

LBL I  COS(N*T - X*sin(T))
RCL N
RCL T
*
RCL X
RCL T
SIN
*
-
COS
RTN

J0(1): 0 STO N 1 STO X XEQ J ENTER -> 0.7652 (correct)

[Thomas Klemm]

It's been a while when I wrote:

(10-25-2015 01:29 PM)Thomas Klemm Wrote:  Since integration and solver can't be mixed a Newton-iteration is used to calculate the inverse.

Looking at the example in the user's guide I'm afraid it is indeed not possible.

[Steve Simpkin]

(04-22-2024 05:46 PM)Rick314 Wrote:  ...
I finally found in the manual under "Restrictions on Solving and Integrating" that "SOLVE and <integral>FN cannot call a routine that contains an FN=label instruction." This surprised me because older calculators such as the HP-15C can SOLVE an expression involving an integral.
...

True but the HP 35s is more comparable with the calculator line that it evolved from, namely: HP-32S --> HP-32SII --> HP-33S --> HP 35s. The HP-32S was more of a replacement for the HP-11C with a few of the features of the HP-15C (a solver, numeric integration and some complex math) plus alpha display of program steps (instead of key codes) and base conversions. In that respect, I think it was a good evolution of that line.

[Thomas Klemm]

(04-22-2024 05:46 PM)Rick314 Wrote:  Is there any way to SOLVE J0(x) = 0 (at x = 2.4048) using an HP-35s?

Here's a program for the HP-15C that uses regula falsi instead of the built-in solver:

Code:

   001 { 42 21 11 } f LBL A
   002 {    44  3 } STO 3
   003 {       33 } Rv
   004 {    44  2 } STO 2
   005 {    32 12 } GSB B
   006 {    44  4 } STO 4
   007 {    45  3 } RCL 3
   008 {    32 12 } GSB B
   009 {    44  5 } STO 5
   010 { 42 21  0 } f LBL 0
   011 {    45  4 } RCL 4
   012 { 45 20  3 } RCL * 3
   013 {    45  5 } RCL 5
   014 { 45 20  2 } RCL * 2
   015 {       30 } -
   016 {    45  4 } RCL 4
   017 { 45 30  5 } RCL - 5
   018 {       10 } /
   019 {    32 12 } GSB B
   020 {    43 34 } g RND
   021 {    43 20 } g x==0
   022 {    22  2 } GTO 2
   023 {    43 36 } g LSTx
   024 { 45 20  4 } RCL * 4
   025 { 43 30  1 } g TEST 1
   026 {    22  1 } GTO 1
   027 {    43 36 } g LSTx
   028 {    44  5 } STO 5
   029 {    45  1 } RCL 1
   030 {    44  3 } STO 3
   031 {    22  0 } GTO 0
   032 { 42 21  1 } f LBL 1
   033 {    43 36 } g LSTx
   034 {    44  4 } STO 4
   035 {    45  1 } RCL 1
   036 {    44  2 } STO 2
   037 {    22  0 } GTO 0
   038 { 42 21  2 } f LBL 2
   039 {    45  1 } RCL 1
   040 {    43 32 } g RTN
   041 { 42 21 12 } f LBL B
   042 {    44  1 } STO 1
   043 {        0 } 0
   044 {    43 26 } g PI
   045 { 42 20 13 } f integrate C
   046 {    43 32 } g RTN
   047 { 42 21 13 } f LBL C
   048 { 45 20  0 } RCL * 0
   049 {       34 } X<=>Y
   050 {       23 } SIN
   051 { 45 20  1 } RCL * 1
   052 {       30 } -
   053 {       24 } COS
   054 {    43 32 } g RTN

Formula

\(
\begin{align}
x = \frac{b \cdot f(a) - a \cdot f(b)}{f(a) - f(b)}
\end{align}
\)

Registers

0: n
1: x
2: a
3: b
4: f(a)
5: f(b)

Example

FIX 7
RAD

0 STO 0

2 ENTER 3
GSB A

2.4048256

But I'm sure that you are able to rewrite it for the HP-35S.

[Rick314]

I have a program on my HP-35s to solve the problem in the original post. Thank you to those that replied, especially Thomas Klemm -- I recognized your HP-15C program as being compatible with the JRPN 15C simulator, ran it, and enjoyed commenting it (below) to understand it. Notably:

• Solving f(x)*constant = 0 gives the same result as solving f(x) = 0.
• Register arithmetic.
• X intercept of a line through (x1, y1) and (x2, y2).
• Regula falsi, that I have used but did not know its name or history.
  

Code:

# Thomas Klemm, https://www.hpmuseum.org/forum/thread-21638.html
# Registers:
# 0: n     Bessel order
# 1: x
# 2: a     = x1
# 3: b     = x2
# 4: f(a)  = y1 = Jn(x1)*pi
# 5: f(b)  = y2 = Jn(x2)*pi
# x0 = (b*f(a)-a*f(b))/(f(a)-f(b)) = (x2*y1 - x1*y2)/(y1 - y2)
#    = where a line through (x1, y1), (x2, y2) crosses the x axis.
# A null in Jn(x) is also a null in Jn(x)*pi, what is found below.

   001 { 42 21 11 } f LBL A   # main
   002 {    44  3 } STO 3     # b
   003 {       33 } Rv
   004 {    44  2 } STO 2     # a
   005 {    32 12 } GSB B
   006 {    44  4 } STO 4     # f(a)
   007 {    45  3 } RCL 3
   008 {    32 12 } GSB B
   009 {    44  5 } STO 5     # f(b)

   010 { 42 21  0 } f LBL 0   # find x0
   011 {    45  4 } RCL 4     # f(a)
   012 { 45 20  3 } RCL * 3   # b*f(a)
   013 {    45  5 } RCL 5     # f(b)
   014 { 45 20  2 } RCL * 2   # a*f(b)
   015 {       30 } -         # b*f(a)-a*f(b)
   016 {    45  4 } RCL 4     # f(a)
   017 { 45 30  5 } RCL - 5   # f(a)-f(b)
   018 {       10 } /         # (b*f(a)-a*f(b))/(f(a)-f(b)) = x0
   019 {    32 12 } GSB B     # find Jn(x0)*pi
   020 {    43 34 } g RND     # round to display resolution
   021 {    43 20 } g x==0
   022 {    22  2 } GTO 2     # if f(x0) == 0, DONE
   023 {    43 36 } g LSTx    # full-resolution f(x0), nearing 0
   024 { 45 20  4 } RCL * 4   # f(a)*f(x0), x0 going to a or b?
   025 { 43 30  1 } g TEST 1  # x>0?
   026 {    22  1 } GTO 1     # if x0->a, branch ahead,
   027 {    43 36 } g LSTx    # else x0->b...
   028 {    44  5 } STO 5     # f(x0) -> f(b)
   029 {    45  1 } RCL 1     # x0
   030 {    44  3 } STO 3     # x0 -> b
   031 {    22  0 } GTO 0     # branch up for another iteration

   032 { 42 21  1 } f LBL 1
   033 {    43 36 } g LSTx    # f(x0)
   034 {    44  4 } STO 4     # f(x0) -> f(a)
   035 {    45  1 } RCL 1     # x0
   036 {    44  2 } STO 2     # x0 -> a
   037 {    22  0 } GTO 0     # branch up for another iteration

   038 { 42 21  2 } f LBL 2   # DONE
   039 {    45  1 } RCL 1     # x0 final value
   040 {    43 32 } g RTN     # DONE

   041 { 42 21 12 } f LBL B   # integrate f(t) 0 to pi dt
   042 {    44  1 } STO 1     # x
   043 {        0 } 0
   044 {    43 26 } g PI
   045 { 42 20 13 } f integrate C
   046 {    43 32 } g RTN

   047 { 42 21 13 } f LBL C   # f(t), Ry=t, Rx=t
   048 { 45 20  0 } RCL * 0   # n*t
   049 {       34 } X<=>Y     # t
   050 {       23 } SIN       # sin(t)
   051 { 45 20  1 } RCL * 1   # x*sin(t)
   052 {       30 } -         # n*t - sin(t)
   053 {       24 } COS       # cos(n*t - sin(t))
   054 {    43 32 } g RTN

# End.

[Thomas Klemm]

(04-25-2024 06:41 PM)Rick314 Wrote:  enjoyed commenting it (below) to understand it.

I'm pleased that you now have a working program for your HP-35s.
And thank you for your excellent comments, which make it much easier to understand.
I often do this myself, but in this case it worked to your advantage that I was lazy.