Eddie's Math and Calculator Blog: November 2022

Saturday, November 26, 2022

HP Prime CAS: Curvature

HP Prime CAS: Curvature

Introduction

The following CAS functions calculates the curvature of:

functions, y(x)

polar functions, r(t) (t: Θ)

parametric functions, x(t), y(t)

Let Δα be the angle of rotation angle and Δs is the slight change of distance. Then the radius of curvature is:

K = abs(Δα ÷ Δs) as Δs → 0

And the radius of curvature is the reciprocal of K.

For circles, the radius of curvature is constant. Wankel engines and rotary engines have their pistons traveling in a circle.

Calculating the curvature depends on the form of the function.

Function: y(x)

K = abs( y''(x) ) ÷ (1 + (y'(x))^2) ^(3/2)

Polar: r(t) (t replaces Θ)

K = abs( r(t)^2 + 2 * (r'(t))^2 - r(t) * r''(t) ) ÷ ( r(t)^2 + r'(t)^2 )^(3/2)

Parametric: x(t), y(t)

K = abs( x'(t) * y''(t) - y'(t) * x''(t) ) ÷ ( x'(t)^2 + y'(t)^2 )^(3/2)

Radius of Curvature:

r = 1 ÷ K

For the CAS functions, they take the form:

#cas

name(arguments):=

BEGIN

...

END;

#end

Clicking on the CAS checkbox will not put the #cas and #end delimiters. And these programs will work in CAS mode only.

HP Prime CAS Program: crvfunc

#cas

crvfunc(y,x):=

BEGIN

// curvature

// function

// radius = 1/curvature

LOCAL a,b;

a:=diff(y,x,2);

b:=diff(y,x,1);

RETURN ABS(a)/(1+b^2)^(3/2);

END;

#end

HP Prime CAS Program: crvpol

#cas

crvpol(r,t):=

BEGIN

// curvature

// polar (t: θ)

// radius = 1/curvature

LOCAL a,b,n,d;

a:=diff(r,t,2);

b:=diff(r,t,1);

n:=simplify(r^2+2*b^2-r*a);

d:=r^2+b^2;

RETURN ABS(n)/(d)^(3/2);

END;

#end

HP Prime CAS Program: crvpar

#cas

crvpar(y,x,t):=

BEGIN

// curvature

// parametric

// radius = 1/curvature

LOCAL y1,y2,x1,x2,n,d;

y2:=diff(y,t,2);

y1:=diff(y,t,1);

x2:=diff(x,t,2);

x1:=diff(x,t,1);

n:=simplify(x1*y2-y1*x2);

d:=simplify(x1^2+y1^2);

RETURN ABS(n)/(d)^(3/2);

END;

#end

Until next time and have a great day,

Eddie

Source:

Svirin, Alex Ph.D. "Curvature and Radius of Curvature" Math24 https://math24.net/curvature-radius.html 2022. Last Updated September 12, 2022.

Gratitude to Arno K. and rombio for helping me with derivatives and CAS programs.

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.

Sunday, November 20, 2022

HP 32SII: Some Algorithms For RPN Calculators

HP 32SII: Some Algorithms For RPN Calculators

Four programs ported to the HP 32SII calculator from algorithms designated for the 1973 HP 45 calculator.

HP 32SII: Euclid Algorithm - Greatest Common Divisor (GCD)

The HP 45 algorithm is found on page 228 in the Algorithms For RPN Calculators book. (see source below)

This algorithm takes up three labels.

E01 LBL E

E02 INPUT M

E03 ENTER

E04 ENTER

E05 ENTER

E06 INPUT N

E07 x<>y

K01 LBL 1

K02 ÷

K03 FP

K04 ×

K05 1

K06 x>y?

K07 GTO L

K08 R↓

K09 ENTER

K10 ENTER

K11 R↓

K12 R↓

K13 GTO K

L01 LBL L

L02 R↓

L03 R↓

L04 RTN

Sizes and Checksums:

E: 10.5 bytes, 9D4D

K: 19.5 bytes, F8AD

L: 6.0 bytes, C304

Total: 36.0 bytes

Instructions:

Press [ XEQ ] E, enter M and N.

Examples:

Input: M = 36, N = 28. Result: 4

Input: M = 48, N = 126. Result: 6

Input: M = 115, N = 300. Result: 5

HP 32SII: GCD Using One Label - John Kenney

The program was provided by Ross Barnes, and the algorithm is from the book ENTER by J. Daniel Dodlin and Keith Jarrett (ISBN 0-9615174-2-1, pg. 84). This is smart, one label program.

Enter both numbers in the stack before running the program.

G01 LBL G

G02 ENTER

G03 ENTER

G04 -

G05 R↓

G06 x<>y

G07 LASTx

G08 /

G09 LASTx

G10 RDN

G11 IP

G12 x

G13 -

G14 x≠0?

G15 GTO G

G16 +

G17 RTN

Size and Checksum: 25.5 bytes, 4E39

Posted with permission.

HP 32SII: Tetens Equation

The HP 45 algorithm is found on page 290 in the Algorithms For RPN Calculators book. The original algorithm took the temperature in Celsius.

Find the saturation of water vapor (e_m) in mmHg (millimeters of Mercury) given the temperature in °F.

Determined Formulas:

T (in °C) = (T°F - 32) * 5/9

α = T/(236.87 + T)

e_m = 4.579 * 10^(7.49 * α)

T01 LBL T

T02 INPUT T

T03 →°C

T04 ENTER

T05 ENTER

T06 236.87

T07 +

T08 ÷

T09 7.49

T10 ×

T11 10^x

T12 4.579

T13 ×

T14 RTN

Size and Checksum:

45.0 bytes, 404A

Examples:

T = 68 °F, Result: 17.53658 mmHg

T = 99 °F, Result: 47.63501 mmHg

HP 32SII: Dew Point Given Relative Humidity and Air Temperature

The HP 45 algorithm is found on page 290 in the Algorithms For RPN Calculators book. The original algorithm took the temperature in Celsius.

Relativity humidity (F) is to be entered as a decimal. For instance, instead of 20%, enter 0.20.

Determined Formulas:

T (in °C) = (T°F - 32) * 5/9

A = T/(T + 236.87)

B = 1/(log F/7.49 + A)

TD = 236.87/(B - 1)

TD = TD * 9/5 + 32

D01 LBL D

D02 INPUT T

D03 →°C

D04 ENTER

D05 ENTER

D06 236.87

D07 STO A

D08 +

D09 ÷

D10 INPUT F

D11 LOG

D12 7.49

D13 ÷

D14 +

D15 1/x

D16 1

D17 -

D18 RCL÷ A

D19 1/x

D20 →°F

D21 RTN

Size and Checksum:

47.5 bytes, 8677

Examples:

T = 80, F = 0.64, Result: 66.725

T = 95, F = 0.32, Result: 60.50684

HP 32SII: Effective Temperature Due to Wind Velocity

The HP 45 algorithm is found on page 291 in the Algorithms For RPN Calculators book. The original algorithm took the temperature in Fahrenheit.

Wind velocity is in miles per hour (mph).

Determined Formulas:

A = 0.634*(0.634 - log V)

ΔT = A*(T - 90)

Effective T = T - ΔT

E01 LBL E

E02 0.634

E03 ENTER

E04 ENTER

E05 INPUT V

E06 LOG

E07 -

E08 ×

E09 INPUT T

E10 ENTER

E11 90

E12 -

E13 ×

E14 RCL T

E15 x<>y

E16 -

E17 RTN

Size and Checksum:

33.5 bytes, 54F7

Examples:

V = 20 mph, T = 15 °F, Result: -16.71728

V = 15 mph, T = 86 °F, Result: 84.62526

Source:

Ball, John A. Algorithms For RPN Calculators John Wiley & Sons: New York, NY. 1978. ISBN 0-471-03070-8

For the second GCD program:

Dodin, J. Daniel and Keith Jarrett. ENTER: Reverse Polish Notation Made Easy Synthetix: Berkeley, CA ISBN 0-9612174-2-1 1984.

Special thanks and gratitude to Ross Barnes.

Saturday, November 19, 2022

Casio fx-9750GIII: Integrals with Infinite Limits

Casio fx-9750GIII: Integrals with Infinite Limits

A Substitution to Get to Infinity...

It is quite a challenge to calculate numerical integrals with infinite limits such as

∫( f(x) dx, a, ∞)

∫( f(x) dx, -∞, a)

∫( f(x) dx, -∞, ∞)

A trick is to substitute x = tan Θ. Then:

Θ = arctan x

dx = sec^2 Θ dΘ = cos^-2 Θ dΘ

Note that

lim t→∞ arctan t = π/2

lim t→-∞ arctan t = -π/2

In this blog, assume that radian angle mode is used in all calculations.

Let's go over some examples and see how it works. I used a Casio fx-9750GIII, however, this technique should work with all calculators with numerical integral calculations. Furthermore, changing the integral will allow for Simpson's Rule or Trapezoid Rule approximation.

For the upper limit, I approximate π/2.

A = 1.5708 (π/2 to 4 places)

B = 1.570796 (π/2 to 6 places)

C = 1.57079633 (π/2 to 8 places)

Example 1:

∫( e^(-x^2) dx, 0, ∞)

transforms to

∫( e^(-tan^2 Θ)/cos^2 Θ dΘ, 0, ≈π/2)

Calculations (fx-9750GIII):

Upper Limit: A (see above), Result: 0.8862269255

Upper Limit: B, Result: 0.8862269255

Upper Limit: C, Result: 0.8862269255

Example 2:

∫(x^0.5 * e^(-x) dx, 0, ∞)

transforms to

∫((tan Θ)^0.5 * e^(-tan Θ))/cos^2 Θ dΘ, 0, π/2)

Calculations:

Upper limits A, B, C: 0.862269255

Coincidently, ∫( e^(-x^2) dx, 0, ∞) = ∫(x^0.5 * e^(-x) dx, 0, ∞) = Γ(1.5) = √(π)/2

Example 3:

∫( e^x*(x^2 + 1) dx, -∞, 0)

transforms to

∫( e^(tan Θ) * (tan^2 Θ + 1)/cos^2 Θ dΘ, -π/2, 0)

= ∫( e^(tan Θ) * 1/cos^2 Θ * 1/cos^2 Θ dΘ, -π/2, 0)

∫( e^(tan Θ)/cos^4 Θ dΘ, -π/2, 0)

For upper limits A, B, C, the answer returned is 3, which is the exact answer.

For integrals like this all is needed is a four digit approximation of π/2 = 1.5708.

Let's keep going:

Example 4:

∫( (x^2 + 1)/(x^4 - 1) dx, 0, ∞)

transforms to

∫( 1/cos^Θ * (tan^2 Θ + 1)/(tan^4 Θ - 1) dΘ, 0, π/2)

= ∫( 1/(sin^4 Θ - cos^4 Θ) dΘ, 0, π/2)

On the fx-9750GIII get MA Error.

Example 5:

∫( 1/√(x^2 + 3*x + 1) dx, 0, ∞)

transforms to

∫( 1/(cos^2 Θ * √(tan^2 Θ + 3 * tan Θ + 1)) dΘ, 0, π/2)

like example 4, I get the MA Error.

Observation: the transformation works best if the integral involves some form of either of the following:

f(x) * e^(g(x))

f(x) * e^(-g(x))

Gamma

The Gamma Function is an excellent candidate for this transformation.

Γ(t) = ∫( x^(t-1) * e^-x dx, 0, ∞)

transforms to

∫( tan^(t-1) Θ * e^(-tan Θ)/cos^2 Θ dΘ, 0, π/2)

Source

Mier-Jedrzejuwicz, W.A.C. Ph.D. Tips And Programs for the HP 32S Synthetix Publication. Berkeley, CA. September 1988. ISBN 0-937637-05-X

Download the book here: https://literature.hpcalc.org/items/1756

Happy calculating,

Eddie

Sunday, November 13, 2022

HP 15C: Weibull Distribution Calculations

HP 15C: Weibull Distribution Calculations

Introduction

The Weibull probability density distribution function is:

f(x) = (b / Θ) * (x / Θ)^(b-1) * exp(-(x / Θ)^b)

with the lower tail cumulative distribution of (-∞ to x):

Area = 1 - exp(-(x / Θ)^b)

The area function tells us what is the probability a device lasts no more than x time units.

Area = 1 - Survival

The survival function is the probability a device lasts more than x time units.

Survival = exp(-(x / Θ)^b)

Generally, the higher Θ is, the flatter the Weibull Distribution curve.

What follows are four calculations regarding the Weibull Distribution. In the following programs, store the following values first prior to running the programs:

R0 = x

R1 = b

R2 = Θ

Use whatever labels you like.

HP 15C Program: Lower Tail Probability - Weibull Distribution

CDF = 1 - exp(-(x/Θ)^b)

Keys:

LBL B

RCL 0

RCL÷ 2

RCL 1

y^x

CHS

e^x

RTN

Key Codes:

42, 21,12

45, 0

45, 10, 2

45, 1

43, 32

Example:

b = 1.96, Θ = 420

x = 300, result: 0.4038

x = 400, result: 0.5970

x = 500, result: 0.7552

HP 15C Program: Failure Rate - Weibull Distribution

FR = b/Θ * (x/Θ)^(b-1)

Keys:

LBL C

RCL 1

RCL÷ 2

RCL 0

RCL÷ 2

RCL 1

y^x

RTN

Key Codes:

42, 21, 13

45, 0

45, 10, 2

45, 0

45, 10, 2

45, 1

43, 32

Example:

b = 1.96, Θ = 420

x = 300, result: 0.0034

x = 400, result: 0.0045

x = 500, result: 0.0055

HP 15C Program: Mean of a Weibull Distribution

µ = (1/b)! * Θ

Keys:

LBL D

RCL 1

1/x

RCL× 2

RTN

Key Codes:

42, 21, 14

45, 1

42, 0

45, 20, 2

43, 32

Example:

b = 1.96, Θ = 420

Result: 373.3720

HP 15C Program: Standard Deviation of a Weibull Distribution

σ = Θ * √((2/b)! - (1/b)!^2)

Keys:

LBL E

RCL 1

1/x

x^2

√

RCL× 2

RTN

Key Codes:

42, 21, 15

45, 1

42, 0

45, 1

42, 0

43, 11

45, 20, 2

43, 32

Example:

b = 1.96, Θ = 420

Result: 198.2208

Sources:

HP55 Statistics Programs Hewlett Packard Company. Cupertino, CA. 1975

Ma, Dan. "The Weibull distribution" Topics in Actuarial Modeling. September 28, 2016. https://actuarialmodelingtopics.wordpress.com/2016/09/28/the-weibull-distribution/ Last Retrieved September 20, 2022.

Eddie

Saturday, November 12, 2022

HP Prime: Reversing an Integer's Digits

HP Prime: Reversing an Integer's Digits

(Inspired by the HHC 2022 programming contest)

What Should I Add To Reverse the Digits?

Let A, B, C, D, and E be individual digits (0-9) of an integer. AB would represent a two digit integer with the value of 10 * A + B. ABC would represent a three digit integer with the value of 100 * A + 10 * B + C.

Reversing a Two Digit Integer

AB + # = BA

10 * A + B + # = 10 * B + A

# = 9 * (B - A)

Example: Let AB = 76.

A = 7, B = 6

# = 9 * (6 - 7) = -9

76 - 9 = 67

Reversing a Three Digit Integer

ABC + # = CBA

100 * A + 10* B + C + # = 100 * C + 10 * B + A

# = 99 * (C - A)

Example: ABC = 469

# = 99 * (9 - 4) = 495

469 + 495 = 964

Reversing a Four Digit Integer

ABCD + # = DCBA

1000 * A + 100 * B + 10 * C + D + # = 1000 * D + 100 * C + 10 * B + A

# = 999 * (D - A) + 90 * (C - B)

Example: ABCD = 7219

# = 999 * (9 - 7) + 90 * (1 - 2) = 1908

7219 + 1908 = 9127

Reversing a Five Digit Integer

ABCDE + # = EDBCA

10000 * A + 1000 * B + 100 * C + 10 * D + E + # =

10000 * E + 1000 * D + 100 * C + 10 * B + A

# = 9999 * (E - A) + 990 * (D - B)

Example: ABCDE = 52693

# = 9999 * (3 - 5) + 990 * (9 - 2) = -13068

52693 - 13068 = 39625

Having the Calculator Do It

The program REVINT reverses the digits of an integer, up to 11 digits. The program does not allow numbers that have non-zero fractional parts or integers more than 11 digits. Instead of solving for # (see above), the program splits the integers into a list in reverse order, and uses list processing to get the final answer.

HP Prime Program: REVINT

Caution: Integers that end or begin with zero may not return accurate results. My suggestion is not use 0s with this program. See examples below for more details.

EXPORT REVINT(N)

BEGIN

// 2022-09-18 EWS

// reverse the integer N

// up to 12 digits

LOCAL D,P,A,I,M,L;

L:={};

P:=XPON(N);

// check size

IF P>11 THEN

RETURN "TOO BIG";

KILL;

END;

// check type

IF FP(N) THEN

RETURN "NOT AN INTEGER";

KILL;

END;

D:=N;

// loop

FOR I FROM P DOWNTO 0 DO

A:=D/ALOG(I);

L:=CONCAT({IP(A)},L);

D:=D-IP(A)*ALOG(I);

END;

// rebuild

M:=ΣLIST(MAKELIST(ALOG(X),X,P,0,−1)*L);

RETURN M;

END;

Examples:

REVINT(4321) returns 1234

REVINT(56765) returns 56765 (56765 is a palindrome, reversing the digits results in the same number)

REVINT(42910) returns 1924 (01924 - be aware about integers ending or beginning with 0)

REVINT(67.28) returns "NOT AN INTEGER" (error)

Eddie

Monday, November 7, 2022

Retro Review: Texas Instruments TI SR-51-II

Retro Review: Texas Instruments TI SR-51-II

Quick Facts

Model: SR-51-II

Company: Texas Instruments

Years: 1976

Type: Scientific

Batteries: 1 × BP6, with AC-Adapter**

Operating Modes: AOS

Number of Registers: 3

Display: 10 digit - red seven segment digits

Leather Case

** I purchased the calculator on eBay not to long ago from night_owl_items, who made it so the SR-51 II can run on two AA batteries, making it one of the kind and extend it's life.

Features

* Trigonometric functions and inverses

* Hyperbolic functions and inverses

* Logarithms and exponentials

* Factorials of positive integers

* Percent and percent change

* Arithmetic, powers, and roots

Conversions

There are eight conversions:

* in · mm

* gal · l

* lb · kg

* °F · °C

* D · R (degrees - radians)

* G · R (grads - radians)

* DMS · DD (degrees-minutes-seconds - decimal degrees)

* P → R (polar - rectangular conversions)

The [ 2nd ] key sequence converts left to right

The [ 2nd ] [ INV ] key sequence converts right to left

Angle Modes

The SR-51-II has the three customary angle modes: degrees, radians, and grads. Unfortunately there is no angle mode indicator. The way to:

Degrees Mode: 0.5 sin^-1 returns 30

Radians Mode: 0.5 sin^-1 returns .5235987756

Grads Mode: 0.5 sin^-1 returns 33.33333333

Storage and Storage Registers

The SR-51-II has three memory registers. The calculator also has storage arithmetic:

SUM: STO+

INV SUM: STO-

PROD: STO×

INV PROD: STO÷

The CONST Function

The CONST is an access to common scientific functions, but allows the user to do a short arithmetic calculation, includes roots, powers, and percent change.

Example: Subtract 100 from 290, 370, and 480 respectively.

100 [ - ] [ 2nd ] (CONST)

290 [ = ] returns 190

370 [ = ] returns 270

480 [ = ] returns 380

Example 2: Multiply 172, 288, and 382 by 6 respectively.

6 [ × ] [ 2nd ] (CONST)

172 [ = ] returns 1032

288 [ = ] returns 1728

382 [ = ] returns 2292

Statistics

The statistics modes have one and two variable statistics. What is unique of the SR-51-II is that the y-variable is primary. However, entering two variable statistics remain the usual: x [ x<>y ] y [ Σ+ ].

[ 2nd ] (MEAN): mean of y data; [ INV ] [ 2nd ] (MEAN): mean of x data

[ 2nd ] (S.DEV): standard deviation of y data;

[ INV ] [ 2nd ] (S.DEV): standard deviation of x data

[ 2nd ] (VAR): population variance of y data;

[ INV ] [ 2nd ] (VAR): population variance of x data

Solid calculator. I wish the SR-51-II had a random number function and angle mode indicator.

Source:

Woerner, Joerg. "Texas Instruments SR-51-II" Datamath. December 5, 2001. Last Retrieved November 1, 2022. http://www.datamath.org/Sci/MAJESTIC/sr-51-II.htm

Eddie

Sunday, November 6, 2022

hp 9g: Sample Programs

hp 9g: Sample Programs

Introduction

As I mentioned in yesterday's blog, today's blog will have some sample programs for the hp 9g.

The hp 9g's programming language is a mix of C and Basic, and is unique from the rest of the Hewlett Packard calculators. The hp 9g has 400 steps of programming memory to be allocated among 10 program files, P0 through P9.

The 9g has 26 variables, but can be extended up to 53 variables at the cost of program steps. The square brackets indicate indirect addressing. The number in the square brackets is the number of steps away from the variable. Any variable beyond Z (Z[0], A[25]) must be addressed in this fashion.

Example:

A[0] calls up A

A[1] calls up B

...

A[25] calls up Z.

We can use any letter, but I recommend using either A or Z.

hp 9g Program: Integer Division

This program gives the quotient and remainder of N ÷ D. The program here assumes you enter positive integers, and some code can be entered for checking inputs.

INPUT N

INPUT D

Q=INT(N/D)

R=N-Q*D

PRINT Q," R",R

END

Examples:

N = 4379, D = 8; Result: 547 R3

N = 5884, D = 29; Result: 202 R26

Note: The PRINT command converts variable values to strings automatically, without the need for a string or STR$ command.

hp 9g Program: Arithmetic-Geometric Mean

This program determines the AGM given two numbers A and G.

INPUT A

INPUT B

PRINT "ACC. 10^-8"

Lbl 0:

X=.5(A+G)

Y=√(AG)

If (ABS(X-Y)<10^(-8))

Then {Goto 1}

A=X

G=Y

Goto 0

Lbl 1:

Print "ANS= ", X

Examples:

A = 1.3, G = 1.5, Result: ANS ≈ 1.39821143

A = 20, G = 45, Result: ANS ≈ 31.23749374

Notes:

10^ is from [ 2nd ] (10^x) the antilog function.

The label command must have the colon character after the label number, or a syntax error occurs.

If you run a program from the Main mode by the [ PROG ] key, you will be transferred to Program mode when the program finishes.

hp 9G Program: The Maximum Value of A through X

Y = maximum value

Z = counter variable

PRINT "MAX(A:X)"; ◢

FOR(Z=2;Z<23;Z++){

IF(A[Z]>Y)

THEN{Y=A[Z]}};

PRINT "MAX=",Y;

END

Example:

Clear all variable in the MAIN mode by CL-VAR

Set up variables:

A = 12, D = 24, E = 37, G = 40, S = 19, T = 16, U = 7

Result: MAX= 40

Notes:

++ is from the Instruction menu, it is a smaller-jointed double plus.

When the program encounters a run/stop instruction (◢), the last message is displayed. Continue execution by pressing the equals key [ = ].

hp 9g Program: Plotting y = A * sin(B * x + C) + D

The programming mode does not offer a lot of support when it comes to graphs.

PRINT "SINE"

SLEEP(.5)

INPUT A,B,C,D

RANGE(-2π,2π,π/4,-A,A,1)

Graph Y=A*sin(BX+C)+D

PRINT "PRESS G<>T"

END

A = 1, B = 1, C = 0, D = 0

Notes:

Prior to running the program, set the angle mode to Radians through the [DRG] key. Also, press [ 2nd ] [ → ] (CLS) to clear the graph screen. Unfortunately, neither of these commands can be programmed.

Plots from programs cannot be traced.

When the program ends, you will need to switch to Main mode, press [G<>T] to show the graph. Not a very efficient way of showing a graph.

That is a sample of the programs for the hp 9g. Good for number crunching, kind of a little bit to be desired for graphics.

Until next time,

Eddie