# hp41programs

Curvature &Torsion

# Curvature & Torsion for the HP-41

Overview

1°)  Derivatives

a)  First, Second & Third Derivatives
b)  Functions of N variables

2°)  Gaussian & Mean Curvatures of a Surface

a)  Surface
b)  Hypersurface

3°)  Curvature & Torsion of a Curve

a)  3-Dimensional Problem - Program#1
b)  3-Dimensional Problem - Program#2
c)  3-Dimensional Problem - Program#3
d)  3-Dimensional Problem - Program#4
e)  N-Dimensional Problem

-In paragraph 1 are listed several programs to calculate various derivatives
-Some of them are already listed elsewhere

Latest Update:   §3-c)&d)

1°)  Derivatives

a) First, Second & Third Derivatives

-This routine uses formulae that give exact results for polynomials of degree < 11
-The formulae for  f '(x) & f ''(x)  are given in "Numerical Differentiation for the HP-41" paragraph 1-b)
-The third derivative is calculated by

f '''(x) ~ [ -e f(x-5h) - d f(x-4h) - c f(x-3h) - b f(x-2h) - a f(x-h) + a f(x+h) + b f(x+2h) + c f(x+3h) + d f(x+4h) + e f(x+5h) ] / h3

with   a = -1669/720 ,  b = 8738/5040 ,  c = -4869/10080 ,  d = 2522/30240 ,  e = -205/30240

Data Registers:           •  R00 = f name                              ( Register R00 is to be initialized before executing "dF" )

R01 = x              R03 = f '(x)              R05 = f '''(x)        R07-R08:  temp
R02 = h              R04 = f ''(x)             R06 = f(x)
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

 01 LBL "dF"  02 STO 01  03 X<>Y  04 STO 02  05 STO 08  06 X<>Y  07 XEQ IND 00  08 STO 06   09 6  10 ST* 08  11 XEQ 01  12 8  13 *  14 STO 04  15 RCL 07  16 STO 03  17 ST+ 03  18 205  19 *  20 CHS  21 STO 05  22 XEQ 01 23 125  24 *  25 ST- 04  26 25  27 RCL 07   28 *  29 ST- 03  30 2522  31 LASTX   32 *  33 ST+ 05  34 XEQ 01          35  E3  36 *  37 ST+ 04  38 150  39 RCL 07   40 *  41 ST+ 03  42 14607  43 LASTX  44 * 45 ST- 05  46 XEQ 01  47 6 E3  48 *  49 ST- 04  50 600  51 RCL 07   52 *  53 ST- 03  54 52428  55 LASTX          56 *  57 ST+ 05  58 XEQ 01   59 42 E3  60 *  61 ST+ 04  62 2100  63 RCL 07   64 *  65 ST+ 03  66 70098 67 LASTX  68 *  69 ST- 05  70 RCL 02   71 ST/ 04  72 2520  73 *  74 ST/ 03  75 10  76 *  77 ST/ 04  78 1.2  79 *  80 RCL 02          81 X^2  82 *  83 ST/ 05  84 GTO 02  85 LBL 01  86 RCL 02  87 ST- 08  88 RCL 01 89 RCL 08  90 +  91 XEQ IND 00  92 STO 07   93 RCL 01  94 RCL 08  95 -  96 XEQ IND 00  97 RCL 07  98 X<>Y  99 ST- 07 100 + 101 RCL 06 102 ST+ X 103 - 104 RTN 105 LBL 02 106 RCL 05 107 RCL 04 108 RCL 03 109 END

( 190 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Z / f '''(x) Y h f ''(x) X x f '(x)

Example:      f(x) = exp(x) + ln(x)     x = 2

 01  LBL "T"  02  E^X  03  LASTX  04  LN  05  +  06  RTN

"T"  ASTO 10

-With  h = 0.1

0.1   ENTER^
2     XEQ "dF"  >>>>   f '(2) = 7.889056111                                            ---Execution time = 14s---
RDN   f "(2) = 7.139056635
RDN  f '''(2) = 7.639050926

h = 0.1 is "often" ( but not always ) a good choice.

b) Functions of N Variables ( N < 10 )

-Formulas of order 10 are again used to estimate the partial derivatives of order 1 & 2 of a function f of n variables

Data Registers:   •  R00 = Function name                               ( Registers R00 thru Rnn are to be initialized before executing "FD" or "SD"  )

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h                                            R11 & R13 to R18 :  temp
R12 =  f / xi  or  2f / xixj
Flags: F10
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 01 LBL "FD"  02 LBL 01  03 STO 13  04 CLX  05 STO 14  06 RCL IND 13  07 STO 11  08 X<> Z  09 STO 10  10 XEQ 03  11 14  12 *  13 STO 12  14 XEQ 03  15 4  16 *  17 ST- 12  18 XEQ 03  19 ST+ 12  20 6  21 ST* 12  22 XEQ 03 23 ST- 12  24 25  25 ST* 12  26 XEQ 03  27 ST+ X  28 ST+ 12  29 2520  30 RCL 10  31 *  32 ST/ 12  33 RCL 11  34 STO IND 13  35 RCL 12  36 RTN  37 GTO 01  38 LBL 03  39 RCL 10  40 ST+ 14  41 RCL 11  42 RCL 14  43 -  44 STO IND 13 45 XEQ IND 00  46 STO 15  47 RCL 11  48 RCL 14  49 +  50 STO IND 13  51 XEQ IND 00  52 RCL 15  53 -  54 RTN  55 LBL "SD"  56 LBL 02  57 CF 10  58 X=Y?  59 SF 10  60 STO 16  61 RDN  62 STO 17  63 CLX  64 STO 15  65 RCL IND T  66 STO 13 67 RCL IND 17  68 STO 14  69 R^  70 STO 10  71 XEQ IND 00  72 STO 11  73 XEQ 03  74 42  75 *  76 STO 12  77 XEQ 03  78 6  79 *  80 ST- 12  81 XEQ 03  82 ST+ 12  83 XEQ 03  84 8  85 /  86 ST- 12  87  E3  88 ST* 12 89 XEQ 03  90 8  91 *  92 ST+ 12  93 RCL 13  94 STO IND 16  95 25200  96 FC? 10  97 ST+ X  98 FC?C 10  99 CHS 100 RCL 10 101 X^2 102 * 103 ST/ 12 104 RCL 12 105 RTN 106 GTO 02 107 LBL 03 108 RCL 10 109 ST+ 15 110 RCL 15 111 RCL 13 112 + 113 STO IND 16 114 XEQ IND 00 115 STO 18 116 RCL 13 117 RCL 15 118 - 119 STO IND 16 120 XEQ IND 00 121 RCL 11 122 ST+ X 123 - 124 ST+ 18 125 RCL 18 126 FS? 10 127 RTN 128 RCL 14 129 RCL 15 130 - 131 STO IND 17 132 XEQ IND 00 133 ST- 18 134 RCL 13 135 STO IND 16 136 XEQ IND 00 137 ST+ 18 138 RCL 14 139 RCL 15 140 + 141 STO IND 17 142 XEQ IND 00 143 ST+ 18 144 RCL 15 145 ST+ IND 16 146 XEQ IND 00 147 ST- 18 148 RCL 14 149 STO IND 17 150 RCL 18 151 RTN 152 END

( 254 bytes / SIZE 019 )

 STACK FD INPUTS OUTPUTS Y h / X i f 'xi = ¶f / ¶xi
 STACK SD  INPUTS OUTPUTS Z h / Y j / X i f "xixj = ¶2f / ¶xi¶xj

Example1:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

LBL "T"        E^X               +                    "T"    ASTO 00
RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
X^2              X^2               *
CHS            RCL 03          RTN

-With  h = 0.1

0.1   ENTER^
1    XEQ "FD"   >>>>   f 'x1 = f / x1 = -0.509989198   ( The last decimal should be a 5 )

0.1   ENTER^
2    XEQ "FD"   >>>>   f 'x2 = f / x2 = 0.367879440     ( The last decimal should be a 1 )

0.1   ENTER^
3    XEQ "FD"   >>>>   f 'x3 = f / x3 = 0.183939720     ( The last decimal should be a 1 )

Example2:       f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

LBL "T"        E^X               +                    "T"    ASTO 00
RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
X^2              X^2               *
CHS            RCL 03          RTN

-With  h = 0.1

0.1  ENTER^
1    ENTER^
XEQ "SD"  >>>>   f "xx = 2f / x2 = 0.509989206   ( exact derivative = 0.509989195 )

0.1   ENTER^
1    ENTER^
2    XEQ "SD"  >>>>   f "xy = 2f / xy = -0.735758889  ( the last decimal should be a 2 )

c) Hessian Matrix + Gradient  ( n < 10 )

"HESSGR" calculates and sores the coefficients of the Hessian matrix H and the coefficients of the gradient G into registers R25  R26 ..... in the following way:

W =   H  Gt      where  Gt  = transposed G              So, we get a square matrix W of order (n+1)
G  0

with        H = [ hi,j ]   with  hi,j = 2f / xixj
G = [ gk ]    with  gk = f / xi

-This form is useful to calculate the Gaussian curvature of a hypersurface.

Data Registers:   •  R00 = Function name                                      ( Registers R00 thru Rnn are to be initialized before executing "HESSGR"  )

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h       R11 to R24 :  temp    R25 ...... = W

Flag:  F10
Subroutines:   "FD"  "SD"  and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 01 LBL "HESSGR"  02 STO 19  03 STO 20   04 1  05 +  06 STO 11  07 X^2  08 23  09 +  10 STO 21  11 STO 22  12 X<>Y  13 STO 10  14 CLX  15 SIGN 16 +  17 RCL 11   18  E-5  19 *  20 +  21 STO 23  22 CLX  23 STO IND 23      24 LBL 01  25 RCL 10   26 RCL 20  27 XEQ "FD"  28 STO IND 21  29 DSE 23  30 STO IND 23 31 DSE 21  32 DSE 20  33 GTO 01  34 DSE 21  35 RCL 21   36 RCL 23   37 FRC  38 +  39 LBL 02              40 STO 23   41 STO 24  42 INT  43 STO 21  44 RCL 19  45 STO 20 46 LBL 03  47 RCL 10  48 RCL 19  49 RCL 20   50 XEQ "SD"  51 STO IND 21   52 STO IND 23      53 DSE 21   54 DSE 23  55 DSE 20  56 GTO 03  57 RCL 24  58 DSE X  59 1  60 - 61 DSE 19  62 GTO 02   63 RCL 22             64 1  65 +  66  E3  67 /  68 25  69 +  70 RCL 24  71 FRC  72 +  73 END

( 130 bytes / SIZE 25+(n+1)^2 )

 STACK INPUTS OUTPUTS Y h / X n < 10 25.eeemm

Where  m = n+1

Example:    f(x,y,z) = x4 y3 z2 - 1    x = y = z = 1

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 03  05  *  06  X^2  07  RCL 02  08  ST* Y  09  X^2  10  *  11  1  12  -  13  RTN

"T"  ASTO 00

1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

0.1   ENTER^
3     XEQ "HESSGR"  >>>>    25.04003                                    ---Execution time = 124s---

-And we find the matrix W in registers  R25 to R40

12   12   8   4
12    6    6   3
8     6    2   2
4     3    2   0

12  12  8
-The Hessian matrix is   12   6   6  and the gradient is  [ 4  3  2 ]
8    6   2

2°)  Gaussian & Mean Curvatures of a Surface

a) Surface

-"KGM" calculates the Gaussian curvature KG and the mean curvature Km of a surface defined by a function  z = z(x,y)

-The Gaussian Curvature KG is be computed by:

KG = [ ( 2z / x2 ) ( 2z / y2 ) - ( 2z / xy )2 ] / [ 1 + ( z / x )2 + ( z / y )2 ]2

-The mean curvature Km is given by:

Km = (1/2) [ t ( 1+p2 ) + r ( 1+q2 ) - 2 p q s ] / ( 1+p2+q2 )3/2

where  p = z / x , q = z / y,  r = 2z / x2,  s = 2z / xy,  t = 2z / y2

Data Registers:           •  R00 = f name                              ( Register R00 is to be initialized before executing "KGM" )

R01 to R15:  temp
Flags: /
Subroutines:   "dF2"  ( cf "Differential Geometry for the HP41" )
and a program that takes x in X-register, y in Y-register and returns z(x,y) in X-register

 01 LBL "KGM"  02 XEQ "dF2"  03 X^2  04 X<>Y   05 X^2  06 +  07 1 08 +  09 STO 10        10 X^2  11 X<> Z   12 *  13 RCL 09  14 X^2 15 -  16 X<>Y   17 /  18 RCL 03        19 X^2  20 RCL 07  21 ST* Y 22 +  23 RCL 06   24 X^2  25 RCL 04         26 ST* Y   27 +  28 + 29 RCL 03  30 RCL 06  31 *  32 RCL 09        33 *  34 ST+ X   35 - 36 RCL 10         37 ST/ Y  38 SQRT  39 ST+ X   40 /  41 X<>Y   42 END

( 60 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS Z h / Y y Km X x KG

Example:   A surface is defined by  z(x,y) = Ln ( 1 + x2 y3 )    Calculate the Gaussian & mean curvatures at the point  (x,y) = (1,1)

 01  LBL "T"  02  X^2  03  X<>Y  04  ST* Y  05  X^2  06  *  07  LN1+X  08  RTN

"T"  ASTO 00   and if you choose h = 0.1

0.1  ENTER^
1    ENTER^    XEQ "KGM"  >>>>   KG = -0.124572624                             ---Execution time = 58s---
X<>Y  Km = -0.171206978

b) Hypersurface ( n < 10 )

-Here, we assume that the hypersurface is defined implicitly by   f( x1 , ........... , xn ) = 0

Formulae:

KG = - ( det W ) / ( - | Grad f | )n+1      where  W is defined in paragraph 1°) c)

Km = [ Grad f x H x TGrad f  -  | Grad f |2 Trace (H) ] / (n-1) / | Grad f |3     where  H = Hessian matrix

Data Registers:   •  R00 = Function name                                      ( Registers R00 thru Rnn are to be initialized before executing "KGMN"  )

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h       R11 to R24 :  temp    R25 ...... = W

Flag:  F10
Subroutines:   "HESSGR"  "DET"  ( the latest version listed in "Determinants for the HP41" )
and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 01 LBL "KGMN"  02 XEQ "HESSGR"  03 STO 11   04 FRC  05 ENTER  06 ISG X  07 INT  08 STO 12   09 STO 13  10 CLX  11  E3  12 *  13 INT  14 DSE X  15 STO 15  16 0  17 DSE 12  18 LBL 01  19 RCL IND Y  20 X^2  21 + 22 DSE Y  23 DSE 12  24 GTO 01  25 STO 16  26 RCL 15   27 RCL 15   28 RCL 13  29 STO 12  30 STO 14  31 -  32 CLA  33 DSE 12  34 DSE 14  35 LBL 02  36 RCL IND Y        37 RCL IND Y  38 *  39 ST+ M  40 RDN  41 DSE Y  42 DSE X 43 DSE 14  44 GTO 02  45 0  46 X<> M  47 RCL IND 15       48 *  49 ST+ N  50 CLX  51 RCL 13   52 DSE X  53 STO 14   54 ST+ Z  55 SIGN  56 -  57 DSE 15  58 DSE 12  59 GTO 02  60 RCL N  61 STO 12  62 RCL 11  63  E-5 64 +  65 0  66 LBL 03  67 RCL IND Y        68 +  69 ISG Y  70 GTO 03  71 RCL 16   72 *  73 RCL 12   74 X<>Y  75 -  76 RCL 16  77 ST/ Y  78 SQRT  79 RCL 13  80 2  81 -  82 *  83 /  84 STO 12 85 RCL 00  86 STO 14  87 RCL 11  88 XEQ "DET"        89 RCL 00  90 RCL 14   91 STO 00  92 CLX  93 RCL 16   94 SQRT  95 CHS  96 RCL 13  97 Y^X  98 /  99 CHS 100 STO 13 101 RCL 12 102 X<>Y 103 CLA 104 END

( 160 bytes / SIZE 25+(n+1)^2 )

 STACK INPUTS OUTPUTS Z / det W Y h Km X n < 10 KG

Example:   The surface ( n = 3 )   x4 y3 z2 - 1 = 0   with  x = y = z = 1

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 03  05  *  06  X^2  07  RCL 02  08  ST* Y  09  X^2  10  *  11  1  12  -  13  RTN

"T"  ASTO 00

1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

0.1   ENTER^
3     XEQ "KGMN"  >>>>      KG  =  0.256837099                             ---Execution time = 151s---
RDN      Km  =  0.518666290
X<>Y  det W = - 216

Warning:

-This program does not check that  f(x1, ..... ,xn) = 0
-The matrix W is not saved.

3°)  Curvature & Torsion of a Curve

a) 3-Dimensional Problem - Program#1

-We assume that the curve is parameterized by 3 functions  x(t) , y(t) , z(t)

-The curvature is calculated by   khi = | r' x r'' | / | r' |3                                      where  the vector  r = [ x(t) , y(t) , z(t) ]
and the torsion by   tau = ( r' x r'' ) .r''' / | r' x r'' |2                        x = cross-product  and  . = dot-product

Data Registers:              R00 & R05 to R15:  temp
R01 = t  R02 = h   R03 = khi   R04 = tau
Flags: /
Subroutines:    "dF"  ( cf paragraph 1 )
and  3 functions named  "X"  "Y"  "Z"   that compute x(t) , y(t) , z(t)  assuming t is in X-register upon entry

 01 LBL "KHT"   02 "Z"  03 ASTO 00   04 XEQ "dF"  05 STO 12  06 RDN  07 STO 13   08 X<>Y  09 STO 14  10 "Y"  11 ASTO 00  12 RCL 02   13 RCL 01  14 XEQ "dF"  15 STO 09  16 RDN  17 STO 10 18 X<>Y  19 STO 11  20 "X"  21 ASTO 00   22 RCL 02  23 RCL 01   24 XEQ "dF"   25 RCL 12  26 RCL 09  27 R-P  28 RCL 03       29 R-P  30 STO 00  31 RCL 13  32 RCL 09  33 *  34 RCL 12 35 RCL 10  36 *  37 -  38 X^2  39 STO 06   40 LASTX  41 RCL 05         42 STO 08   43 *  44 RCL 12   45 RCL 04  46 STO 07  47 *  48 RCL 13  49 RCL 03  50 *  51 - 52 ENTER  53 X^2  54 ST+ 06  55 CLX  56 RCL 11   57 *  58 +  59 RCL 10   60 RCL 03         61 *  62 RCL 09  63 RCL 07  64 *  65 -  66 ENTER  67 X^2  68 ST+ 06 69 CLX  70 RCL 14  71 *  72 +  73 RCL 06   74 ST/ Y  75 SQRT  76 RCL 00          77 3  78 Y^X  79 /  80 X<>Y  81 STO 04  82 X<>Y  83 STO 03  84 END

( 110 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS Y h tau(t) X t khi(t)

Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1.3

-Load the following routines in main memory

 01  LBL "X"  02  E^X  03  LASTX  04  COS  05  *  06  RTN  07  LBL "Y"  08  E^X  09  LASTX  10  SIN  11  *  12  RTN  13  LBL "Z"  14  LN  15  RTN

XEQ "RAD"   and if you choose  h = 0.1

0.1   ENTER^
1.3   XEQ "KHT"  >>>>  khi (1.3)  = 0.194807823                       ---Execution time = 46s---
X<>Y  tau (1.3)  = 0.123665257

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529

b) 3-Dimensional Problem - Program#2

-Instead of 3 subroutines, the curve may be computed by a single routine.
-The main program becomes longer but faster:

Data Registers:           •  R00 = Function name           R05 to R18:  temp           ( Register R00 is to be initialized before executing "KHT"  )

R01 = t    R03 = khi
R02 = h   R04 = tau
Flags: /
Subroutine:   A program that takes t in X-register and returns x(t)  y(t)  z(t)  in registers  X  Y  Z  respectively

 01 LBL "KHT"      02 STO 01  03 X<>Y  04 STO 02  05 6  06 *  07 STO 18  08 XEQ 01   09 RCL 12  10 STO 03  11 ST+ 03  12 205  13 CHS  14 *  15 STO 09  16 RCL 13  17 STO 04  18 ST+ 04  19 LASTX  20 *  21 STO 10  22 RCL 14  23 STO 05  24 ST+ 05  25 LASTX  26 *  27 STO 11  28 RCL 15  29 8  30 ST* 16  31 ST* 17  32 *  33 STO 06  34 RCL 16  35 STO 07  36 RCL 17  37 STO 08  38 XEQ 01 39 RCL 12  40 25  41 *  42 ST- 03  43 RCL 13  44 LASTX  45 *  46 ST- 04  47 RCL 14           48 LASTX  49 *  50 ST- 05  51 RCL 15  52 125  53 ST* 16  54 ST* 17  55 *  56 ST- 06  57 RCL 16  58 ST- 07  59 RCL 17  60 ST- 08  61 RCL 12  62 2522  63 ST* 13  64 ST* 14  65 *  66 ST+ 09  67 RCL 13  68 ST+ 10  69 RCL 14  70 ST+ 11  71 XEQ 01  72 RCL 12  73 150  74 *  75 ST+ 03  76 RCL 13 77 LASTX  78 *  79 ST+ 04  80 RCL 14  81 LASTX  82 *  83 ST+ 05  84 RCL 15   85  E3  86 ST* 16  87 ST* 17  88 *  89 ST+ 06  90 RCL 16           91 ST+ 07  92 RCL 17  93 ST+ 08  94 RCL 12  95 14607  96 ST* 13  97 ST* 14  98 *  99 ST- 09 100 RCL 13 101 ST- 10 102 RCL 14 103 ST- 11 104 XEQ 01 105 RCL 12 106 600 107 * 108 ST- 03 109 RCL 13 110 LASTX 111 * 112 ST- 04 113 RCL 14 114 LASTX 115 * 116 ST- 05 117 RCL 15 118 6 E3 119 ST* 16 120 ST* 17 121 * 122 ST- 06 123 RCL 16  124 ST- 07 125 RCL 17  126 ST- 08 127 RCL 12 128 52428 129 ST* 13 130 ST* 14 131 * 132 ST+ 09 133 RCL 13         134 ST+ 10 135 RCL 14 136 ST+ 11 137 XEQ 01 138 RCL 12 139 2100 140 * 141 ST+ 03 142 RCL 13 143 LASTX 144 * 145 ST+ 04 146 RCL 14 147 LASTX 148 * 149 ST+ 05 150 RCL 15 151 42 E3 152 ST* 16 153 ST* 17 154 * 155 ST+ 06 156 RCL 16 157 ST+ 07 158 RCL 17 159 ST+ 08 160 RCL 12 161 70098 162 ST* 13 163 ST* 14 164 * 165 ST- 09 166 RCL 13 167 ST- 10 168 RCL 14 169 ST- 11 170 RCL 01 171 XEQ IND 00  172 73766 173 ST* T 174 ST* Z 175 * 176 ST- 06 177 RDN 178 ST- 07 179 X<>Y 180 ST- 08 181 GTO 00 182 LBL 01 183 RCL 02 184 ST- 18 185 RCL 01 186 RCL 18 187 + 188 XEQ IND 00 189 STO 12 190 STO 15 191 RDN 192 STO 13 193 STO 16 194 X<>Y 195 STO 14 196 STO 17 197 RCL 01 198 RCL 18 199 - 200 XEQ IND 00  201 ST- 12 202 ST+ 15 203 RDN 204 ST- 13 205 ST+ 16 206 X<>Y 207 ST- 14 208 ST+ 17 209 RTN 210 LBL 00 211 RCL 04 212 RCL 08 213 * 214 RCL 05 215 RCL 07 216 * 217 - 218 ST* 09 219 X^2 220 RCL 05 221 RCL 06 222 * 223 RCL 03 224 RCL 08 225 * 226 - 227 ST* 10 228 X^2 229 + 230 RCL 03 231 RCL 07 232 * 233 RCL 04 234 RCL 06 235 * 236 - 237 ST* 11 238 X^2 239 + 240 RCL 09          241 RCL 10 242 + 243 RCL 11  244 + 245 2100 246 * 247 X<>Y 248 STO 12 249 / 250 X<> 04 251 X^2 252 RCL 03 253 X^2 254 + 255 RCL 05 256 X^2 257 + 258 ENTER 259 SQRT 260 * 261 RCL 12 262 SQRT 263 252 264 * 265 X<>Y 266 / 267 STO 03 268 RCL 04 269 X<>Y 270 END

( 422 bytes / SIZE 019 )

 STACK INPUTS OUTPUTS Y h tau(t) X t khi(t)

Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1.3

-Load the following routines in main memory

 01  LBL "T"  02  LN  03  LASTX  04  ENTER  05  E^X  06  P-R  07  RTN

"T"  ASTO 00     XEQ "RAD"   and if you choose  h = 0.1

0.1   ENTER^
1.3   XEQ "KHT"  >>>>  khi (1.3)  = 0.194807828                       ---Execution time = 25s---
X<>Y  tau (1.3)  = 0.123665240

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529

Note:

-About 21 seconds have been saved

c) 3-Dimensional Problem - Program#3

-This variant employs other formulae to compute the derivatives:

y'  ~ (1/(60.h)) ( -y-3 + 9.y-2 - 45.y-1 + 45.y1 - 9.y2 + y3 )
y''  ~ (1/(90.h2)) [ y-3 - (27/2).y-2 + 135.y-1 - 245.y0 + 135.y1 - (27/2).y2 + y3 ]
y'''  ~ (1/(240.h3)) [ -7.y-4 + 72.y-3 - 338.y-2 + 488.y-1 - 488.y1 + 338.y2 - 72.y3 + 7.y4 ]

-These formulas are of order 6 ( except y''' of order 8 )

Data Registers:           •  R00 = Function name           R03 to R15:  temp           ( Register R00 is to be initialized before executing "KHT"  )

R01 = t    R02 = h

Flags: /
Subroutine:   A program that takes t in X-register and returns x(t)  y(t)  z(t)  in registers  X  Y  Z  respectively

 01 LBL "KHT"  02 STO 01  03 X<>Y  04 STO 02  05 4  06 *  07 STO 15  08 XEQ 02  09 RCL 12  10 7  11 ST* 13  12 ST* 14  13 *  14 STO 09  15 RCL 13  16 STO 10  17 RCL 14  18 STO 11  19 XEQ 01  20 STO 08  21 RDN  22 STO 07  23 X<>Y  24 STO 06  25 RCL 12  26 STO 03            27 72  28 *  29 ST- 09 30 RCL 13  31 STO 04  32 LASTX  33 *  34 ST- 10  35 RCL 14  36 STO 05  37 LASTX  38 *  39 ST- 11  40 XEQ 01  41 13.5  42 ST* T  43 ST* Z  44 *  45 ST- 08  46 RDN  47 ST- 07  48 X<>Y  49 ST- 06  50 RCL 12  51 9  52 *  53 ST- 03  54 RCL 13  55 LASTX  56 *  57 ST- 04  58 RCL 14 59 LASTX  60 *  61 ST- 05  62 RCL 12  63 338  64 ST* 13  65 ST* 14  66 *  67 ST+ 09  68 RCL 13  69 ST+ 10  70 RCL 14  71 ST+ 11  72 XEQ 01  73 135  74 ST* T  75 ST* Z  76 *  77 ST+ 08  78 RDN  79 ST+ 07  80 X<>Y  81 ST+ 06  82 RCL 12  83 45  84 *  85 ST+ 03  86 RCL 13            87 LASTX 88 *  89 ST+ 04  90 RCL 14  91 LASTX  92 *  93 ST+ 05  94 RCL 12  95 488  96 ST* 13  97 ST* 14  98 *  99 ST- 09 100 RCL 13 101 ST- 10 102 RCL 14 103 ST- 11 104 GTO 03 105 LBL 01 106 RCL 02 107 ST- 15 108 LBL 02 109 RCL 01 110 RCL 15 111 + 112 XEQ IND 00 113 STO 12 114 RDN 115 STO 13           116 X<>Y 117 STO 14 118 RCL 01 119 RCL 15 120 - 121 XEQ IND 00 122 RCL 12 123 X<>Y 124 ST- 12 125 + 126 X<>Y 127 RCL 13 128 X<>Y 129 ST- 13 130 + 131 RCL 14 132 R^ 133 ST- 14 134 + 135 RTN 136 LBL 03 137 RCL 01           138 XEQ IND 00 139 245 140 ST* T 141 ST* Z 142 * 143 ST- 06 144 RDN 145 ST- 07 146 X<>Y 147 ST- 08 148 RCL 04 149 RCL 08 150 * 151 RCL 05 152 RCL 07 153 * 154 - 155 ST* 09 156 X^2 157 RCL 05 158 RCL 06 159 * 160 RCL 03 161 RCL 08 162 * 163 - 164 ST* 10 165 X^2 166 + 167 RCL 03           168 RCL 07 169 * 170 RCL 04 171 RCL 06 172 * 173 - 174 ST* 11 175 X^2 176 + 177 RCL 09 178 RCL 10 179 RCL 11 180 + 181 + 182 22.5 183 * 184 X<>Y 185 / 186 RCL 03 187 X^2 188 RCL 04 189 X^2 190 + 191 RCL 05           192 X^2 193 + 194 ENTER 195 SQRT 196 * 197 R^ 198 SQRT 199 40 200 * 201 X<>Y 202 / 203 END

( 283 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS Y h tau(t) X t khi(t)

Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1.3

-Load the following routines in main memory

 01  LBL "T"  02  LN  03  LASTX  04  ENTER  05  E^X  06  P-R  07  RTN

"T"  ASTO 00     XEQ "RAD"   and if you choose  h = 0.05

0.05   ENTER^
1.3    XEQ "KHT"  >>>>  khi (1.3)  = 0.194807823                       ---Execution time = 19s---
X<>Y  tau (1.3)  = 0.123665666

-The exact values are   khi (1.3) = 0.194807823...  &  tau (1.3) = 0.123665529...

Note:

-About 139 bytes & 6 extra seconds have been saved and the results are even more accurate !

d) 3-Dimensional Problem - Program#4

-If we use more data registers, the number of bytes may be decreased:
-The same formulas of order 6 ( or 8 ) are employed.

Data Registers:           •  R00 = Function name           R03 to R18:  temp                                   ( Register R00 is to be initialized before executing "KHT"  )

R01 = t    R02 = h

Flags: /
Subroutine:   A program that takes t in X-register and returns x(t)  y(t)  z(t)  in registers  X  Y  Z  respectively

-If you don't have an HP41CX, add  STO 03  STO 04  STO 05  STO 06  STO 07  STO 08  STO 09  STO 10  STO 11  after line 10 and delete lines 08-09.

 01 LBL "KHT"  02 STO 01  03 X<>Y  04 STO 02  05 5  06 *  07 STO 18  08 3.011  09 CLRGX  10 CLX  11 STO 15  12 7  13 XEQ 01  14 1  15 STO 15  16 72  17 CHS  18 XEQ 01  19 9  20 CHS  21 STO 15            22 13.5  23 CHS 24 338  25 XEQ 01  26 45  27 STO 15  28 3  29 *  30 488  31 CHS  32 XEQ 01  33 GTO 02  34 LBL 01  35 STO 17  36 X<>Y  37 STO 16  38 RCL 02  39 ST- 18  40 RCL 01  41 RCL 18            42 +  43 XEQ IND 00  44 STO 12  45 RDN  46 STO 13 47 X<>Y  48 STO 14  49 RCL 01  50 RCL 18  51 -  52 XEQ IND 00  53 RCL 14  54 R^  55 ST- 14  56 +  57 RCL 13  58 R^  59 ST- 13  60 +  61 RCL 12  62 R^  63 ST- 12  64 +  65 RCL 16            66 ST* T  67 ST* Z  68 *  69 ST+ 06 70 RDN  71 ST+ 07  72 X<>Y  73 ST+ 08  74 RCL 12  75 RCL 15  76 *  77 ST+ 03  78 RCL 13  79 LASTX  80 *  81 ST+ 04  82 RCL 14  83 LASTX  84 *  85 ST+ 05  86 RCL 12            87 RCL 17  88 ST* 13  89 ST* 14  90 *  91 ST+ 09  92 RCL 13 93 ST+ 10  94 RCL 14  95 ST+ 11  96 RTN  97 LBL 02  98 RCL 01  99 XEQ IND 00 100 245 101 ST* T 102 ST* Z 103 * 104 ST- 06 105 RDN 106 ST- 07 107 X<>Y 108 ST- 08 109 RCL 04 110 RCL 08           111 * 112 RCL 05 113 RCL 07 114 * 115 - 116 ST* 09 117 X^2 118 RCL 05 119 RCL 06 120 * 121 RCL 03 122 RCL 08 123 * 124 - 125 ST* 10 126 X^2 127 + 128 RCL 03 129 RCL 07 130 * 131 RCL 04           132 RCL 06 133 * 134 - 135 ST* 11 136 X^2 137 + 138 RCL 09 139 RCL 10 140 RCL 11 141 + 142 + 143 22.5 144 * 145 X<>Y 146 / 147 RCL 03 148 X^2 149 RCL 04 150 X^2 151 + 152 RCL 05           153 X^2 154 + 155 ENTER 156 SQRT 157 * 158 R^ 159 SQRT 160 40 161 * 162 X<>Y 163 / 164 END

( 236 bytes / SIZE 019 )

 STACK INPUTS OUTPUTS Y h tau(t) X t khi(t)

Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1.3

-Load the following routines in main memory

 01  LBL "T"  02  LN  03  LASTX  04  ENTER  05  E^X  06  P-R  07  RTN

"T"  ASTO 00     XEQ "RAD"   and if you choose  h = 0.05

0.05   ENTER^
1.3    XEQ "KHT"  >>>>  khi (1.3)  = 0.194807823                       ---Execution time = 21s---
X<>Y  tau (1.3)  = 0.123665666

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529

Note:

47 bytes are saved but the program is slightly slower ( 2 extra seconds )

e) N-Dimensional Problem

-The method follows the Gram-Schmidt process ( cf references [1] & [2] ). After some calculations, we get:

-With    r = [ x1(t) , x2(t) , ........ , xn(t) ]

•   khi = | e'1 | / | r' |

where   e1 = r' / | r' |    whence   e'1 = r'' / | r' | - r' < r' , r'' > / | r' |3                         < , > = dot product.

•   tau = < e'2 , ê3 > / ( | r' | | ê3 | )

where   e'2 = e''1 / | e'1 | - e'1 < e'1 , e''1 > / | e'1 |3
e''1 = r''' / | r' | - 2 r'' < r' , r'' > / | r' |3 + 3 r' < r' , r'' >2 / | r' |5 - r' | r'' |2 / | r' |3 - r' < r' , r''' >2 / | r' |3
and     ê3  = r''' - r' < r' , r''' > / | r' |2 - e'1 < e'1 , r''' > / | e'1 |2

Data Registers:              R00 & R05 to R09 & R11 to R5n+12:  temp

R01 = t    R03 = khi
R02 = h   R04 = tau
Flags: /
Subroutines:    "dF"  ( cf paragraph 1 )
and  n functions named  "X1"  "X2" ........... "XN"   that compute x1(t) , x2(t) , ........ , xn(t)  assuming t is in X-register upon entry

 01 LBL "KHTN"  02 STO 01  03 RDN  04 STO 02   05 X<>Y  06 STO 09  07 STO 10  08 FIX 00  09 CF 29  10 LBL 01  11 "X"  12 ARCL 09  13 ASTO 00  14 RCL 02  15 RCL 01  16 XEQ "dF"  17 RCL 09  18 12  19 +  20 X<>Y  21 STO IND Y   22 CLX  23 RCL 10  24 +  25 RCL 04  26 STO IND Y    27 CLX  28 RCL 10  29 +  30 RCL 05  31 STO IND Y  32 DSE 09  33 GTO 01  34 FIX 09  35 SF 29  36 RCL 10  37 12.012  38 +  39 STO 05 40 INT  41 RCL 10  42 +  43 STO 06   44 LASTX  45 +  46 STO 07  47 LASTX  48 +  49 STO 08  50 LASTX  51 +  52 STO 09  53 RCL 05  54 INT  55 XEQ 11  56 STO 03  57 RCL 06  58 RCL 05  59 INT  60 XEQ 13  61 STO 04  62 RCL 06  63 RCL 05  64 LBL 02  65 RCL IND X   66 RCL 04  67 *  68 RCL 03  69 X^2  70 /  71 RCL IND Z    72 X<>Y  73 -  74 RCL 03  75 /  76 STO IND 08  77 RDN  78 DSE 08 79 DSE Y  80 DSE X  81 GTO 02  82 RCL 10  83 ST+ 08  84 RCL 08  85 XEQ 11  86 STO 11  87 RCL 06  88 XEQ 11  89 STO 00  90 RCL 07  91 RCL 05  92 INT  93 XEQ 13  94 STO 12  95 LBL 03  96 RCL 04  97 RCL 03  98 /  99 X^2 100 3 101 * 102 RCL 00 103 X^2 104 - 105 RCL 12 106 - 107 RCL IND 05  108 * 109 RCL 04 110 ST+ X 111 RCL IND 06 112 * 113 - 114 RCL 03 115 X^2 116 / 117 RCL IND 07 118 + 119 RCL 03  120 / 121 STO IND 09 122 DSE 09 123 DSE 07 124 DSE 06 125 DSE 05 126 GTO 03 127 RCL 10 128 ST+ 05 129 ST+ 06 130 ST+ 07 131 ST+ 09 132 RCL 08 133 RCL 09 134 XEQ 13 135 STO 00 136 RCL 10 137 STO M 138 LBL 04 139 RCL IND 09 140 RCL IND 08 141 RCL 00 142 * 143 RCL 11 144 X^2 145 / 146 - 147 RCL 11 148 / 149 STO IND 09 150 DSE 09 151 DSE 08 152 DSE M 153 GTO 04 154 RCL 10 155 ST+ 08 156 ST+ 09 157 RCL 07 158 RCL 08  159 XEQ 13 160 STO 00 161 LBL 05 162 RCL IND 07 163 RCL IND 05 164 RCL 12 165 * 166 RCL 03 167 X^2 168 / 169 - 170 RCL IND 08 171 RCL 00 172 * 173 RCL 11 174 X^2 175 / 176 - 177 STO IND 08 178 DSE 08 179 DSE 07 180 DSE 05 181 GTO 05 182 RCL 10 183 ST+ 05 184 ST+ 07 185 ST+ 08 186 RCL 08 187 XEQ 11 188 STO 00 189 RCL 08 190 RCL 09 191 XEQ 13 192 GTO 00 193 LBL 11 194 RCL 10 195 STO M 196 CLX 197 LBL 12 198 RCL IND Y 199 X^2 200 + 201 DSE Y 202 DSE M 203 GTO 12 204 SQRT 205 RTN 206 LBL 13 207 CLA 208 RCL 10  209 STO N 210 RDN 211 LBL 14 212 RCL IND Y 213 RCL IND Y  214 * 215 ST+ M 216 RDN 217 DSE Y 218 DSE X 219 DSE N 220 GTO 14 221 CLX 222 X<> M 223 RTN 224 LBL 00 225 RCL 03 226 RCL 00 227 * 228 / 229 STO 04 230 RCL 11 231 RCL 03 232 / 233 STO 03 234 RCL 10 235 3 236 X#Y? 237 GTO 00 238 RCL 14  239 RCL 18 240 * 241 RCL 15 242 RCL 17 243 * 244 - 245 RCL 19  246 * 247 RCL 13 248 RCL 18 249 * 250 RCL 15 251 RCL 16  252 * 253 - 254 RCL 20 255 * 256 - 257 RCL 13 258 RCL 17 259 * 260 RCL 14 261 RCL 16 262 * 263 - 264 RCL 21          265 * 266 + 267 SIGN 268 ST* 04 269 LBL 00 270 RCL 04 271 RCL 03  272 END

( 389 bytes / SIZE 5n+13 )

 STACK INPUTS OUTPUTS Z n < 62 / Y h tau(t) X t khi(t)

Example:   In a 4-dimensional space, the curve defined by

x1(t) = t4       x3(t) = t3      at the point  t = 1
x2(t) = t2       x4(t) = et

 01  LBL "X1"  02  X^2  03  X^2  04  RTN  05  LBL "X2"  06  X^2  07  RTN  08  LBL "X3"  09  ENTER^  10  X^2  11  *  12  RTN  13  LBL "X4"  14  E^X  15  RTN

-If you choose  h = 0.1

4   ENTER^
0.1  ENTER^
1   XEQ "KHTN"   >>>>   khi (1)  = 0.142277239                       ---Execution time = 71s---
X<>Y  tau (1)  = 0.121469324

Note:

-In the Gram-Schmidt process, the orthogonal basis does not always satisfy the right-hand rule.
-So, if n = 3 , the torsion given by the formulas above is multiplied by sign [ det ( r' , r'' , r''' ) ] ( lines 238 to 268 )
-Otherwise, the sign could be different from the result obtained by "KHT"

References:

[1]  https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_curves&oldid=740874421
[2]  https://en.wikipedia.org/w/index.php?title=Frenet%E2%80%93Serret_formulas&oldid=731811089
[3]  Ron Goldman - "Curvature formulas for implicit curves and surfaces"