Division Formulae are Wrong : Dual Numbers

[rakeshvar]

Thank you for this extremely wonderful tutorial. I learnt a lot. Also please mention that they need to read backandforth.ipynb after the intro.ipynb. It makes most sense that way. And it took me a while to figure that out.

And now the bug...

struct Dual{T<:Real} <: Real
    x::T
    y::T
end
Base.show(io::IO, d::Dual) = println(io, d.x, d.y≥0 ? "+" : "", d.y, "ϵ")

import Base: +, -, *, /
a::Dual + b::Dual = Dual(a.x + b.x, a.y + b.y)
a::Dual * b::Dual = Dual(a.x * b.x, b.x * a.y + a.x * b.y)
a::Dual / b::Dual = Dual(a.x / b.x,(b.x * a.y - a.x * b.y)/b.x^2)

Base.convert(::Type{Dual{T}}, a::Dual) where T = Dual(convert(T, a.x), convert(T, a.y))
Base.convert(::Type{Dual{T}}, a::Real) where T = Dual(convert(T, a), zero(T))
Base.promote_rule(::Type{Dual{T}}, ::Type{R}) where {T,R} = Dual{promote_type(T,R)}

Now let us try a simple function. I am also giving a numerical derivative for reference.

f(x) = x/(1+x*x)
ndf(x, ϵ=√eps()) = (f(x+ϵ)-f(x-ϵ))/2ϵ
f(7), ndf(7)

(0.14, -0.01919999998062849)

Now let us try with Dual numbers

f(Dual(7., 1.))

0.14-0.0192ϵ

Works!

But with your definition of division,

a::Dual / b::Dual = Dual(a.x * b.x, b.x * a.y - a.x * b.y)
f(Dual(7., 1.))

350.0-48.0ϵ

It is obviously wrong. Also you can see that my formula for division matches with your high school text book.