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) + h² f''(x) / 2! + ..... + hⁿ f⁽ⁿ⁾(x) / n! + ................

-Given a function f(x), we seek approximations of a₁ = f'(x) , a₂ = f''(x) / 2! , .... , a_n = f⁽ⁿ⁾(x) / n!
-The following programs calculate a_k 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 a₁₀ 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
h² f"(x) ~ ( -73766 F + 42000 G - 6000 H + 1000 I - 125 J + 8 K ) / 25200
h³ f"'(x) ~ ( -70098 A + 52428 B - 14607 C + 2522 D - 205 E ) / 30240
h⁴ f⁽⁴⁾(x) ~ ( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / 15120
h⁵ f⁽⁵⁾(x) ~ ( 1938 A - 1872 B + 783 C - 152 D + 13 E ) / 288
h⁶ f⁽⁶⁾(x) ~ ( -233244 F + 184110 G - 88920 H + 24795 I -3610 J + 247 K ) / 4560
h⁷ f⁽⁷⁾(x) ~ ( -378 A + 408 B - 207 C + 52 D - 5 E ) / 24
h⁸ f⁽⁸⁾(x) ~ ( 462 F - 378 G + 204 H - 69 I + 13 J - K ) / 3
h⁹ f⁽⁹⁾(x) ~ ( 42 A - 48 B + 27 C - 8 D + E ) / 2
h¹⁰ f⁽¹⁰⁾(x) ~ ( -252 F + 210 G - 120 H + 45 I - 10 J + K )

Data Registers: R10 thru R24: temp

>>>> When the program stops: R00 = a₀ = f(x) , R01 = a₁ , .......................... , R10 = a₁₀

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	/	a₃
Z	/	a₂
Y	h	a₁
X	x	a₀

R00 = a₀ , R01 = a₁ , ................... , R10 = a₁₀

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" >>>> a₀ = 2.718281828 = R00                                    ---Execution time = 24s---
                                   RDN   a₁ = 2.718281832 = R01
                                   RDN   a₂ = 1.359140923 = R02
                                   RDN   a₃ = 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 a_k = 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
h² f"(x) ~ ( -2480478 M + 1425600 G - 222750 H + 44000 I - 7425 J + 864 K - 50 L ) / 302400
h³ f"'(x) ~ ( -764208 A + 603315 B - 198760 C + 46296 D - 6840 E + 479 F ) / 30240
h⁴ f⁽⁴⁾(x) ~ ( 6222216 M - 4585248 G + 1809945 H - 397520 I + 69444 J - 8208 K + 479 L ) / 453600
h⁵ f⁽⁵⁾(x) ~ ( 99744 A - 101559 B + 48176 C - 12500 D + 1936 E - 139 F ) / 12096
h⁶ f⁽⁶⁾(x) ~ ( -3735732 M + 2992320 G - 1523385 H + 481760 I - 93750 J + 11616 K - 695 L ) / 60480
h⁷ f⁽⁷⁾(x) ~ ( -11652 A + 13275 B - 7550 C + 2404 D - 410 E + 31 F ) / 480
h⁸ f⁽⁸⁾(x) ~ ( 84084 M - 69912 G + 39825 H - 15100 I + 3606 J - 492 K + 31 L ) / 360
h⁹ f⁽⁹⁾(x) ~ ( 216 A - 261 B + 164 C - 60 D + 12 E - F ) / 4
h¹⁰ f⁽¹⁰⁾(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 = a₀ = f(x) , R01 = a₁ , .......................... , R10 = a₁₀

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	/	a₃
Z	/	a₂
Y	h	a₁
X	x	a₀

R00 = a₀ , R01 = a₁ , ................... , R10 = a₁₀

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" >>>> a₀ = 2.718281828 = R00                                    ---Execution time = 30s---
                                   RDN   a₁ = 2.718281831 = R01
                                   RDN   a₂ = 1.359140933 = R02
                                   RDN   a₃ = 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.8² h , .............
until the sequence | d_i+1 - d_i | is no more decreasing, where d_i 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 = a₀ = f(x) , R01 = a₁ , .......................... , R10 = a₁₀

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	/	a₃
Z	/	a₂
Y	h	a₁
X	x	a₀

R00 = a₀ , R01 = a₁ , ................... , R10 = a₁₀

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+" >>>> a₀ = 2.718281828 = R00                                    ---Execution time = 2mn20s---
                                RDN   a₁ = 2.718281829 = R01
                                RDN   a₂ = 1.359140917 = R02
                                RDN   a₃ = 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 <= ( d_i+1 - d_i ) / ( d_i - d_i-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"

hp41programs

Taylor Series for the HP-41