mathr / blog / #

The essence of perturbation is to find the difference between the high precision values of a function at two nearby points, while using only the low precision value of the difference between the points. In this post I'll write the high precision points in CAPITALS and the low precision deltas in lowercase. There are two auxiliary operations needed to define the perturbation \(P\), \(B\) replaces all variables by their high precision version, and \(W\) replaces all variables by the sum of the high precision version and the low precision delta. Then \(P = W - B\):

\[\begin{aligned} B(f) &= f(X) &\text{ (emBiggen)}\\ W(f) &= f(X + x) &\text{ (Widen)}\\ P(f) &= W(f) - B(f) \\ &= f(X + x) - f(X) &\text{ (Perturb)} \end{aligned}\]

For example, perturbation of \(f(z, c) = z^2 + c\), ie, \(P(f)\), works out like this:

\[\begin{aligned} & P(f) \\ \to & f(Z + z, C + c) - f(z, c) \\ \to & (Z + z)^2 + (C + c) - (Z^2 + C) \\ \to & Z^2 + 2 Z z + z^2 + C + c - Z^2 - C \\ \to & 2 Z z + z^2 + c \end{aligned}\]

where in the final result the additions of \(Z\) and \(z\) have mostly cancelled out and all the terms are "small".

For polynomials, regular algebraic manipulation can lead to successful outcomes, but for other functions it seems some "tricks" are needed. For example, \(|x|\) (over \(\mathbb{R}\)) can be perturbed with a "diffabs" function proceeding via case analysis:

// evaluate |X + x| - |X| without catastrophic cancellation
function diffabs(X, x) {
  if (X >= 0) {
    if (X + x >= 0) { return x; }
    else { return -(2 * X + x); }
  } else {
    if (X + x > 0) { return 2 * X + x; }
    else { return -x; }
  }
}

This formulation was developed by laser blaster at fractalforums.com.

For transcendental functions, other tricks are needed. Here for example is a derivation of \(P(\sin)\):

\[\begin{aligned} & P(\sin) \\ \to & \sin(X + x) - \sin(X) \\ \to & \sin(X) \cos(x) + \cos(X) \sin(x) - \sin(X) \\ \to & \sin(X) (\cos(x) - 1) + \cos(X) \sin(x) \\ \to & \sin(X) \left(-2\sin^2\left(\frac{x}{2}\right)\right) + \cos(X) \sin(x) \\ \to & \sin(X) \left(-2\sin^2\left(\frac{x}{2}\right)\right) + \cos(X) \left(2 \cos\left(\frac{x}{2}\right) \sin\left(\frac{x}{2}\right)\right) \\ \to & 2 \sin\left(\frac{x}{2}\right) \left(-\sin(X) \sin\left(\frac{x}{2}\right) + \cos(X) \cos\left(\frac{x}{2}\right)\right) \\ \to & 2 \sin\left(\frac{x}{2}\right) \cos\left(X + \frac{x}{2}\right) \end{aligned}\]

Knowing when to apply the sum- and double-angle-formulae, is a bit of a mystery, especially if the end goal is not known beforehand. This makes implementing a symbolic algebra program that can perform these derivations quite a challenge.

In lieu of a complete symbolic algebra program that does it all on demand, here are a few formulae that I calculated, some by hand, some using Wolfram Alpha:

\[\begin{aligned} P(a) &= 0 \\ P(a f) &= a P(f) \\ P(f + g) &= P(f) + P(g) \\ P(f g) &= P(f) W(g) + B(f) P(g) \\ P(f^{n+1}) &= P(f) \sum_{k=0}^n W(f^k) B(f)^{n-k} \\ P\left(\frac{1}{f}\right) &= -\frac{P(f)}{B(f)W(f)} \\ P\left(\frac{f}{g}\right) &= \frac{P(f) B(g) - B(f) P(g)}{B(g) W(g)} \\ P(|f|) &= \operatorname{diffabs}(B(f), P(f)) \\ P(\exp) &= \exp(X) \operatorname{expm1}(x) \\ P(\exp \circ f) &= 2 \sinh\left(\frac{P(f)}{2}\right)\exp\left(\frac{W(f)+B(f)}{2}\right) \\ P(\log) &= \operatorname{log1p}\left(\frac{x}{X}\right) \\ P(\sin \circ f) &= \phantom{-}2 \sin\left(\frac{P(f)}{2}\right)\cos\left(\frac{W(f)+B(f)}{2}\right) \\ P(\cos \circ f) &= -2 \sin\left(\frac{P(f)}{2}\right)\sin\left(\frac{W(f)+B(f)}{2}\right) \\ P(\tan \circ f) &= \frac{\sin(P(f))}{\cos(B(f))\cos(W(f))} \\ P(\sinh \circ f) &= 2 \sinh\left(\frac{P(f)}{2}\right)\cosh\left(\frac{W(f)+B(f)}{2}\right) \\ P(\cosh \circ f) &= 2 \sinh\left(\frac{P(f)}{2}\right)\sinh\left(\frac{W(f)+B(f)}{2}\right) \\ P(\tanh \circ f) &= \frac{\sinh(P(f))}{\cosh(B(f))\cosh(W(f))} \\ P(\sec \circ f) &= \frac{2\sin(\frac{P(f)}{2})\sin(\frac{W(f)+B(f)}{2})}{\cos(B(f))\cos(W(f))} \\ P(\csc \circ f) &= \frac{-2\sin(\frac{P(f)}{2})\cos(\frac{W(f)+B(f)}{2})}{\sin(B(f))\sin(W(f))} \\ P(\cot \circ f) &= \frac{\cos(P(f))}{\sin(B(f))\sin(W(f))} \\ \end{aligned}\]

I hope to find time to add these to et soon.

EDIT there is a simpler and more general way to derive \(P(\sin)\) and so on, using \(\sin(a) \pm \sin(b)\) formulae...

EDIT 2019-12-23 added \(P(\exp \circ f)\) via \(\exp = \sinh + \cosh\)

EDIT 2020-05-28 added \(P(f^n)\) via induction

EDIT 2021-09-29 added \(\sec, \csc, \cot\) thanks to FractalAlex

EDIT 2021-10-05 added \(\frac{f}{g}\) by combining the product and reciprocal rules