Re: Interval Arithmetic project

For light-duty work, try something like:

infix("@+");
infix("@*");
infix("@-");
infix("@/");
intervalp(e) := not(mapatom(e)) and op(e) = 'interval;

makeinterval(x) := if intervalp(x) then x else interval(x,x);

"@+"(a,b) :=  (
   a : makeinterval(a),
   b : makeinterval(b),
   interval(first(a) + first(b), second(a) + second(b)));

intervalnegate(a) := (
  a : makeinterval(a),
  interval(-second(a), -first(a)));

"@-"(a,b) :=  a @+ intervalnegate(b);

"@*"(a,b) := block([r],
  a : makeinterval(a),
  b : makeinterval(b),
  r : [first(a) * first(b), first(a) * second(b), second(a) * first(b), 
second(a) * second(b)],
  interval(lmin(r),lmax(r)));

/* The interval(minf,inf) case could be a union of intervals.  But this 
code doesn't support that. */

intervalreciprocal(a) := block([r],
  a : makeinterval(a),
  if first(a) <= 0 and second(a) >= 0 then interval(minf,inf) else (
    r : [1/first(a), 1/second(a)],
    interval(lmin(r),lmax(r))));

"@/"(a,b) := a @* intervalreciprocal(b);

(%i1) load("E:/interval.mac")$
(%i2) x : interval(11/10,12/10)$

(%i3) (x @+ 10) @/ (x @+ 20);
(%o3) interval(111/212,112/211)

(%i4) r1 : interval(10-p,10+p)$
(%i5) r2 : interval(20-p,20+p)$
(%i6) assume(p > 0, p < 1)$
(%i7) 1 @/ ((1 @/ r1 ) @+ (1 @/ r2));
(%o7) 
interval(min(max(1/(1/(p+20)+1/(p+10)),1/(1/(20-p)+1/(10-p))),1/(1/(p+20)+1/(p+10)),1/(1/(20-p)+1/(10-p))),
max(min(1/(1/(p+20)+1/(p+10)),1/(1/(20-p)+1/(10-p))),1/(1/(p+20)+1/(p+10)),1/(1/(20-p)+1/(10-p))))

(%i8) subst(p=1/10,%);
(%o8) interval(19701/2980,20301/3020)

Notice the difference between @- and -:

(%i9) interval(-1,1) @- interval(-1,1);
(%o9) interval(-2,2)

(%i10) interval(-1,1) - interval(-1,1);
(%o10) 0


Barton