finding the minimmum/maximum of a function using Newton's method.
Memory Map
R0 = tolerance
R1 = x
R2 = h
R3 = f(x), f''(x)
R4 = f(x+h)
R5 = f(x-h)
HP-67 Implementation
Code:
1 LBL A
2 STO 1 # store guess for optimum x
3 X<>Y
4 STO 0 # store tolerance
5 LBL 0
6 RCL 1 # display current guess for optimum x
7 PAUSE
8 ABS
9 1
10 +
11 EEX
12 CHS
13 3
14 *
15 STO 2 # calculate and store h
16 RCL 1
17 GSB E
18 STO 3 # Calculate and store f(x)
19 RCL 1
20 RCL 2
21 +
22 GSB E
23 STO 4 # Calculate and store f(x+h)
24 RCL 1
25 RCL 2
26 -
27 GSB E
28 STO 5 # Calculate and store f(x-h)
29 RCL 4
30 +
31 RCL 3
32 2
33 *
34 -
35 RCL 2
36 X^2
37 /
38 STO 3 # calculate and store f''(x)
39 RCL 4
40 RCL 5
41 -
42 RCL 2
43 /
44 2
45 / # calculate f'(x)
46 RCL 3
47 / # calculate diff
48 STO- 1 # x = x - diff
49 ABS
50 RCL 0
51 X<=Y?
52 GTO 0
53 RCL 1
54 RTN
55 LBL E
56 EXP
57 LASTX
58 X^2
59 3
60 *
61 -
62 RTN
Usage
LBL E has the code for the function whose minimum or maximum you seek.
1) Enter the tolerance and press [ENTER].
2) Enter a guess for the minimum/maximum and pres [A].
3) The program pauses showing intermediate values for the refined min/max values. When the program has reached the tolerance level it displays the value of the X for the min/max value.
Thank you for sharing.
A pity the listings don't come with key codes. This would allow me to copy/paste the programs right into my HP-67/97 emulators for the iPad...
(05-26-2014 12:14 PM)Willy R. Kunz Wrote: [ -> ]A pity the listings don't come with key codes.
Code:
001 31 25 11 LBL A
002 33 01 STO 1 # store guess for optimum x
003 35 52 X<>Y
004 33 00 STO 0 # store tolerance
005 31 25 00 LBL 0
006 34 01 RCL 1 # display current guess for optimum x
007 35 72 PAUSE
008 35 64 ABS
009 01 1
010 61 +
011 43 EEX
012 42 CHS
013 03 3
014 71 *
015 33 02 STO 2 # calculate and store h
016 34 01 RCL 1
017 31 22 15 GSB E
018 33 03 STO 3 # Calculate and store f(x)
019 34 01 RCL 1
020 34 02 RCL 2
021 61 +
022 31 22 15 GSB E
023 33 04 STO 4 # Calculate and store f(x+h)
024 34 01 RCL 1
025 34 02 RCL 2
026 51 -
027 31 22 15 GSB E
028 33 05 STO 5 # Calculate and store f(x-h)
029 34 04 RCL 4
030 61 +
031 34 03 RCL 3
032 02 2
033 71 *
034 51 -
035 34 02 RCL 2
036 32 54 X^2
037 81 /
038 33 03 STO 3 # calculate and store f''(x)
039 34 04 RCL 4
040 34 05 RCL 5
041 51 -
042 34 02 RCL 2
043 81 /
044 02 2
045 81 / # calculate f'(x)
046 34 03 RCL 3
047 81 / # calculate diff
048 33 51 01 STO- 1 # x = x - diff
049 35 64 ABS
050 34 00 RCL 0
051 32 71 X<=Y?
052 22 00 GTO 0
053 34 01 RCL 1
054 35 22 RTN
055 31 25 15 LBL E
056 32 52 EXP
057 35 82 LASTX
058 32 54 X^2
059 03 3
060 71 *
061 51 -
062 35 22 RTN
Cheers
Thomas
Namir - Nice program
Any reason to take out and use my HP-67 is a good one.
Thanks for that!
(05-26-2014 10:13 PM)Thomas Klemm Wrote: [ -> ]
Code:
001 31 25 11 LBL A
002 33 01 STO 1 # store guess for optimum x
003 35 52 X<>Y
004 33 00 STO 0 # store tolerance
...
Cheers
Thomas
Thanks, Thomas. Very kind of you.
This also prompted me to have a hard look at RPN-67/97's import capabilities. The result being that in the next update the emulators will accept listings without key codes. In fact, you can now enter and edit programs by simply typing calculator function names in the listing area.
All of Namir's programs shown here run without any editing after pasting.
UPDATE:
RPN-67 Pro v2.5 and RPN-97 Pro v1.5 are now available, featuring the new import capabilities.
Also in these updates: for the first time ever, real breakpoints to help you debug your program.
(05-26-2014 10:13 PM)Thomas Klemm Wrote: [ -> ][quote='Willy R. Kunz' pid='12139' dateline='1401106474']
A pity the listings don't come with key codes.
Dear Thomas:
I am using HP 67 SD, I got for Tolerance=5x10-4 and a guessing initial number of 9, a final figure of 2,833144107. Is this correct ?
Now some basic questions please:
1- What equation represents steps 55-62? How I entry an equation to be evaluated?
2- From my example 2,833 is a Max. or a min.? How to obtain both if corresponds?
Thank you in advance, Pedro
(05-21-2016 10:50 PM)PedroLeiva Wrote: [ -> ]1- What equation represents steps 55-62? How I entry an equation to be evaluated?
As usual, you enter f(x) at LBL E. The argument x is expected in the X-register. Here f(x) obviously is e
x–3x².
(05-21-2016 10:50 PM)PedroLeiva Wrote: [ -> ]2- From my example 2,833 is a Max. or a min.? How to obtain both if corresponds?
Since f"(x) is stored in R3, you could simply inpect this value (RCL 3) an see whether it's positive (minimum) or negative (maximum).
Here the function has a minimum at x=2,833147892... (solution of f'(x) = e
x–6x = 0).
At ths point f"(x) = e
x–6 is 10,998887... > 0, so this is a minimum. The content of R3 should be close to the latter value.
Dieter
(05-22-2016 01:11 AM)Dieter Wrote: [ -> ] (05-21-2016 10:50 PM)PedroLeiva Wrote: [ -> ]1- What equation represents steps 55-62? How I entry an equation to be evaluated?
As usual, you enter f(x) at LBL E. The argument x is expected in the X-register. Here f(x) obviously is ex–3x².
(05-21-2016 10:50 PM)PedroLeiva Wrote: [ -> ]2- From my example 2,833 is a Max. or a min.? How to obtain both if corresponds?
Since f"(x) is stored in R3, you could simply inpect this value (RCL 3) an see whether it's positive (minimum) or negative (maximum).
Here the function has a minimum at x=2,833147892... (solution of f'(x) = ex–6x = 0).
At ths point f"(x) = ex–6 is 10,998887... > 0, so this is a minimum. The content of R3 should be close to the latter value.
Dieter
Dieter
It's more complicated than I imagined. Luckily I made the query. Thank you very much, Peter
(05-22-2016 01:26 AM)PedroLeiva Wrote: [ -> ]It's more complicated than I imagined. Luckily I made the query.
Hmmm... honestly, I do not quite understand the problem.
You simply enter f(x) as a short routine just as you do with any other "classic" HP, be it the 65, the 41, the 67/97 or most others. The argument x can be expected in the X-register, so you simply type what you'd also do in a manual calculation. In this example f(x) is e
x–3x², so it's [e
x] [LstX] [x²] [3] [x] [–]. This is what you enter at LBL E so that f(x) can be calculated by a simple [E] or GSB E.
Checking for a minimum or maximum is easily done by RCL 3 after the program has finished. If the result is positive it's a minimum, and if it's negative there's a maximum. A value of (or very close to) zero would indicate a possible inflection point, but this should not happen as the program divides by R3 so that it would have stopped with an error before.
You could add these lines:
Code:
... ...
053 RCL 3
054 ENTER
055 ABS
056 /
057 CHS
058 RCL 1
059 RTN
This calculates –sign(f"(x)). So after the program has finished, a simple X<>Y will show –1 (minimum) or +1 (maximum).
Dieter
This calculates –sign(f"(x)). So after the program has finished, a simple X<>Y will show –1 (minimum) or +1 (maximum).
Dieter
[/quote]
The point is that in the original post of this program was no explanation about how to determine Max. or min. (or I did not realize).
From the two options for cheking, I prefer: x<>y (+1, -1), for Max. & min., respectively. This one is more intuitive (+ is related to Max.), the other procedure is inverse.
Thank you for clarifying this math issue, and modifying the program
Pedro
(05-22-2016 10:35 PM)PedroLeiva Wrote: [ -> ]The point is that in the original post of this program was no explanation about how to determine Max. or min. (or I did not realize).
The program implements a simple Newton iteration to find the zeroes
of the derivative f'(x). In a minimum or maximum f'(x) must be zero, and this is what the program finds. It evaluates x_new = x – f'(x) / f"(x).
(05-22-2016 10:35 PM)PedroLeiva Wrote: [ -> ]From the two options for cheking, I prefer: x<>y (+1, -1), for Max. & min., respectively.
This is what the proposed code does. It returns +1 for a maximum and –1 for a minimum.
Dieter
Examples:
Equation (began step 061) ex–3x²
Tolerance= 1x10-5
Initial guess: 9……………X= 2,83314411………..min…………..Y= -7,08129358
Initial guess: 0,1…………X= 0,20448151……….Max…………..Y= 1,10145071
For X=0……………………Y= 1
For Y=0……………………X= 0,91
For the original equation, input TOLERANCE [ENTER] 9 [A] you get 2,8331, [x<>y] -1 so you know that is a minimun. Later, [x<>y] [E] you get Y=-7.081 wich is the Y value for the min.
When X=0 [E] you get Y=1, that is to say intercept
The value of X when Y=0 was solved by iteraction
Follow the same procedure for the Maximun
Pedro
(05-26-2014 10:13 PM)Thomas Klemm Wrote: [ -> ] (05-26-2014 12:14 PM)Willy R. Kunz Wrote: [ -> ]A pity the listings don't come with key codes.
Code:
001 31 25 11 LBL A
002 33 01 STO 1 # store guess for optimum x
003 35 52 X<>Y
004 33 00 STO 0 # store tolerance
005 31 25 00 LBL 0
006 34 01 RCL 1 # display current guess for optimum x
007 35 72 PAUSE
008 35 64 ABS
009 01 1
010 61 +
011 43 EEX
012 42 CHS
013 03 3
014 71 *
015 33 02 STO 2 # calculate and store h
016 34 01 RCL 1
017 31 22 15 GSB E
018 33 03 STO 3 # Calculate and store f(x)
019 34 01 RCL 1
020 34 02 RCL 2
021 61 +
022 31 22 15 GSB E
023 33 04 STO 4 # Calculate and store f(x+h)
024 34 01 RCL 1
025 34 02 RCL 2
026 51 -
027 31 22 15 GSB E
028 33 05 STO 5 # Calculate and store f(x-h)
029 34 04 RCL 4
030 61 +
031 34 03 RCL 3
032 02 2
033 71 *
034 51 -
035 34 02 RCL 2
036 32 54 X^2
037 81 /
038 33 03 STO 3 # calculate and store f''(x)
039 34 04 RCL 4
040 34 05 RCL 5
041 51 -
042 34 02 RCL 2
043 81 /
044 02 2
045 81 / # calculate f'(x)
046 34 03 RCL 3
047 81 / # calculate diff
048 33 51 01 STO- 1 # x = x - diff
049 35 64 ABS
050 34 00 RCL 0
051 32 71 X<=Y?
052 22 00 GTO 0
053 34 01 RCL 1
054 35 22 RTN
055 31 25 15 LBL E
056 32 52 EXP
057 35 82 LASTX
058 32 54 X^2
059 03 3
060 71 *
061 51 -
062 35 22 RTN
Cheers
Thomas
Thanks Thomas! I was able to copy and paste to my iPhone's RPN-67 SD and it works like a charm!
(05-28-2016 08:09 PM)bshoring Wrote: [ -> ]I found a pretty good article about the concept at
http://www.themathpage.com/acalc/max.htm
It helped me to understand what's happening with this program.
Excellent webpage, also for the others 14th. lessons. Thank´s for sharing
Pedro