Automatic differentiation using dual numbers

[Albert Chan]

Apply f to dual number {g, g'} produce {f(g), (f(g))' = f'(g) * g'}

To emphasize its "symbolic" nature, I use f and g for variables.

r = x^(1/n)
r' = x^(1/n-1)/n = r / (x*n)

Code:

D = {}  -- note: f1, f2 = f, f'
function D.new(f1,f2) return setmetatable({f1,f2 or 0}, D) end
function D.surd(f,n) local r=surd(f[1],n); return D.new(r, f[2]*r/(f[1]*n)) end
function D.__add(f,g) return D.new(f[1]+g[1], f[2]+g[2]) end
function D.__sub(f,g) return D.new(f[1]-g[1], f[2]-g[2]) end
function D.__mul(f,g) return D.new(f[1]*g[1], (f[2]*g[1]+f[1]*g[2])) end
function D.__div(f,g) return D.new(f[1]/g[1], (f[2]*g[1]-f[1]*g[2]) / g[1]^2) end
function D.__unm(f)   return D.new(-f[1], -f[2]) end -- negate
function D.log1p(f) return D.new(log1p(f[1]), f[2]/(f[1]+1)) end
function D.expm1(f) local r=expm1(f[1]); return D.new(r, (r+1)*f[2]) end
function D.newton(x,f) local r=-f[1]/f[2]; x[1]=x[1]+r; return x[1],f[1] end

Example, solve (³√x + x² - 3/5*x = 0), guess x = -1

lua> x = D.new(-1,1)
lua> for k=1,6 do print(D.newton(x, D.surd(x,3) + x*x - D.new(3/5)*x)) end
-0.7352941176470589     0.6
-0.6875931861095272     0.07925115704554603
-0.6858395921931165     0.0027133429258336395
-0.6858371853881947     3.7138631565070135e-006
-0.6858371853836587     6.999345547598068e-012
-0.6858371853836587     0

Note: outputs are (extrapolated x, f(x)); Shift up 2nd column to get points (x, f(x))