• Initialization
• Specification of Quantities with Error
  • A Note about the Display of Maple Floating-Point Numbers
• Simple Error Analysis - Propagation of Errors
• The Formulas of Error Propagation
• Working with Rounding Rules
  • Introduction
  • Another Note about Floating-Point Numbers
  • Adding a New Rounding Rule
• Error Analysis with Physical Constants
  • The amplitude of a simple pendulum
  • Ideal Gas Law
  • Working at higher Digits
• Error in a Derived Physical Constant: electron mass

The ScientificErrorAnalysis Package in Maple 9

by Maplesoft

The new ScientificErrorAnalysis  package provides representation and construction of numerical quantities that have a value  and associated error , which is a measure of the degree of precision to which the quantity's value is known. Given the errors of base quantities, you can automatically compute the errors of quantities derived from them, as well as correlations and covariances between them.

Initialization

>

restart:

>	with( ScientificConstants  ); with( ScientificErrorAnalysis  );

Warning, the names GetError and GetValue have been rebound

Specification of Quantities with Error

The same quantity-with-error can be entered three different ways.

>	Quantity( 20.0, 0.1 );

>	Quantity( 20.0, 0.005, 'relative' );

>	Quantity( 20.0, 1, 'uld' );

Other examples:

>	Quantity( 3.4e15, 0.2e15 );

>	Quantity( 5.1e-20, 0.1, 'relative' );

>	Quantity( 4.2012883, 12, 'uld' );

Important: These quantities do not  represent intervals  in which possible values must be contained - they represent unknown values with central tendency.

A Note about the Display of Maple Floating-Point Numbers

If a floating-point number is not  displayed in scientific notation, it may not be apparent whether trailing zeros are significant figures.

>	Quantity( 4.0e2, 1, 'uld');

>	SFloatMantissa( GetValue( % ) ), SFloatMantissa( GetError( % ) );

Simple Error Analysis - Propagation of Errors

Errors can be propagated through calculations by using the combine/errors  command.

>	q1  := Quantity( 10.0, 1.0 );

>	q2  := Quantity( 20.0, 3.0 );

Error in the sum of q1 and q2:

>	SumQ1Q2  := combine( q1+q2 , 'errors' );

Relative error in their sum:

>	GetError( SumQ1Q2 ) / GetValue( SumQ1Q2 );

Error in their difference:

>	DifferenceQ1Q2  := combine( q2-q1 , 'errors' );

Relative error in their difference:

>	GetError( DifferenceQ1Q2  ) / GetValue( DifferenceQ1Q2  );

As expected, the relative error in the difference is greater than that in the sum.

Error in the product of the above two quantities:

>	combine( q1*q2 , 'errors' );

Relative error in their product:

>	GetError(%)/GetValue(%);

Error in their quotient:

>	combine( q2/q1 , 'errors' );

Relative error in their quotient:

>	GetError(%)/GetValue(%);

As expected, the relative errors in the product and quotient are the same.

Error in a power:

>	combine( q1^2 , 'errors' );

Error in differentiable functions:

>	combine( sin( q1  ), 'errors' );

>	combine( exp( q1/q2  ), 'errors' );

Important: The above calculations are not   interval arithmetic.  That is, the quantities involved do not  represent intervals  in which values must be contained.

The Formulas of Error Propagation

The usual expressions for error propagation through simple functions are easily obtained.

>	e1  := Quantity( x, delx );

>	e2  := Quantity( y, dely );

Error a sum:

>	combine( e1+e2 , 'errors' );

Error a product:

>	combine( e1*e2 , 'errors' );

Relative error^2 in the product:

>	expand(GetError(%)^2/GetValue(%)^2);

Error in a square:

>	combine( e1^2 , 'errors' );

General Formula

The general formula for the error   in y , where y  is a function of variables  is:

where  is the error in  and the partials are evaluated at the central values of the  .

With covariances  taken into account, the formula is:

The covariance  may be expressed in terms of the correlation  and errors ,  as follows:

Working with Rounding Rules

Introduction

The initial default rounding rule is digits , for which no rounding is performed.

>	combine( q1+q2 , 'errors' );

We can explicitly round a Quantity( )  object:

>	ApplyRule( %, round[2] );

specify a rounding rule  in combine/errors :

>	combine( q1+q2 , 'errors', 'rule'=round[2] );

>	combine( q1*q2 , 'errors', 'rule'=round[2] );

>	combine( q1^2 , 'errors', 'rule'=round[2] );

or change the default rounding rule using the UseRule  function.

>	UseRule( round[2] );

>	combine( q1+q2 , 'errors' );

>	combine( q1*q2 , 'errors' );

>	combine( q1^2 , 'errors' );

>	UseRule( digits );

Another Note about Floating-Point Numbers

If a floating-point number is not  displayed in scientific notation, it may not be apparent whether trailing zeros are significant figures.

>	Quantity( 23900.0, 1854.0 );

>	ApplyRule( %, round[2] );

>	SFloatMantissa( GetValue( % ) ), SFloatMantissa( GetError( % ) );

Adding a New Rounding Rule

The rounding rule round3g [n]  rounds the error to n  figures for a significant figure of 3  or greater and n+1  figures for a significant figure of   1  or 2 . Here, we add a similar rule round2g[n] , which rounds the error to n  figures for a significant figure of 2  or greater and n+1  figures for a significant figure of 1 .

>	AddRule( round2g  =

>	proc( x, y )

>	local xp, yp, n;

>	if not type( procname, 'indexed' ) then

>	error "`round2g must be indexed";

>

end if;

>	n := op( procname );

>	if not type( n, 'posint' ) then

>	error "invalid index to `round2g`";

>

end if;

>	xp, yp := evalf( x ), evalf( y );

>	if iquo( SFloatMantissa( yp ), 10^ilog10( SFloatMantissa( yp ) ) ) < 2 then

>

n := n+1;

>

end if;

>	`if`( xp<>0, evalf[n + ilog10( SFloatMantissa( xp ) ) - ilog10( SFloatMantissa( yp ) ) + SFloatExponent( xp ) - SFloatExponent( yp )]( xp ), 0.0 ), evalf[n]( yp );

>

end proc

>

);

We can use the new rule round2g[n] .

>	HasRule( round2g  );

>	Quantity( 11.22, 2.3 );

>	ApplyRule( %, round2g [1] );

>	Quantity( 11.22, 1.3 );

>	ApplyRule( %, round2g [1] );

Error Analysis with Physical Constants

The amplitude of a simple pendulum

The period of oscillation T  of a simple pendulum for small amplitude is given by:

where g  is the standard acceleration of gravity and L  is the length of the string of the pendulum (or, more precisely, the distance from the point of suspension to the center of the pendulum bob).

A student constructs a pendulum to determine a value for g .  The student measures L  to be 153.2 centimeters (cm).  It is decided that the uncertainty in this measurement is primarily due to the difficulty in measuring L  to the center of the bob, and that this uncertainty is about plus or minus 2 millimeters (mm).  (In other words, it is estimated that the measurement is probably within 2 millimeters of the correct value.)  Thus, L  is taken to be 153.2 +/- 0.2 centimeters.  In standard units:

>	Lq  := Quantity( 153.2, 2, 'uld' );

Now, the student sets the pendulum swinging with a small amplitude (a few centimeters).  Using a stopwatch, the student measures the time for 30 oscillations to be 75.3 seconds (s).  The uncertainty in this measurement is due to the human reaction time in operating the stopwatch, and so the time measurement is taken to be 75.3 +/- 0.5 seconds.

>	t30  := Quantity( 75.3, 0.5 );

Because the number of oscillations (30) is exact, the period of a single oscillation is:

>	Tq  := combine( t30 /30, 'errors', rule=round[2] );

The previous formula rearranges to:

and so, the student calculates the experimentally determined value of g  to be:

>	combine( 4*3.1416^2* Lq/Tq^2 , 'errors', rule=round[2] );

The accepted value of g  is:

>	evalf( Constant( g ) );

Assuming the error in the experimental value represents a standard deviation, the accepted value thus lies farther than one standard deviation from the experimental central value.  The student then must decide whether to repeat the experiment, or consider possibly overlooked sources of measurement error.

Ideal Gas Law

A sample of a gas is measured to have a volume of 20.3 liters (L), pressure of 105 kilopascals (kPa), and temperature of 297 kelvins (K).  If all the measurements have uncertainty of 1%, what is the value and uncertainty of the number of moles in the sample?

The ideal gas law is:

where P  is pressure, V  is volume, T  is temperature, and R  is the ideal gas constant.

Constructing quantities-with-error for P, V, and T in standard units:

>	Pq  := Quantity( 105e3, 0.01, 'relative' );

>	Vq  := Quantity( 20.3e-3, 0.01, 'relative' );

>	Tq  := Quantity( 297.0, 0.01, 'relative' );

Thus, the number of moles n  is:

>	combine( Pq*Vq /(Constant( R )* Tq ), 'errors', 'rule'=round[2] );

Thus, the sample has 0.863 moles (mol) of gas molecules, with a relative error of:

>	GetError(%)/GetValue(%);

or about 2%.

Working at higher Digits

Some physical constants have their values determined to more than 10 significant figures, hence calculations involving these objects at the default setting of Digits may result in a loss of precision.  If more precision is desired, Digits may be set to a higher value, 15 for example.

>	GetConstant( R[infinity] );

>	combine( Constant( R[infinity] )*Constant( c ), 'errors', 'rule=round[2]' );

>	Digits:=15:

>	combine( Constant( R[infinity] )*Constant( c ), 'errors', 'rule=round[2]' );

>	Digits:=10:

Error in a Derived Physical Constant: electron mass

Let us examine the formula for the electron mass.

>	GetConstant( m[e] );

As a demonstration, we manually construct the same expression using individual physical constant objects.

>	e3  := 2*Constant( R[infinity] )*Constant( h )/( Constant( c )*Constant( alpha )^2 );

Propagating the errors gives us:

>	combine( e3 , 'errors', 'rule=round[2]' );

The above calculation is similar to the method ScientificConstants  and ScientificErrorAnalysis  use internally.

>	(GetValue,GetError)( Constant( m[e] ) );

>	ApplyRule( Quantity( % ), 'round[2]' );