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 )
c)  Biharmonic Operator ( 3-Dim )
d)  Biharmonic Operator ( N-Dim )
e)  Biharmonic Operator ( N-Dim ) - Radial Functions

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

-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

 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

 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

 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

 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)xxyy4f / x2y2 =  - 0.038395989                                             ---Execution time = 33s---

-Exact result = - 0.0384

b)  Biharmonic Operator - 2 Dim

-The 1st 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

-This 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 ?