hp41programs

Curvature &Torsion

Curvature & Torsion for the HP-41


Overview
 

 1°)  Derivatives

  a)  First, Second & Third Derivatives
  b)  Functions of N variables
  c)  Hessian Matrix+Gradient

 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 / xixj

 
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"