The Mathematica Journal
Departments
Download This Issue
Home
Feature Articles
Graphics Gallery
Tricks of the Trade
In and Out
Columns
The Mathematica Programmer
New Products
New Publications
Classifieds
Calendar
News Bulletins
Editor's Pick
Mailbox
Letters
Write Us
About the Journal
Staff and Contributors
Submissions
Subscriptions
Advertising
Back Issues

Interval Plots

If we compute [Graphics:../Images/index_gr_83.gif] for an interval [Graphics:../Images/index_gr_84.gif], we know that [Graphics:../Images/index_gr_85.gif] for [Graphics:../Images/index_gr_86.gif]. Therefore, the points of the graph of [Graphics:../Images/index_gr_87.gif] (that is, [Graphics:../Images/index_gr_88.gif]) lie in the rectangle with lower left corner [Graphics:../Images/index_gr_89.gif] and upper right corner [Graphics:../Images/index_gr_90.gif]. This observation allows us to produce interval plots of functions. By subdividing the domain interval i into many small pieces, we can hope to close in on the values of [Graphics:../Images/index_gr_91.gif].

Here is a simple function with interesting behavior near [Graphics:../Images/index_gr_92.gif].

[Graphics:../Images/index_gr_93.gif]

This is our first domain interval.

[Graphics:../Images/index_gr_94.gif]

Here we see an ordinary plot of [Graphics:../Images/index_gr_95.gif] in this domain.

[Graphics:../Images/index_gr_96.gif]

[Graphics:../Images/index_gr_97.gif]

Applying [Graphics:../Images/index_gr_98.gif] to the interval gives us a rather crude idea of where the points of the graph of [Graphics:../Images/index_gr_99.gif] lie.

[Graphics:../Images/index_gr_100.gif]
[Graphics:../Images/index_gr_101.gif]

This definition produces a rectangle from a pair of intervals.

[Graphics:../Images/index_gr_102.gif]

Here we superimpose the rectangle [Graphics:../Images/index_gr_103.gif] on the graph of [Graphics:../Images/index_gr_104.gif].

[Graphics:../Images/index_gr_105.gif]

[Graphics:../Images/index_gr_106.gif]

We define a function split so that split[Interval[[Graphics:../Images/index_gr_107.gif]]] returns the two intervals {Interval[[Graphics:../Images/index_gr_108.gif]], Interval[[Graphics:../Images/index_gr_109.gif]]}.

[Graphics:../Images/index_gr_110.gif]

We split the domain into two intervals.

[Graphics:../Images/index_gr_111.gif]
[Graphics:../Images/index_gr_112.gif]

Now, we apply [Graphics:../Images/index_gr_113.gif] to these two intervals.

[Graphics:../Images/index_gr_114.gif]
[Graphics:../Images/index_gr_115.gif]

A transposition gives us the two rectangles.

[Graphics:../Images/index_gr_116.gif]
[Graphics:../Images/index_gr_117.gif]

Here we see that we have successfully narrowed the possible values of [Graphics:../Images/index_gr_118.gif] on the second interval.

[Graphics:../Images/index_gr_119.gif]

[Graphics:../Images/index_gr_120.gif]

Here we split the domain five more times to obtain [Graphics:../Images/index_gr_121.gif] intervals.

[Graphics:../Images/index_gr_122.gif]

We apply [Graphics:../Images/index_gr_123.gif] to the intervals by repeating the previous computation, which was set up to be independent of the number of intervals present.

[Graphics:../Images/index_gr_124.gif]

Here we see the result. In areas where [Graphics:../Images/index_gr_125.gif] varies only little, the rectangles are oblong and we could stop the subdivisions. In other areas, the rectangles are still rather tall and we should proceed with further subdivision.

[Graphics:../Images/index_gr_126.gif]

[Graphics:../Images/index_gr_127.gif]

The observations from our example suggest that we should subdivide the domain intervals [Graphics:../Images/index_gr_128.gif] until the length of [Graphics:../Images/index_gr_129.gif] is small enough (within the resolution of the image). Because the function being plotted could have singularities (such as 0 for [Graphics:../Images/index_gr_130.gif]) we should also stop if the length of [Graphics:../Images/index_gr_131.gif] gets small.

This is an adaptive technique, because not all intervals will be subdivided the same number of times. Such a method avoids unnecessary work and stops when at least one of several criteria is fulfilled.

The command IntervalPlot[expr,{x,[Graphics:../Images/index_gr_132.gif],[Graphics:../Images/index_gr_133.gif]}] implements these ideas. It is shown, together with most of the auxiliary definitions, in Listing 1. The auxiliary procedure refine[f,init,tol] contains the code for the adaptive subdivision, starting with an initial interval init. In the main loop, the nested list new contains the queue of intervals that have not yet been dealt with. We take the first of these intervals [Graphics:../Images/index_gr_134.gif] and compute [Graphics:../Images/index_gr_135.gif]. Note that we put the domain interval [Graphics:../Images/index_gr_136.gif] into a list. The reason for this will become clear in the next section when we use the same code for three-dimensional plots. If either [Graphics:../Images/index_gr_137.gif] or [Graphics:../Images/index_gr_138.gif] is small enough (smaller than tol, the value of the Tolerance option), we add the rectangle [Graphics:../Images/index_gr_139.gif] to the list of results; otherwise we split [Graphics:../Images/index_gr_140.gif] and put the resulting intervals back into the queue.

The queue is maintained as a nested list, [Graphics:../Images/index_gr_141.gif], instead of a simple list because append and drop operations are relatively slow compared to the manipulations of a nested list. The language LISP uses the same data structure. Table 1 compares the operations of adding an element to the queue and removing one for ordinary lists and nested lists. Note that the combined operation of obtaining the first element and removing it from the list is particularly simple for nested lists (because of Mathematica's parallel assignments of list elements). We use a similar nested structure for the list of final rectangles. Because the elements of this nested list are lists themselves, we use a head h, different from List, to hold the elements together to avoid problems when flattening the list. The nested list is converted into an ordinary list with the help of Flatten.

The function interval[arg] deals with cases where Mathematica's interval arithmetic does not return a proper interval, either because it returns a single number (which is turned into an interval) or for some other reason, in which case we return the interval [Graphics:../Images/index_gr_142.gif]. If an interval boundary is infinite, we replace it by a huge finite value and color the rectangle red to alert us to the singularity present.

[Graphics:../Images/index_gr_143.gif] [Graphics:../Images/index_gr_144.gif] [Graphics:../Images/index_gr_145.gif]
[Graphics:../Images/index_gr_146.gif] [Graphics:../Images/index_gr_147.gif] [Graphics:../Images/index_gr_148.gif]
[Graphics:../Images/index_gr_149.gif] [Graphics:../Images/index_gr_150.gif] [Graphics:../Images/index_gr_151.gif]
[Graphics:../Images/index_gr_152.gif] [Graphics:../Images/index_gr_153.gif] [Graphics:../Images/index_gr_154.gif]
[Graphics:../Images/index_gr_155.gif] [Graphics:../Images/index_gr_156.gif] [Graphics:../Images/index_gr_157.gif]
[Graphics:../Images/index_gr_158.gif] -- [Graphics:../Images/index_gr_159.gif]
Table 1. Recursive and iterative queue operations.


Options[IntervalPlot] = { Tolerance -> 0.01 }

length[i_] := Max[i] - Min[i]

split[i_Interval] :=
  With[{m = (Min[i] + Max[i])/2},
       {Interval[{Min[i], m}],Interval[{m, Max[i]}]} ]

SetAttributes[interval, Listable]

interval[i_Interval] := i
interval[i_?NumericQ] := Interval[i]
interval[_] := Interval[{-Infinity, Infinity}] (* junk *)

infcolor = RGBColor[1,0,0]; huge = 10.0^10;

rect[{xi__,yi_}] /; Min[yi] == -Infinity || Max[yi] == Infinity :=
 {infcolor, rect[{xi, yi /. DirectedInfinity[d:(-1|1)] :> d huge}]}

rect[l:{_,_}] := Rectangle[Min /@ l, Max /@ l]

RectangleGraph[l:{{_,_}..}] := Graphics[{Thickness[0], rect /@ l}]

IntervalPlot[expr_, {x_Symbol, x0_, x1_}, opts___?OptionQ]:=
  With[{tol = (x1-x0)Tolerance /. {opts} /. Options[IntervalPlot]},
    Module[{finals},
        finals = refine[Function[x, expr],
                        Interval /@ N[{{x0, x1}}], tol];
        Show[RectangleGraph[finals],
             Evaluate[FilterOptions[Graphics, opts]],
             PlotRange -> {{x0, x1}, Automatic}, Axes -> True]
  ] ]

refine[f_, init_, tol_] :=
  Module[{new = {init, {}}, finals = h[], i, j},
    While[new =!= {},
      {i, new} = new; (* dequeue *)
      j = interval[f @@ i];
      If[ length[j]<tol || Max[length /@ i] < tol,
         finals = h[Append[i,j], finals]
        ,(*else*)
         Scan[ (new = {#, new})&, split[i] ] (* queue new ones *)
      ]
    ];
    List @@ Flatten[finals, Infinity, h]
  ]
Listing 1. Interval plots.

Our code is in the Intervals.m package which we load after restarting the kernel.

[Graphics:../Images/index_gr_160.gif]

The difficult function [Graphics:../Images/index_gr_161.gif] gives a rather good result, including the singularity at zero.

[Graphics:../Images/index_gr_162.gif]

[Graphics:../Images/index_gr_163.gif]

For the function [Graphics:../Images/index_gr_164.gif], the simple pole is clearly marked in red.

[Graphics:../Images/index_gr_165.gif]

[Graphics:../Images/index_gr_166.gif]

The function[Graphics:../Images/index_gr_167.gif] is not treated in the best possible way at zero. The reason is that the numerator and denominator of the fraction are treated as independent values, as explained in the previous section.

[Graphics:../Images/index_gr_168.gif]

[Graphics:../Images/index_gr_169.gif]

A function with periodic features that are smaller than the tolerance shows the effects of aliasing. (Aliasing denotes features that occur in a sampled function when the sampling rate is too low for the true shape of the function to be reconstructed; it is observed frequently in digitized images and sounds.) In an ordinary plot, we might miss essential features; here, we are alerted to the problems present.

[Graphics:../Images/index_gr_170.gif]

[Graphics:../Images/index_gr_171.gif]

If the tolerance is smaller than the period of the oscillation, a satisfactory picture results.

[Graphics:../Images/index_gr_172.gif]

[Graphics:../Images/index_gr_173.gif]

This next example is often used to justify the use of interval plots. The small singularity is clearly detected.

[Graphics:../Images/index_gr_174.gif]

[Graphics:../Images/index_gr_175.gif]

An ordinary plot may miss the singularity and show only a small spike.

[Graphics:../Images/index_gr_176.gif]

[Graphics:../Images/index_gr_177.gif]


Converted by Mathematica      October 15, 1999 [Prev Page][Next Page]