home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / INTERVA1.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  7.1 KB  |  239 lines
  1.  
  2. (* :Name: NumericalMath`IntervalArithmetic` *)
  3.  
  4. (* :Title: Approximate Interval Arithmetic *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Implements a simple form of interval arithmetic in which the
  10.     endpoints of the intervals are calculated correctly, but with
  11.     standard machine-precision arithmetic.  In particular, no effort
  12.     is made to round results outward.  It is not intended for
  13.     rigorous calculation, but rather for educational purposes.
  14. *)
  15.  
  16. (* :Context: NumericalMath`IntervalArithmetic` *) 
  17.  
  18. (* :Package Version: Mathematica 2.0 *)
  19.  
  20. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  21.     Permission is hereby granted to modify and/or make copies of
  22.     this file for any purpose other than direct profit, or as part
  23.     of a commercial product, provided this copyright notice is left
  24.     intact.  Sale, other than for the cost of media, is prohibited.
  25.  
  26.     Permission is hereby granted to reproduce part or all of
  27.     this file, provided that the source is acknowledged.
  28. *)
  29.  
  30. (* :History:
  31.     Revised by Jerry B. Keiper, December, 1990.
  32.     Originally by Jerry B. Keiper, May, 1990.
  33. *)
  34.  
  35. (* :Keywords: interval arithmetic *)
  36.  
  37. (* :Source:
  38.     R. E. Moore, Interval Analysis, Prentice-Hall, Englewood Cliffs,
  39.     NJ, 1966
  40.  
  41.     Germund Dahlquist and Ake Bjorck (Translated by Ted Anderson),
  42.     Numerical Methods, Prentice-Hall, Englewood Cliffs, NJ, 1974,
  43.     pp. 37-39
  44. *)
  45.  
  46. (* :Mathematica Version: 2.0 *)
  47.  
  48. (* :Limitation:
  49.     This package cannot be used to provide rigorous bounds on a
  50.     result, since the rounding is provided by the hardware and
  51.     occasionally may give smaller intervals than are justified.
  52. *)
  53.  
  54. (* :Discussion:
  55.     A simple form of interval arithmetic is implemented to provide
  56.     students with some (albeit limited) experience with the ideas
  57.     of interval arithmetic.  It is assumed that the intervals
  58.     overestimate the error sufficiently that incorrect rounding
  59.     in calculation of the endpoints of the intervals will not
  60.     significantly affect the result.  It is undoubtedly the case
  61.     that examples can be found where this asumption is violated.
  62.     Basic arithmetic on intervals is automatic."
  63.  
  64.     The basic data object used is Interval[a, x, b], where a and b
  65.     are the left and right endpoints of the interval, respectively,
  66.     and x is the value of the number that the interval represents.
  67.     Of course x is not correct either; it is simply the result of
  68.     machine arithmetic used in evaluating x.
  69.  
  70.     Interval[a, b] is converted to Interval[a, x, b] where x is
  71.     the midpoint of the interval (a, b).  Interval[x] is converted
  72.     to Interval[x - fuzz, x, x + fuzz] where fuzz is determined
  73.     by the magnitude of x and a constant, which can be set with
  74.     the function SetEpsilon[ ].
  75.  
  76.     SetEpsilon[dig] sets a constant that is used in converting
  77.     simple numbers to intervals.  When Interval[x]is converted
  78.     to Interval[x - fuzz, x, x + fuzz], the value of fuzz used is
  79.     Abs[x] * epsilon, where epsilon is 10^(-dig).  The value of
  80.     dig must be between 1 and 10, inclusive, and the default value
  81.     is 6."
  82. *)
  83.  
  84. (* :Examples:
  85.  
  86. In[1]:= << NumericalMath`IntervalArithmetic`    (* read in the package *)
  87.  
  88. Out[1]= NumericalMath`IntervalArithmetic`
  89.  
  90. In[2]:= a = Interval[-1,3]
  91.  
  92. Out[2]= Interval[-1, 1, 3]
  93.  
  94. In[3]:= b = Interval[1,5]
  95.  
  96. Out[3]= Interval[1, 3, 5]
  97.  
  98. In[4]:= a b
  99.  
  100. Out[4]= Interval[-5, 3, 15]
  101.  
  102. In[5]:= a a
  103.  
  104. Out[5]= Interval[-3, 1, 9]
  105.  
  106. In[6]:= a^2
  107.  
  108. Out[6]= Interval[0, 1, 9]
  109.  
  110. In[7]:= a-a
  111.  
  112. Out[7]= Interval[-4, 0, 4]
  113.  
  114. In[8]:= 3b
  115.  
  116. Out[8]= Interval[3, 9, 15]
  117.  
  118. In[9]:= Interval[Pi]        (* inexact arithmetic is used when necessary *)
  119.  
  120. Out[9]= Interval[3.14159, 3.14159, 3.1416]
  121.  
  122.     (* Note that the interval has positive length: *)
  123.  
  124. In[10]:= SetPrecision[SetPrecision[%,20],10]
  125.  
  126. Out[10]= Interval[3.141589512, 3.141592654, 3.141595795]
  127.  
  128.     (* We want a larger epsilon (i.e., more fuzz, fewer digits): *)
  129.  
  130. In[11]:= SetEpsilon[3]
  131.  
  132. In[12]:= sum = 0; Do[sum += 1/Interval[i], {i, 1., 500.}]; sum
  133.  
  134. Out[12]= Interval[6.78604, 6.79282, 6.79962]
  135.  
  136.     (* exact interval arithmetic is not affected by the order of a sum.
  137.        inexact interval is slightly affected by order, but not
  138.        significantly if the intervals are much bigger than
  139.        rounding errors. *)
  140.  
  141. In[13]:= sum = 0; Do[sum += 1/Interval[i], {i, 500., 1., -1}]; sum
  142.  
  143. Out[13]= Interval[6.78604, 6.79282, 6.79962] 
  144. *)
  145.  
  146. BeginPackage["NumericalMath`IntervalArithmetic`"];
  147.  
  148. SetEpsilon::usage =
  149. "SetEpsilon[dig] sets a constant that is used in converting simple numbers
  150. to intervals.  When Interval[x] gives Interval[x - fuzz, x, x + fuzz], the
  151. fuzz used is Abs[x] * epsilon, where epsilon is 10^(-dig).  The value of
  152. dig must be between 1 and 10, inclusive, and the default value is 6."
  153.  
  154. Interval::usage =
  155. "Interval[lowerbound, x, upperbound] is a data object that represents a
  156. number in interval arithmetic.  The second argument is the value that the
  157. number would have in ordinary machine arithmetic, and lowerbound and
  158. upperbound are the bounds of the interval.  Interval[x] gives
  159. Interval[x - fuzz, x, x + fuzz] where fuzz is Abs[x] * epsilon and epsilon
  160. is set using SetEpsilon[].  Interval[a,b] gives Interval[a,(a+b)/2,b].
  161. Basic arithmetic on intervals is automatic."
  162.  
  163. Begin["NumericalMath`IntervalArithmetic`Private`"]; 
  164.  
  165. SetEpsilon::digs = "The argument `1` to SetEpsilon is not between 1 and 10."
  166.  
  167. SetEpsilon[n_Integer] := 
  168.     (If[1 <= n <= 10, $epsilon = 10^(-n), Message[SetEpsilon::digs, n]];)
  169.  
  170. If[!NumberQ[$epsilon], SetEpsilon[6]];
  171.  
  172. Interval[x_] :=
  173.     Module[{nx = N[x]}, If[NumberQ[nx], Interval[nx], x]] /; !NumberQ[x]
  174.  
  175. Interval[x_] :=
  176.     Module[{fuzz}, fuzz = Abs[x] $epsilon; Interval[x - fuzz, x, x + fuzz]] /;
  177.     (NumberQ[x] && (Head[x] =!= Complex))
  178.  
  179. Interval[a_, b_] := 
  180.     Module[{na = a, nb = b},
  181.     Interval[na, (na+nb)/2, nb] /;
  182.         (If[!NumberQ[na], na = N[a]];
  183.          If[!NumberQ[nb], nb = N[b]];
  184.          (NumberQ[na] && NumberQ[nb] && na <= nb))];
  185.  
  186. Interval /:
  187. Plus[x_Interval, y_Interval] :=
  188.     Interval[x[[1]]+y[[1]], x[[2]]+y[[2]], x[[3]]+y[[3]]] /;
  189.         ((Length[x] == 3) && (Length[y] == 3))
  190.  
  191. Interval /:
  192. Plus[x_Interval, y_] :=
  193.     Interval[x[[1]] + y, x[[2]] + y, x[[3]] + y] /;
  194.         (NumberQ[y] && (Head[y] =!= Complex))
  195.  
  196. Interval /:
  197. Times[s_, Interval[a_,x_,b_]] :=
  198.     If[s >= 0,
  199.         Interval[a s, x s, b s],
  200.         (* else *)
  201.         Interval[b s, x s, a s]
  202.     ] /; (NumberQ[s] && (Head[s] =!= Complex))
  203.  
  204. Interval /:
  205. Times[x_Interval, y_Interval] :=
  206.     Module[{a1, a2, a3, a4},
  207.         a1 = x[[1]] y[[1]];
  208.         a2 = x[[1]] y[[3]];
  209.         a3 = x[[3]] y[[1]];
  210.         a4 = x[[3]] y[[3]];
  211.         Interval[Min[a1,a2,a3,a4], x[[2]] y[[2]], Max[a1,a2,a3,a4]]
  212.     ] /; (Length[x] == 3 && Length[y] == 3)
  213.  
  214. Interval /: Power[x_Interval, 0] := Interval[1, 1, 1]
  215.  
  216. Interval /:
  217. Power[x_Interval, n_Integer] :=
  218.     Module[{min, max, a1, a3},
  219.         If[x[[1]] <= 0 <= x[[3]],
  220.             If[n <= 0, Return[Indeterminate]];
  221.             a1 = x[[1]]^n;
  222.             a3 = x[[3]]^n;
  223.             min = Min[a1, 0, a3];
  224.             max = Max[a1, 0, a3],
  225.             (* else *)
  226.             a1 = x[[1]]^n;
  227.             a3 = x[[3]]^n;
  228.             min = Min[a1, a3];
  229.             max = Max[a1, a3]
  230.         ];
  231.         Interval[min, x[[2]]^n, max]
  232.     ] /; (Length[x] == 3)
  233.  
  234. End[ ] (* "NumericalMath`IntervalArithmetic`Private`" *)
  235.  
  236. Protect[SetEpsilon, Interval];
  237.  
  238. EndPackage[ ] (* "NumericalMath`IntervalArithmetic`" *)
  239.