Numeric error propagation: Difference between revisions

=={{header|Fortran}}==

===Direct calculation===

Following the propagation of error estimates through a computation is a nightmare of convoluted formulae wherein mistakes are easily made. The basic method is to derive the formulae according to the standard rules while carrying forward the calculation by hand with the utmost caution. A computer can perform the calculations, but the real problem is in ensuring that it performs the correct calculations, not some incorrect calculation... So, rather than attempt to "optimise" the calculation, the objective is to reduce brain strain by producing code whose ''checkability'' is optimised instead, somewhat as follows: <lang Fortran> PROGRAM CALCULATE !A distance, with error propagation.

REAL X1, Y1, X2, Y2 !The co-ordinates.

REAL X1E,Y1E,X2E,Y2E !Their standard deviation.

DATA X1, Y1 ,X2, Y2 /100., 50., 200.,100./ !Specified

DATA X1E,Y1E,X2E,Y2E/ 1.1, 1.2, 2.2, 2.3/ !Values.

REAL DX,DY,D2,D,DXE,DYE,E !Assistants.

CHARACTER*1 C !I'm stuck with code page 437 instead of 850.

PARAMETER (C = CHAR(241)) !Thus ± does not yield this glyph on a "console" screen. CHAR(241) does.

REAL SD !This is an arithmetic statement function.

SD(X,P,S) = P*ABS(X)**(P - 1)*S !SD for X**P where SD of X is S

WRITE (6,1) X1,C,X1E,Y1,C,Y1E, !Reveal the points

1 X2,C,X2E,Y2,C,Y2E !Though one could have used an array...

1 FORMAT ("Euclidean distance between two points:"/ !A heading.

1 ("(",F5.1,A1,F3.1,",",F5.1,A1,F3.1,")")) !Thus, One point per line.

DX = (X1 - X2) !X difference.

DXE = SQRT(X1E**2 + X2E**2) !SD for DX, a simple difference.

DY = (Y1 - Y2) !Y difference.

DYE = SQRT(Y1E**2 + Y2E**2) !SD for DY, (Y1 - Y2)

D2 = DX**2 + DY**2 !The distance, squared.

DXE = SD(DX,2,DXE) !SD for DX**2

DYE = SD(DY,2,DYE) !SD for DY**2

E = SQRT(DXE**2 + DYE**2) !SD for their sum

D = SQRT(D2) !The distance!

E = SD(D2,0.5,E) !SD after the SQRT.

WRITE (6,2) D,C,E !Ahh, the relief.

2 FORMAT ("Distance",F6.1,A1,F4.2) !Sizes to fit the example.

END !Enough.</lang>

This is old-style Fortran, except for the CHARACTER variable caused by problems with character codes and their screen glyphs. As can be seen, the formulae invite mistakes which is why there is no attempt to produce a single arithmetic expression for the result and its error estimate. Further, rather than attempt to emplace appropriate instances of the formula for a value raised to some power (squaring, and square root), risking all manner of misthinks, a function to do so is prepared, here using Fortran's "arithmetic statement function" protocol, expressly intended for such situations. And the results are...

<pre>

Euclidean distance between two points:

(100.0±1.1, 50.0±1.2)

(200.0±2.2,100.0±2.3)

Distance 111.8±2.49

</pre>

All outputs are formatted for the specific values of the test data so as to avoid excessive spacing. Since the given errors are of the order of one, showing more than one decimal digit would be silly. Similarly, standard precision suffices though for more elaborate calculations, double precision could well be preferable for the working out.

===More general===

Rather than agonise over devising adjoint formulae for the error propagation through some calculation, one can perform the desired calculation via routines that will carry along the error term.

===Fortran 90 ''et seq''.===

A latter-day expansion of Fortran makes it possible to define a compound entity such as a value and its associated error, for instance,<lang Fortran> TYPE DATUM

REAL VALUE

REAL SD

END TYPE DATUM

TYPE POINT

TYPE(DATUM) X

TYPE(DATUM) Y

END TYPE POINT

TYPE(POINT) P1,P2

</lang>

Whereupon, instead of a swarm of separate variables named according to some scheme, you can have a collection of variables with subcomponents named systematically. Further, via a great deal of syntax one can devise functions dealing with those compound types and moreover, prepare procedures that will perform operations such as addition and subtraction, etc. and merge these with the ordinary usages of addition, etc. of ordinary variables...