hp41programs

Derive2

Numerical Differentiation(2) for the HP-41


Overview
 

 1°)  Third Derivatives

  a)  MDV3 - 3 variables
  b)  MDV21 - 2 variables
  c)  Functions of N-variables

 2°)  Fourth Derivatives

  a)  MDV4 - 2 variables
  b)  Biharmonic Operator ( 2-Dim )  ( 3 programs )
  c)  Biharmonic Operator ( 3-Dim )  ( 2 programs )
  d)  Biharmonic Operator ( N-Dim )
  e)  Biharmonic Operator ( N-Dim ) - Radial Functions

 3°)  Triharmonic Operator ( 3-Dim ) ( 2 programs )

 
Latest Update:    §2°)b) & 2°)c)

-This page completes the programs listed in "Numerical Differentiation" & "Curvature and Torsion"
-Please take also a look at "Taylor Series for the HP-41"
-The programs hereunder deal with partial derivatives of order 3 & 4 and with partial derivatives of order 6 in paragraph 3.
 

1°)  Third Derivatives
 

     a)  MDV3 - 3 Variables
 

 "MDV3" employs a method of order 11 to estimate  f'''xyz = 3f / xyz

    h3 f'''xyz ~  SUMk=1,2,...,5  Ak [ fk00 - f-k00 + f0k0 - f0-k0 + f00k - f00-k + fkkk - f-k-k-k - ( fkk0 - f-k-k0 + fk0k - f-k0-k + f0kk - f0-k-k ) ]

 With    A1 = 5/6 ,  A2 = -5/84  ,  A3 = +5/756  ,  A4 = -5/8064  ,  A5 = +1/31500
 

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

                                         R01 = x     R04 = h
                                         R02 = y     R05 = f'''xyz = 3f / xyz
                                         R03 = z     R06-R07:  temp
Flags: /
Subroutine:  A program that takes x , y , z in registers X , Y , Z respectively and returns  f(x,y,z)  in X-register.
 
 

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

 
       ( 228 bytes / SIZE 008 )
 
 

      STACK        INPUTS      OUTPUTS
           T             h             /
           Z             z             /
           Y             y             /
           X             x          f'''xyz

 
Example:     f(x,y,z) = Ln ( 1 + x2 + 2.y + z3 )    x = y = z = 1
 

-Load this routine:
 
 

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

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

    0.1   ENTER^
     1     ENTER^  ENTER^   XEQ "MDV3"  >>>>   f'''xyz = 3f / xyz  =  0.191999728                ---Execution time = 62s---

-The exact value is  0.192
 

     b)  MDV21 - 2 Variables
 

 "MDV21" employs a similar method of order 11 to estimate  f'''xxy = 3f / x2y

    h3 f'''xxy ~  SUMk=1,2,...,5  Ak [ - 2 f0k + 2 f0-k  + ( fkk - f-k-k + f-kk - fk-k ) ]

 With    A1 = 5/6 ,  A2 = -5/84  ,  A3 = +5/756  ,  A4 = -5/8064  ,  A5 = +1/31500
 

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

                                         R01 = x     R03 = h
                                         R02 = y     R04 = f'''xxy = 3f / x2y  or  f'''xyy = 3f / y2x     R05-R06  temp

Flag:  F02            CF 02 ->  f'''xxy = 3f / x2y
                             SF 02  ->  f'''xyy = 3f / y2x

Subroutine:  A program that takes x , y  in registers X , Y respectively and returns  f(x,y)  in X-register.
 
 
 

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

 
    ( 162 bytes / SIZE 007 )
 
 

      STACK        INPUTS      OUTPUTS
           T             h             /
           Z             z             /
           Y             y             /
           X             x    f'''xxy or f'''xyy

     CF 02 -> f'''xxy
     SF 02 -> f'''xyy

Example:       f(x,y,z) = Ln ( 1 + x2 y )    x = 1  ,  y = 2
 

-Load this routine:
 
 

 01  LBL "T"
 02  X^2
 03  *
 04  LN1+X
 05  RTN

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

    CF 02

    0.1   ENTER^
     2     ENTER^
     1     XEQ "MDV21"  >>>>   f'''xxy = 3f / x2y  =  - 0.370369497                ---Execution time = 25s---

-The exact value is  -10/27 =  - 0.370370370.....

-With   SF 02  it yields   f'''xyy = 3f / y2x  =  - 0.148148208

-The exact result is  -4/27 =  - 0.148148148....
 

     c)  Functions of N-Variables
 

 "TD"  combines the 2 formulas above and the formula of order 10 given in "Curvature & Torsion" to evaluate
  the 3rd derivatives of a function of n variables ( n < 10 )
 

Data Registers:           •  R00 = f name                                 ( Registers R00 & R01 thru Rnn are to be initialized before executing "TD" )

                                      •  R01 = x1  ,  ..............  , •  Rnn = xn    ( n < 10 )

                                          R10 = h   R11 = D3f      R12 to R19:  temp
Flags: F09-F10
Subroutine:  A program that takes x1, .........,xn   in registers R01 .....  Rnn  respectively and returns  f(x1, .........,xn)  in X-register.
 

-Lines 48 and 75 are 3-byte GTOs
 
 

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

 
       ( 380 bytes / SIZE 020 )
 
 

      STACK        INPUTS      OUTPUTS
           T             h             /
           Z            k             /
           Y             j             /
           X             i       f(3)XiXjXk

 
Example:             f(x,y,z,t) = exp(-x2) Ln( y2 + z + t3 )     x = 1 , y = 2 , z = 3 , t = 1
 

-Load this short routine:
 
 

 01  LBL "T"
 02  RCL 01
 03  X^2
 04  CHS
 05  E^X
 06  RCL 02
 07  X^2
 08  RCL 03
 09  +
 10  RCL 04
 11  ENTER^
 12  X^2
 13  *
 14  +
 15  LN
 16  *
 17  RTN
 18  END

 
    "T"  ASTO 00

     1   STO 01
     2   STO 02
     3   STO 03
     1   STO 04

  •  f'''xyz = 3f / xyz  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      2    ENTER^
      3    XEQ "TD"   >>>>    f'''xyz = 3f / xyz  = 0.045984769                                          ---Execution time = 88s---

-Exact result = 0.045984930....

  •  f'''xtt = 3f / xt2  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      4    ENTER^    XEQ "TD"   >>>>    f'''xtt = 3f / xt2  = - 0.448353739                                    ---Execution time = 42s---

-Exact result = - 0.448353069....

  •  f'''yyy = 3f / y3  = ?

-With   h = 0.1

     0.1  ENTER^
      2    ENTER^   ENTER^    XEQ "TD"   >>>>    f'''yyy = 3f / y3  = - 0.045984325                               ---Execution time = 15s---

-Exact result = - 0.0459849301...

Note:

-Lines 02 to 09 sort  i j k  in increasing order.
 

2°)  Fourth Derivatives
 

     a)  MDV4 - 2 Variables
 

 "MDV4" evaluates   f'''xxyy = 4f / x2y2  with a formula of order 12

  h4 f'''xxyy ~  20881861 f00 / 3240000  + SUMk=1,....,5  Ak [ fkk + f-k-k + fk-k + f-kk - 2 ( fk0 + f-k0 + f0k + f0-k ) ]

  where   A1 = 5/3 ,  A2 = -5/84  ,  A3 = +5/1134  ,  A4 = -5/16128  ,  A5 = +1/78750
 

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

                                         R01 = x     R03 = h
                                         R02 = y     R04 = f""xxyy = 4f / x2y2     R05-R06-R07:  temp

Flags:  /
Subroutine:  A program that takes x , y  in registers X , Y respectively and returns  f(x,y)  in X-register.
 
 
 

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

 
    ( 164 bytes SIZE 008 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y             /
           X             x         f(4)xxyy

 
Example:    f(x,y) = Ln ( x2 + y3 )     x = 2 , y = 1
 

-Load for example this routine:
 
 

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

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

    0.1  ENTER^
     1    ENTER^
     2    XEQ "MDV4"   >>>>  f(4)xxyy =  4f / x2y2 =  - 0.038395989                                             ---Execution time = 33s---

-Exact result = - 0.0384
 

     b)  Biharmonic Operator - 2 Dim  ( 3 programs )

 
-This program employs the formula in reference1 to compute:    D2 f = 4f / x4 + 4f / y4 + 2 4f / x2y2

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

                                         R01 = x        R03 = h                  R05-R06-R07:  temp
                                         R02 = y        R04 = D2 f

Flags:  /
Subroutine:         A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.
 
 
 

 01 LBL "BHRM2"
 02 STO 01
 03 RDN
 04 STO 02
 05 X<>Y
 06 STO 03
 07 X<> T
 08 XEQ IND 00
 09 184
 10 *
 11 STO 04          
 12 RCL 03
 13 3
 14 *
 15 STO 05
 16 1
 17 CHS
 18 XEQ 01
 19 RCL 03		  20 ST+ X
 21 STO 05
 22 14
 23 XEQ 01
 24 RCL 03
 25 STO 05
 26 STO 06
 27 77
 28 CHS
 29 XEQ 01
 30 20
 31 STO 07          
 32 XEQ 02
 33 RCL 03
 34 ST+ X
 35 STO 05
 36 1
 37 CHS
 38 STO 07		  39 XEQ 02
 40 RCL 05
 41 X<> 06
 42 STO 05
 43 XEQ 02
 44 GTO 04
 45 LBL 01
 46 STO 07
 47 XEQ 01
 48 RCL 05
 49 CHS
 50 STO 05          
 51 LBL 01
 52 RCL 02
 53 RCL 01
 54 RCL 05
 55 +
 56 XEQ 03
 57 RCL 02		  58 RCL 05
 59 +
 60 RCL 01
 61 GTO 03
 62 LBL 02
 63 XEQ 02
 64 RCL 05
 65 CHS
 66 STO 05
 67 XEQ 02
 68 RCL 06          
 69 CHS
 70 STO 06
 71 XEQ 02
 72 RCL 05
 73 CHS
 74 STO 05
 75 LBL 02
 76 RCL 02		  77 RCL 05
 78 +
 79 RCL 01
 80 RCL 06
 81 +
 82 LBL 03
 83 XEQ IND 00
 84 RCL 07
 85 *
 86 ST+ 04
 87 RTN
 88 LBL 04
 89 RCL 03          
 90 X^2
 91 X^2
 92 6
 93 *
 94 ST/ 04
 95 RCL 04
 96 END

 
     ( 143 bytes / SIZE 008 )
 
 

           STACK           INPUTS         OUTPUTS
               Z                h                /
               Y                y                /
               X                x          D2 f(x,y)

 
Example:      f(x) =.Ln(2+x2+y)    x = 1 , y = 1
 
 

 01  LBL "T"
 02  X^2
 03  +
 04  2
 05  +
 06  LN
 07  RTN

 
      T  ASTO 00

-If you choose h = 0.1

      0.1  ENTER^
        1   ENTER^   XEQ "BHRM2"  >>>>   D2 f(1,1) = 0.289065                                      ---Execution time = 24s---

-Exact result = 37 / 128 = 0.2890625


-The 2nd program combines the method of order 12 given in paragraph 2°) a)
  and the method of order 10 in "Taylor Series for the HP-41"  to evaluate
 

            D2 f = 4f / x4 + 4f / y4 + 2 4f / x2y2
 

   ( like in "Differential Geometry for the HP-41" )
 

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

                                         R01 = x        R03 = h                  R05-R06-R07-R08:  temp
                                         R02 = y        R04 = D2 f

Flags:  /
Subroutine:         A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.
 
 
 

 01 LBL "BHRM2"
 02 STO 01
 03 RDN
 04 STO 02
 05 X<>Y
 06 STO 03
 07 STO 08
 08 X<>Y
 09 R^
 10 XEQ IND 00 
 11 4
 12 *
 13 STO 05
 14 6
 15 ST* 08
 16 XEQ 01
 17 24
 18 *
 19 RCL 06
 20 5173
 21 *
 22 -
 23 125
 24 /
 25 STO 04
 26 XEQ 01
 27 75		  28 *
 29 RCL 06
 30 10238
 31 *
 32 -
 33 16
 34 /
 35 ST- 04
 36 XEQ 01
 37 200
 38 *
 39 RCL 06        
 40 15007
 41 *
 42 -
 43 3
 44 /
 45 ST+ 04
 46 6
 47 ST/ 04
 48 XEQ 01
 49 150
 50 *
 51 RCL 06
 52 4669
 53 *
 54 -		  55 ST- 04
 56 7
 57 ST/ 04
 58 XEQ 01
 59 600
 60 *
 61 RCL 06
 62 2869
 63 *
 64 -
 65 RCL 04
 66 +
 67 180
 68 /
 69 GTO 02
 70 LBL 01
 71 RCL 03
 72 ST- 08
 73 RCL 02
 74 RCL 08
 75 +
 76 RCL 01
 77 XEQ IND 00 
 78 STO 06 
 79 RCL 02
 80 RCL 08
 81 -		  82 RCL 01
 83 XEQ IND 00
 84 ST+ 06
 85 RCL 02
 86 RCL 01
 87 RCL 08
 88 +
 89 XEQ IND 00
 90 ST+ 06
 91 RCL 02
 92 RCL 01
 93 RCL 08
 94 -
 95 XEQ IND 00
 96 ST+ 06
 97 RCL 02
 98 RCL 01
 99 RCL 08
100 ST+ Z
101 +
102 XEQ IND 00 
103 STO 07
104 RCL 02
105 RCL 01
106 RCL 08
107 ST+ Z
108 -		 109 XEQ IND 00
110 ST+ 07
111 RCL 02
112 RCL 01
113 RCL 08
114 ST- Z
115 +
116 XEQ IND 00 
117 ST+ 07
118 RCL 02
119 RCL 01
120 RCL 08
121 ST- Z
122 -
123 XEQ IND 00 
124 RCL 07
125 +
126 RCL 05
127 ST- 06
128 -
129 RTN
130 LBL 02
131 RCL 03
132 X^2
133 X^2
134 /
135 STO 04
136 END

 
     ( 213 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y             /
           X             x       D2 f(x,y)

 
Example:      f(x) =.Ln(2+x2+y)    x = 1 , y = 1
 
 

 01  LBL "T"
 02  X^2
 03  +
 04  2
 05  +
 06  LN
 07  RTN

 
      T  ASTO 00

-If you choose h = 0.1

      0.1  ENTER^
        1   ENTER^   XEQ "BHRM2"  >>>>   D2 f(1,1) = 0.288987222                                      ---Execution time = 35s---

-Exact result = 37 / 128 = 0.2890625

Note:

-In the following variant, all the partial derivatives are approximated by formulas of order 12
 
 

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

                                         R01 = x        R03 = h                  R05-R06-R07-R08:  temp
                                         R02 = y        R04 = D2 f

Flag:  F10
Subroutine:         A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.
 
 
 

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

 
     ( 233 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y             /
           X             x       D2 f(x,y)

 
Example:   the same one:   f(x) =.Ln(2+x2+y)    x = 1 , y = 1
 
 

 01  LBL "T"
 02  X^2
 03  +
 04  2
 05  +
 06  LN
 07  RTN

 
      T  ASTO 00

-If you choose h = 0.1

      0.1  ENTER^
        1   ENTER^   XEQ "BHRM2"  >>>>   D2 f(1,1) = 0.288981300                                      ---Execution time = 39s---

-Exact result = 37 / 128 = 0.2890625

Note:

-Difficult to know which method is better...
 

     c)  Biharmonic Operator - 3 Dim   ( 2 programs )

 
-The 1st program uses a similar formula employed in the 1st program in §2°)b) to evaluate:  

     D2 f = 4f / x4 + 4f / y4 + 4f / z4 + 2 ( 4f / x2y2 + 4f / x2z2 + 4f / y2z2 )
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D2 f
                                         R03 = z        R06-R07-R08:  temp
Flags: /
Subroutine:         A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.
 
   
 

 01 LBL "BHRM3"
 02 STO 01
 03 RDN
 04 STO 02
 05 RDN
 06 STO 03
 07 RDN
 08 STO 04
 09 RDN
 10 XEQ IND 00
 11 384
 12 *
 13 STO 05
 14 RCL 04
 15 3
 16 *
 17 STO 06
 18 1
 19 CHS
 20 XEQ 01
 21 RCL 04          
 22 ST+ X
 23 STO 06
 24 16
 25 XEQ 01		  26 RCL 04
 27 STO 06
 28 STO 07
 29 115
 30 CHS
 31 XEQ 01
 32 20
 33 STO 08
 34 XEQ 02
 35 RCL 04
 36 ST+ X
 37 STO 06
 38 1
 39 CHS
 40 STO 08
 41 XEQ 02
 42 RCL 06
 43 X<> 07
 44 STO 06          
 45 XEQ 02
 46 GTO 04
 47 LBL 01
 48 STO 08
 49 XEQ 01
 50 RCL 06		  51 CHS
 52 STO 06
 53 LBL 01
 54 RCL 03
 55 RCL 06
 56 +
 57 RCL 02
 58 RCL 01
 59 XEQ 03
 60 RCL 03
 61 RCL 02
 62 RCL 06
 63 +
 64 RCL 01
 65 XEQ 03
 66 RCL 03
 67 RCL 02
 68 RCL 01
 69 RCL 06
 70 +
 71 GTO 03
 72 LBL 02
 73 XEQ 02
 74 RCL 06          
 75 CHS		  76 STO 06
 77 XEQ 02
 78 RCL 07
 79 CHS
 80 STO 07
 81 XEQ 02
 82 RCL 06
 83 CHS
 84 STO 06
 85 LBL 02
 86 RCL 03
 87 RCL 07
 88 +
 89 RCL 02
 90 RCL 06
 91 +
 92 RCL 01
 93 XEQ 03
 94 RCL 03          
 95 RCL 02
 96 RCL 07
 97 +
 98 RCL 01
 99 RCL 06		 100 +
101 XEQ 03
102 RCL 03
103 RCL 07
104 +
105 RCL 02
106 RCL 01
107 RCL 06
108 +
109 LBL 03
110 XEQ IND 00
111 RCL 08
112 *
113 ST+ 05
114 RTN
115 LBL 04
116 RCL 04
117 X^2
118 X^2
119 6
120 *
121 ST/ 05
122 RCL 05          
123 END

 
       ( 176 bytes / SIZE 009 )
 
 

           STACK            INPUTS          OUTPUTS
               T                 h                 /
               Z                 z                /
               Y                 y                /
               X                 x         D2 f(x,y,z)

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

 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN

 
      T  ASTO 00

-If you choose h = 0.1

   0.1   ENTER^
    3     ENTER^
    2     ENTER^
    1     XEQ "BHRM3"  >>>> D2 f (1,2,3) = -14.33893267 = R05             ---Execution time = 67s---

-With  h = 0.075  it yields:  -14.34215664

-The exact result is  -14.34264116

-The 2nd program combines the methods of order 12 given in paragraph 2°) a)  and in "Taylor Series for the HP-41"  to evaluate
 

            D2 f = 4f / x4 + 4f / y4 + 4f / z4 + 2 ( 4f / x2y2 + 4f / x2z2 + 4f / y2z2 )
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D2 f
                                         R03 = z        R06-R07-R08-R09:  temp
Flag:  F10
Subroutine:         A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.
 
 
 
 

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

 
       ( 254 bytes / SIZE 010 )
 
 

      STACK        INPUTS      OUTPUTS
           T          h > 0             /
           Z             z             /
           Y             y             /
           X             x     D2 f(x,y,z)

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

 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN

 
      T  ASTO 00

-If you choose h = 0.1

   0.1   ENTER^
    3     ENTER^
    2     ENTER^
    1     XEQ "BHRM3"  >>>> D2 f (1,2,3) = -14.34283200 = R05             ---Execution time = 1m51s---

-The exact result is  -14.34264116
 

     d)  Biharmonic Operator - N Dim
 

-The same formulas of order 12 are now used to estimate the n-dimensional biharmonic operator
 

         D2 f = SUMk=1,....,n  4f / xk4 + 2  SUMj<k  4f / xj2yk2
 

Data Registers:           •  R00 = f name                         ( Registers R00 & R01 thru Rnn are to be initialized before executing "BHRMN" )

                                      •  R01 = x1  ,  ..............  , •  Rnn = xn    ( n < 10 )

                                          R10 = h   R11 = D2 f      R12 to R19:  temp
Flags: F09-F10
Subroutine:  A program that takes x1, .........,xn   in registers R01 .....  Rnn  respectively and returns  f(x1, .........,xn)  in X-register.
 

-Lines 71-108 are three-byte GTOs
 
 

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

 
     ( 316 bytes / SIZE 020 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             /
           X             n    D2 f(x1,...,xn)

  With  n < 10

Example1:    f(x,y,z,t) = exp(-x2) Ln( y2 + z + t3 )     x = 1 , y = 2 , z = 3 , t = 1
 

-With this short routine in main memory:
 
 

 01  LBL "T"
 02  RCL 01
 03  X^2
 04  CHS
 05  E^X
 06  RCL 02
 07  X^2
 08  RCL 03
 09  +
 10  RCL 04
 11  ENTER^
 12  X^2
 13  *
 14  +
 15  LN
 16  *
 17  RTN
 18  END

 
    "T"  ASTO 00

     1   STO 01  STO 04
     2   STO 02
     3   STO 03

    .1  ENTER^
    4   XEQ "BHRMN"  >>>>  D2 f (1,2,3,1) = -14.93977513 = R11                         ---Execution time = 3m53s---
 

Note:

-The function f is evaluated  10 n2 + 2 n + 1  times
 
 

      N   EVAL
      1     13
      2     45
      3     97
      4    169
      5    261
      6    373
      7    505
      8    657
      9    829

 
Example2:    9 variables,     f (x,y,z,t,u,v,w,r,s)  = [ exp(-x y z) + t u v w ] / Ln ( 1 + w r s )    x = y = ....... = s = 1
 
 
 

 01  LBL "T"
 02  RCL 01
 03  RCL 02
 04  RCL 03
 05  *
 06  *
 07  CHS
 08  E^X
 09  RCL 04
 10  RCL 05
 11  *
 12  RCL 06
 13  *
 14  RCL 07
 15  *
 16  +
 17  RCL 07
 18  RCL 08
 19  *
 20  RCL 09
 21  *
 22  LN1+X
 23  /
 24  RTN

 
    "T"  ASTO 00

     1   STO 01  STO 02  ..................  STO 09

    .1  ENTER^
     9  XEQ "BHRMN"  >>>>  D2 f(1,.....,1) = 103.2553329 = R11                                ---Execution time = 21m03s---
 

Notes:

-The exact value is  103.2389124...
-The result is returned in less than 3 seconds with V41.

-We can also combine a formula of order 10 for 4f / xk4 and a formula of order 12 for 4f / xj2yk2  ( like in "Differential Geometry for the HP-41" )

-Lines 60-97 are three-byte GTOs
 
 

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

 
          ( 296 bytes / SIZE 020 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             /
           X             n    D2 f(x1,...,xn)

     With  n < 10

Example1:    f(x,y,z,t) = exp(-x2) Ln( y2 + z + t3 )     x = 1 , y = 2 , z = 3 , t = 1

-Here, we get, with h = 0.1     D2 f (1,2,3,1) = -14.93977249 = R11                         ---Execution time = 3m43s---

-Since the exact result is  -14.9397000647...  the precision is slightly better than with the 1st version

Note:

-The function f is evaluated  10 n2  + 1  times
 
 

      N   EVAL
      1     11
      2     41
      3     91
      4    161
      5    251
      6    361
      7    491
      8    641
      9    811

 
Example2:        f (x,y,z,t,u,v,w,r,s)  = [ exp(-x y z) + t u v w ] / Ln ( 1 + w r s )    x = y = ....... = s = 1

-Here we get, again with h = 0.1,      D2 f(1,.....,1) = 103.2103466 = R11                                                     ---Execution time = 20m38s---

-Ulike the 1st example, the precision is smaller than the 1st version ( exact value =  103.2389124... )
-But it's difficult to decide which one is the best on average
 

      e)  Biharmonic Operator - N Dim - Radial Functions
 

-The formula becomes much simpler if the function f only depends on  r = ( x12 + ............ + xn2 )1/2

-The Laplacian is  D f = f " + 2 f ' ( n - 1 ) / r  ,  so the bilaplacian operator is:

      D2 f = f "" + 2 f ''' ( n - 1 ) / r  + f '' ( n - 1 ) ( n - 3 ) / r2 -  f ' ( n - 1 ) ( n - 3 ) / r3
 

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

                                         R01 = x        R04 = D2 f       R07 = f "(r)      R10-R11: temp
                                         R02 = h        R05 = 2 f(r)     R08 = f '''(r)
                                         R03 = n        R06 =  f '(r)     R09 = f ""(r)
Flags:  /
Subroutine:         A program which computes f(r) assuming r is in X-register.
 
 
 

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

 
     ( 264 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             n             /
           X             r         D2 f(r)

 
Example:    n = 7  ,   f(r) = Ln ( 1 + r4 )   ,   r = 2
 
 

 01  LBL "T"
 02  X^2
 03  X^2
 04  LN1+X
 05  RTN
 06  END

 
    "T"   ASTO 00    and with  h = 0.1

    0.1   ENTER^
     7     ENTER^
     2     XEQ "BHRMR"   >>>>   D2 f(2) = -7,695049202                                          ---Execution time = 14s---

Notes:

-Exact result = - 642696 / 83521 = - 7.695022808.....

-"BHRMR" employs formulas of order 10 to evaluate the derivatives.
-As usual, 4th derivatives are not easy to get with a great precision...
 

3°)  Triharmonic Operator ( 3-Dim )
 

-The triharmonic operator is the trilaplacian operator = the Laplacian of the Laplacian of the Laplacian of a function f

-"THRM" employs a method of order 10  to evaluate
 

   D3 f = 6f / x6 + 6f / y6 + 6f / z6 + 3 ( 6f / x4y2 + 6f / x4z2 + 6f / y4z2 + 6f / x2y4 + 6f / x2z4 + 6f / y2z4 ) + 6 6f / x2y2z2
 

   h6 D3 f  ~  SUMm=1....5  Am ( fm00 + f-m00 + f0m0 + f0-m0 + f00m + f00-m - 6 f000 )
                                     + Bm ( fmm0 + f-m-m0 + fm-m0 + f-mm0 + fm0m + f-m0-m + fm0-m + f-m0m + f0mm + f0-m-m + f0m-m + f0-mm - 12 f000 )
                                     + Cm ( fmmm + f-m-m-m + fmm-m + f-m-mm + fm-mm + f-mm-m + fm-m-m + f-mmm - 8 f000 )

     where  fijk = f ( x+i.h , y+j.h , z+k.h )    and
 

     A1 = +22997 / 120                  B1 = - 2869 / 60                C1 = + 10
     A2 = - 12709 / 420                  B2 = + 667 / 240               C2 = - 5 / 56
     A3 = + 858391 / 136080         B3 = - 15007 / 68040        C3 = + 5 / 1701
     A4 = - 137843 / 161280          B4 = + 5119 / 322560       C4 = - 5 / 43008
     A5 = + 894317 / 15750000     B5 = - 739 / 1125000        C5 = + 1 / 328125
 
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D3 f
                                         R03 = z        R06 to R12:  temp   ( R07 = 6 f(x,y,z) )
Flags:  /
Subroutine:         A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.
 

-Line 89 is a three-byte GTO 05
 
 

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

 
       ( 330 bytes / SIZE 013 )
 
 

      STACK        INPUTS      OUTPUTS
           T          h > 0             /
           Z             z             /
           Y             y             /
           X             x     D3 f(x,y,z)

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

 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN

 
      T  ASTO 00

-If you choose h = 0.1

   0.1   ENTER^
    3     ENTER^
    2     ENTER^
    1     XEQ "THRM"  >>>> D3 f (1,2,3) = 133.675100 = R05                                                 ---Execution time = 2m42s---

-With a similar program, the HP48 - which works with 12 digits - gives  133.533179
 

Notes:

-Fourth derivatives are already difficult to evaluate numerically but here,
 there will be seldom more than 3 or 4 significant digits !
-The function f is evaluated 131 times.

-Actually, the formulas are of order 12 except for the first 3 terms.

-In the variant hereunder, formulas of order 12 are also used to estimate  6f / x6 + 6f / y6 + 6f / z6
-The coefficients B & C are unchanged, but Ai  become ( there are now 6 coefficients A ):
 

     A1 = +7026 / 35                  B1 = - 2869 / 60                C1 = + 10
     A2 = - 80523 / 2240            B2 = + 667 / 240               C2 = - 5 / 56
     A3 = + 75151 / 8505           B3 = - 15007 / 68040        C3 = + 5 / 1701
     A4 = - 28907 / 17920          B4 = + 5119 / 322560       C4 = - 5 / 43008
     A5 = + 21293 / 109375       B5 = - 739 / 1125000        C5 = + 1 / 328125
     A6 = - 695 / 60480
 
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D3 f
                                         R03 = z        R06 to R12:  temp   ( R07 = 6 f(x,y,z) )
Flags:  /
Subroutine:         A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.
 

-Line 99 is a three-byte GTO 05
 
 

 01 LBL "THRM"
 02 STO 01
 03 RDN
 04 STO 02
 05 RDN
 06 STO 03
 07 X<>Y
 08 STO 04
 09 STO 06
 10 X<>Y
 11 R^
 12 R^
 13 XEQ IND 00
 14 6
 15 *
 16 STO 07
 17 7
 18 ST* 06
 19 SF 10
 20 XEQ 01
 21 RCL 08
 22 695
 23 *
 24 60480
 25 /
 26 CHS
 27 STO 05
 28 XEQ 01
 29 24
 30 *
 31 RCL 09
 32 5173		  33 *
 34 -
 35 RCL 08
 36 1533096
 37 *
 38 +
 39 7875 E3
 40 /
 41 ST+ 05
 42 XEQ 01
 43 75
 44 *
 45 RCL 09         
 46 10238
 47 *
 48 -
 49 RCL 08
 50 1040652
 51 *
 52 +
 53 645120
 54 /
 55 ST- 05
 56 XEQ 01
 57 200
 58 *
 59 RCL 09
 60 15007
 61 *
 62 -
 63 RCL 08
 64 601208		  65 *
 66 +
 67 68040
 68 /
 69 ST+ 05
 70 XEQ 01
 71 600
 72 *
 73 RCL 09
 74 18676
 75 *
 76 -
 77 RCL 08         
 78 241569
 79 *
 80 +
 81 6720
 82 /
 83 ST- 05
 84 XEQ 01
 85 10
 86 *
 87 RCL 09
 88 20083
 89 *
 90 RCL 08
 91 84312
 92 *
 93 -
 94 420
 95 /
 96 -		  97 RCL 05
 98 +
 99 GTO 05
100 LBL 01
101 RCL 04
102 ST- 06
103 CLX
104 STO 08
105 STO 09
106 STO 10
107 LBL 02
108 RCL 03
109 RCL 06
110 +
111 RCL 02
112 RCL 01
113 XEQ IND 00
114 ST+ 08
115 RCL 03
116 RCL 02
117 RCL 06
118 +
119 RCL 01
120 XEQ IND 00
121 ST+ 08
122 RCL 03
123 RCL 02
124 RCL 01
125 RCL 06
126 +
127 XEQ IND 00
128 ST+ 08		 129 RCL 06
130 CHS
131 STO 06
132 X<0?
133 GTO 02
134 STO 11
135 RCL 07
136 ST- 08
137 FS?C 10
138 RTN
139 LBL 03
140 RCL 03
141 RCL 06
142 +
143 RCL 02
144 RCL 11
145 +
146 RCL 01
147 XEQ IND 00
148 ST+ 09
149 RCL 03
150 RCL 06
151 +
152 RCL 02
153 RCL 01
154 RCL 11
155 +
156 XEQ IND 00
157 ST+ 09
158 RCL 03
159 RCL 02
160 RCL 06		 161 +
162 RCL 01
163 RCL 11
164 +
165 XEQ IND 00
166 ST+ 09
167 RCL 06
168 CHS
169 STO 06
170 X<0?
171 GTO 03
172 RCL 11
173 CHS
174 STO 11
175 X<0?
176 GTO 03
177 STO 12
178 LBL 04
179 RCL 03
180 RCL 06
181 +
182 RCL 02
183 RCL 11
184 +
185 RCL 01
186 RCL 12
187 +
188 XEQ IND 00
189 ST+ 10
190 RCL 06
191 CHS
192 STO 06		 193 X<0?
194 GTO 04
195 RCL 11
196 CHS
197 STO 11         
198 X<0?
199 GTO 04
200 RCL 12
201 CHS
202 STO 12
203 X<0?
204 GTO 04
205 RCL 10
206 RCL 07
207 ST+ X
208 ST- 09
209 1.5
210 /
211 -
212 RTN
213 LBL 05
214 RCL 04
215 X^2
216 ST/ Y
217 X^2
218 /
219 STO 05
220 END

 
       ( 355 bytes / SIZE 013 )
 
 

      STACK        INPUTS      OUTPUTS
           T          h > 0             /
           Z             z             /
           Y             y             /
           X             x     D3 f(x,y,z)

 
Example:    the same one:    f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3
 
 

 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN

 
      T  ASTO 00

-If you choose h = 0.1

   0.1   ENTER^
    3     ENTER^
    2     ENTER^
    1     XEQ "THRM"  >>>> D3 f (1,2,3) = 133.688400 = R05                                                 ---Execution time = 2m49s---
 

Notes:

-The HP48 - which works with 12 digits - gives  133.533179 and the previous version returned 133.6751

-The exact result, given by  Wolframalpha  is    ( 21647416 Ln 7 + 579720 ) / ( 117649 e ) = 133.53104241.....

-So, the precision is slightly lower !
-But it's perhaps not significant ?
 
 
Reference:
 
[1]  Abramowitz & Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4