HP 75C Program Collection

HP 75C Program Collection

Merry Christmas everyone!   I hope everyone is safe, happy, and blessed.

LAWCOS:  Law of Cosines

335 bytes

100 DISP "Law of Cosines" @ WAIT .25

105 OPTION ANGLE DEGREES

180 DISP "Θ is opp. side a" @ WAIT .25

110 INPUT "1) Side, 2) Angle? "; H

115 ON H GOTO 1000, 2000

1000 INPUT "b?, c?, Θ? "; B,C,T

1010 A=SQR(B^2+C^2-2*B*C*COS(T))

1020 DISP "a: "; STR$(A) @ END

2000 INPUT "a?, b?, c? "; A,B,C

2010 T=ACOS((B^2+C^2-A^2)/(2*B*C))

2020 DISP "Θ:"; STR$(T); "°" @ END

Example:

1) Side:  1

B = 11.78, C = 32.12, T = 60

Result:  a = 28.1440793063

2)  Angle: 2

A = 14.55, B = 12.65, C = 17.00

Result:  Θ = 56.5108780171°

SPECPOLY:  Special Polynomials

526 bytes  (10/27/2022)

1.Cheb:  Chebyshev Polynomial of the 1st order

2.Herm:  Hermite Polynomial 

3.Lagur:  Laguerre Polynomial (with α=0)

4.Legen:  Legendre Polynomial of the 1st Order

100 DISP "Polynomials:  Type?" @ WAIT .5

105 INPUT "1.Cheb 2.Herm 3.Lagur 4.Legen:"; H

110 INPUT "x? "; X

115 INPUT "order? ";N

120 A=1 @ B=X*(H=1 OR H=4)+2*X*(H=2)+(1-X)*(H=3)

300 IF N=0 THEN DISP A @ END

302 IF N=1 THEN DISP B @ END

305 FOR I=2 TO N

310 ON H GOSUB 500,505,510,515

315 A=B

320 B=C

325 NEXT I

330 DISP C

335 END

500 C=2*X*B-A @ RETURN

505 C=2*X*B-2*(I-1)*A @ RETURN

510 C=((2*(I-1)+1-X)*B-(I-1)*A)/I @ RETURN

515 C=((2*I-1)*X*B-(I-1)*A)/I @ RETURN

Source:

Davidson, James J.  "Chebyshev Polynomials, First & Second Kind, Tn(x) & Un(x)", "Hermite Polynomials, H(n); n≥n", "Laguerre Polynomials, Generalized, Ln^(α)(x)", "Legendre Functions, First & Second Kinds, Pn(x) & Qn(x)"  ENTER 65 NOTES  Vol. 3 No. 8 September 1976  pp. 11-12

Examples:

x = 11, order = 6

1.  Cheb:  55989361

2.  Herm: 106439224

3.  Lagur: -35.5902777778

4.  Legen: 25289461

ENGINE:  Automotive Engine Mathematics

630 bytes (10/27/2022)

100 DISP "Engine Math" @ WAIT .5

110 INPUT "# cylinders? "; N

120 DISP "Solve for?" @ WAIT .25

125 INPUT "1)DSP, 2)STROKE, 3)BORE "; H

130 IF H#1 THEN INPUT "DSP? "; D

135 IF H#2 THEN INPUT "STROKE? "; S

140 IF H#3 THEN INPUT "BORE? "; B

145 IF H=1 THEN LET D=PI/4*B^2*S*N

150 IF H=2 THEN LET S=D/(PI/4*B^2*N)

155 IF H=3 THEN LET B=SQR(D/(PI/4*S*N))

160 DISP "DPS= "; D @ STOP

165 DISP "STORKE= "; S @ STOP

170 DISP "BORE= "; B @ STOP

205 INPUT "RPM? (y/n) "; I$

210 IF UPRC(I$)="N" THEN 400

300 INPUT "RPM? "; R

310 P=S*R/6

320 E=R*D*.85/3456

350 DISP "Piston Speed="; P @ STOP

360 DISP "St. Carb CFM= "; E @ STOP

370 INPUT "Again? (y/n) "; I$ 

375 IF UPRC(I$)="Y" THEN 300

400 DISP "Done."

Source:

Lawlor, John.  Auto Math Handbook.  Calculations, Formulas, Equations, and Theory for Automotive Enthusiasts.  HP Books:  Berkeley Publishing Group:  New York, NY.  1991.  ISBN 1-55788-020-4

LINREG:  Linear Regression

568 bytes (10/31/2022)

100 DISP "y=mx+b" @ WAIT .5

110 S0=0 @ S1=0 @ S2=0 @ S3=0 @ S4=0 @ S5=0

120 DISP 'Σ+ Σ- Calc'

130 K$=KEY$ @ IF K$='' THEN 130

140 IF K$="+" THEN H=1 @ GOTO 200

150 IF K$="-" THEN H=-1 @ GOTO 200

160 IF UPRC$(K$)='C' THEN 400

170 GOTO 130

200 INPUT 'x,y? '; X,Y

210 S0=S0+H @ S1=S1+H*X @ S2=S2+H*Y

220 S3=S3+H*X^2 @ S4=S4+H*Y^2 @ S5=S5+H*X*Y

230 DISP S0 @ WAIT .5 @ GOTO 120

400 B=(S5-S1*S2/S0)/(S3-S1^2/S0)

440 A=(S2-B*S1)/S0

450 R=B*SQR((S0*S3-S1^2)/(S0*S4-S2^2))

470 DISP 'type cont...' @ WAIT .5

500 DISP 'SLP=';B @ STOP

510 DISP 'ITC=';A @ STOP

520 DISP 'CORR=';R

Instructions:

Press + to add data.  Separate  x and y with a comma.

Press - to delete data.  Separate x and y with a comma.

Press C to determine the slope, intercept, and correlation.

Example:

Data:

5, 4.066

8, 4.125

11, 4.202

18, 4.293

21, 4.349

ITC ≈ 3.99046

SLP ≈ 0.01719

CORR ≈ 0.99376

Source:

Hewlett Packard.  HP 12C User Guide  Edition 4.  San Diego, CA  2004.

SIMPSON:  Integral Using Simpson’s Rule Approximation

size varies - around 300 bytes

∫ FNF(X) dX from X = a to X = b

Line 100 is where you have to define your function (of X). 

100 DEF FNF(X) = [define f(x) here]

115 OPTION ANGLE RADIANS

120 INPUT "a? ";A

125 INPUT "b? ";B

130 INPUT "# parts (even)? ";N

135 T=FNF(A)+FNF(B)

140 H=(B-A)/N

145 FOR I TO N-1

150 IF FP(I/2)=0 THEN 155 ELSE 160

155 T=T+2*FNF(A+I*H) @ GOTO 165

160 T=T+4*FNF(A+I*H)

165 NEXT I

170 T=T*H/3

175 DISP "Integral = ";T

Examples:

For both examples, n = 24

FNF(X)= 2*X^2+3,  a = 1, b = 6, result:  158.333333333

FNF(X)=2*COS(X), a = 0, b = .785398163398, result:  1.41421357139  (exact √2)

Happy Holidays everyone!

Next post will be on New Years' Eve:  December 31, 2022.

Eddie 

All original content copyright, © 2011-2022.  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. 

Python - Lambda Week: Integration by Simpson's Rule

Python - Lambda Week: Integration by Simpson's Rule

Welcome to Python Week!  This we we're going to cover calculus and the keyword lambda.

Note:  All Python scripts presented this week were created using a TI-NSpire CX II CAS.   As of June 2022, the lambda keyword is available on all calculators (in the United States) that have Python.   If you are not sure, please check your calculator manual. 

Simpson's Rule

The Simpson's Rule estimates numeric integrals by:

∫( f(x) dx, x = a to b) ≈

(b - a) /(3 * n) * (f(a) + 4 * f1 + 2 * f2 + 4 * f3 + .... + 2 * f_n-2 + 4 * f_n-1 + f(b))

n must be an even number of partitions.  The more partitions, the higher the accuracy and the higher computation time.

integrallam.py:  Numeric Integer

from math import *

print("The math module is imported.")

print("Integra of f(x), 6 places")

f=eval("lambda x:"+input("f(x)? "))

# input parameters

a=eval(input("lower = "))

b=eval(input("upper = "))

n=int(input("even parts: "))

# checksafe, add 1 if n is odd

if n/2-int(n/2)==0:

  n=n+1

# integral calculus

s=f(a)+f(b)

w=1

# 1 to n-1

for i in range(1,n):

  w=f(a+i*(b-a)/n)

  s+=(2*w) if (i/2-int(i/2)==0) else (4*w)

s*=(b-a)/(3*n)

print("Integral: "+str(round(s,6)))

All original content copyright, © 2011-2022.  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. 

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. 

Retro Review: Radio Shack EC-4027

Retro Review:  Radio Shack EC-4027

Quick Facts:

Company:  Radio Shack

Year Introduced:  1993

Type:  Scientific, Formula Programming

Memory:  40 steps

Operating System:  Algebraic

Memory Registers:  7, A through F, M

Batteries:  2 x LR44 (357)

Equivalent Of:  Sharp EL-5020

Background

From the 1970s to about the 1990s (maybe early 2000s), Radio Shack sold OEM equivalents of calculators from Sharp, Casio, and Texas Instruments, and the EC-4027 is no exception.  One advantage to buying Radio Shack is that the price tends to be cheaper but with equivalent quality of the original.

An Update to the Sharp EL-5100

The EC-4027 (and the Sharp EL-5020) is a scientific calculator which holds multi-step formulas.   There are no switches, everything is keys this time.  The following modes offered are:

*  Normal:  Computational mode

*  AER:  Program and store formulas

* ∫dx:  Integral Mode

* STAT 1:  Single-variable statistics mode

* STAT 2:  Linear Regression, y = a+bx

Base conversions and fractions are added.  

Like the EL-5100, expressions are entered algebraically as you would write it down.   A frustrating part is that the entire expression is not shown on the screen, regulating functions to the sides of the screen, which makes editing expressions difficult if not impossible.

Integration and Loops but With Less Memory

The memory of the EC-4027 is half of the EL-5100:

* 40 steps instead of 80

* Only two formulas can be stored instead of five

We do get comparisons against the M register, all six comparisons (<, ≤, =, >, ≥, ≠).  If the condition is met, the formula is starts over from the beginning, creating a loop.   

There is a separate mode for integration.  Either formula 1 or 2 can be called, and Simpson's Rule is used.  Any function that can be integrated is in terms of M.  

I wish the memory wasn't halved but instead increased.  Also, expand the ALPHA characters to include the entire alphabet.  But, the EC-4027 serves it's purpose.

Sources:

Toth, Viktor T. "Radio Shack EC-4027"  2000-2018.  https://www.rskey.org/ec4027  Last accessed April 9, 2021

Voidware.  "Sharp EL-5020" http://www.voidware.com/calcs/el5020.htm  Last accessed April 9, 2021. 

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. 

Numworks/Casio MicroPython/Python: Calculus

Numworks/Casio MicroPython/Python:  Calculus 

Introduction

The following scripts creates the user functions for calculus:

f(x):  define your function in terms of x here.   This needs to be loaded into the script before running it.   Each of the functions that follow will use f(x).  You can call f(x) to evaluate the function at any value.

deriv(x):  The approximate derivative at point x.  The Five Stencil approximation is used. 

sigma(a,b):  Calculate the sum  (Σ f(x)) from x = a to x = b. 

integral(a,b,n):  Calculates the definite integral ( ∫ f(x) dx) from x = a to x = b.  The Simpson's rule is used with n divisions (n needs to be even)

solve(x0):  Uses Newton's Rule to find roots for f(x). 

Example

f(x) = -2x*^2 + 3x + 5
In Python:  -2*x**2+3*x+5

f(0): 5
f(10): -165
f(-10): -225

deriv(10): -37.00004450971999

Σ f(x): x = 1 to 25:  sigma(1,25): -8780

∫ f(x) dx: x = -3 to 1, n = 20:  integral(-3,1,20):  -10.666666666667

Solve f(x)=0, initial condition x0 = 2.5:  solve(5):  2.5 

Python Script: calculus.py

from math import *

# 2020-04-15 EWS

# define f(x) here
def f(x):
  return -2*x**2+3*x+5
  
# derivative
def deriv(x):
  # uses f(x), 5 stencil
  # h is tolerance
  h=1e-10
  d=(12*h)**-1*(f(x-2*h)-8*f(x-h)+8*f(x+h)-f(x+2*h))
  return d

# sum/sigma
def sigma(a,b):
  t=0
  n=b-a
  for i in range(n):
    t=t+f(i+1)
  return t

# integral by simpsons rule
def integral(a,b,n):
  t=f(a)+f(b)
  h=(b-a)/n
  for i in range(n-1):
    w=(i+1)/2
    if (w-int(w))==0:
      t=t+2*f(a+(i+1)*h)
    else:
      t=t+4*f(a+(i+1)*h)
  t=t*h/3
  return t

# solver
def solve(x0):
  tol=1e-14
  x1=x0-f(x0)/deriv(x0)
  while abs(x1-x0)>tol:
    x0=x1
    x1=x0-f(x0)/deriv(x0)
  return x1


Eddie

All original content copyright, © 2011-2020.  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.

Fun with the Sharp EL-5500 III

Fun with the Sharp EL-5500 III

See my review for the EL-5500 III here:  http://edspi31415.blogspot.com/2017/01/retro-review-sharp-scientific-computer.html

Even though the EL-5500 III has one programming space, we can fit many programs in it.  One of the great features is that how programs take relatively little space.  Obviously you can change the line numbers, the sequence of lines are most important.

Euclid Algorithm – RUN 10

10 PAUSE “GCD: A>B”

18 INPUT “A: “; A, “B: “; B

20 C = A – INT(A/B)*B

30 IF C = 0 THEN 40

32 A = B

34 B = C

36 GOTO 20

40 PRINT “GCD = “; B

46 END

My first program on the EL-5500 III.  We can probably make the coding more efficient.  PRINT with a proceeding WAIT pauses execution until [ ENTER ] is pressed.

The command END ends the program execution and allows the space to have more than one routines.

Examples:

A: 100, B: 20.  Result:  GCD = 20

A: 63, B: 60.  Result:  GCD = 3

Binomial Expansion – RUN 50

50 PAUSE “(Ax + B)^n”

55 INPUT “A: “;A, “B: “;B, “N: “,N

60 FOR I=0 TO N

65 C = FACT(N) / (FACT(I) * FACT(N-I))

70 T = C * A^(N-I) * B^I

75 PRINT T; “*x^”; N-I

80 NEXT I

85 END

The command PAUSE shows the text for about 0.85 seconds.

This program displays each coefficient one at time.

Example:  A = 3, B = -2, N = 3.  Expand (3 – 2x)^3.

Result:  27, -54, 36, -8  (27x^3 – 54x^2 + 36x – 8)

Days Between Dates – RUN 100

100 PAUSE “Days between Dates”

112 INPUT “M1:”; M1, “D1:”; D1, “Y1:”; Y1

113 INPUT “M2:”; M2, “D2:”; D2, “Y2:”; Y2

115 M = M1: D = D1: Y = Y1

117 GOSUB 140

119 F1 = F

121 M = M2: D = D2: Y = Y2

123 GOSUB 140

125 F2 = F

127 N = F2 – F1

129 PRINT “# DAYS: “; N

131 END

140 IF M > 2 THEN 150

142 X = 0: Z = Y – 1

144 GOTO 160

150 X = INT(.4 * M + 2.3): Z = Y

160 F = 365 * Y + 31* (M – 1) + D + INT(Z/4) – X

165 RETURN

If I had it my way, every calculator that isn’t a basic 4-function calculator will have a Days Between Dates function.

Example:  November 4, 1986 to January 12, 2017

M1:  11, D1: 4, Y1: 1986 

M2: 1, D2: 12, Y2: 2017

Result:  11,027

Power of a Complex Number – RUN 200

This calculates (a + bi)^n

200 PAUSE “(A + Bi)^n”

203 INPUT “A: “;A, “B: “;B, “n: “;N

205 IF A<>0 THEN 210

207 A = A^N: GOTO 220

210 R = √(A^2 + B^2)^N

212 T = ATN(B/A)*N

214 A = R * COS T

216 B = R * SIN T

220 PRINT A; “+”; B; “i”

225 END

The command <> means not equals (≠) and ATN is the arctangent function.

Examples:

(2 + 3i)^3:  A = 2, Bi = 3.  Result: -46 + 9i

(-5 + i)^2:  A = -5, Bi = 1.  Result: 24 – 10i

Error Function Approximation – RUN 250

250 PAUSE “APPROX ERF(X)”

252 CLEAR

254 DIM A(3)

256 E = 0

260 INPUT “X > 0 : “; X

262 T = RCP(1 + .47047 * X)

270 FOR I = 1 TO 3

272 READ A(I)

274 E = E + A(I) * T^I

276 NEXT I

280 DATA .3480242, -.0958798, .7478556

290 E = 1 – E * EXP(-(X^2))

292 PRINT “ERF : “, F

294 CLEAR

296 END

Use the command CLEAR to clear all the variables.  At this point, I will use this command both at the beginning and end of calculations.  The command RCP is the reciprocal.  Finally, EXP is from e^x function, not the [EXP] key.

Example:

x = 0.53, erf ≈ 0.546455764

x = 0.88, erf ≈ 7.86708E-01 = 0.786708

Keep in mind this approximation is accurate 2.5 parts in 10^-5.

Source:  Smith, Jon M.  Scientific Analysis on the Pocket Calculator  John Wiley & Sons: New York. 1975.  ISBN 0-471-79997-1

Simpson’s Rule – RUN 300

Calculates the numeric integral ∫ f(x) dx from x = L to x = U.

You can edit the function at line 342.  In BASIC-PRO mode, call up line 342 by LIST 342.  Edit as desired.

300 PAUSE “Simpsons Rule”

302 CLEAR

304 INPUT “LOWER: “; L, “UPPER: “; U, “PARTS (even): “; N

306 RADIAN

308 X = L: GOSUB 340: T = F

310 X = U: GOSUB 340: T = T + F

314 H = (U – L)/N

320 FOR I = 1 TO N-1

322 X = L + I * H: GOSUB 340

324 R = I/2 – INT(I/2)

326 IF R = 0 THEN LET T = T + 2*F

328 IF R <> 0 THEN LET T = T + 4*F

330 NEXT I

332 T = T * H/3

334 PRINT “Integral: “, T

336 CLEAR

338 END

340 REM f(x)

342 F = [insert f(X) here]

344 RETURN

Note:  I finally learn why LET is included in BASIC.  In an IF-THEN statement, if you want to assign a value or calculation to a variable, you will need a LET command.

I use a REM (remark) command on line 340.  REM is for comments and nothing said after REM is executed.

Examples:

342 F = X^2 – 1

L = -2, U = 2, N = 14.  Result:  1.333333334 (Actual: 4/3 ≈ 1.333333333)

342 F = √(X – 1)

L = 2, U = 3, N = 14.  Result: 1.218951372 (Actual:  ≈ 1.2158951416)

Dancing Star Demo – RUN 400

400 PAUSE “Dancing Star Demo”

405 CLEAR

410 S$ = “         “  (9 spaces)

415 FOR I = 1 TO 48

420 N = RND 7 + 1

425 T$ = MID$(S$, 1, N-1) + “*” + MID$(S$, N+1, 9-(N+1))

430 WAIT 15: PRINT T$

440 NEXT I

445 BEEP 3: CLEAR: END

This is just show a dancing asterisk (*) on the screen.

WAIT 15: PRINT T$:   This causes the string T$ after a quarter of a second.  WAIT operates in 1/60 seconds.

RND n:  This is the random command.  If n = 1, the result is between 0 and 1. 

BEEP n:  This makes the EL-5500 III beep n times. 

Eddie

This blog is property of Edward Shore, 2017.

HP 71B: Simpson’s Rule Approximation for f(x,y)

HP 71B: Simpson’s Rule Approximation for f(x,y)

Happy New Year everyone!  Hope everyone had a great New Years Day!  J

For the TI-84 Plus version, please visit here:  http://edspi31415.blogspot.com/2016/02/ti-84-plus-simpsons-rule-for-fxy.html

Today, I will present apply the Simpson’s Rule to functions of two variables like f(x,y) for the HP 71B. Here is the general algorithm:

Start with the function f(x,y) and your integration limits A, B, C, and D.

Then determine Δx and Δy (labeled E and F in the programs below) by:

Δx = (B – A)/(N – 1)

Δy = (D – C)/(N – 1)

Where N is the number of partitions.  Unlike the Simpson’s Rule for one variable, in this case N must be odd.  Generally, the higher N is, the more accurate the approximation is, with the expense of additional computational time.

Next, build a matrix, let’s say [ I ].   This is your Simpson’s Matrix.  The Simpson’s Matrix is built by the expression

[ I ] = [1, 4, 2, 4, 2, 4, 2, 4, 2, …, 4, 1]^T * [1, 4, 2, 4, 2, 4, 2, 4, 2, …, 4, 1]

The length of the vector used to determine [ I ] depends on N.  A way to build it by the routine:

Store 1 into the element 1 of [ I ]  (first element)

Store 1 into the element N of [ I ]  (last element)

For J from 2 to N – 1

If J is divisible by 2, then store 4 in the jth element of [ I ],

Else store 2 in the jth element of [ I ]

For N = 5, the vector would be built is [1, 4, 2, 4, 1]

And

[ I ] = [1,4,2,4,1]^T * [1,4,2,4,1] =

[[1, 4, 2, 4, 1]

[4, 16, 8, 16, 4]

[2, 8, 4, 8, 2]

[4, 16, 8, 16, 4]

[1, 4, 2, 4, 1]]

Build another matrix [ J ].  The elements are determined by the following formulas:

For row j and column k, the element is f(A + Δx*(j – 1),  C + Δy*(k – 1))

Once finished, multiply every element of [ I ] by [ J ].  This is NOT matrix multiplication.   Then sum all of the elements of the results.  In essence:

S = ∑ (j = 1 to N) ∑ (k = 1 to N)   [I](j,k)* [J](j,k)

Determine the final integral approximation as:

Integral = Δx * Δy * 1/9 * S

The program DBLSIMP uses N = 5.  This program works best for f(x,y) where they are polynomials.  On the HP 71B, matrices cannot be typed directly, elements have to be stored and recalled on element at a time.  The program presented does not use modules. 

HP 71B Program DBLSIMP

At least 580 Bytes

Edit f(x,y) at line 10. Use variables X and Y.

5 DESTROY I,J,A,B,C,D
8 DESTROY E,F,K,S,X,Y

10 DEF FNF(X,Y)= [ enter f(X,Y) here ]

12 RADIANS

14 DIM I(5,5)

20 I(1,1) = 1

21 I(1,2) = 4

22 I(1,3) = 2

23 I(1,4) = 4

24 I(1,5) = 1

25 I(2,1) = 4

26 I(2,2) = 16

27 I(2,3) = 8

28 I(2,4) = 16

29 I(2,5) = 4

30 I(3,1) = 2

31 I(3,2) = 8

32 I(3,3) = 4

33 I(3,4) = 8

34 I(3,5) = 2

35 I(4,1) = 4

36 I(4,2) = 16

37 I(4,3) = 8

38 I(4,4) = 16

39 I(4,5) = 4

40 I(5,1) = 1

41 I(5,2) = 4

42 I(5,3) = 2

43 I(5,4) = 4

44 I(5,5) = 1

50 DISP “X: from a to b” @ WAIT 1

52 INPUT “a = “; A

54 INPUT “b = “; B

56 DISP “Y: from c to d” @ WAIT 1

58 INPUT “c = “; C

60 INPUT “d = “; D

62 E = .25 * (B – A)

64 F = .25 * (D – C)

66 S = 0

70 FOR J = 1 TO 5

72 FOR K = 1 TO 5

74 X = A + E * (K – 1)

76 Y = C + F * (J – 1)

78 S = S + FNF(X,Y) * I(J,K)

80 NEXT K

82 NEXT J

90 S = S * E * F/9

95 DISP “INTEGRAL = “ @ WAIT 1

97 DISP S

Examples:

F(X,Y) = 2*Y – 3*X

A = 1, B = 2, C = 2, D = 5

Result:  Integral = 7.5

(Actual answer: 7.5)

F(X,Y) = X^2/Y^2

A = 1, B = 2, C = 2, D = 5

Result:  Integral ≈ 0.70212579101

(Actual answer:  0.7)

F(X,Y) = 0.5*X*EXP(Y)

A = 1, B = 2, C = 2, D = 5

Result:  Integral ≈ 105.942243008

(Actual answer is about 105.768077253)

Source:

Cooper, Ian.  “Doing Physics With Matlab:  Mathematical Routines”  School of Physics, University of Sydney

http://www.physics.usyd.edu.au/teach_res/mp/doc/math_integration_2D.pdf

Retrieved January 30, 2016

Thank you, and wishing you a happy, healthy, and successful 2017!

Eddie

This blog is property of Edward Shore, 2017

Fun with the HP 71B II

Fun with the HP 71B II

Four years ago, I had fun with the HP 71B programming:

http://edspi31415.blogspot.com/2012/06/fun-with-hp-71b.html

http://edspi31415.blogspot.com/2012/06/cubic-formula-basic-program-hp-71b.html

That has been too long ago.  Time to pull out the 71B again.

Digital Root of an Integers

Add all the numbers of an integer and repeat until you have a single digit (1-9).

Program DROOT (52 bytes)

10 DESTROY N

20 INPUT “INTEGER:”;N

22 N=IP(N)

30 D=1+MOD(N-1,9)

40 DISP D

Easy Traverse Calculation

Calculates the new point knowing the original coordinates, direction, and angle of travel.  The angle 0° comes from due east and rotates counterclockwise (see diagram below). 

Program TRAVEZ (216 bytes)

15 DEGREES

20 INPUT “INIT. EASTING:”;E

25 INPUT “INIT. NORHTING:”;N

30 D=0

40 INPUT “DISTANCE:”;I

45 INPUT “ANGLE:”;A

50 E=I*COS(A)+E

55 N=I*SIN(A)+N

57 D=D+I

60 DISP N;’,’;E

65 PAUSE

70 INPUT “DONE? (Y=1,N=0)”;X

75 IF X=0 THEN 40

90 DISP “TOTAL DIST:”;D

Integral Using Simpson’s Rule

∫ FNF(X) dX from X = a to X = b

Line 10 is where you have to define your function (of X). 

Program SIMPSON  (200+ bytes)

10 DEF FNF(X)=insert function of X here

30 INPUT “LOWER:”;A

32 INPUT “UPPER:”;B

34 INPUT “# PART’NS (EVEN):”;N

36 RADIANS

38 T=FNF(A)+FNF(B)

40 H=(B-A)/N

50 FOR I=1 TO N-1

52 IF FP(I/2)=0 THEN 54 ELSE 56

54 T=2*FNF(A+I*H)+T @ GOTO 58

56 T=4*FNF(A+I*H)+T

58 NEXT I

60 T=T*H/3

62 DISP “INTEGRAL:”;T

Expanding the Binomial (ax+b)^n

BINOMEXP gives the coefficients of the expansion of (ax+b)^n. 

Program BINOMEXP (about 189 bytes)

10 DESTROY A,B,N,L,I,C

12 OPTION BASE 0

14 DISP ‘exapand(Ax+B)^N’ @ WAIT 1

16 INPUT “A,B,N”;A,B,N

18 DIM L(N)

20 FOR I=0 TO N

22 C=FACT(N)/(FACT(I)*FACT(N-I))  \\ FACT is the factorial function

24 L(I)=C*A^(N-I)*B^I

26 DISP “L(“; I; ”):”; L(I); “x^”; N-I  @ PAUSE

28 NEXT I

40 DISP “DONE, CHECK L.”

Synthetic Division

Divide the polynomial p(x) by (x-R).  P is the array of p(x), Q is the array representing the quotient q(x), and E is the remainder.  Hence:  P(x)/(x-R) = Q(x) + E/(x-R)

Program SYNTH (304 bytes)

10 DESTROY P,I,Q,R,E

12 OPTION BASE 0

14 DIPS “P(X)/(X-R)” @ WAIT 1

16 INPUT “DEGREE OF P(X):”;N

18 DIM P(N),Q(N)

20 FOR I=0 TO N

22 DISP “COEF OF x^”;N-1 @ PAUSE

24 INPUT P(I) @ P(I)=Q(I)

26 NEXT I

30 INPUT “R:’;R

40 FOR I=0 TO N-1

42 Q(I+1)=R*Q(I)+P(I+1)

44 NEXT I

50 E=Q(N)

60 DIM Q(N-1)

70 DISP “COEF OF Q(X)”

72 FOR I=0 TO N-1

74 DISP Q(I); “x^”; N-I-1 @ PAUSE

76 NEXT I

80 DISP “REMAIN:”; E; “/(x-“; R; “)”

Dew Point Measurement (in °F)

 (source:  http://ag.arizona.edu/azmet/dewpoint.html )

Program DEWPOINT (159 bytes)

10 DESTROY T,H,V,W,C,D

20 INPUT “TEMP °F”;T

22 C=(T-32)*5/9

26 INPUT “REL HUMIDITY (%):”;H

30 V=(LOG(H%1)+17.27*C/(237.3+C))/17.27

32 W=237.3*V/(1-V)

40 D=9/5*W+32

42 DISP “DEW POINT °F:”;D

// Degree symbol:  [ g ] [RUN] (CTRL) [ A ]

 This blog is property of Edward Shore, 2016