HP-35s SOLVE with integration - Rick314 - 04-22-2024 05:46 PM
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)
RE: HP-35s SOLVE with integration - Thomas Klemm - 04-22-2024 10:48 PM
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.
RE: HP-35s SOLVE with integration - Steve Simpkin - 04-23-2024 12:06 AM
(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.
RE: HP-35s SOLVE with integration - Thomas Klemm - 04-23-2024 12:44 AM
(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.
RE: HP-35s SOLVE with integration - Rick314 - 04-25-2024 06:41 PM
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.
RE: HP-35s SOLVE with integration - Thomas Klemm - 04-26-2024 05:19 PM
(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.
|