Showing posts with label numerical derivative. Show all posts
Showing posts with label numerical derivative. Show all posts

Saturday, June 26, 2021

7000G Retro Month - June 26 Edition

7000G Retro Month - June 26 Edition





Introduction


Welcome to the 7000G Retro Month, which features programming for the classic Casio calculators from the mid/late 1980s:  primarily fx-7000G and fx-7500G.  Since the programming language stays similar throughout the years, programs can be translated to the fx-6300G and later graphing calculators with little to no adjustments.  Non graphic programs should be ported to the fx-4000P, fx-4500P (A), fx-3650p (II), fx-50F Plus (II), and fx-5800P with little to no adjustments.  


7000G Retro Month takes place every Saturday during June 2021.


To make text easier to type, I can going to use the following text friendly symbols for the following:


->  for →


/I for ⊿


=> for ⇒


What do you think?   Unicode or simple text equivalents?  


- - - - - - -- - -- - -


Today's subject revolves around Calculus.  Enjoy!


- - - -- - - -- - -- -- -


The three programs listed here call another program as a subroutine.  I use Prog 0 as a subroutine.  


Prog 0 


[insert f(x) here]



Example:  Prog 0 contains X^2+4.   The results gets stored in in Ans.  


Sum


S = ∑( f(x), X = A to B)


"A"? -> A

"B"? -> B

A -> X

0 -> S

Lbl 1

Prog 0

Ans + S -> S

X + 1 -> X

X≤B => Goto 1

S


Numeric Derivative - Simple Approximation


f'(x) ≈ (f(x+h) - f(x-h))/(2h),  h = tolerance (default to 10^-5)


The derivative is stored in the variable D. 


"X0"? -> Z

Z+10^-5 -> X : Prog 0 : Ans -> A

Z-10^-5 -> X : Prog 0 : Ans -> B

(A-B)÷(2 10^-5) -> D


Definite Integral - Simpson's Rule


∫ f(x) dx ≈ h/3 * (y_a + ∑(4*y_odd + 2*y_even) + y_b)

where h=(b-a)/n   (n is even)


The integral is stored in the variable I.


"A"? -> A 

"B"? -> B

"N"? -> N

A -> X

Prog 0

Ans -> I

(B-A)÷N -> D

N÷2 -> K

Lbl 2

X+D -> X : Prog 0 : I + 4 Ans -> I

X+D -> X : Prog 0 : I + 2 Ans -> I

K - 1-> K

K≠0 => Goto 2

B -> X : Prog 0 : (I - Ans)D÷3 -> I



Eddie


All original content copyright, © 2011-2021.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Monday, May 20, 2019

Intensity/Illumination, Days Since Jan. 1, Derivatives: HP 20S and 21S

Intensity/Illumination, Days Since Jan. 1, Derivatives: HP 20S and 21S

Table of contents

1.  Intensity and Illumination
2.  Days Since January 1
3.  Numerical Derivative

Disclaimer: I believe the key codes for the programs in this blog entry are the all the same even though the HP 20S and HP 21S have slightly different codes.  Format:   Step  Key: Key Code.  I took turns programming the HP 20S and HP 21S.

1.  Intensity and Illumination

The follow equation relates the luminous intensity (measured in candelas, cd) and illuminance (measured in lux) of a light source.  The equation assumes the light source radiates a spherical matter.

E = I / R^2

E = illuminance
I = luminous intensity
R = radius of the sphere's light (meters)

LBL A:  Solve for E
LBL B:  Solve for I
LBL C:  Solve for R

Registers:
R0 = E
R1 = I
R2 = R

Store the following values in the register and execute the appropriate label.

Program:

01  LBL A:  61,41,A
02 RCL 1:  22, 1
03 ÷:  45
04 RCL 2:  22, 2
05 x^2:  51, 11
06 =: 74
07 STO 0:  21, 0
08 R/S:  26
09 LBL B:  61,41,B
10  RCL 0: 22,0
11  *:  55
12  RCL 2: 22,2
13  x^2:  51,11
14  =:  74
15  STO 1: 21, 1
16  R/S:  26
17  LBL C: 61,41,C
18  RCL 1: 22,1
19 ÷: 45
20 RCL 0: 22,0
21  =:  74
22 √:  11
23 STO 2: 21, 2
24 R/S:  26

Example 1:
R1 = I = 400, R2 = R = 2.
Solve for E,  XEQ A returns 100

Example 2:
R0 = E = 180, R2 = R = 3
Solve for I, XEQ B returns 1620

Example 3:
R1 = I =420, R0 = E = 195
Solve for R, XEQ C returns 1.467598771

2.  Days Since January 1

Calculate the number of days since January 1.  For more information, please see:  http://edspi31415.blogspot.com/2019/03/ti-84-plus-and-hp-41c-number-of-days.html

Input:
R1:  day
R2: month
R3: 0 if we are working in a non-leap year, 1 if we are working in a leap year

Output:
R4:  number of days since January 1

Program:

01 LBL A: 61,41,A
02 RCL 1:  22,1
03 STO 4: 21, 4
04 3:  3
05 5:  5
06 STO - 4:  21,65,4
07 RCL 2: 22,2
08  INPUT:  31
09  2:  2
10  x ≤ y?:  61,42
11  GTO 2:  51,41,2
12  RCL 2: 22,2
13  *:  55
14  3:  3
15  0:  0
16  . : 73
17 6:  6
18  +:  75
19  1:  1
20  .  : 73
21  6:  6
22  =:  74
23  IP:  51, 45
24  STO + 4:  21,75,4
25  RCL 3:  22,3
26  STO + 4:  21,75,4
27  RCL 4: 22,4
28  RTN:  61, 26
29  LBL 2:  61,41,2
30  RCL 2:  22,2
31   *:  55
32  3:  3
33  0:  0
34  .  : 73
35  6:  6
36  +:  75
37  3:  3
38  6:  6
39  8:  8
40:  .  : 73
41  8:  8
42 =: 74
43  IP:  51,45
44  STO + 4:  21,75,4
45  3:  3
46  6:  6
47  5:  5
48  STO - 4:  21,65,4
49  RCL 4:  22,4
50  RTN:  61,26

Example 1:
1/1/2019 - 5/7/2019  (non-leap year)
R1:  7,  R2:  5,  R3:  0
Result:  R4 = 126

Example 2:
1/1/2020 - 11/14/2020  (leap year)
R1:  14,  R2:  11, R3:  1
Result:  R4 = 318

3.  Numerical Derivative

f'(x0) ≈ ( f(x0 + h) - f(x0 - h) ) / ( 2*h )

x = point
h = small change of x, example h = 0.0001

LBL A:  Main Progam
LBL F:  f(X), where R0 acts as X

Input variables:
R1 = h
R2 = point x0

Used variables:
R0 = x   (use R0 for f(x), LBL F)

Calculated Variables:
R3 = f'(x)

Radians mode will be set.

Program:

01  LBL A:  61,41A
02  RAD:  61,24
03  RCL 2:  22,2
04  +:  75
05  RCL 1:  22, 1
06  =:  74
07  STO 0:  21,0
08 XEQ F:  41,F
09  STO 3:  21, 3
10  RCL 2: 22,2
11 -:  65
12 RCL 1: 22,1
13 =:  74
14 STO 0: 21,0
15 XEQ F:  41,F
16 STO - 3:  21,65,3
17 2:  2
18 STO ÷ 3:  21,45,3
19  RCL 1:  22,1
20 STO ÷ 3: 21,45,3
21  RCL 3:  22,3
22 R/S:  26
23 LBL F:  61,41,F
...
xx  RTN:  61,26  (end f(X) with RTN)

Example:  e^x * sin x

LBL F
RCL 0
e^x
*
RCL 0
SIN 

RTN

R1 = 0.0001
R2 = x0 = 0.03
Result:  1.060899867

R1 = 0.0001
R2 = x0 = 1.47
Result:  4.7648049


Note:  I am going on vacation this week and I have jury duty in June. So far, I have blog entries scheduled to be posted throughout June 22.  I plan to have a weekly post every Saturday in June.  - E

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Sunday, July 10, 2016

HP 42S Programming Part III: Numerical Derivative, Recurring Sequences, Fractions

HP 42S Programming Part III: Numerical Derivative, Recurring Sequences, Fractions

Click here for Part I:  Matrix Column Sum, GCD, Error Function


HP 42S:  Numerical Derivative

The program DERVX is based off my numerical derivative program for the HP 15C.  To see the HP 15C program, click here:  http://edspi31415.blogspot.com/2011/11/numerical-derivatives-in-part-9-we-will.html

DERVX calculates numerical derivatives of f(x), given a small increment h. Generally, the smaller h is, the better the calculation. However with certain methods, if h is too small, the final result may be unexpected.

This program uses a five-point formula:

f'(x) ≈ (f(x - 2h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2h)) / (12h)

The error is of the order of h^4. 

Source: Burden, Richard L. and J. Douglas Faires. "Numerical Analysis 8th Edition" Thomson Brooks/Cole Belton, CA 2005

Labels Used:
Label DERVX: Main routine
Label 00: f(X).  This function starts with X on the X register. 

Numerical memory registers will be used, and they are the following:
R00 = the numerical derivative
R01 = X
R02 = h

Radians mode is automatically set.

HP 42S:  Program DERVX
00 {“70+”-Byte Prgm}  \\ number varies depends on what is f(X)
01 LBL “DERVX”
02 RAD
03 STO 02  \\ store h
04 R↓ \\ roll down
05 STO 01 \\ store x
06 2
07 RCL* 02
08 –
09 XEQ 00 \\ subroutine
10 STO 00
11 RCL 01
12 RCL- 02
13 XEQ 00
14 8
15 +/-
16 *
17 STO+ 00
18 RCL 01
19 RCL+ 02
20 XEQ 00
21 8
22 *
23 STO+ 00
24 RCL 01
25 2
26 RCL* 02
27 +
28 XEQ 00
29 +/-
30 STO+ 00
31 RCL 00
32 RCL÷ 02
33 12
34 ÷
35 STO 00
36 “d/dX =”   \\ lower case d:  [shift] { D } in ALPHA
37 ARCL ST X
38 RTN \\ main program ends

39 LBL 00
NN RTN \\ end f(X) with RTN
NN+1 END

How to run DERVX:

1.  If desired, edit the function by [shift] [XEQ] (GTO) 00, [shift] [R/S] (PRGM).  Remember you start with x in the X stack.  End f(x) with RTN. 
2.  Input x, press [ENTER].
3.  Input h, press [XEQ] {DERVX}.


Test 1:  f(x) = x^2 + x*sin x, where x = Ï€/4, h = 1E-5

The code for f(x) would look like this:

39 LBL 00
40 X^2
41 LASTX
42 ENTER
43 SIN
44 *
45 +
46 RTN
47 END

Ï€ [ENTER] 4 [÷]  1 [E] -5 [XEQ] {DERVX}
Result:  2.8333

Test 2:  f(x) = x*e^x, where x = 1.75, h = 1E-5

The code for f(x) would look like this:

39 LBL 00
40 ENTER
41 E^X
42 *
43 RTN
44 END

1.75 [ENTER] 1 [E] -5 [XEQ] {DERVX}
Result:  15.8252

HP 42S: Recurring Sequence Matrix

The program RECUR builds a table for the sequence f(n, u(n-1)).  The table is stored in the matrix MAT.  At the completion of the program, you will be to see the elements of MAT in matrix edit mode.

LBL 00 holds the sequence function. 
R00 = u(n-1)

Specify as many entries as you like, up to 999.

To use the counter in the sequence function, use RCL 02, IP. 

Other registers that are used are:
R01 = n
R02 = counter, in the form of 2.nnn

HP 42S Program RECUR
00 {“95+” Byte Prgm}  \\ varies depending on the sequence function
01 LBL “RECUR”
02 “U(1)=”
03 PROMPT
04 STO 00  \\ store the initial condition
05 “# TERMS:”
06 PROMPT
07 STO 01  \\ store n
08 2
09 DIM “MAT” \\ create the matrix MAT (n x 2)
10 1E-3
11 RCL* 01
12 2
13 +
14 STO 02 \\ set up counter
15 INDEX “MAT” \\ index the matrix
16 1
17 ENTER
18 STOIJ  \\ set index counters at I=1 (row), J=1 (col)
19 1
20 STOEL  \\  store 1 into MAT(1,1)
21 J+  \\  XEQ “J+” 
22 RCL 00
23 STOEL  \\ store U(1) into MAT(1,2)
24 J-  \\ XEQ “J-“, set column pointer back to 1
25 LBL 10  \\ main loop
26 I+ \\ XEQ “I+”, go to the next row
27 RCL 02
28 IP
29 STOEL
30 J+  \\ XEQ “J+”
31 RCL 00  \\ recall u(n-1)
32 XEQ 00 \\ calculate f(n, u(n-1))
33 STO 00 \\ store u(n) for the next loop
34 STOEL
35 J-  \\ XEQ “J-“
36 ISG 02  \\ increase counter
37 GTO 10
38 RCL “MAT”
39 EDIT \\ put MAT in viewing/editing mode
40 RTN  \\ main program ends

41  LBL 00
NN RTN \\ end sequence function with RTN
NN+1 END

How to run RECUR:

1.  If desired, edit the function by [shift] [XEQ] (GTO) 00, [shift] [R/S] (PRGM).  Remember you start with u(n-1) in the X stack.  End f(x) with RTN. 
2.  Press [XEQ] {RECUR}.  You will be prompted for U(1) and the number of rows desired.

Test 1:  u(n) = 2*u(n-1),  u(1) = 1, 5 rows wanted

41 LBL 00
42 2
43 *
44 RTN
45 END

[XEQ] {RECUR} 1 [XEQ] 5 [XEQ]

MAT =
1
1
2
2
3
4
4
8
5
16

Test 2:  u(n) = u(n-1)^2 – u(n-1)^3/6,  u(1) = 1.5,  5 rows wanted

41 LBL 00
42 ENTER
43 X^2
44 X<>Y
45 3
46 Y^X
47 6
48 ÷
49 –
50 RTN
51 END

MAT =
1
1.5000
2
1.6875
3
2.0468
4
2.7602
5
4.1138


 HP 42S: Fractions

The program FRACT provides three calculations involving fractions:

ADD:  add two fractions
W/X + Y/Z = (W*Z + X*Y)/(X*Z)

MULT: multiply two fractions
W/X * Y/Z = (W*Y)/(X*Z)

POWER: take a fraction to a power
(W/X)^Y = W^Y/X^Y

The end of the program stores the resulting numerator in N and the resulting denominator in D.  FRACT simplifies the fraction after calculation.  Note:  FRACT sets the HP 42S to FIX 0 mode. Mixed fractions are not included.

HP 42S:  Program FRACT
00 {255-Byte Prgm}
01 LBL “FRACT”
02 FIX 00
03 LBL 99 \\ start menu
04 CLMENU
05 LBL 00
06 “ADD”
07 KEY 1 GTO 01  \\ KEYG
08 “MULT”
09 KEY 2 GTO 02
10 “POWER”
11 KEY 3 GTO 03
12 MENU
13 STOP
14 GTO 99  \\ end menu
15 LBL 01  \\ addition routine
16 CLMENU
17 “W/X + Y/Z”
18 AVIEW
19 PSE
20 INPUT “W”
21 INPUT “X”
22 INPUT “Y”
23 INPUT “Z”
24 RCL “W”
25 RCL* “Z”
26 RCL “X”
27 RCL* “Y”
28 +
29 STO “N”
30 RCL “X”
31 RCL* “Z”
32 STO “D”
33 GTO 10 \\ go to GCD routine
34 LBL 02 \\ multiply routine
35 CLMENU
36 “W/X x Y/Z”
37 AVIEW
38 PSE
39 INPUT “W”
40 INPUT “X”
41 INPUT “Y”
42 INPUT “Z”
43 RCL “W”
44 RCL* “Y”
45 STO “N”
46 RCL “X”
47 RCL* “Z”
48 STO “D”
49 GTO 10 \\ go to GCD routine
50 LBL 03 \\ power routine
51 CLMENU
52 “(W/X)/\Y”  \\ the power symbol is made of two symbols: / and \
53 AVIEW
54 PSE
55 INPUT “W”
56 INPUT “X”
57 INPUT “Y”
58 RCL “W”
59 RCL “Y”
60 Y^X
61 STO “N”
62 RCL “X”
63 RCL “Y”
64 Y^X
65 STO “D”
66 GTO 10  \\ go to GCD routine
67 LBL 10 \\ GCD routine
68 RCL “N”
69 RCL “D”
70 X<Y?
71 X<>Y
72 STO 01
73 X<>Y
74 STO 00
75 LBL 11
76 RCL 01
77 ENTER
78 RCL÷ 00
79 IP
80 RCL* 00
81 –
82 X=0?
83 GTO 12
84 X<>Y
85 STO 01
86 X<>Y
87 STO 00
88 GTO 11
89 LBL 12 \\ ending and display
90 CLA \\ clear ALPHA register
91 RCL “N”
92 RCL÷ 00 \\ GCD = R00
93 STO “N”
94 ARCL ST X
95 ├ “/”   \\ ├ is the attach symbol, press [shift] [ENTER] (ALPHA) [ENTER]
96 RCL “D”
97 RCL÷ 00
98 STO “D”
99 ARCL ST X
100 AVIEW
101 END

Test 1:  Addition (ADD)
12/19 + 8/17 = 356/323
W = 12, X = 19, Y = 8, Z = 17
Result:  “356./323.”   Y: 356, X: 323

Test 2: Multiplication (MULT)
12/19 * 8/17 = 96/323
W = 12, X = 19, Y = 8, Z = 17
Result:  “96./323.”   Y: 96, X: 323

Test:  Power (POWER)
(2/7)^3 = 8/343
W = 2, X = 7, Y = 3
Result:  “8./343.”,  Y: 8, X: 343

That is Part III, I am happy I was able to post this a lot sooner than I thought.  

Eddie

This blog is property of Edward Shore, 2016.






RPN HP 12C: Fibonacci and Lucas Sequences

  RPN HP 12C: Fibonacci and Lucas Sequences Golden Ratio, Formulas, and Sequences Let φ be the Golden Ratio: φ = (1 + √5) ÷ 2...