Fun with the 71B '20

Fun with the 71B '20

Differential Equations:  Runge Kutta Method 4th Order

Find a numerical solution to the differential equation:

dy/dx = f(x, y)

x is the independent variable, y is the dependent variable.  You define f(x,y) on line 10.   For example, dy/dx = sin(x * y) should have this as line 10:

10 DEF FNF(X,Y) = SIN(X*Y)

Line 5 is a remark line.  Remarks are followed by exclamation points on the HP 71B.

HP 71B Program: RK4 

Size: About 300 - 330 bytes

5 ! FNF(X,Y) = dY/dX
10 DEF FNF(X,Y) = [ insert f(x,y) here ]
15 DESTROY X,Y,H,K1,K2,K3,K4
20 INPUT "X0? "; X
25 INPUT "Y0? "; Y
30 INPUT "STEP? "; H
40 K1 = FNF(X,Y)
45 K2 = FNF(X+H/2,Y+H*K1/2)
50 K3 = FNF(X+H/2,Y+H*K2/2)
55 K4 = FNF(X+H,Y+H*K3)
60 X=X+H
65 Y=Y+H*(K1+2*K2+2*K3+K4)/6
80 DISP "(";X;",";Y;")" @ PAUSE
85 DISP "NEXT? Y/N"
90 A$=KEY$
95 IF A$="Y" THEN 40
99 IF A$="N" THEN DISP "DONE" @ END ELSE 85

Example:
dy/dx = sin(x * y) with inital condition y(0) = 0.5,  h = 0.2 * π

First three results:
( .628318530718 , .607747199386 )   ( [ f ] [ +] (CONT), [ Y ] )
( 1.25663706144, 1.02432288082 )
( 1.88495559216, 1.51038862362 ) 

Hyperbolic Functions

[ S ] = sinh(x)
[ C ] = cosh(x)
[ A ] = asinh(x)
[ H ] = acosh(x)
[ X ]  to exit

HP 71B Program: HYP

Size:  401 bytes
acosh(x) requires that | x | ≥ 1

100 DESTROY A,X
115 DISP "sinh S/A, cosh C/H, X"
120 A$=KEY$
125 IF A$="S" THEN INPUT "X? ";X @ CALL SINH(X)
130 IF A$="C" THEN INPUT "X? ";X @ CALL COSH(X)
135 IF A$="A" THEN INPUT "X? ";X @ CALL ASINH(X)
140 IF A$="H" THEN INPUT "X? ";X @ CALL ACOSH(X)
145 IF A$="X" THEN 150 ELSE 115
150 DISP "DONE" @ END
200 SUB SINH(X)
205 DISP (EXP(X)-EXP(-X))/2 @ PAUSE
210 END SUB
300 SUB COSH(X)
305 DISP (EXP(X)+EXP(-X))/2 @ PAUSE
310 END SUB
400 SUB ASINH(X)
405 DISP LOG(X+SQR(X^2+1)) @ PAUSE
410 END SUB
500 SUB ACOSH(X)
510 DISP LOG(X+SQR(X^2-1)) @ PAUSE
515 END SUB

Example:
X = 2.86

sinh(2.86) returns 8.70212908815
cosh (2.86) returns 8.75939784845
asinh(2.86) returns 1.77321957441
acosh(2.86) returns  1.71190019325

Arithmetic-Geometric Mean

The arithmetic-geometric mean (AGM) is found by the iterative process:

a = 0.5 * (x + y)
g = √(x * y)

The values of a and g are stored into x and y, respectively.  The process repeats until the values of a and g converge.   A tolerance of 10^(-9) is used to display an 8-digit approximation.

HP 71B Program: AGM

Size:  148 Bytes

10 DESTROY X,Y,A,B
15 DISP "AGM(X,Y)" @ WAIT .5
20 INPUT "X? ";X
25 INPUT "Y? ";Y
30 A=.5*(X+Y)
35 G=SQR(X*Y)
40 X=A
45 Y=G
50 IF ABS(X-Y)>1E-9 THEN 30
55 DISP USING 60;X
60 IMAGE 10D.8D    // (10 digit integer parts with rounding to 8 decimal places)
65 END

Example:
AGM(178, 136)

Result:  156.29380544

Pythagorean Triple Generator

Given two positive integers m, n; where m > n, a Pythagorean triple is generated with the following calculations:

a = 2*m*n
b = m^2 - n^2
c = m^2 + n^2

Properties:

a^2 + b^2 = c^2
Perimeter: p = a + b + c
Area: r = a * b / 2

HP 71B Program: PYTHTRI

Size: 217 bytes

10 DESTROY M,N,A,B,C,R,P
20 DISP "M>N, INTEGERS" @ WAIT .5
25 INPUT "M? "; M
30 INPUT "N? "; N
35 A=2*M*N
40 B=M^2-N^2
45 C=M^2+N^2
50 P=A+B+C
55 R=A*B/2
60 DISP 'A = ';A @ PAUSE
65 DISP 'B = ';B @ PAUSE
70 DISP 'C = ';C @ PAUSE
75 DISP 'PERIM.=';P @ PAUSE
80 DISP 'AREA =';R
85 END

Example:
M = 16, N = 11

Results:
A = 352, B = 153, C = 377, P = 864, R = 23760

Impedance of An Alternating Current

The program ALTCURR calculates the impedance (magnitude and phase angle) of a sinusoidal alternating current consisting of one resistor, one capacitor, and one inductor in a series.

HP 71B Program: ALTCURR

Size: 210 bytes

10 DESTROY F,L,C, R,W,Z,T
15 DEGREES
20 INPUT "FREQUENCY? ";F
25 INPUT "INDUCTANCE? ";L
30 INPUT "CAPACITANCE? ";C
35 INPUT "RESISTANCE? ";R
40 W=2*PI*F
45 Z=SQR(R^2+(W*L-1/(W*C))^2)
50 T=ATAN((W*L-1/(W*C))/R)
55 DISP "MAGNITUDE= "; Z @ PAUSE
60 DISP "PHASE ANGLE = "; T

Example:
F = 152 Hz
L = 4.75E-3 H  (4.75 mH)
C = 8E-6 F  (8 μF)
R = 6400 Ω

Results: 
Magnitude:  6401.24704262
Phase Angle: -1.1309750812°

Source:
Rosenstein, Morton.  Computing With the Scientific Calculator Casio.  Japan. 1986.  ISBN 1124161430

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.

TI-89 and HP Prime: Approximation Digits of π to Lots of Digits

TI-89 and HP Prime:  Approximation Digits of π to Lots of Digits 

Ways to find digits of π. 

Four steps with a calculator with CAS, which allows computation with large integers:

1.  Take a series and multiply each term by 10^n, where n is the number of digits desired. The exact form of each calculation 

2. Take the integer portion of the sum. Make any adjustments when necessary. 

3. Convert the sum to a string.  Insert a decimal point when appropriate. 

4. Return the answer as a string. 

There are three methods that are demonstrated, one routine is presented for the TI-89 and the HP Prime.  All programs presented takes two arguments: n number of digits, m is the number of iterations. Remember that m must be sufficiently high to get an accurate computation of π. 

In the TI-89 programs, (C) marks a comment.  All HP Prime programs on this blog entry are to be run in CAS mode (make sure CAS is checked when creating the program).  

Keep in mind the HP Prime is significantly faster than the TI-89.  

Newton/Euler - newpi

π = 2 * Σ ( 2^k * (k!)^2 / (2*k + 1)!, k, 0, infinity)

TI-89 Program newpi

newpi(n,m)

(C) Newton/Euler π 

(C) 2016-11-13

Local s,k,t

(C) initialize

0 → t

(C) loop

For k,0,m

t + exact( 10^n * 2^k * (k!)^2 / ( (2*k + 1)! ) ) → t

EndFor

iPart(2*t) → t

string(t) → s

left(s,1) & "." & right(s,n) → s

Disp s

EndPrgm

HP Prime CAS Program newpi

#cas

newpi(n,m):=

BEGIN

LOCAL t,k,s;

// Netwon/Euler π 

// 2016-11-13 EWS

t:=0;

FOR k FROM 0 TO m DO

t:= t + exact( 10^n * 2^k * (k!)^2 /  ( (2*k + 1)! ));

END;

t:= IP(2*t);

// string, left, right 

// must be in lower case

s:= string(t);

s:= left(s,1) + "." + right(s,n);

return s;

END;

#end

newpi(20, 175) returns "3.14159265358979323846"

Arctangent - arctanpi

π = Σ ( 2^(k+1) / ( COMB(2*k, k) * (2*k +1) ), k, 0, infinity)

TI-89 Program arctanpi

arctanpi(n,m)

(C) Arctan π 

(C) 2016-11-13

Local s,k,t

(C) initialize

0 → t

(C) loop

For k,0,m

t + exact( 10^n * 2^(k + 1) / ( nCr(2*k,k) * (2*k+1) ) ) → t

EndFor

iPart(t) → t

string(t) → s

left(s,1) & "." & right(s,n) → s

Disp s

EndPrgm

HP Prime CAS Program arctanpi

#cas

arctanpi(n,m):=

BEGIN

LOCAL t,k,s;

// Arctan π 

// 2016-11-13 EWS

t:=0;

FOR k FROM 0 TO m DO

t:= t + exact( 10^n * 2^(k+1)/(COMB(2*k,k) * (2*k + 1)));

END;

t:= IP(t);

// string, left, right 

// must be in lower case

s:= string(t);

s:= left(s,1) + "." + right(s,n);

return s;

END;

#end

arctanpi(25, 200) returns "3.1415926535897932384626433"

Bailey-Borwein-Plouffe Formula (BBP) - bbppi

π = Σ ( 1/16^k * ( 4/(8*k +1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k + 6) ), k, 0, infinity)

TI-89 Program bbppi

bbppi(n,m)

(C) BBP π 

(C) 2016-11-13

Local s,k,t

(C) initialize

0 → t

(C) loop

For k,0,m

t + exact( 10^n/16^k * ( 4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) - 1/(8*k + 6) ) ) → t

EndFor

iPart(t) → t

string(t) → s

left(s,1) & "." & right(s,n) → s

Disp s

EndPrgm

HP Prime CAS Program bbppi

#cas

bbppi(n,m):=

BEGIN

LOCAL t,k,s;

// Netwon/Euler π 

// 2016-11-13 EWS

t:=0;

FOR k FROM 0 TO m DO

t:= t + exact( 10^n/16^k * ( 4/(8*k + 1) - 2/(8*k + 4) - 1/(8*k + 5) 

- 1/(8*k + 6) ) );

END;

t:= IP(t);

// string, left, right 

// must be in lower case

s:= string(t);

s:= left(s,1) + "." + right(s,n);

return s;

END;

#end

bbppi(25,200) returns "3.1415925635897932846"

Source:  Wikipeida.  "Approximations of π"  https://en.m.wikipedia.org/wiki/Approximations_of_π. Retrieved November 12, 2016

This blog is property of Edward Shore, 2016.