Numeric error propagation: Difference between revisions
→More general: Add the multiply and divide routines.
m (→{{header|REXX}}: moved a statement.) |
(→More general: Add the multiply and divide routines.) |
||
Line 614:
REAL STACKV(VMAX) !Holds the values.
REAL STACKE(VMAX) !And the corresponding error estimate.
LOGICAL VTRACE !Perhaps progress is to be followed in detail.
DATA VSP,VOUT,VTRACE/0,0,.FALSE./ !Start with nothing.
PARAMETER (
CONTAINS !The servants.
INTEGER OUT !Unit number for output.
VSP = 0 !My stack is empty.
VOUT = OUT !Save this rather than have extra parameters.
VTRACE = VOUT .GT. 0 !By implication.
END SUBROUTINE VINIT !Ready.
SUBROUTINE VSHOW(WOT) !Show the topmost element.
CHARACTER*(*) WOT !The caller identifies itself.
IF (VSP.LE.0) THEN !Just in case
WRITE (VOUT,1) "Empty",VSP !My stack may be empty.
ELSE !But normally, it is not.
WRITE (VOUT,1) WOT,VSP,STACKV(VSP),PM,STACKE(VSP) !Topmost.
1 FORMAT (A8,": Vstack(",I2,") =",F8.1,A1,F6.2) !Suits the example.
END IF !So much for protection.
END SUBROUTINE VSHOW !Could reveal all the stack...
SUBROUTINE VLOAD(V,E) !Load the stack.
REAL V,E !The value and its error.
Line 622 ⟶ 642:
STACKV(VSP) = V !Place the value.
STACKE(VSP) = E !And the error.
IF (VTRACE) CALL VSHOW("vLoad")
END SUBROUTINE VLOAD !That was easy!
SUBROUTINE VADD !Add the top two elements.
IF (VSP.LE.1) STOP "VADD: underflow!" !Maybe not.
STACKV(VSP - 1) = STACKV(VSP - 1) + STACKV(VSP) !Do the deed.
STACKE(VSP - 1) = SQRT(STACKE(VSP - 1)**2 + STACKE(VSP)**2) !The errors follow.
VSP = VSP - 1 !Two values have become one.
IF (VTRACE) CALL VSHOW("vAdd")!The result.
END SUBROUTINE VADD !The variance of the sum is the sum of the variances.
SUBROUTINE VSUB !Subtract the topmost element from the one below.
IF (VSP.LE.1) STOP "VSUB: underflow!" !Perhaps not.
STACKV(VSP - 1) = STACKV(VSP - 1) - STACKV(VSP) !The topmost was the second loaded.
STACKE(VSP - 1) = SQRT(STACKE(VSP - 1)**2 + STACKE(VSP)**2) !Add the variances also.
VSP = VSP - 1 !Two values have become one.
IF (VTRACE) CALL VSHOW("vSub")!The result.
END SUBROUTINE VSUB !Could alternatively play with the signs and add...
SUBROUTINE VMUL !Multiply the top two elements.
REAL R1,R2 !Use relative errors in place of plain SD.
IF (VSP.LE.1) STOP "VMUL: underflow!" !Perhaps not.
R1 = STACKE(VSP - 1)/STACKV(VSP - 1) !The relative errors for multiply
R2 = STACKE(VSP) /STACKV(VSP) !Are treated as are variances in addition.
STACKV(VSP - 1) = STACKV(VSP - 1)*STACKV(VSP) !Perform the multiply.
VSP = VSP - 1 !Unstack, but not quite finished.
STACKE(VSP) = SQRT((R1**2 + R2**2)*STACKV(VSP)**2) ![SD/xy]² = [SD/x]² + [SD/y]²
IF (VTRACE) CALL VSHOW("vMul") !Thus SD² = [[SD/x]² + [SD/y]²]xy²
END SUBROUTINE VMUL !The square means that the error's sign is not altered.
SUBROUTINE VDIV !Divide the penultimate element by the top elements.
REAL R1,R2 !Use relative errors in place of plain SD.
IF (VSP.LE.1) STOP "VDIV: underflow!" !Perhaps not.
R1 = STACKE(VSP - 1)/STACKV(VSP - 1) !The relative errors for divide
R2 = STACKE(VSP) /STACKV(VSP) !Are treated as are variances in subtraction.
STACKV(VSP - 1) = STACKV(VSP - 1)/STACKV(VSP) !Perform the divide.
VSP = VSP - 1 !X/Y is Load X, Load Y, Divide; Y is topmost.
STACKE(VSP) = SQRT((R1**2 + R2**2)*STACKV(VSP)**2) ![SD/(x/y)]² = [SD/x]² + [SD/y]²
IF (VTRACE) CALL VSHOW("vDiv") !Thus SD² = [[SD/x]² + [SD/y]²](x/y)²
END SUBROUTINE VDIV !Worry over y ± SD spanning zero...
SUBROUTINE VSQRT !Now for some fun with the topmost element.
IF (VSP.LE.0) STOP "VSQRT: underflow!" !Maybe not.
STACKV(VSP) = SQRT(STACKV(VSP)) !Negative? Let the system complain.
STACKE(VSP) = 0.5/STACKV(VSP)*STACKE(VSP) !F(x ± s) = F(x) ± F'(x).s
IF (VTRACE) CALL VSHOW("vSqrt") !Here, F' can't be negative.
END SUBROUTINE VSQRT !No change to the pointer.
SUBROUTINE VSQUARE !Another raise-to-a-power.
STACKE(VSP) = 2*ABS(STACKV(VSP))*STACKE(VSP) !The error's sign is not to be messed with.
STACKV(VSP) = STACKV(VSP)**2 !This will never be negative.
IF (VTRACE) CALL VSHOW("vSquare") !Keep away from zero though.
▲ END SUBROUTINE VSQUARE !Same formula as VSQRT.
END SUBROUTINE VSQUARE !Same formula as VSQRT, just a different power.
SUBROUTINE VPOW(P) !Now for the more general.
INTEGER P !Though only integer powers for this routine, so no EXP(P*LN(x)).
Line 650 ⟶ 700:
STACKE(VSP) = P*ABS(STACKV(VSP))**(P - 1)*STACKE(VSP) !Negative values a worry.
STACKV(VSP) = ABS(STACKV(VSP))**P !I only want the magnitude.
IF (VTRACE) CALL VSHOW("vPower") !So, what happened?
END SUBROUTINE VPOW !Powers with fractional parts are troublesome.
END MODULE ERRORFLOW !That will do for the test problem.
PROGRAM CALCULATE !A distance, with error propagation.
USE ERRORFLOW !For the details.
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.
▲ 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.
▲ 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.
Calculate SQRT
CALL VLOAD(X2,X2E)
CALL VLOAD(X1,X1E)
Line 678 ⟶ 728:
CALL VADD !(X2 - X1)**2 + (Y2 - Y1)**2
CALL VSQRT !SQRT((X2 - X1)**2 + (Y2 - Y1)**2)
WRITE (6,2) STACKV(1),
2 FORMAT ("Distance",F6.1,A1,F4.2) !Sizes to fit the example.
END !Enough.</lang>
This is closer to the idea of extending the language to supply additional facilities. At the cost of hand-compiling the arithmetic expression into a sequence of pseudo-machine code subroutines, it is apparent that the mind-tangling associated error formulae need no longer be worried over. The various arithmetic subroutines have to be coded correctly with careful attention to the V or E statements in fact being for the V and E terms (cut&paste followed by inadequate adjustment the culprit here: such mistakes are less likely when using a card punch because copying is more troublesome), but this is a straightforward matter of checking. And indeed the VPOW(2) routine has the same effect as VSQUARE.
<pre>
Euclidean distance between two points:
(100.0±1.1, 50.0±1.2)
(200.0±2.2,100.0±2.3)
vLoad: Vstack( 1) = 200.0± 2.20
vLoad: Vstack( 2) = 100.0± 1.10
vSub: Vstack( 1) = 100.0± 2.46
vSquare: Vstack( 1) = 10000.0±491.93
vLoad: Vstack( 2) = 100.0± 2.30
vLoad: Vstack( 3) = 50.0± 1.20
vSub: Vstack( 2) = 50.0± 2.59
vSquare: Vstack( 2) = 2500.0±259.42
vAdd: Vstack( 1) = 12500.0±556.15
vSqrt: Vstack( 1) = 111.8± 2.49
Distance 111.8±2.49
</pre>
Although the example source file uses the F90 MODULE protocol, this is merely for convenience. Otherwise, there would be the same collection of routines, all sharing a COMMON work area.
The trace as the calculation progresses shows a disconcerting surge in the intermediate values of the error estimate. One ought not "compute from the hip" and not worry over intermediate results that are not shown! More generally, the progression of error should be watched carefully lest assumptions prove invalid. For instance, that x ± e is a ''small'' span, within which F(x) and F'(x) do not change greatly, still less pass a discontinuity. Seemingly helpful events such as F'(x) = 0 should be thought about... An alternative approach to F(x ± e) = F(x) ± |F'(x)|.e is to compute F(x - e) and F(x + e) as well as F(x). Some systems offer "interval arithmetic" that might assist with such a procedure, but they are not widely available.
▲Although the example source file uses the F90 MODULE protocol, this is merely for convenience. Otherwise, there would be the same collection of routines, all sharing a COMMON work area. In more general use one would extend the routines to provide more checking and some sort of optional trace facility. Indeed, it would be easy enough to prepare an interactive calculator using this scheme so that different calculations (and data) could be more easily experimented with.
===Fortran 90 ''et seq''.===
| |||