Automatic differentiation using dual numbers

[Thomas Klemm]

As already mentioned we can extend the idea of automatic differentiation to higher derivatives.
If we represent the trial numbers by matrices we again can use the HP-42S.

Basis

We use the following basis:

1 0 0
0 1 0
0 0 1
STO "1"

0 1 0
0 0 1
0 0 0
STO "e"

0 0 1
0 0 0
0 0 0
STO "e↑2"

Constructors

We can use "TRIAL" to enter a trial number of the form \(x +\varepsilon\):

7
XEQ "TRIAL"

7 1 0
0 7 1
0 0 7

Or then, if we want to specify the 3 coefficients we can use "NEW":

1
ENTER
2
ENTER
3
XEQ "NEW"

1 2 3
0 1 2
0 0 1

Display

We can use "SHOW" to display the 3 coefficients.
I recommend to use FIX 3:

XEQ "SHOW"

2.000 3.000 5.000
x: [ 3×3 Matrix ]

This function is also called at the end of "TRIAL" and "NEW":

Code:

00 { 71-Byte Prgm }
01▸LBL "TRIAL"
02 RCL× "1"
03 RCL+ "e"
04 GTO 00
05▸LBL "NEW"
06▸LBL 01
07 RCL× "e↑2"
08 X<>Y
09 RCL× "e"
10 +
11 X<>Y
12 RCL× "1"
13 +
14▸LBL "SHOW"
15▸LBL 00
16 CLA
17 EDIT
18 ARCL ST X
19 ├" "
20 →
21 ARCL ST X
22 ├" "
23 →
24 ARCL ST X
25 EXITALL
26 AVIEW
27 END

Addition and Subtraction

Since adding a number to a matrix is done on each element, we can not simply add constants.
Instead we have to multiply the constant by the "1" matrix.

Multiplication and Division

However we can use the ordinary multiplication and division both for constants and trial numbers.

Thus without much programming we can already calculate first and second derivatives of polynomials or even rational functions.

Example

\(
\begin{align}
p(x) = 2x^3 + 3x^2 + 5x + 7
\end{align}
\)

Code:

00 { 38-Byte Prgm }
01▸LBL "p(x)"
02 RCL "x"
03 2
04 ×
05 3
06 RCL× "1"
07 +
08 RCL× "x"
09 5
10 RCL× "1"
11 +
12 RCL× "x"
13 7
14 RCL× "1"
15 +
16 END

Let's assume that we want to evaluate the polynomial at \(x = -1\) and calculate its 1st and 2nd derivative:

-1
XEQ "TRIAL"
STO "x"
XEQ "p(x)"
XEQ "SHOW"

3.000 5.000 -3.000
x: [ 3×3 Matrix ]

Thus we get:

\(
\begin{align}
p(x) &= 3 \\
{p}'(x) &= 5 \\
\tfrac{1}{2}{p}''(x) &= -3 \\
\end{align}
\)

Halley's Method

It's straightforward to implement Halley's method:

\(
\begin{align}
x_{{n+1}}=x_{n}-{\frac{2f(x_{n})f'(x_{n})}{2{[f'(x_{n})]}^{2}-f(x_{n})f''(x_{n})}}
\end{align}
\)

However we rather use this formula:

\(
\begin{align}
x_{n+1}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})-{\frac {f(x_{n})}{f'(x_{n})}}{\frac {f''(x_{n})}{2}}}} \\
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1-{\frac {f(x_{n})}{f'(x_{n})}}\cdot {\frac {f''(x_{n})}{2f'(x_{n})}}\right]^{-1} \\
\end{align}
\)

Code:

00 { 56-Byte Prgm }#
01▸LBL "HALLEY"    #
02 STO 03          # x -> R03
03 RCL× "1"        # x*1
04 RCL+ "e"        # x*1+e
05 STO "x"         # x*1+e -> x
06 XEQ "p(x)"      #
07 EDIT            # matrix -> register
08 STO 00          # f -> R00
09 →               # f'
10 STO 01          # f' -> R01
11 →               # f"/2
12 STO 02          # f"/2 -> R02
13 EXITALL         #
14 RCL 00          # f
15 RCL÷ 01         # f/f'
16 1               # 1 f/f'
17 RCL ST Y        # f/f' 1 f/f'
18 RCL× 02         # f/f'*f"/2 1 f/f'
19 RCL÷ 01         # f/f'*f"/2/f' 1 f/f'
20 -               # 1-f/f'*f"/2/f' f/f'
21 ÷               # f/f'/[1-f/f'*f"/2/f']
22 STO- 03         # x-dx -> x
23 RCL 03          # x
24 END             #

Now we can use it to find the root of the aforementioned polynomial:

-1 XEQ "HALLEY" R/S R/S …

-1.441176470588235294117647058823529
-1.445528443543314459809602861040273
-1.445528458679522302862795272237467
-1.445528458679522302862795910270395
-1.445528458679522302862795910270395

We notice that the rate of convergence is cubic compared to the quadratic convergence of Newton's method.

Dottie Number

Code:

00 { 19-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 XEQ "COS"
04 RCL- "x"
05 END

1 XEQ "HALLEY"

0.7408739950803435700746289353295154
0.7390851338775818843562183300927446
0.7390851332151606416553120877075439
0.7390851332151606416553120876738734
0.7390851332151606416553120876738734


Time Value of Money

Consider the following parameters:

Quote:n, pv, pmt, fv = T.new(10), T.new(50), T.new(-30), T.new(100)

Code:

00 { 46-Byte Prgm }
01▸LBL "TVM"
02 10
03 STO "n"
04 50
05 RCL× "1"
06 STO "pv"
07 -30
08 RCL× "1"
09 STO "pmt"
10 100
11 RCL× "1"
12 STO "fv"
13 END

We want to find the root \(x\) of the following function:

Quote:function f(x) return ((pv+fv)/T.expm1(T.log1p(x)*n) + pv)*x + pmt end

Code:

00 { 50-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 XEQ "LN1+X"
04 RCL× "n"
05 XEQ "E↑X-1"
06 RCL "pv"
07 RCL+ "fv"
08 X<>Y
09 ÷
10 RCL+ "pv"
11 RCL× "x"
12 RCL+ "pmt"
13 END

-0.3 XEQ "HALLEY" R/S R/S …
-0.284444698106904440553487332229463
-0.2844359988802575605180815267548689
-0.2844359988802559255723584032186149
-0.284435998880255925572358403218615
-0.284435998880255925572358403218615


Implementation

I haven't implemented all standard functions (yet) but I hope it is still useful to get an idea on how this can be done.

Code:

00 { 397-Byte Prgm }#
01▸LBL "HALLEY"     #
02 STO 03           # x -> R03
03 RCL× "1"         #
04 RCL+ "e"         #
05 STO "x"          #
06 XEQ "f(x)"       #
07 XEQ 02           # matrix -> register
08 RCL 00           # f
09 RCL÷ 01          # f/f'
10 1                # 1 f/f'
11 RCL ST Y         # f/f' 1 f/f'
12 RCL× 02          # f/f'*f"/2 1 f/f'
13 RCL÷ 01          # f/f'*f"/2/f' 1 f/f'
14 -                # 1-f/f'*f"/2/f' f/f'
15 ÷                # f/f'/[1-f/f'*f"/2/f']
16 STO- 03          # x-dx -> x
17 RCL 03           # x
18 RTN              #
19▸LBL "X↑2"        #
20 ENTER            #
21 ×                #
22 RTN              #
23▸LBL "E↑X"        # e^f(x) + y e^f(x) f'(x) + 1/2 y^2 e^f(x) (f''(x) + f'(x)^2)
24 XEQ 02           # matrix -> register
25 RCL 02           # x"/2
26 RCL 01           # x' x"/2
27 X↑2              # x'^2 x"/2
28 2                # 2 x'^2 x"/2
29 ÷                # x'^2/2 x"/2
30 +                # x'^2/2+x"/2
31 RCL× "e↑2"       # (x'^2/2+x"/2)*e^2
32 RCL 01           # x' (x'^2/2+x"/2)*e^2
33 RCL× "e"         # x'*e (x'^2/2+x"/2)*e^2
34 +                # x'*e+(x'^2/2+x"/2)*e^2
35 RCL "1"          # 1 x'*e+(x'^2/2+x"/2)*e^2
36 +                # 1+x'*e+(x'^2/2+x"/2)*e^2
37 RCL 00           # x 1+x'*e+(x'^2/2+x"/2)*e^2
38 E↑X              # exp(x) 1+x'*e+(x'^2/2+x"/2)*e^2
39 ×                # exp(x)[1+x'*e+(x'^2/2+x"/2)*e^2]
40 RTN              #
41▸LBL "LN"         # log(f(x)) + (y f'(x))/f(x) - (y^2 (f'(x)^2 - f(x) f''(x)))/(2 f(x)^2)
42 XEQ 02           # matrix -> register
43 RCL 01           # x'
44 RCL÷ 00          # x'/x
45 RCL× "e"         # x'/x*e
46 LASTX            # x'/x x'/x*e
47 X↑2              # [x'/x]^2 x'/x*e
48 2                # 2 [x'/x]^2 x'/x*e
49 ÷                # [x'/x]^2/2 x'/x*e
50 RCL 02           # x"/2 [x'/x]^2/2 x'/x*e
51 RCL÷ 00          # x"/x/2 [x'/x]^2/2 x'/x*e
52 -                # -x"/x/2+[x'/x]^2/2 x'/x*e
53 RCL× "e↑2"       # (-x"/x/2+[x'/x]^2/2)*e^2 x'/x*e
54 -                # x'/x*e+(x"/x/2-[x'/x]^2/2)*e^2
55 RCL 00           # x x'/x*e+(x"/x/2-[x'/x]^2/2)*e^2
56 LN               # log(x) x'/x*e+(x"/x/2-[x'/x]^2/2)*e^2
57 RCL× "1"         # log(x)*1 x'/x*e+(x"/x/2-[x'/x]^2/2)*e^2
58 +                # log(x)*1+x'/x*e+(x"/x/2-[x'/x]^2/2)*e^2
59 RTN              #
60▸LBL "E↑X-1"      # (e^f(x) - 1) + y e^f(x) f'(x) + 1/2 y^2 e^f(x) (f''(x) + f'(x)^2)
61 XEQ 02           # matrix -> register
62 RCL 02           # x"/2
63 RCL 01           # x' x"/2
64 X↑2              # x'^2 x"/2
65 2                # 2 x'^2 x"/2
66 ÷                # x'^2/2 x"/2
67 +                # x'^2/2+x"/2
68 RCL× "e↑2"       # (x'^2/2+x"/2)*e^2
69 RCL 01           # x' (x'^2/2+x"/2)*e^2
70 RCL× "e"         # x'*e (x'^2/2+x"/2)*e^2
71 +                # x'*e+(x'^2/2+x"/2)*e^2
72 RCL 00           # x x'*e+(x'^2/2+x"/2)*e^2
73 E↑X              # exp(x) x'*e+(x'^2/2+x"/2)*e^2
74 ×                # exp(x)[x'*e+(x'^2/2+x"/2)*e^2]
75 RCL 00           # x exp(x)[x'*e+(x'^2/2+x"/2)*e^2]
76 E↑X-1            # exp(x)-1 exp(x)[x'*e+(x'^2/2+x"/2)*e^2]
77 RCL× "1"         # [exp(x)-1]*1 exp(x)[x'*e+(x'^2/2+x"/2)*e^2]
78 +                # [exp(x)-1]*1+exp(x)[x'*e+(x'^2/2+x"/2)*e^2]
79 RTN              #
80▸LBL "LN1+X"      # log(f(x) + 1) + (y f'(x))/(f(x) + 1) + (y^2 ((f(x) + 1) f''(x) - f'(x)^2))/(2 (f(x) + 1)^2)
81 XEQ 02           # matrix -> register
82 1                # 1
83 RCL+ 00          # 1+x
84 RCL 01           # x' 1+x
85 RCL÷ ST Y        # x'/(1+x) 1+x
86 RCL× "e"         # x'/(1+x)*e 1+x
87 LASTX            # x'/(1+x) x'/(1+x)*e 1+x
88 X↑2              # [x'/(1+x)]^2 x'/(1+x)*e 1+x
89 2                # 2 [x'/(1+x)]^2 x'/(1+x)*e 1+x
90 ÷                # [x'/(1+x)]^2/2 x'/(1+x)*e 1+x 1+x
91 RCL 02           # x"/2 [x'/(1+x)]^2/2 x'/(1+x)*e 1+x
92 RCL÷ ST T        # x"/(1+x)/2 [x'/(1+x)]^2/2 x'/(1+x)*e
93 -                # -x"/(1+x)/2+[x'/(1+x)]^2/2 x'/(1+x)*e
94 RCL× "e↑2"       # (-x"/(1+x)/2+[x'/(1+x)]^2/2)*e^2 x'/(1+x)*e
95 -                # x'/(1+x)*e+(x"/(1+x)/2-[x'/(1+x)]^2/2)*e^2
96 RCL 00           # x x'/(1+x)*e+(x"/(1+x)/2-[x'/(1+x)]^2/2)*e^2
97 LN1+X            # log(1+x) x'/(1+x)*e+(x"/(1+x)/2-[x'/(1+x)]^2/2)*e^2
98 RCL× "1"         # log(1+x)*1 x'/(1+x)*e+(x"/(1+x)/2-[x'/(1+x)]^2/2)*e^2
99 +                # log(1+x)*1+x'/(1+x)*e+(x"/(1+x)/2-[x'/(1+x)]^2/2)*e^2
100 RTN             #
101▸LBL "SIN"       # sin(f(x)) + y f'(x) cos(f(x)) + 1/2 y^2 (f''(x) cos(f(x)) - f'(x)^2 sin(f(x)))
102 XEQ 02          # matrix -> register
103 RCL 00          # x
104 1               # 1 x
105 →REC            # cos(x) sin(x)
106 RCL 01          # x' cos(x) sin(x)
107 X↑2             # x'^2 cos(x) sin(x)
108 RCL× ST Z       # sin(x)*x'^2 cos(x) sin(x)
109 -2              # -2 sin(x)*x'^2 cos(x) sin(x)
110 ÷               # -sin(x)*x'^2/2 cos(x) sin(x)
111 RCL 02          # x"/2 -sin(x)*x'^2/2 cos(x) sin(x)
112 RCL× ST Z       # cos(x)*x"/2 -sin(x)*x'^2/2 cos(x) sin(x)
113 +               # cos(x)*x"/2-sin(x)*x'^2/2 cos(x) sin(x)
114 RCL× "e↑2"      # (cos(x)*x"/2-sin(x)*x'^2/2)*e^2 cos(x) sin(x)
115 X<>Y            # cos(x) (cos(x)*x"/2-sin(x)*x'^2/2)*e^2 sin(x)
116 RCL× 01         # cos(x)*x' (cos(x)*x"/2-sin(x)*x'^2/2)*e^2 sin(x)
117 RCL× "e"        # cos(x)*x'*e (cos(x)*x"/2-sin(x)*x'^2/2)*e^2 sin(x)
118 +               # cos(x)*x'*e+(cos(x)*x"/2-sin(x)*x'^2/2)*e^2 sin(x)
119 X<>Y            # sin(x) cos(x)*x'*e+(cos(x)*x"/2-sin(x)*x'^2/2)*e^2
120 RCL× "1"        # sin(x)*1 cos(x)*x'*e+(cos(x)*x"/2-sin(x)*x'^2/2)*e^2
121 +               # sin(x)*1+cos(x)*x'*e+(cos(x)*x"/2-sin(x)*x'^2/2)*e^2
122 RTN             #
123▸LBL "COS"       # cos(f(x)) - y f'(x) sin(f(x)) + 1/2 y^2 (f'(x)^2 (-cos(f(x))) - f''(x) sin(f(x)))
124 XEQ 02          # matrix -> register
125 RCL 00          # x
126 1               # 1 x
127 →REC            # cos(x) sin(x)
128 RCL 01          # x' cos(x) sin(x)
129 X↑2             # x'^2 cos(x) sin(x)
130 RCL× ST Y       # cos(x)*x'^2 cos(x) sin(x)
131 2               # 2 cos(x)*x'^2 cos(x) sin(x)
132 ÷               # cos(x)*x'^2/2 cos(x) sin(x)
133 RCL 02          # x"/2 cos(x)*x'^2/2 cos(x) sin(x)
134 RCL× ST T       # sin(x)*x"/2 cos(x)*x'^2/2 cos(x) sin(x)
135 +               # sin(x)*x"/2+cos(x)*x'^2/2 cos(x) sin(x)
136 RCL× "e↑2"      # (sin(x)*x"/2+cos(x)*x'^2/2)*e^2 cos(x) sin(x)
137 +/-             # -(sin(x)*x"/2+cos(x)*x'^2/2)*e^2 cos(x) sin(x)
138 X<>Y            # cos(x) -(sin(x)*x"/2+cos(x)*x'^2/2)*e^2 sin(x)
139 RCL× "1"        # cos(x)*1 -(sin(x)*x"/2+cos(x)*x'^2/2)*e^2 sin(x)
140 +               # cos(x)*1-(sin(x)*x"/2+cos(x)*x'^2/2)*e^2 sin(x)
141 X<>Y            # sin(x) cos(x)*1-(sin(x)*x"/2+cos(x)*x'^2/2)*e^2
142 RCL× 01         # sin(x)*x' cos(x)*1-(sin(x)*x"/2+cos(x)*x'^2/2)*e^2
143 RCL× "e"        # sin(x)*x'*e cos(x)*1-(sin(x)*x"/2+cos(x)*x'^2/2)*e^2
144 -               # cos(x)*1-sin(x)*x'*e-(sin(x)*x"/2+cos(x)*x'^2/2)*e^2
145 RTN             #
146▸LBL 02          # matrix -> register
147 EDIT            # f
148 STO 00          # f -> R00
149 →               # f'
150 STO 01          # f' -> R01
151 →               # f"/2
152 STO 02          # f"/2 -> R02
153 EXITALL         #
154 RTN             #
155▸LBL "TRIAL"     #
156 RCL× "1"        #
157 RCL+ "e"        #
158 GTO 00          #
159▸LBL "NEW"       #
160▸LBL 01          #
161 RCL× "e↑2"      #
162 X<>Y            #
163 RCL× "e"        #
164 +               #
165 X<>Y            #
166 RCL× "1"        #
167 +               #
168▸LBL "SHOW"      #
169▸LBL 00          #
170 CLA             #
171 EDIT            #
172 ARCL ST X       #
173 ├" "            #
174 →               #
175 ARCL ST X       #
176 ├" "            #
177 →               #
178 ARCL ST X       #
179 EXITALL         #
180 AVIEW           #
181 END             #