(* :Copyright: © 1994 by Dmitri Linde Permission is hereby granted to make copies of this file for any purpose other than direct profit, or as part of a commercial product, provided this copyright notice is left intact. *) (* :Title: ErrorPropagation *) (* :Author: Dmitri Linde dmitri@physics.stanford.edu *) (* :Summary: This package provides tools to do calculations on numbers with uncertainties. Calculations like: (2.0 +/- 0.1)*(3.0 +/- 0.01) are performed automatically once the package is loaded. *) (* :Discussion: Once this package is loaded, a number together with its uncertainty can be represented in Mathematica using the +/-. For example: In[2]:= 2.5535345 +/- 0.01242 -2 Out[2]= (2.554 +/- 1.2 10 ) or equivalently by using the function error[x,y] In[3]:= error[2.5535345,0.01242] -2 Out[3]= (2.554 +/- 1.2 10 ) or by using the function percenterror[x,y] (the uncertainty of x is y percent) In[4]:= percenterror[50,10] Out[4]= (50 +/- 5.) +/- in the output stands for plus-minus. Although full accuracy is retained, notice that for printing the number is automatically rounded within its uncertainty. Simple calculations with errors will work automatically. For example In[8]:= (2.0 +/- 0.1)*(3.0 +/- 0.01) -1 Out[8]= (6. +/- 3. 10 ) In[9]:= Exp[%] 2 2 Out[9]= (4. 10 +/- 1.2 10 ) In[10]:= Log[% + 32.4 +/- 0.5] -1 Out[10]= (6.08 +/- 2.8 10 ) Uncertainty can be printed in percentage form: In[11]:= PercentForm[%] Out[11]= (6.08 +/- 4.6%) However you cannot do calculations in percent form. Yet. Should you find a function which is not automatically supported by the package, you can use SubstituteError to plug in the number with an uncertainty into that function. For example In[13]:= SubstituteError[BesselI[3,x],x->5.0+/-0.3] 1 Out[13]= (1. 10 +/- 3.4) More complicated functions, with several variables, will work as well: In[14]:= centerOfMass= SubstituteError[ (m1*x1+m2*x2)/(m1+m2), {m1->3.42+/-0.01, m2->4.42+/-0.01, x1->10.0+/-0.4, x2->5.0+/-1}] -1 Out[14]= (7.18 +/- 5.9 10 ) All calculations are performed symbolically, so errors don't have to be numbers: In[16]:= 1/error[x0,epsilon] 2 1 epsilon Out[16]= (-- +/- Sqrt[--------]) x0 4 x0 In[17]:= PowerExpand[%] 1 epsilon Out[17]= (-- +/- -------) x0 2 x0 *) (* :Limitations: When a number with uncertainty is plugged into a function, the function should be approximately linear in the neighborhood of the number. For very fast changing functions or large uncertainties this package will not yield accurate results for uncertainties. *) (* :Package Version: 1.0 *) (* :Warning: This package modifies the $PreRead command *) BeginPackage["ErrorCalculations`"] SubstituteError::usage = "SubstituteError[function,substitutions] plugs in numbers with uncertainties into the function. e.g. SubstituteError[x+Sin[x*y],{x->error[2.0,0.1],y->error[1.4,0.6]}]" PercentForm::usage = "PercentForm[error[A,B]] will calculate and display a number A with percentage uncertainty (B is the absolute uncertainty)" error::usage = "error[x,dx] defines a number x with an uncertainty dx" percenterror::usage "percenterror[x,y] defines a number x with uncertainty y percent" Format[error[x_,y_]]:=DisplayMyError[x,y] Format[error[x_?NumberQ,y_?NumberQ]]:=Module[{myy,myx,precis}, myy=ScientificForm[y,2]; precis=Ceiling[N[x/y]]; If[precis>0,precis=Ceiling[Log[10.0,precis]]]; myx=ScientificForm[x,precis+1]; DisplayMyError[myx,myy] ] /; y!=0 Begin["`Private`"] (* set the private context *) percenterror[x_,y_]:=error[x,x*y/100.0] DisplayMyError[x_,y_]:=StringForm["(`` +/- ``)", x, y] DisplayMyErrorPercent[x_,y_]:=StringForm["(`` +/- ``%)", x, y] PercentForm[error[x_,0]]:=DisplayMyError[x,"0%"] PercentForm[error[x_,y_]]:=DisplayMyErrorPercent[x,Simplify[y/x*100]] PercentForm[error[x_?NumberQ,y_?NumberQ]]:= Module[{myy,myx,precis}, precis=Ceiling[N[x/y]]; If[precis>0,precis=Ceiling[Log[10.0,precis]]]; myx=ScientificForm[x,precis+1]; myy=N[100.0*y/x,2]; DisplayMyErrorPercent[myx,myy] ] /; y!=0 error /: Plus[error[x1_,y1_],error[x2_,y2_]]:=error[x1+x2,Sqrt[y1^2+y2^2]] error /: Plus[error[x1_,y1_],x2_]:=error[x1+x2,y1] error /: Times[error[x1_,y1_],x2_]:=error[x1*x2,y1*x2] error /: Times[error[x1_,y1_],error[x2_,y2_]]:= error[x1*x2 , (x2^2*y1^2 + x1^2*y2^2)^(1/2)] error /: Power[error[x1_,y1_],x2_]:=error[x1^x2, (x1^(-2 + 2*x2)*x2^2*y1^2)^(1/2)] error /: Power[error[x1_,y1_],error[x2_,y2_]]:= error[x1^x2 , (x1^(-2 + 2*x2)*x2^2*y1^2 + x1^(2*x2)*y2^2*Log[x1]^2)^(1/2)] error /: Exp[error[x_,y_]]:=error[E^x, Abs[E^x*y]] error /: Log[error[x_,y_]]:=error[Log[x], Abs[y/x]] error /: Sin[error[x_,y_]]:=error[Sin[x], Abs[y*Cos[x]]] error /: Cos[error[x_,y_]]:=error[Cos[x], Abs[-(y*Sin[x])]] error /: Tan[error[x_,y_]]:=error[Tan[x], Abs[y*Sec[x]^2]] SubstituteError[f_,subs_List]:= Module[{n,vars,realstuff,err,i,myhead,symberr}, n=Length[subs]; realstuff=Map[(#[[2,1]])&,subs]; err=Map[(#[[2,2]])&,subs]; vars=Map[(#[[1]])&,subs]; symberr=error[f,Sqrt[Sum[(D[f,vars[[i]]]*err[[i]])^2,{i,1,n}]]]; symberr/.MapThread[Rule,{vars,realstuff}] ] SubstituteError[f_,subs_Rule]:=SubstituteError[f,{subs}] $PreRead = StringReplace[#, "+/-" -> "~error~"]& (* Provides shorthand input form. Thanks J. Keiper and D. Withoff and Michael Ibrahim *) End[] (* end the private context *) Attributes[SubstituteError]={ReadProtected} EndPackage[] (* end the package context *)