# Convenient error analysis

Is there any convenient way to find out the error in evaluating an expression numerically, providing all errors of the arguments are given? Like

Sin[DataError[123, .1] + 3*DataError[456, .2]/DataError[789, .5]]^2


DataError[xxx, xx]

• You might consider using Interval[] for this application. – J. M. will be back soon Dec 12 '16 at 14:51
• What would be the form of the result you have in mind? – Daniel Lichtblau Dec 12 '16 at 15:51
• @DanielLichtblau , like what's written above, plug in all tuples of data and it's error, giving out a tuple of result and calculated error (according to error analysis and pde stuff) – ZisIsNotZis Dec 12 '16 at 15:54

It was not stated how DataError is to be defined, and in particular whether the second argument is to be taken as an absolute or relative error. I am assuming absolute.

Easiest is, as suggested by @J.M., to use Interval.

dataError1[a_, b_] := Interval[{a - b, a + b}]
val1 = Sin[
dataError1[123, .1] + 3*dataError1[456, .2]/dataError1[789, .5]]^2

(* Out= Interval[{0.542555735765, 0.736481835669}] *)


Could use explicit optimization to avoid the "dependency problem" (which does not arise in this example).

dataError2[a_, b_, eps_] := a + b*eps
val2 = Sin[
dataError2[123, .1, eps1] +
3*dataError2[456, .2, eps2]/dataError2[789, .5, eps3]]^2;

evars = {eps1, eps2, eps3};
NMinimize[{val2, Thread[-1 <= evars <= 1]}, evars]
NMaximize[{val2, Thread[-1 <= evars <= 1]}, evars]

(* Out= {0.542555735765, {eps1 -> 1., eps2 -> 1., eps3 -> -1.}}

Out= {0.736481835669, {eps1 -> -1., eps2 -> -1., eps3 -> 1.}} *)


If the error offsets are fairly tight, could do this as a first-order estimate. Much better for speed, in complicated cases.

firstorder =
Normal[Series[val2 /. Thread[evars -> t*evars], {t, 0, 1}]] /. t -> 1
NMinimize[{firstorder, Thread[-1 <= evars <= 1]}, evars]
NMaximize[{firstorder, Thread[-1 <= evars <= 1]}, evars]

(* Out= -0.095854777084 eps1 - 0.000728933666038 eps2 +
0.00105321214104 eps3 + Sin[32805/263]^2

Out= {0.544829012526, {eps1 -> 1., eps2 -> 1., eps3 -> -1.}}

Out= {0.740102858308, {eps1 -> -1., eps2 -> -1., eps3 -> 1.}} *)


One could use the built-in arbitrary precision numerics, I suppose. I don't know the accuracy of Accuracy[], but it seems to be close if an approximate error is what is desired and not hard bounds. The smaller the relative error, the better it will be. Accuracy is measure on a log-10 scale. See also SetAccuracy[].

Sin[SetAccuracy[123, 10^-.1] +
3*SetAccuracy[456, 10^-.2]/SetAccuracy[789, 10^-.5]]^2
% + {1, -1} Log10[Accuracy[%]]
(*
0.*10^-1
{0.549593, 0.735339}
*)


We can see from the 0.*10^-1, whose FullForm is 0.6424659354173337650.6153213200244898`, is that one drawback is that the error can overwhelm the computed value. If it does, Mathematica will return a value of zero and this method will fail.