This is the concept: https://en.wikipedia.org/wiki/Automatic_differentiation It's nicely done here: http://codon.com/automatic-differentiation-in-ruby, where I was looking for help. The implementation is based on a concept of Dual Numbers: https://en.wikipedia.org/wiki/Dual_number. The basic idea is to (not like in numeric differentiation) compute change and value together (Dual Numbers!). The first thing - implement a dual numbers class!
# dual numbers class
import math
class Dual_number:
"""Dual numbers API"""
def __init__(self, _real, _dual=0):
self.real = _real
self.dual = _dual
# present a number to the world:
def __str__(self):
def sign_help(x):
return '+' if self.dual >= 0 else '-'
return str(self.real) + ' ' + sign_help(self.dual) + ' ' + str(abs(self.dual)) + 'e'
def __add__(self, other):
if not isinstance(other, Dual_number):
other = dual_number(other)
return Dual_number(self.real + other.real, self.dual + other.dual)
else:
return Dual_number(self.real + other.real, self.dual + other.dual)
def __radd__(self, other):
return Dual_number(self.real + other, self.dual)
def __sub__(self, other):
if not isinstance(other, Dual_number):
other = dual_number(other)
return Dual_number(self.real - other.real, self.dual - other.dual)
else:
return Dual_number(self.real - other.real, self.dual - other.dual)
def __rsub__(self, other):
return Dual_number(-(self.real - other), -self.dual)
def __mul__(self, other):
if not isinstance(other, Dual_number):
other = dual_number(other)
return Dual_number(self.real * other.real, self.real * other.dual +
self.dual * other.real)
else:
return Dual_number(self.real * other.real, self.real * other.dual +
self.dual * other.real)
def __rmul__(self, other):
return Dual_number(self.real * other, self.dual * other)
def __truediv__(self, other):
if not isinstance(other, Dual_number):
other = dual_number(other)
return Dual_number(self.real / other.real , (self.dual * other.real -
self.real * other.dual) / (other.real * other.real))
else:
return Dual_number(self.real / other.real , (self.dual * other.real -
self.real * other.dual) / (other.real * other.real))
def __rtruediv__(self, other):
return Dual_number((1 / self.real) * other, ((-1 * self.dual) / self.real * self.real) * other )
def dual_number(x):
return Dual_number(x)
# dual number with dual part eq to zero is a real number (similar to complex), constant value
As seen arithmetic operations are constructed in similar way like in complex class, just have in mind, that e is infinitesimal so e*e is zero; add and substract like normal algebraic expressions, multiplication is based on a matrix representation[of dual number]. Let's check if everything's all right.
d3 = Dual_number(1, 1)
d4 = Dual_number(1, 2)
print(d3 * d4)
print(d3 - d4)
print(1 + d3)
print(1 - d3)
print(3 * d3 * 3)
print((3 - d4) * (5 * d3))
print(d3 / d4)
print(d3 + 1.09)
print(d4 - 3)
print(d3 / 2)
print(2 / d3)
1 + 3e 0 - 1e 2 + 1e 0 - 1e 9 + 9e 10 + 0e 1.0 - 1.0e 2.09 + 1e -2 + 2e 0.5 + 0.5e 2.0 - 2.0e
All seem to be fine; see how it works. The function is 2t^2 + 3 and we have:
# distance function:
def s(t):
return 2 * (t * t) + 3
def ds(t):
"""numeric differentiation"""
dt = 0.001
s0 = s(t)
s1 = s(t + dt)
return (s1 - s0)/dt
def ds_aut(x):
"automatic differnetiation"
return 2 * x * x + 3
print("Time = 0:")
t0 = Dual_number(0, 1)
print(ds(0))
print(ds_aut(t0))
print("Time = 4:")
t0 = Dual_number(4, 1)
print(ds(4))
print(ds_aut(t0))
Time = 0: 0.0019999999998354667 3 + 0e Time = 4: 16.0020000000074 35 + 16e
# Another function - tangent
def s(t):
return tan(t)
def ds(t):
"""numeric differentiation"""
dt = 0.001
s0 = s(t)
s1 = s(t + dt)
return (s1 - s0)/dt
def ds_aut(x):
"automatic differnetiation"
return tan(x)
print("Time = 0:")
t0 = Dual_number(0, 1)
print(ds(0))
print(ds_aut(t0))
print("Time = PI / 3:")
t0 = Dual_number(math.pi / 3, 1)
print(ds(math.pi / 3))
print(ds_aut(t0))
Time = 0: 1.0000003333334668 0.0 + 1.0e Time = PI / 3: 4.006941562015198 1.7320508075688767 + 3.9999999999999982e
As seen our dual numbers do the job, we set dual part of time to one: this maybe the change for one second, the real part is the distance travelled or whatever - the point where we calculate the derivative. The cool thing is, that the result is accurate - like in symbolic calculation; but numeric is as acurate as the step (dt) is.
All what is left is define functions, I've done a few here:
def sin(x):
if isinstance(x, Dual_number):
return Dual_number(math.sin(x.real), x.dual * math.cos(x.real))
else:
return math.sin(x)
def cos(x):
if isinstance(x, Dual_number):
return Dual_number(math.cos(x.real), -x.dual * math.sin(x.real))
else:
return math.sin(x)
def tan(x):
"""tangent x for |x| < pi / 2 """
if isinstance(x, Dual_number):
return sin(x) / cos(x)
else:
return math.tan(x)
def cot(x):
if isinstance(x, Dual_number):
return cos(x) / sin(x)
else:
return math.cos(x) / math.sin(x)
def sinh(x):
if isinstance(x, Dual_number):
return Dual_number(math.sinh(x.real), x.dual * math.cosh(x.real))
else:
return math.sinh(x)
def cosh(x):
if isinstance(x, Dual_number):
return Dual_number(math.cosh(x.real), -x.dual * math.sinh(x.real))
else:
return math.cosh(x)
def tanh(x):
if isinstance(x, Dual_number):
return sinh(x) / cosh(x)
else:
return math.tanh(x)
def coth(x):
"""coth if x or x.real not eq to 0"""
if isinstance(x, Dual_number):
return cosh(x) / sinh(x)
else:
return math.cosh(x) / math.sinh(x)
def sqrt(x):
if isinstance(x, Dual_number):
return Dual_number(math.sqrt(x.real), (x.dual) / (2 * sqrt(x.real)))
else:
return math.sqrt(x)
def cube_root(x):
if isinstance(x, Dual_number):
return Dual_number(x.real ** (1/3), (x.dual) / (3 * (x.real ** 2)**(1 / 3)))
else:
return x ** (1 / 3)
def exp(x):
if isinstance(x, Dual_number):
return Dual_number(math.exp(x.real), math.exp(x.real) * x.dual)
else:
return math.exp(x)
def ln(x):
if isinstance(x, Dual_number):
if not x.dual is 0:
return Dual_number(math.log(x.real) , (x.real / x.dual))
else:
return Dual_number(math.log(x.real))
else:
return math.log(x)
def pow(a, x):
if isinstance(x, Dual_number):
if not x.dual is 0:
return Dual_number(math.pow(a, x.real), math.pow(a, x.real) * x.dual * ln(x.dual))
else:
return Dual_number(math.pow(a, x.real))
else:
return math.pow(a, x)
To derrive them formulas we can used Maclaurin expansions with all factors with e^2 and higher dropped. That's all.