HP 65 Programs: Triangulation, Percentile, Roots of Unity, Partial Fractions

 HP 65 Programs: Triangulation, Percentile, Roots of Unity, Partial Fractions

Triangulation

(done with HP-65 Emulator for Windows, Bernhard Emese)

Given:

Side angle α

Side angle ß

r = Length of base

Find:

h = height which crosses the base line at right angles

Code:

23: LBL

12: B

35: g

41: DEG

33 01: STO 1

35 08: R↓

41: ENTER

41: ENTER

44: CLX

51: -

31: f

04: SIN

35 09: R↑

31: f

04: SIN

71: *

33: STO

71: *

01: 1

35 08: R↓

61: +

31: f

04: SIN

33: STO

81: ÷

01: 1

34 01: RCL 1

24: RTN

Example: ß = 50°, α = 30°, r = 10

Result: h ≈ 3.89

Source:

“Triangulation (surveying)” Wikipedia. https://en.wikipedia.org/wiki/Triangulation (surveying) (last edited October 24, 2025). Retried March 17, 2026

Percentile in a Range

The program calculates the percentile of x in the range [a, b].

percentile = (x – a) ÷ (b – a) * 100%

Stack:

Z: x

Y: a

X: b

23: LBL

11: B

35 07: x<>y

51: -

35 00: LSTx

35 07: R↓

51: -

35 09: R↑

81: ÷

02: 2

32: f^-1

08: LOG (10^x)

71: ×

24: RTN

Example: The range: [210, 470], x = 376

Stack: Z: 376, Y: 210, X: 470

Result: 63.85

Roots of Unity

(done with HP-65 Emulator for Windows, Bernhard Emese)

w^n = 1

Registers used: R1 = n, R8 = n (counter)

n [ A ]

Cycle:

root # [ R/S ]

real part [ R/S ]

imaginary part [ R/S ] (the cycle starts again)

Repeat the cycle until the display is flashing zeros (forced 1/0 error). Hit [ CLx ] to stop the flashing.

This program switches the mode to radians angle measurement.

Code:

23: LBL

11: A

33 01: STO 1

33 08: STO 8

35: g

42: RAD

23: LBL

00: 0 (LBL 0 - subroutine)

34 08: RCL 8

84: R/S (show root #)

35: g

02: π

71: ×

02: 2

71: ×

34 01: RCL 1

81: ÷

01: 1

32: f^-1

01: R->P (->rect)

84: R/S (show real part)

35 07: x<>y

84: R/S (show imaginary part)

35: g

83: DSZ (decrement R8 by 1)

22: GTO

00: 0 (GOTO LBL 0)

34 08: RCL 8

35: g

04: 1/x

24: RTN

Example (FIX 4):

w^3 = 1

#; R/S; real part; R/S; imaginary part

3; 1.0000; -0.0000

2; -0.5000; -0.8660

1; -0.5000; 0.8660

(flashing zeroes)

Partial Fractions

(w is the variable, and T, Z, Y, and X are the stack contents)

(Tw+Z) ÷ ((w+Y)(w+X)) → B÷(w+Y) + A÷(w+X)

Input Stack: T, Z, Y, X 

Output Stack: B, A

Code: 

23: LBL 

11: A 

35, 07: x<>y 

51: - 

35, 00: LST x 

35, 09: R↑ 

71: × 

35, 00: LST x 

35, 07: x<>y 

35, 09: R↑ 

35, 07: x<>y 

51: - 

35, 09: R↑ 

81: ÷ 

41: ENTER 

35, 08: R↓ 

51: - 

35, 09: R↑ 

35, 07: x<>y 

24: RTN

Example 1: 

(12*w + 28) ÷ ((w+2) * (w+3)) = 4 ÷ (w+2) + 8 ÷ (w+3) 

Input Stack: T: 12, Z: 28, Y: 2, X: 3 press [ A ] 

Output Stack: Y (B): 4, X (A): 8

Example 2: 

(33) ÷ ((w-5) * (w+6)) = 3 ÷ (w-5) + -3 ÷ (w+6) 

Input Stack: T: 0, Z: 33, Y: -5, X: 6 press [ A ] 

Output Stack: Y (B): 3, X (A): -3

Example 3: 

(2*w - 3) ÷ ((w-1) * (w+4)) = -0.2 ÷ (w-1) + 0.2 ÷ (w+4) 

Input Stack: T: 2, Z: -3, Y: -1, X: 4 press [ A ] 

Output Stack: Y (B): -0.2, X (A): 0.2

Eddie

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

TI-84 Plus CE and Numworks: System of Linear Differential Equations (Simple System)

TI-84 Plus CE and Numworks: System of Linear Differential Equations (Simple System) 

Solving  The System

The program presented today will solve the following system of equations:

dx/dt = A * x(t) + B * y(t)

dy/dt = C * x(t) + D * y(t)

with initial conditions x(0) = x0 and y(0) = y0

The formulas used in the program are derived using the Laplace Transform.

dx/dt = A * x(t) + B * y(t)

dy/dt = C * x(t) + D * y(t)

Laplace Transforms to be used:

ℒ{f(t)} = F(s)

ℒ{df/dt} = s * F(s) - f(0)

Apply the Laplace Transform transform:

s * X(s)  - x0 = A * X(s) + B * Y(s)

s * Y(s) - y0 = C * X (s) + B * Y(s)

(s - A) * X(s) + (-B) * Y(s) = x0

(-C) * X(s) + (s - D) * Y(s) = y0

Define three matrices:

Md = [ [ s - A, -B ], [ -C, s - D ] ]

where

det(Md) = s^2 - (A + D) * s + (A * D - B * C)

which can be factored into:

(s - P) * (s - Q)

where:

P = (( A + D ) + √(( A + D)^2 - 4 * (A * D - B * C)) /2

Q = (( A + D ) - √(( A + D)^2 - 4 * (A * D - B * C)) /2

Mx = [ [ x0, -B ], [ y0, s - D ] ] 

where

det(Mx) = x0 * (s - D) + B * y0

My = [ [ s - A, x0 ], [ - C, y0 ] ]

where 

det(My) = (s - A) * y0 + C * x0

Per Cramer's Rule, the solution to X(s) and Y(s) are:

X(s) = det(Mx) / det(Md)

Y(s) = det(My) / det(Md)

What happens next depends of the roots P and Q.  We're assuming that P and Q are real here but the situation applies if P and Q are complex. 

If P ≠ Q:

If P ≠ Q, then after simplifying and finding partial fractions:

X(s) = ( x0 * (s - D) + B * y0 ) / ((s - P) * (s - Q)) = α / (s - P) + β / (s - Q)

where

α = ( B * y0 + (P - D) *x0) / (P - Q)

β = x0 - α

Y(s) = ( (s - A) * y0 + C * x0 ) / ((s - P) * (s - Q)) = γ / (s - P) + δ / (s - Q)

where

γ =(C * x0 + y0 * (P - A)) / (P - Q)

δ = y0 - γ

Applying the inverse Laplace transform:

x(t) = α * e^(P * t) + β * e^(Q * t)

y(t) = γ * e^(P * t) + δ * e^(Q * t)

If P = Q:

If P = Q, then after simplifying and finding partial fractions:

X(s) = ( x0 * (s - D) + B * y0 ) / (s - P)^2 = α / (s - P)^2 + β / (s - Q)

where

α = B * y0 + (P - D) * x0 

β = x0

Y(s) = ( (s - A) * y0 + C * x0 ) / (s - P)^2 = γ / (s - P)^2 + δ / (s - Q)

where

γ =C * x0 + y0 * (P - A)

δ = y0

Applying the inverse Laplace transform:

x(t) = α * t * e^(P * t) + β * e^(Q * t)

y(t) = γ * t * e^(P * t) + δ * e^(Q * t)

TI-84 Plus CE Program:  LINSYS1

Download the program here:  https://drive.google.com/file/d/14zYpT4xhLRGVvBNriOacIOGBpqXLWswR/view?usp=share_link

Program Listing:

Radian

a+bi

ClrHome

Disp "DX/DT=AX+BY DY/DT=CX+DY"

Prompt A,B,C,D

Input "X0? ",M

Input "Y0? ",N

((A+D)+√((A+D)²-4*(A*D-B*C)))/2→P

((A+D)-√((A+D)²-4*(A*D-B*C)))/2→P

If P≠Q

Then

(B*N-D*M+P*M)/(P-Q)→R

M-R→S

(C*M-A*N+P*N)/(P-Q)→U

N-U→V

ClrHome

Output(2,1,"X=Re^(PT)+Se^(QT)")

Output(3,1,"R= "+toString(R))

Output(4,1,"P= "+toString(P))

Output(5,1,"S= "+toString(S))

Output(6,1,"Q= "+toString(Q))

Pause

ClrHome

Output(2,1,"Y=Ue^(PT)+Ve^(QT)")

Output(3,1,"U= "+toString(U))

Output(4,1,"P= "+toString(P))

Output(5,1,"V= "+toString(V))

Output(6,1,"Q= "+toString(Q))

Pause

Else

B*N-D*M+M*P→P

M→S

C*M-A*N+N*P→U

N→V

ClrHome

Output(2,1,"X=RTe^(PT)+Se^(PT)")

Output(3,1,"R= "+toString(R))

Output(4,1,"P= "+toString(P))

Output(5,1,"S= "+toString(S))

Pause

ClrHome

Output(2,1,"Y=UTe^(PT)+Ve^(PT)")

Output(3,1,"U= "+toString(U))

Output(4,1,"P= "+toString(P))

Output(5,1,"V= "+toString(V))

Pause

End

ClrHome

Python:  dfsys1.py

Link from my Numworks page:  https://my.numworks.com/python/ews31415/dfsys1

Modules: 

cmath:  complex math

mathplotlib.pyplot:  screen plotting and text placement

Script listing:  

from cmath import *

from matplotlib.pyplot import *

# 2023-01-28 EWS

# parameter input

print("systems of differential\nequations")

print("dx/dt=a*x+b*y")

print("dy/dt=c*x+d*y")

a=float(input("a? "))

b=float(input("b? "))

c=float(input("c? "))

d=float(input("d? "))

m=float(input("x0? "))

n=float(input("y0? "))

# characteristic roots

p=((a+d)+sqrt((a+d)**2-4*(a*d-b*c)))/2

q=((a+d)-sqrt((a+d)**2-4*(a*d-b*c)))/2

# set up result screen

axis((0,10,0,10))

axis("off")

# determine solutions

# best used for real solutions

if p!=q:

  r=(b*n+m*(p-d))/(p-q)

  s=m-r

  u=(c*m+n*(p-a))/(p-q)

  v=n-u

  text(0,9,"x=re**(pt)+se**(qt)")

  text(0,8.5,"r="+str(r))

  text(0,8,"p="+str(p))

  text(0,7.5,"s="+str(s))

  text(0,7,"q="+str(q))

  text(0,6,"y=ue**(pt)+se**(vt)")

  text(0,5.5,"u="+str(u))

  text(0,5,"p="+str(p))

  text(0,4.5,"v="+str(v))

  text(0,4,"q="+str(q))

else:

  r=(b*n+m*(p-d))

  s=m

  u=(c*m+n*(p-a))

  v=n

  text(0,9,"x=rte**(pt)+se**(pt)")

  text(0,8.5,"r="+str(r))

  text(0,8,"p="+str(p))

  text(0,7.5,"s="+str(s))

  text(0,6,"y=ute**(pt)+se**(vt)")

  text(0,5.5,"u="+str(u))

  text(0,5,"p="+str(p))

  text(0,4.5,"v="+str(v))

show()

Examples

Example 1

dx/dt = 3 * x(t) - 2 * y(t)

dy/dt = 2 * x(t) - 2 * y(t)

x(0) = 6, y(0) = 9

Solutions:

x(t) = 2 * e^(2*t) + 4 * e^(-t)

y(t) = e^(2*t) + 8 * e^(-t)

Example 2

dx/dt = 3 * x(t) - 4 * y(t)

dy/dt = x(t) - y(t)

x(0) = 1, y(0) = 0

Solutions:

x(t) = 2 * t * e^(t) + e^(t)

y(t) = t * e^(t) 

Eddie 

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

Casio fx-3650P and fx-CG50: Partial Fraction Decomposition

Casio fx-3650P and fx-CG50:   Partial Fraction Decomposition

Casio fx-3650P (II) Program: Partial Fraction Decomposition

The following program calculates the coefficients decomposition, X and Y, as follows:

(A ∙ t + B) ÷ ((t + C) ∙ (t + D)) = X ÷ (t + C) + Y ÷ (t + D)

Note:  C ≠ D

where:

X = (B - A ∙ C) ÷ (D - C)

Y = A - X

Program:  (35 Steps)

? → A : ? → B : ? → C : ? → D : 

( B - A C ) ÷ ( D - C → X ◢ A - X → Y

Example:  

(2t - 3) ÷ ((t - 1)(t + 4))  = -0.2÷(t-1) + 2.2÷(t + 4)

A = 2, B = -3, C = -1, D = 4

Result:  -0.2, 2.2

The above algorithm should be good for the fx-7000G, fx-50FII, fx-4000P, fx-5800P, and other programming and graphing calculators (modifications may be necessary).

Casio fx-CG50 Program:  PARTFRAC

The code listed will also be good for the fx-9750G/fx-9860G series.

The program PARTFRAC, 720 bytes, will execute partial fraction decomposition of any of the three types:

(1)    (A ∙ x + B)÷((x + C) ∙ (x + D)) = R÷(x + C) + S÷(x + D)

Inputs:  A, B, C, D

Outputs:  R, S

[[ R ] [ S ]] = [ [ D, C ] [ 1, 1 ] ]^(-1) ∙ [ [ B ] [ A ] ]

(2)  (A ∙ x^2 + B ∙ x + C)÷((x + D)^2 ∙ (x + E)) = R÷(x + D)^2 + S÷(X + D) + T÷(x + E)

Inputs:  A, B, C, D, E

Outputs:  R, S, T

[[ R ][ S ][ T ]] = [[ E, D ∙ E, D^2 ][ 1, D + E, 2 ∙ D ][ 0, 1, 1 ]]^-1) ∙ [[ C ][ B ][ A ]] 

(3)  (A ∙ x^2 + B ∙ x + C)÷((x+D)∙(x+E)∙(x+F)) = R÷(x+D) + S÷(x+E) + T÷(x+F)

Inputs:  A, B, C, D, E, F

Outputs:  R, S, T

[[ R ][ S ][ T ]] = 

[[ E ∙F, D ∙ F, D ∙ E ][ E+F, D+F, D+E ][ 1, 1, 1 ]]^-1) ∙ [[ C ][ B ][ A ]] 

Program Listing:

'ProgramMode:RUN

"2021-11-05 EWS"

"PARTIAL FRACTION"

"N(X)/D(X)"◢

Lbl 0

Menu "DENOMINATOR?","(X+C)(X+D)",1,"(X+D)_^<2>_(X+E)",2,"(X+D)(X+E)(X+F)",3,"EXIT",E

Lbl 1

"(AX+B)/((X+C)(X+D))"◢

"A"?->A

"B"?->B

"C"?->C

"D"?->D

[[D,C][1,1]]^<-1>*[[B][A]]->Mat A

Mat A[1,1]->R

Mat A[2,1]->S

"=R/(X+C)+S/(X+D)"◢

"R="

R◢

"S="

S◢

Goto 0

Lbl 2

"(AX_^<2>_+BX+C)/((X+D)_^<2>_(X+E))"◢

"A"?->A

"B"?->B

"C"?->C

"D"?->D

"E"?->E

[[E,D*E,D^<2>][1,D+E,2*D][0,1,1]]^<-1>*[[C][B][A]]->Mat A

Mat A[1,1]->R

Mat A[2,1]->S

Mat A[3,1]->T

"=R/(X+D)_^<2>_+S/(X+D)+T/(X+E)"◢

"R="

R◢

"S="

S◢

"T="

T◢

Goto 0

Lbl 3

"(AX_^<2>_+BX+C)/((X+D)(X+E)(X+F))"◢

"A"?->A

"B"?->B

"C"?->C

"D"?->D

"E"?->E

"F"?->F

[[E*F,D*F,D*E][E+F,D+F,D+E][1,1,1]]^<-1>*[[C][B][A]]->Mat A

Mat A[1,1]->R

Mat A[2,1]->S

Mat A[3,1]->T

"=R/(X+D)+S/(X+E)+T/(X+F)"◢

"R="

R◢

"S="

S◢

"T="

T◢

Goto 0

Lbl E

"DONE."

Download (fx-CG 50):

https://drive.google.com/file/d/1hAIFrlX6xLMYPRiLz_dt7hpiP64CxmKC/view?usp=sharing

Examples To Try:

Type 1: 

(8 ∙ x + 3) ÷((x + 4) ∙ (x - 6))  = (29/10) ÷ (x + 4) + (51/10) ÷ (x - 6)

Type 2:

(x^2 + 4 ∙x + 5) ÷ ((x - 3)^2 ∙ (x + 1)) = 

(13/2) ÷ (x - 3)^2 + (7/8) ÷ (x - 3) + (1/8) ÷ (x + 1)

Type 3:

(2 ∙ x^2 + 3 ∙ x - 1) ÷ ((x - 1) ∙ (x + 2) ∙ (x + 4)) = 

(4/15) ÷ (x - 1) + (-1/6) ÷ (x + 2) + (19/10) ÷ (x + 4)

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.