Taylor Series

# Taylor Series for the HP-41

Overview

1°)  Program#1
2°)  Program#2
3°)  Program#3

-The Taylor expansion of a function f is  f(x+h) = f(x) + h f'(x) + h2 f''(x) / 2! + ..... + hn f(n)(x) / n! + ................

-Given a function f(x), we seek approximations of   a1 = f'(x) , a2 = f''(x) / 2! , .... , an = f(n)(x) / n!
-The following programs calculate  ak  up to k = 10

1°)  Program#1

-To approximate the derivatives, a polynomial of degree 10 is used which gives - theoretically - perfect accuracy for polynomials of degree <= 10
-Practically - due to roundoff-errors - the precision decreases as k increases and the estimation of a10 is often very doubtful.

Formulae:

-Let  A = f(x+h) - f(x-h)  ,  B = f(x+2h) - f(x-2h)  ,  C = f(x+3h) - f(x-3h)  ,  D = f(x+4h) - f(x-4h)  , E = f(x+5h) - f(x-5h)

G = f(x+h) + f(x-h)  ,  H = f(x+2h) + f(x-2h)  ,  I = f(x+3h) + f(x-3h)  ,  J = f(x+4h) + f(x-4h)  , K = f(x+5h) + f(x-5h)

and   F = f(x)

-We have:

h f'(x)  ~  ( 2100 A - 600 B + 150 C - 25 D + 2 E ) / 2520
h2 f"(x) ~  ( -73766 F + 42000 G - 6000 H + 1000 I - 125 J + 8 K ) / 25200
h3 f"'(x) ~  ( -70098 A + 52428 B - 14607 C + 2522 D - 205 E ) / 30240
h4 f(4)(x) ~  ( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / 15120
h5 f(5)(x) ~  ( 1938 A - 1872 B + 783 C - 152 D + 13 E ) / 288
h6 f(6)(x) ~  ( -233244 F + 184110 G - 88920 H + 24795 I -3610 J + 247 K ) / 4560
h7 f(7)(x) ~  ( -378 A + 408 B - 207 C + 52 D - 5 E ) / 24
h8 f(8)(x) ~  ( 462 F - 378 G + 204 H - 69 I + 13 J - K ) / 3
h9 f(9)(x) ~  ( 42 A - 48 B + 27 C - 8 D + E ) / 2
h10 f(10)(x) ~  ( -252 F + 210 G - 120 H + 45 I - 10 J + K )

Data Registers:                              R10 thru R24: temp

>>>>   When the program stops:   R00 = a0 = f(x) ,  R01 = a1 ,  .......................... , R10 = a10

Flags: /
Subroutine:  A program that takes x in X-register and returns f(x) in X-register

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

( 503 bytes / SIZE 025 )

 STACK INPUTS OUTPUTS Alpha function name function name T / a3 Z / a2 Y h a1 X x a0

R00 = a0 , R01 = a1 , ................... ,  R10 = a10

Example:   f(x) = exp(x)  for x = 1

 01  LBL "T"  02  E^X  03  RTN

-Place  "T"  in the alpha "register"

-If you choose h = 0.2

0.2  ENTER^
1    XEQ "TAYLR"  >>>>  a0 = 2.718281828 = R00                                    ---Execution time = 24s---
RDN   a1 = 2.718281832 = R01
RDN   a2 = 1.359140923 = R02
RDN   a3 = 0.453046944 = R03  and in R04 thru R10:

 Approximations Exact values R00 = 2.718281828    R01 = 2.718281832    R02 = 1.359140923    R03 = 0.453046944    R04 = 0.113261615    R05 = 0.022652271    R06 = 0.003775851    R07 = 0.000539241    R08 = 0.000064587    R09 = 0.000007643    R10 = 0.000001884 2.718281828        2.718281828        1.359140914        0.453046971        0.113261743        0.022652349        0.003775391        0.000539342        0.000067418        0.000007491        0.000000749

Notes:

-The exact values are  ak =  exp(1) / k!
-Registers R00 thru R09 may be used by the subroutine.

>>>  h between  0.1 to 0.2 is "often" a good choice.

2°)  Program#2

-This variant uses a polynomial of degree 12 to approximate the first 10 derivatives.
-So, we can hope a better precision.

Formulae:

-Let  A = f(x+h) - f(x-h)  ,  B = f(x+2h) - f(x-2h)  ,  C = f(x+3h) - f(x-3h)  ,  D = f(x+4h) - f(x-4h)  , E = f(x+5h) - f(x-5h) , F = f(x+6h) - f(x-6h)

G = f(x+h) + f(x-h)  ,  H = f(x+2h) + f(x-2h)  ,  I = f(x+3h) + f(x-3h)  ,  J = f(x+4h) + f(x-4h)  , K = f(x+5h) + f(x-5h) , L = f(x+6h) + f(x-6h)

and   M = f(x)

-We have:

h f'(x)  ~  ( 23760 A - 7425 B + 2200 C - 495 D + 72 E - 5 F ) / 27720
h2 f"(x) ~  ( -2480478 M + 1425600 G - 222750 H + 44000 I - 7425 J + 864 K - 50 L ) / 302400
h3 f"'(x) ~  ( -764208 A + 603315 B - 198760 C + 46296 D - 6840 E + 479 F ) / 30240
h4 f(4)(x) ~  ( 6222216 M - 4585248 G + 1809945 H - 397520 I + 69444 J - 8208 K + 479 L ) / 453600
h5 f(5)(x) ~  ( 99744 A - 101559 B + 48176 C - 12500 D + 1936 E - 139 F ) / 12096
h6 f(6)(x) ~  ( -3735732 M + 2992320 G - 1523385 H + 481760 I - 93750 J + 11616 K - 695 L ) / 60480
h7 f(7)(x) ~  ( -11652 A + 13275 B - 7550 C + 2404 D - 410 E + 31 F ) / 480
h8 f(8)(x) ~  ( 84084 M - 69912 G + 39825 H - 15100 I + 3606 J - 492 K + 31 L ) / 360
h9 f(9)(x) ~  ( 216 A - 261 B + 164 C - 60 D + 12 E - F ) / 4
h10 f(10)(x) ~  ( -7644 M + 6480 G - 3915 H + 1640 I - 450 J + 72 K - 5 L ) / 12

Data Registers:                              R10 thru R26: temp

>>>>   When the program stops:   R00 = a0 = f(x) ,  R01 = a1 ,  .......................... , R10 = a10

Flags: /
Subroutine:  A program that takes x in X-register and returns f(x) in X-register

 01  LBL "TAYLR"   02  ASTO 10   03  STO 11   04  X<>Y   05  STO 12            06  6   07  STO 25   08  24   09  STO 26   10  LBL 00   11  RCL 11            12  RCL 12   13  RCL 25   14  *   15  STO 14   16  +   17  XEQ IND 10   18  STO 13   19  RCL 11   20  RCL 14   21  -   22  XEQ IND 10   23  RCL 13   24  X<>Y   25  +   26  STO IND 26   27  RCL 13   28  LASTX   29  -   30  DSE 26   31  STO IND 26   32  DSE 26    33  DSE 25   34  GTO 00   35  RCL 11   36  XEQ IND 10    37  STO 00   38  CLA   39  ARCL 10    40  RCL 21   41  72   42  *    43  RCL 23   44  5    45  *   46  -   47  RCL 19   48  495   49  *    50  -   51  RCL 17   52  2200   53  *   54  +   55  RCL 15   56  7425   57  STO 02   58  *   59  - 60  RCL 13   61  23760   62  *   63  +   64  27720   65  /   66  STO 01   67  RCL 22            68  864   69  *   70  RCL 24   71  50   72  *   73  -   74  RCL 20   75  RCL 02   76  *   77  -   78  RCL 18   79  44 E3   80  *   81  +   82  RCL 16   83  222750   84  *   85  -   86  RCL 14   87  1425600   88  *   89  +   90  RCL 00   91  2480478   92  *   93  -   94  1663200   95  /   96  STO 02   97  RCL 23   98  479   99  STO 03 100  * 101  RCL 21 102  6840 103  * 104  - 105  RCL 19 106  46296 107  * 108  + 109  RCL 17 110  198760 111  * 112  - 113  RCL 15 114  603315 115  * 116  + 117  RCL 13 118  764208 119  * 120  - 121  10 122  FACT 123  STO 10          124  2 125  / 126  / 127  X<> 03 128  RCL 24 129  * 130  RCL 22 131  8208 132  * 133  - 134  RCL 20 135  69444 136  * 137  + 138  RCL 18 139  397520 140  * 141  - 142  RCL 16 143  1809945 144  * 145  + 146  RCL 14 147  4585248 148  * 149 - 150  RCL 00 151  6222216 152  * 153  + 154  RCL 10 155  3 156  * 157  / 158  STO 04 159  RCL 21 160  1936 161  * 162  RCL 23 163  139 164  * 165  - 166  RCL 19 167  12500 168  * 169  - 170  RCL 17 171  48176 172  * 173  + 174  RCL 15 175  101559 176  * 177  - 178  RCL 13 179  99744 180  * 181  + 182  RCL 10          183  .4 184  * 185  STO 09 186  / 187  STO 05 188  RCL 22 189  11616 190  * 191  RCL 24 192  695 193  * 194  - 195  RCL 20 196  93750 197 * 198  - 199  RCL 18 200  481760 201  * 202  + 203  RCL 16 204  1523385 205  * 206  - 207  RCL 14 208  2992320 209  * 210  + 211  RCL 00 212  3735732 213  * 214  - 215  RCL 10 216  12 217  * 218  STO 10 219  / 220  STO 06 221  RCL 23 222 31 223  * 224  RCL 21 225  410 226  * 227  - 228  RCL 19 229  2404 230  * 231  + 232  RCL 17 233  7550 234  * 235  - 236  RCL 15 237  13275 238  * 239  + 240  RCL 13          241  11652 242  * 243  - 244  RCL 10 245  18 246  / 247  / 248  STO 07 249  RCL 24 250  31 251  * 252  RCL 22 253  492 254  * 255  - 256  RCL 20 257  3606 258  * 259  + 260  RCL 18 261  15100 262  * 263  - 264  RCL 16 265  39825 266  * 267  + 268  RCL 14 269  69912 270  * 271  - 272  RCL 00 273  84084 274  * 275  + 276  RCL 10 277  3 278  / 279  / 280  STO 08 281  RCL 21 282  12 283  * 284  RCL 23 285  - 286  RCL 19 287  60 288  * 289  - 290  RCL 17 291  164 292  * 293  + 294  RCL 15 295  261 296  * 297  - 298  RCL 13          299  216 300  * 301  + 302  RCL 09 303  / 304  STO 09 305  RCL 22 306  72 307  * 308  RCL 24 309  5 310  * 311  - 312  RCL 20 313  450 314  * 315  - 316  RCL 18 317  1640 318  * 319  + 320  RCL 16 321  3915 322  * 323  - 324  RCL 14 325  6480 326  * 327  + 328  RCL 00 329  7644 330  * 331  - 332  RCL 10 333  / 334  STO 10 335  1.01 336  STO 25 337  RCL 12 338  ENTER^ 339  ENTER^ 340  ENTER^ 341  LBL 01 342  ST/ IND 25 343  * 344  ISG 25 345  GTO 01 346  RCL 03 347  RCL 02 348  RCL 01 349  RCL 00 350  END

( 649 bytes / SIZE 027 )

 STACK INPUTS OUTPUTS Alpha function name function name T / a3 Z / a2 Y h a1 X x a0

R00 = a0 , R01 = a1 , ................... ,  R10 = a10

Example:   f(x) = exp(x)  for x = 1

 01  LBL "T"  02  E^X  03  RTN

-Place  "T"  in the alpha "register"

-If you choose h = 0.2

0.2  ENTER^
1    XEQ "TAYLR"  >>>>  a0 = 2.718281828 = R00                                    ---Execution time = 30s---
RDN   a1 = 2.718281831 = R01
RDN   a2 = 1.359140933 = R02
RDN   a3 = 0.453046944 = R03  and in R04 thru R10:

 Approximations Exact values R00 = 2.718281828    R01 = 2.718281831    R02 = 1.359140933    R03 = 0.453046944    R04 = 0.113261730    R05 = 0.022652332    R06 = 0.003774779    R07 = 0.000539305    R08 = 0.000064587    R09 = 0.000007266    R10 = 0.000002243 2.718281828        2.718281828        1.359140914        0.453046971        0.113261743        0.022652349        0.003775391        0.000539342        0.000067418        0.000007491        0.000000749

Notes:

-Registers R00 thru R09 may be used by the subroutine.
-Several results remain disappointing !

3°)  Program#3

-The 2 routines above employ the same h-value for all the derivatives.
-But a good choice for f'(x) may be a bad choice for f'''(x)

-The program below calls "TAYLR" with a starting h-value and calls "TAYLR" again with 0.8 h , 0.82 h , .............
until the sequence  | di+1 - di |  is no more decreasing, where di is one of the Taylor coefficients.

-So, "TLR+" can choose a better  h-value than the initial one, which may be different for different derivatives.

Data Registers:                              R10 thru R50: temp

>>>>   When the program stops:   R00 = a0 = f(x) ,  R01 = a1 ,  .......................... , R10 = a10

Flags:  F01 to F10 are cleared by this routine.
Subroutine:  A program that takes x in X-register and returns f(x) in X-register and "TAYLR" ( paragraph 1°) or 2°) above )

 01  LBL "TLR+"  02  XEQ "TAYLR"  03  1.02701  04  REGMOVE  05  RCL 12             06  .8  07  *  08  RCL 11  09  XEQ "TAYLR"  10  1.03701  11  REGMOVE  12  10  13  STO 47 14  LBL 00  15  SF IND X  16  DSE X  17  GTO 00  18  LBL 01  19  RCL 12             20  .8  21  *  22  RCL 11  23  XEQ "TAYLR"  24  10  25  STO 48  26  36 27  STO 49  28  +  29  STO 50             30  LBL 02  31  FC? IND 48  32  GTO 03  33  RCL IND 48   34  ENTER^  35  X<> IND 50  36  ST- Y  37  ENTER^  38  X<> IND 49  39  - 40  ABS  41  X<>Y  42  ABS  43  X>Y?   44  CF IND 48  45  X>Y?  46  DSE 47  47  LBL 03   48  DSE 50  49  DSE 49  50  DSE 48  51  GTO 02  52  RCL 47 53  X#0?  54  GTO 01  55  27.00101  56  REGMOVE  57  RCL 03  58  RCL 02  59  RCL 01  60  RCL 00             61  END

( 136 bytes / SIZE 051 )

 STACK INPUTS OUTPUTS Alpha function name function name T / a3 Z / a2 Y h a1 X x a0

R00 = a0 , R01 = a1 , ................... ,  R10 = a10

Example:   f(x) = exp(x)  for x = 1

 01  LBL "T"  02  E^X  03  RTN

-Place  "T"  in the alpha "register"

-If you choose h = 0.4

0.4  ENTER^
1    XEQ "TLR+"  >>>>  a0 = 2.718281828 = R00                                    ---Execution time = 2mn20s---
RDN   a1 = 2.718281829 = R01
RDN   a2 = 1.359140917 = R02
RDN   a3 = 0.453046975 = R03  and in R04 thru R10:

 Approximations Exact values R00 = 2.718281828    R01 = 2.718281829    R02 = 1.359140917    R03 = 0.453046975    R04 = 0.113261790    R05 = 0.022652341    R06 = 0.003775285    R07 = 0.000539348    R08 = 0.000067420    R09 = 0.000007449    R10 = 0.000000755 2.718281828        2.718281828        1.359140914        0.453046971        0.113261743        0.022652349        0.003775391        0.000539342        0.000067418        0.000007491        0.000000749

Notes:

-Except R06 the results are much more satisfactory, especially R09 & R10 !
-Don't choose a too small h-value.

-Lines 06 & 20 may be replaced by another coefficient: 0.9 or 0.7
-Registers R00 thru R09 may be used by the subroutine.

-Here, we have used the 2nd version of "TAYLR" ( paragraph 2°) )
-If you have used the 1st version listed in paragraph 1°) , you got the following results:

 Approximations Exact values R00 = 2.718281828    R01 = 2.718281830    R02 = 1.359140935    R03 = 0.453046948    R04 = 0.113261703    R05 = 0.022652494    R06 = 0.003775442    R07 = 0.000539309    R08 = 0.000067370    R09 = 0.000007752    R10 = 0.000000786 2.718281828        2.718281828        1.359140914        0.453046971        0.113261743        0.022652349        0.003775391        0.000539342        0.000067418        0.000007491        0.000000749

Notes:

-Though the approximations are slightly less accurate, the precision remains acceptable.
-The calculations may be performed again with another initial h-value.

>>>>  In the PPC ROM user's manual, we find a slightly different termination criterion:

-The iterations stops if the relation    0 <= ( di+1 - di ) / ( di - di-1 ) < 1  is no more satisfied.

-If you prefer this method, replace lines 40 to 46 by

X=0?          GTO 02      LBL 02
GTO 02      1                 CF IND 48
/                  X>Y?          DSE 47
X<0?          GTO 03

-With these modifications, we get, with the 2nd version of "TAYLR":

 Approximations Exact values R00 = 2.718281828    R01 = 2.718281829    R02 = 1.359140917    R03 = 0.453046970    R04 = 0.113261790    R05 = 0.022652341    R06 = 0.003775285    R07 = 0.000539348    R08 = 0.000067420    R09 = 0.000007499    R10 = 0.000000755 2.718281828        2.718281828        1.359140914        0.453046971        0.113261743        0.022652349        0.003775391        0.000539342        0.000067418        0.000007491        0.000000749

-So, the results are almost identical, but slightly better !  ( in R03 & R09 )

Reference:

[1]   "PPC ROM user's manual"