Initializations

Initial Code: Evaluate First

Graphics Settings

<< Miscellaneous`Units`;
Remove[Gray];
<< Graphics`Master`;
Off[General::spell1];
fs = 7;
fsbig = 9;
width  =  Convert[8.64 Centi Meter,  Point]/Point;
th  =  0.0008;

Function Definitions Relevant to Problem 4

Clear[f, fx, fy]; Off[General::spell1];

Needs["Graphics`Graphics`"];

f[x_, y_] := (E^Sin)[50 x] + 
            Sin[60 (E^y)] + Sin[70 Sin[x]] + Sin[Sin[80 y]] - 
            Sin[10 (x + y)] + (1/4) ((x^2) + (y^2));

g[x_, y_, z_]:=(E^Sin)[50 x]+Sin[60 (E^y)] Sin[60 z]+Sin[70 Sin[x]] Cos[10 z]+Sin[Sin[80 y]]-Sin[10 (x+z)]+(1/4) ((x^2)+(y^2)+(z^2));
gx[x_,  y_, z_]  :=  
      Evaluate[D[g[x, y, z],  x]]
gy[x_,  y_, 
        z_]  :=  Evaluate[
        D[g[x, y, z],  y]]
gz[x_,  y_, 
        z_]  :=  Evaluate[
        D[g[x, y, z],  z]]
fy[x_,  y_]  :=  
      Evaluate[D[f[x, y],  y]]
fx[x_,  
        y_]  :=  Evaluate[D[f[x, y],  x]]
fy[
        x_,  y_]  :=  Evaluate[D[f[x, y],  y]]

true  =  
        f[x, y] /. 
          FindRoot[{0 ==  fx[x, y]  ,   0 == 
                fy[x, y]},  {x, -0.01  },  {y, 0.2 },  
            AccuracyGoal ->  150,  
            WorkingPrecision -> 300];
truepos =  {x, y} /. 
          FindRoot[{0 ==  fx [x, y] ,   0 == 
                fy[x, y]},  {x, -0.01  },  {y, 0.2 },  
            AccuracyGoal ->  150,  
            WorkingPrecision -> 300];

Utilities

Off[General::"spell1"];

DigitsAgreeCount[a_, 
              b_] := (prec = 
                Ceiling@
                  Min[Precision /@ {a, b}]; {{ad, ae}, {bd, 
                    be}} = RealDigits[#, 10, prec] & /@ {a, b}; 
              If[ae ≠ be ∨ a b ≤ 0, Return[0]]; 
              If[ad == bd, Return@Length[ad]]; {{com}} = 
                Position[MapThread[Equal, {ad, bd}], False, 1, 1] - 1; 
              com);
DigitsAgreeCount[Interval[{a_, b_}]] := 
            DigitsAgreeCount[a, b];

IntervalForm[
              Interval[{a_, 
                  b_}]] := (If[(com = DigitsAgreeCount[a, b]) == 
                  0, Return@Interval[{a, b}]]; 
              start = Sign[a] N[
                    FromDigits@{ad[[
                          Range@com]], 1}, com]; {low, 
                  up} = SequenceForm @@ 
                      Take[#, {com + 1, prec}] & /@ {ad, bd}; 
              If[ae == 0, start /= 10; ae++]; 
              SequenceForm[
                DisplayForm@SubsuperscriptBox[NumberForm@start, low, up], 
                If[ae ≠ 1, 
                  Sequence @@ {" × ", 
                      DisplayForm@SuperscriptBox[10, ae - 1]}, ""]]);

diam[X_] := Max[X] - Min[X];
diam[{X__Interval}]  :=  
            Max[diam  /@  {X}];

gradient[f_,  
              v_]  :=  D[f,  #] &  /@  v;

ShowIntervals[a_List,  
              opts___] := (Show[
                  Graphics[
({PointSize[ .025], {{Hue[
                                  Random[],  .5,  .7], 
                                tt = Rectangle @@ (tt = 
                                        Transpose[{#1[[1, 
                                        1]], #1[[2, 1]]}])}, 
                              Point[ .5 Plus  @@  tt]}} &) /@ a], 
                  opts]) /; Depth[a] == 5;

containszerovector[{i__Interval}] := 
          And @@ (IntervalMemberQ[#, 0] &  /@  {i})
containszeronumber[z_Interval] := 
          IntervalMemberQ[z, 0]

  

subdivide[{ints__Interval},  wp_] := 
            With[{n = Length[{ints}]}, 
N[
                Flatten[
                  Table  @@ 
                    Prepend[
                      Table[{ind[j], 0, 1}, {j, n}], 
Table[
                        Interval[{ints}[[k, 1, 1]] + {ind[k], 
                                ind[k] + 
                                  1} ({ints}[[k, 1, 2]] - {ints}[[k, 
                                        1, 1]])/2],  {k,  n}]], n - 1], 
                wp]];

subdivide[z_ /; Depth[N[z]] == 5,  
              wp_] := Join  @@  (subdivide[#,  wp]  &  /@  
                  z);

IntervalContains[Interval[{a_,  b_}],  
              Interval[{c_,  d_}]]  :=  a < c < d < b;
IntervalContains[Interval[{a_,  b_}],  
              Interval[
                q__List]]  :=   And  @@ 
( IntervalContains[Interval[{a, b}],  Interval[#]] &  /@  {q}) /;  
              Length[{q}] > 1;

IntervalContains[{i1__Interval},  {i2__Interval}]  :=  
            And  @@  ((IntervalContains  @@  #) &  /@  
                  Transpose[{{i1},  {i2}}]);


combine[s_]  :=  
            Module[{newseeds = s}, 
              Do[
If[
                  FreeQ [
                    intervalintersectionbox [a = s[[i]], 
                      b = s[[j]]],  
                    Interval[]], 
newseeds[[j]] = newseeds[[i]] = 
                      intervalunionbox[a, 
                        b]], 
{i,   Length[s] - 
                    1},  {j,  i + 1,  
                  Length[s]}]; 
Union[newseeds]];

intervalintersectionbox[{i1__Interval},  {i2__Interval}] :=  
Apply[IntervalIntersection, 
              Transpose[{{i1},  {i2}}], {1}];
intervalunionbox[ {i1__Interval},  {i2__Interval}]  :=  
Apply[IntervalUnion, Transpose[{{i1},  {i2}}], {1}];

rectToIteratorList[{i__},  
              vv_] :=   Flatten  /@  Transpose[{vv,  First /@ {i}}];

IteratorListToBox[{i__}] :=   Interval  /@  Rest /@ {i};

intervalmemberQbox[{i__Interval},  
              r_] :=  
And @@ 
              Apply[IntervalMemberQ, Transpose[{{i},  r}], {1}];

Attributes[midpoint] = Listable;

midpoint[X_] := (Min[X] + Max[X])/2;
midpoint[x_?NumericQ]  :=  x;

Attributes[midpoint] = Listable;
Attributes[thirdpoint] = Listable;

thirdpoint[X_] := (2 Min[X] + Max[X])/
              3;
thirdpoint[x_?NumericQ]  :=  
            x;

IntervalToReal[
              Interval[{a_, b_}]] := ((a + b)/2);
IntervalToReal[Interval[]] :=  {};
IntervalToReal[l_List]  :=  IntervalToReal  /@  l;
IntervalToReal[a_] := a;
RealToInterval[r_?NumericQ, d_] := 
            Interval[r + {-1, 1} d];
Attributes[RealToInterval]  =  Listable;
RealToInterval[v_?VectorQ, 
              d_] := RealToInterval[#,  d]  &  /@  v;
MidRad[x_, r_] := 
            x + Interval[{-1, 1} ] r;

Main Functions: GeneticMinimize, IntervalKrawczyk, IntervalMinimize, FindRootsIntervalKrawczyk, ValidateRoots, FindAndValidateRoots

Algorithm 4.1: GeneticMinimize

\(\(Needs["Utilities`FilterOptions`"];\)

 (*\ 
      subtle\ \(point : \ 
          in\ this\ code\ I\ allow\ the\ parents\ to\ go\ outside\ the\ box . \ When\ that\ is\ forcibly\ disallowed\), \ the\ results\ are\ much\ worse; \ 
      there\ is\ an\ option\ to\ set\ \(\(this\)\(.\)\)\ *) 
\)
GeneticMinimize::usage  =  "GeneticMinimize[f, {x,a,b}, \
{y,c,d}, ...] uses a genetic algorithm to try to find the minimum of f on the \
n-dimensional box defined by the iterators.";

    
GenerationSize::usage  =  "GenerationSize is an option to \
genetic minimize that sets the number of points kept to form each generation.\
";

    
StoreGenerations::usage  =  "StoreGenerations is an option to \
GeneticMinimize that causes each genreation to be stored in a global variable \
called store.";

    
ScaleFactor::usage  =  "ScaleFactor is an option to \
GeneticMinimize that sets the factor by which the region containing the \
children is shrunk at each generation.";

    
\(\(\(GenerationSize::usage\  = \ "GenerationSize is an option to GeneticMinimize that sets the number of points in each generation.";\)\(
\)
    \)\)
StayInDomain::usage  =  "StayInDomain is an option to \
GeneticMinimize that forces children to stay in the given domain. The default \
is False.";

    
SwitchToNewtonTolerance::usage  =  "SwitchToNewtonTolerance is \
an option to GeneticMinimize that sets a tolerance on the size of the current \
domain in which the children will live below which Newton's method is used to \
find the nearest critical point.";

    
TraceQ::usage  =  "TraceQ is an option to several functions \
that causes various messages showing the progress of the computation to be \
printed.";

    
Options[GeneticMinimize] = {AccuracyGoal -> 6, 
          StoreGenerations -> False,  MaxIterations ->  100,  
          ScaleFactor -> 1/2, GenerationSize -> 30, 
          TraceQ -> False, 
          StayInDomain -> 
            False, 
SwitchToNewtonTolerance  ->  -∞, 
WorkingPrecision -> 
            MachinePrecision,  
          Compiled  -> True};

    
\(GeneticMinimize[f_, iters__List, opts___Rule] := 
      Module[{func, switch = False, 
          h}, 
{sgQ, ag, \ ml, sf, wp, n, trQ, stayQ, 
            newtol, \ 
            comQ} = \({StoreGenerations, AccuracyGoal, \ MaxIterations, \ 
                ScaleFactor, \ WorkingPrecision, \ GenerationSize, \ 
                TraceQ, \ StayInDomain, \ SwitchToNewtonTolerance, \ 
                Compiled} /. {opts}\) /. Options[GeneticMinimize]; 
tol = 
          10\^\(-ag\); 
start = 
          N[Apply[Plus, \ Rest\  /@ \ {iters}, \ {1}]/2, 
            wp]; 
h = 
          N[\(-Apply[Subtract, \ Rest\  /@ \ {iters}, \ {1}]\)/2, 
            wp]; 
vbles = 
          First /@ {iters}; 
func[wp_]\  := \(func[wp] = 
            If[comQ && \ wp <= $MachinePrecision, \ 
              Compile[Evaluate[vbles], f], 
              Function[Evaluate[vbles], 
                Evaluate[N[f, wp]]]]\); 
parents = {start}; 
        lev = 0; 
fvals = {func[wp] @@ 
              start}; 
store = {{fvals, 
              start}}; 
While[
          Min[h] > tol\  && \ 
            lev\  < \ ml, \(lev++\); 
new = 
            If[stayQ, \ 
Flatten[\((Table[# + 
                          SetPrecision[
                            Table[Random[
                                Real, 
{Min[
                                    h[[
                                      kk]], \ {iters}[\([kk, 2]\)] - #[\([kk]\)]], 
Min[
                                    h[[
                                      kk]], \ {iters}[\([kk, 3]\)] - #[\([kk]\)]]}], {kk, Length[h]}], wp], {n}] &)\) /@ 
                  parents, 
                1], 
Flatten[\((Table[# + 
                          SetPrecision[
                            Table[Random[
                                Real, 
{\(-1\), 1} 
                                  h[[kk]]], {kk, Length[h]}], wp], {n}] &)\) /@ parents, 
                1]]; 
{fvals, parents} = 
            Transpose[
              Take[Sort[
                  Transpose[{Join[fvals, Apply[func[wp], new, {1}]], 
                      Join[parents, new]}]], n]]; 
          h *= sf; 
If[sgQ, 
            AppendTo[
              store, {fvals[[1]], 
                parents}]]; 
If[trQ, \ 
            Print[StringForm["Iteration ``; max spread = ``; best = ``", 
                lev, NumberForm[N[Max[h]], 4], 
                NumberForm[\ fvals[\([1]\)], 12]]]]; 
If[
            Max[h]\  <= \ newtol, switch = True; \ 
            Break[]]]; 

If[switch, \ 
          Print["Switching to Newton's Method on the gradient."]; 
grad = \(D[f, #] &\)\  /@ \ 
              vbles; 
fr = \((FindRoot\  @@ \ 
                Join[Prepend[
                    Transpose[{vbles, \ 
                        parents[[1]]}], 
                    grad\ ], {FilterOptions[FindRoot, 
                      opts]}])\); 
{func[
                wp]\ \  @@ \ \((Last\  /@ \ fr)\), \ 
            fr}, 
{fvals[[1]], 
            Thread[vbles -> 
                parents[[1]]]}]]\)

The Krawczyk Operator

Code for the interval Newton operator is included too; Krawczyk is faster and is used throughout the routines later.

\(\(\(\ \)\(\(Krawczyk::usage\  = \ "Krawczyk[f, {x, a, b},..., jacobian] returns a triple: the first indicates whether the Krawczyk containment holds, the second is True if the Krawczyk image is disjoint from the given box, the third is the Krawczyk operator on the given box.";\)\ 


    \(KrawczykQ::usage\  = \ "KrawczykQ[f, {x, a, b},..., jacobian] returns True if the Krawczyk condition holds on the given n-dimensional interval.";\)\ 


    \(NewtonBox::usage\  = \ "NewtonBox[f, {x,a,b},..., jacobianinv] returns a triple: the first indicates whether the Newton containment holds, the second is True if the Newton image is disjoint from the given box, the third is the Newton operator on the given box.";\)\ 


    \(Krawczyk[f_, \ i__List, jacobian_]\  := \ 
        Block[{Experimental`$EqualTolerance\  = \ 0}, 
          Module[{m, Y, int = IteratorListToBox[{i}], 
              vb\  = \ First\  /@ \ {i}}, 
m\  = \ 
              midpoint[int]; \ 
midjac = 
              midpoint[
                jj = jacobian\  /. \ 
                    Thread[vb -> \ int]]; 
If[
              Abs[Det[midjac]] < 10^\(-6\), \ 
              midjac\  = 
                thirdpoint[
                  jj = jacobian\  /. \ 
                      Thread[vb -> \ int]]]; \ 
Y = 
              Inverse[midjac]; 
kk = 
              m - Y\ \  . \ \((f\  /. \ 
                      Thread[vb\  -> \ m]\ )\) + \ \((IdentityMatrix[
                        Length[f]] - \ Y . \ jj)\) . \((int - 
                      m)\)\ ; 
{IntervalContains[int, \ 
                kk], \(! FreeQ[intervalintersectionbox[int, kk], \ 
                  Interval[]]\), \ 
              kk}]];\)


    \(KrawczykQ[f_, \ i__List, jacobian_]\  := \ 
        Block[{Experimental`$EqualTolerance\  = \ 0}, 
          Module[{m, Y, int = IteratorListToBox[{i}], 
              vb\  = \ First\  /@ \ {i}}, 
m\  = \ 
              midpoint[int]; \ 
midjac = 
              midpoint[
                jacobian\  /. \ 
                  Thread[vb -> \ m]]; 
If[
              Abs[Det[midjac]] < 10^\(-6\), \ 
              midjac\  = 
                thirdpoint[
                  jj = jacobian\  /. \ 
                      Thread[vb -> \ int]]]; 
Y = 
              Inverse[midjac]; 
IntervalContains[int, \ 
              m - Y\ \  . \ \((f\  /. \ 
                      Thread[vb\  -> \ m]\ )\) + \ \((IdentityMatrix[
                        Length[f]] - \ 
                      Y . \ \((jacobian\  /. \ 
                            Thread[vb\  -> \ int])\))\) . \((int - 
                      m)\)]]];\)


    \(NewtonBox[f_, \ i__List, jacobianinv_]\  := \ 
        Block[{Experimental`$EqualTolerance\  = \ 0}, 
          Module[{m, int = IteratorListToBox[{i}], 
              vb\  = \ First\  /@ \ {i}}, 
m\  = \ 
              midpoint[int]; \ 
kk = 
              m - \((jacobianinv\  /. \ 
                      Thread[vb\  -> \ int])\) . \((f\  /. \ 
                      Thread[vb\  -> \ 
                          m]\ )\); 
{IntervalContains[
                int, \ kk], \(! FreeQ[intervalintersectionbox[int, kk], \ 
                  Interval[]]\), \ 
              kk}]];\)


    \(NewtonBoxQ[f_, \ i__List, jacobianinv_]\  := \ 
        Block[{Experimental`$EqualTolerance\  = \ 0}, 
          Module[{m, int = IteratorListToBox[{i}], 
              vb\  = \ First\  /@ \ {i}}, 
m\  = \ 
              midpoint[int]; \ 
IntervalContains[int, \ 
              m - \((jacobianinv\  /. \ 
                      Thread[vb\  -> \ int])\) . \((f\  /. \ 
                      Thread[vb\  -> \ m]\ )\)]]];\)\)\)\)

IntervalKrawczyk, a needed subroutine for the rootfinding part of Algorithm 4.2

\(\(\(IntervalKrawczyk::usage\  = \ "IntervalKrawczyk[f, iters (,opts)] iterates the Krawczyk operator on the given box."\)\(
\)
    \)\)
Options[
          IntervalKrawczyk] = {WorkingPrecision -> MachinePrecision, 
          AccuracyGoal -> 6, MaxIterations -> 40, 
          CheckUniqueness -> False, Jacobian  ->  Automatic,  
          EvaluationCount -> 0};

    
\(\(IntervalKrawczyk[f_, iters__List, opts___Rule] := 
        Block[{Experimental`$EqualTolerance\  = \ 0}, 
          Module[{t, count = 0, ag, wp, maxiter, cuQ, tol, det, der, invder, 
              invderinterval, temp, Xans, pp, fail, start, uvbles, vbles, 
              kk}, {ag, wp, maxiter, cuQ, \ der, \ 
                ec} = \({AccuracyGoal, WorkingPrecision, MaxIterations, 
                    CheckUniqueness, \ Jacobian, \ 
                    EvaluationCount} /. {opts}\) /. Options[IntervalKrawczyk]; 
fcount = 
              ec; 
fff = 
              If[\(! ListQ[f]\), {f}, f]; 
Xans = 
              SetPrecision[IteratorListToBox[{iters}], 
                wp]; 
uvbles = 
              First /@ {iters}; 
dim = Length[vbles]; 
            tol = 10\^\(-ag\); 
fail = False; 
            count = 0; 
deriv\ \  = \ 
              If[\((der\  === \ Automatic\ )\) || \ 
                  cuQ\ , \ 
\ \ \ \ Transpose[\((∂\_# fff &)\) /@ uvbles], \ 
                der\ ]; 
If[cuQ, det = Det[deriv]]; 
            While[\(count++\); 
count < maxiter && 
                FreeQ[Xans, 
                  Interval[]] && \((\((Max[
                          fff\  /. \ Thread[uvbles -> Xans]]\  > \ 
                        tol\ )\) || \ \ \ \ \ \((Or @@ \(\((\(-Subtract @@ #1[[1]]\) > tol &)\) /@ Xans\))\) || 
                    If[cuQ, 
                      containszeronumber[
                        det /. Thread[ubles -> Xans]], 
                      False])\), 
kk = 
                Krawczyk[fff, \ 
                  Sequence @@ rectToIteratorList[Xans, \ uvbles], \ 
                  deriv]; 
fcount += \ 
                dim^2; 
Xans = 
                intervalintersectionbox[Xans, 
                  kk[\([3]\)]]]; 
If[fail, $Failed, 
              Xans]]];\)\)

Algorithm 4.2: IntervalMinimize

\(OpportunisticEvaluations::usage\  = \ "OpportunisticEvaluations is an option to IntervalMinimize that checks f-values on the center of each box and lowers the bound accordingly. When this technique is used, the result will be an interval containing the absolute minimum on the domain, instead of the lowest critical point, since these evaluations can find values below the critical points.";\)

ShowRectangles::usage  =  "ShowRectangles is an option to \
some interval algorithm functions that works only in dimension 2 and causes \
the rectangles at each stage of the algorithm to be shown, and stored as \
im[level]. The storage makes for easy generation of art as in Figure 4.5 in \
the book.";

\(UseIntervalKrawczykWhenOneBoxLeft::usage\  = \ "UseIntervalKrawczykWhenOneBoxLeft is an option to IntervalMinimize that causes an interval root-finding method on the gradient to be used when only one interval remains. This speeds convergence.";\)

TraceQ::usage  =  "TraceQ is an option to several functions \
that causes various messages showing the progress of the computation to be \
printed.";

IntervalMinimize::nocp =  "There is no critical point in the \
domain. Try setting CheckGradient to False to turn off this search.";

IntervalMinimize::maxlev =  "The MaxIterations bound has \
been exceeded, and so the accuracy goal is likely not met.";

MergeTolerance::usage  =  "MergeTolerance is an option to \
certain interval-algorithm functions that merges neighboring intervals when \
their diameter is less than the tolerance.";

FormatIntervalsNicely::usage  =  "FormatIntervalsNicely is \
an option to certain interval-algorithm functions that formats the intervals \
in the output in a sub-super-script form.";

CheckGradient::usage  =  "CheckGradient is an option to \
IntervalMinimize that causes intervals that cannot contain a critical point \
to be discarded.";

Options[IntervalMinimize]  =  {MaxIterations ->  30,  
            WorkingPrecision  ->  MachinePrecision,  
            AccuracyGoal  ->  6,  TraceQ  ->  False,  
            MergeTolerance  ->  0,  
            OpportunisticEvaluations  -> False,  
 
            CheckGradient -> True,  ShowRectangles  -> False,  
            UseIntervalKrawczykWhenOneBoxLeft -> 
              True,    FontSize -> 9,  
            FormatIntervalsNicely -> True};

IntervalMinimize[f_, iters__List, 
          givenmin_?((# == ∞  || NumericQ[#]) &), 
          opts___]:=Block[{Experimental`$EqualTolerance  =  0},Module[{uvbles,  vbles, dim , kk, wp},
{wp,  maxlev, ag, trQ,  mtol,  oppQ, inQ,  
                      cgQ,  artQ, fs,  
                      formQ}  = {WorkingPrecision, 
                          MaxIterations,   AccuracyGoal, TraceQ,  
                          MergeTolerance,  
                          OpportunisticEvaluations,  
UseIntervalKrawczykWhenOneBoxLeft, CheckGradient,  
                          ShowRectangles, 
FontSize,  
                          FormatIntervalsNicely}  /.  {opts}  /.  
                      Options[IntervalMinimize];


hf  =  
                    HoldForm[IntervalMinimize[f, iters, givenmin, opts]];
rects = 
                    N[ {Interval  /@  Rest  /@  {iters}}, wp];fcount = 0;
ends  =  
                    rects[[1]]  /.  
                      Interval[{i_,  
                            j_}] :> {Interval[{-∞,  i}],  Interval[{j, ∞}]};
borderQ[{ii__Interval}]  :=  
                    Or  @@  Table[
                        Max[{ii}[[j]]]  >=  
                            Min[ ends[[j, 
                                  2]]]  ||  
Min[{ii}[[j]]]  <=  Max[ ends[[j, 1]]],  {j,  
                          Length[ends]}] ;

exterior  =  rects[[1]];
smallstage  =  False;
vbles  =  
                    Unique  /@  First  /@  {iters};
uvbles  =  
                    First  /@  {iters};
dim = Length[uvbles];
If[dim  ≤  2,  
                    ccQ = False];
Clear[fcn, gradf];
fcn[m_:  ∞] := 
                    Function[
                      zz, (N[f /.  Thread[uvbles ->  zz], m])];
gradf[Pattern[#, 
                            Blank[]] &  /@  vbles]  := 
                    Evaluate[(
  (D[f, #] & /@  
                              uvbles) /. 
                          Thread[uvbles ->  vbles])];

If[inQ,  
                    hessian = 
                      Transpose[(D[gradf [uvbles], #] & /@  
                            uvbles)]; 
invder = 
                      Inverse[hessian]];

ansmax = N[givenmin, wp]; ansmin  =  -∞; lev = 0;
rectcount = {};
\(While[
                    rects\  ≠ \ {}\  && 
                      ansmax - ansmin > 10^\(-ag\)\  && \ 
                      lev\  ≤ \ 
                        maxlev, \(lev++\); 

If[dim == 2\  && \ artQ, 
                      im[lev] = 
                        ShowIntervals[rects, FilterOptions[Graphics, opts], \ 
                          AspectRatio -> Automatic, \ 
                          Frame -> True, 
                          PlotRange\  -> \ {Rest[{iters}[\([1]\)]], \ 
                              Rest[{iters}[\([2]\)]]}, 
TextStyle -> \(FontSize -> fs\), \ 
                          PlotLabel -> 
                            StringForm["Bound = ``", \ 
                              NumberForm[ansmax, 
                                5]]]]; 
rects = 
                      Join\  @@ \ 
                        Map[subdivide[#, wp] &, 
                          rects]; 

fvals\  = Map[fcn[wp], rects]; \ 
                    fcount\  += \ Length[rects]; 
If[trQ, \ 
                      rectcount\  = \ {rectcount, \ 
                          Length[rects]}]; 
ansmax\  = \ 
                      If[oppQ, \ 
 (*\ 
                          if\ requested\ do\ a\ filtering\ evaluation\ to\ try\ to\ improve\ ansmax\ *) 
ansmax\  = \ 
                          Min[ansmax, 
\ 
                            Min[Map[\(\(fcn[
                                        wp]\)[#]\)[\([1, \(-1\)]\)]\  &, \ 
                                rects\  /. \ 
                                  Interval[{a_, b_}] :> \ 
                                    Interval[{\((a + b)\)/2, \((a + b)\)/
                                        2}]]]], 
Min[
                          ansmax, 
                          Min[fvals /. \ 
                              Interval[{_, \ m_}] :> \ 
                                m]]]; 

 (*\ Delete\ based\ on\ size\ *) 
pos = 
                      Flatten[
                        Position[fvals, 
                          Interval[{m_, _}] /; \ 
                            m\  ≤ \ 
                              ansmax]]; 
rects = \ 
                      rects[[pos]]; \ 
                    fvals = \ 
                      fvals[[pos]]; 

If[
                      cgQ, \  (*\ 
                        delete\ based\ on\ whether\ gradient\ has\ a\ zero; \ 
                        keep\ border\ boxes\ *) 
pos = 
                        Flatten\ [
                          Position[
rects, {ii__Interval}\  /; \ \((borderQ[{ii}]\  || \ \ \((fcount += \ \ dim; 
                                    containszerovector[
                                      gradf[{ii}]])\))\)]]; 
rects = rects[[pos]]; 
                      fvals = 
                        fvals[[
                          pos]]]; 

\  (*\ improve\ ansmin\ *) \ \  (*\ 
                      pos\ can\ NEVER\ be\ \(empty!\)\ *) 
If[pos ≠ \ {}, 
                      ansmin\  = \ 
                        Min[fvals]]; 

 (*\ if\ there\ is\ one\ box\ left, \ 
                      try\ to\ use\ rootfinding\ *) 
If[
                      inQ\  && \ Length[rects] == 1\  && \ 
                        lev\  > \ 1, 
test = \ 
                        Krawczyk[
                          gradf[uvbles], 
Evaluate[
                            Sequence\  @@ \ \ rectToIteratorList[
                                rects[\([1]\)], \ uvbles]], 
                          hessian]; 
If[
                        test[\([1]\)], 
ans = 
IntervalKrawczyk[
                            gradf[uvbles], \ 
Evaluate[
                              Sequence @@ 
                                rectToIteratorList[test[\([3]\)], \ uvbles]], 
                            Jacobian -> 
                              hessian, 
WorkingPrecision -> wp, \ AccuracyGoal\  -> \ ag, \ 
                            EvaluationCount -> 
                              fcount]; 
Return[{\(fcn[wp]\)[
                                ans], \(Thread[
                                    uvbles\  -> \ #] &\)\  /@ \ \ {\ 
                                  ans}}\  /. \ 
                            z_Interval\  :> \ \(If[formQ, 
                                  IntervalForm, \ Identity]\)[
                                z]]]]; 

If[trQ, \ 
                      Print[StringForm["Finished round ``; `` < answer < ``; `` rectangles to be subdivided", lev, NumberForm[ansmin, ag], \ 
                          NumberForm[\ ansmax, ag], \ \ Length[rects]]]]]\);

If[
                    NumericQ[mtol]  &&  Length[rects]  >  1  &&  
                      Max[maxdim /@  rects]  <  mtol,  
                    rects  =  FixedPoint[combine,  rects]];
If[trQ,  
                    Print[
StringForm["Finished; no. of \
current boxes = ``, total rects = ``, ``; total function evals = ``", 
                        Length[rects] + Length[smallrects], 
                        Flatten[ rectcount],  
                        Plus  @@  Flatten[rectcount],  fcount]]];
If[lev > maxlev,  
                    Message[IntervalMinimize::maxlev]];
{Interval[{ansmin,  
                          ansmax}] ,  Thread[
                            uvbles  ->  #] &  /@  rects} /.  
                    z_Interval  :>  If[formQ, IntervalForm,  
                          Identity][z]]];

Algorithm 4.3: FindRootsContourMethod

Needs["Utilities`FilterOptions`"];  Off[General::spell1];  
    Off[General::spell];

\(\(\(FindRootsContourMethod::usage\  = \ "FindRootsFindRootsContourMethod[{f, g}, {x, a, b}, {y, c, d}] finds all the  roots of f = g = 0 in the rectangle determined by a ≤ x ≤ b and c ≤ y ≤ d. It also works with arguments (f, {z, z1, z2}) where f is a complex function and z1 and z2 define a rectangle in the complex numbers.";\)\(
\)
    \)\)
\(\(\(NewtonsMethod::usage\  = \ "NewtonMethod is an option to FindRootsContourMethod that sets whether Newtons method or a secant method will be used to refine the seeds into accurate roots.";\)\(
\)
    \)\)
\(\(\(SeedsOnly::usage\  = \ "SeedsOnly is an option to FindRootsContourMethod that, when true, causes only the seeds to be returned.";\)\(
\)
    \)\)
\(\(\(CountRoots::usage\  = \ "CountRoots is an option to FindRootsContourMethod that, when true, prints a message counting the number of seeds and the number of high-precision roots. If these numbers do not agree, then the seeds have not led to distinct roots, which indicates that a higher Resolution setting should be used.";\)\(
\)
    \)\)
\(\(\(ContourResolution::usage\  = \ "ContourResolution is an option to FindRootsFindRootsContourMethod that controls the resolution of the contour plo.";\)\(
\)
    \)\)
\(\ \(DuplicateTolerance::usage\  = \ "DuplicateTolerance is an option to FindRootsContourMethod that sets the tolerance used when considering two roots the same. Setting it to 0 causes no elimination of duplicates.";\)
\)
\(\(\(FindRootsContourMethod::notfine\  = \ "Newton's method does not contract. Refining the subdivision.";\)\(
\)
    \)\)
EliminateDuplicates[l_?VectorQ,  tol_]  :=  
        Union[l,  
          SameTest  ->  (Abs[#1 - #2] < tol  &)];
\(\(\(EliminateDuplicates[l_, \ tol_]\  := \ 
        Union[l, \ 
          SameTest\  -> \ \((normsq[#1 - #2] < tol^2\ \  &)\)];\)\(
\)
    \)\)
Options[FindRootsContourMethod] = {ContourResolution -> 59,  
          NewtonsMethod -> True,  
          SeedsOnly -> False,  
CountRoots -> False, 
          WorkingPrecision -> $MachinePrecision, 
          AdaptiveQ  -> 
            False, 
DuplicateTolerance -> 10^-5,  
          Trace -> False}; 

normsq[a_]  :=  a . a;

    
\(\(FindRootsContourMethod[\ {f1_, f2_}, {x_, x1_, x2_}, {y_, y1_, y2_}, 
          opts___] := 
Module[{ (*seeds, scount, met, tt, 
              cd, *) e1 = \((x1 - x2)\)/
                50, 
e2 = \((y1 - y2)\)/
                50, 
fy = Compile[{x, y}, Evaluate[f2]], 
            ans, tol}, 
{met, \ seedQ, \ tol, \ ctQ, \ 
              pp}\  = \({NewtonsMethod, \ SeedsOnly, \ DuplicateTolerance, \ 
                  CountRoots, \ ContourResolution}\  /. \ {opts}\)\  /. \ 
              Options[FindRootsContourMethod]; 
cd = 
            Map[First, 
Cases[
                Graphics[
                  ContourPlot[
\ 
                    f1, {x, x1, x2}, {y, y1, 
                      y2}, 
\ \ \ Evaluate[
                      FilterOptions[ContourPlot, 
                        opts]], 
\ \ \ Contours -> {0}, ContourShading -> False, 
		\ PlotPoints -> 
                      pp, \ \ DisplayFunction -> 
                      Identity]], _Line, ∞]]; 
seeds = 
            Flatten[Map[#[[1 + 
                      Flatten[
		Position[
                          Rest[tt = Sign[Apply[fy, #, 2]]]*
		Rest[
                              RotateRight[
                                tt]], \(-1\)]]]] &, cd], 
              1]; 
scount = 
            Length[seeds]; 
If[ctQ, \ 
            Global`seedcount = scount]; 

\ 
          If[seedQ\  || \ seeds == {}, \ 
            Return[seeds]]; 
ans = 
            Select[Union[
                mmm = Map[
{x, y} /. 
                        FindRoot[{f1 == 0, 
                            f2 == 0}, 
Evaluate@
                            If[met, {x, #[[1]]}, 
                              Prepend[\ #[[1]] + {\(-1\), 1} e1, x]], 
Evaluate@
                            If[met, {y, #[[2]]}, 
                              Prepend[\ #[[2]] + {\(-1\), 1} e2, y]], 
	
                          Evaluate[FilterOptions[FindRoot, opts]]] &, 
                    seeds]], 
\ \ \ \ x1 < #[[1]] < x2 ∧ 
                  y1 < #[[2]] < 
                    y2\  &]; 

If[\ 
            tol > 0, \ \ \ ans = 
              EliminateDuplicates[ans, \ tol]]; 
If[ctQ, \ 
            Print[StringForm["`` seeds were found, yielding `` roots.", \ 
                scount, \ 
                Length[ans]]]]; 
ans];\)\)

Algorithm 4.4: FindRootsIntervalKrawczyk

\(FindRootsIntervalKrawczyk::usage\  = \ "FindRootsIntervalKrawczyk[f, {x, a, b}, {y, c, d}, ..., opts] returns a pair. The first contains trapping intervals for roots, with one in each. The second contains intervals that might contain roots.";\)

FindRootsIntervalKrawczyk::badinp  =  "The number of \
functions does not match the number of unknowns.";

Options[
            FindRootsIntervalKrawczyk]  =  {WorkingPrecision ->  
              MachinePrecision,  TraceQ -> False,  
            AccuracyGoal -> 6,  
            MaxIterations -> 40,   MergeTolerance  ->  0, 
            ShowIntervals -> False,  
            FormatIntervalsNicely -> True};

FindRootsIntervalKrawczyk[f_, 
            iters__List ,   opts___Rule] := Block[{Experimental`$EqualTolerance  =  0},Module[{wp, fff, reg, grad, int,  vals,  pos, len2,  
                    maxiter, invder, det, count, trQ, t},
{trQ,  wp,  maxiter,  artQ,  mtol,  
                        formQ}  = { TraceQ,  WorkingPrecision,  
                            MaxIterations,  ShowIntervals,  
                            MergeTolerance,  
                            FormatIntervalsNicely}  /.  {opts}  /.  
                        Options[FindRootsIntervalKrawczyk];
  
hf  =  
                      HoldForm[
                        FindRootsIntervalKrawczyk[f, iters ,   opts]];
reg =  {Interval  /@  Rest  /@  {iters}};
uvbles  =  
                      First  /@  {iters};
vbles  = 
                      Unique  /@   uvbles;
dim = Length[uvbles];
fff = 
                      If[! ListQ[f],  {f},  f];
If[
                      Length[fff]  ≠  dim,  
                      Message[FindRootsIntervalKrawczyk::badinp];  
                      Return[hf]];
h[m_:  ∞] := 
                      Function[Evaluate[uvbles], Evaluate[N[fff, m]]];
der = 
                      Transpose[D[fff,  #]  &  /@  uvbles];

det  =  Det[der];
count = 0;   answer = {};

While[reg  !=  {}  &&  count  <  maxiter,count++;
reg = 
                            subdivide[reg,  wp];
len1 = Length[reg];
If[
                            dim == 2  &&  artQ, 
                            Show[Graphics[
                                reg  /. {Interval[{aa_,  bb_}], 
                                      Interval[{cc_,  
                                        dd_}]} :> {Hue[
                                        Random[]],  
                                      Rectangle[{aa, cc}, {bb  , 
                                        dd}]}], 
PlotRange -> Rest /@  Take[{iters},  2],  AspectRatio -> Automatic,  
                              Frame -> True]];
  
pos = 
                            Flatten [ 
                              Position[
                                Apply[h[wp], reg, 1],  _?
                                  containszerovector]] ;
reg = 
                            reg[[pos ]];    
len2 = Length[reg];

  
If[
                            NumericQ[mtol]  &&  Length[reg]  >  1  &&  
                              Max[diam /@  reg]  <  mtol,  
                            reg  =  FixedPoint[combine,  reg]];

If[ trQ,  
                            Print[StringForm[
"Round ``. \
Out of `` regions searched, `` might contain zeros and need to be searched. \
`` answers so far", count, len1, len2,  Length[answer]]]];
  
pos = {};posbad = {};
newroots =   MapIndexed[(
kk  =  
                                    Krawczyk[fff, 
                                      Evaluate[
                                        Sequence  @@   
rectToIteratorList[#,  uvbles]], 
                                      der]; 
Which[
                                    kk[[2]],  
                                    posbad  = {posbad,  #2[[1]]};  {}, 
kk[[1]],   pos =  {pos,  #2[[1]]}; 
IntervalKrawczyk[fff, 
                                      Evaluate[
                                        Sequence  @@   
rectToIteratorList[#,  uvbles]], opts], 
                                    True,  {}]) &,  reg];

newroots = 
                            DeleteCases[newroots,  {}];

  

If[NumericQ[mtol] && Length[newroots]  >  1 && Max[diam/@ newroots] < mtol, newroots  =  
                                FixedPoint[combine,  
                                  newroots]; 

newroots = MapIndexed[(
kk  =  
                                        Krawczyk[fff, 
                                        Evaluate[
                                        Sequence  @@   
rectToIteratorList[#,  uvbles]], 
                                        der]; 
Which[
                                        kk[[2]],  
       posbad  = {posbad,  #2[[1]]};  {}, 
kk[[1]],   pos =  {pos,  #2[[1]]}; 
IntervalKrawczyk[fff, 
                                        Evaluate[
                                        Sequence  @@   
rectToIteratorList[#,  uvbles]], opts], 
                                        True,  {}]) &,  newroots]];

If[
                            newroots ≠ {}, 
                            answer  =  
                              Join[answer, 
                                newroots]; 
reg  =  
                              reg[[Complement[Range[Length[reg]],  
                                    Flatten[Join[pos,  posbad]]]]]]];

If[ 
                      trQ  &&  reg  == {} &&  
                        Length[answer]  >  1,  
                      Print[StringForm[
"Finished last \
round (``) successfully. `` roots have been found and the set is \
complete.", count, Length[answer]]]];

\(If[\ 
                      trQ\  && \ reg\  == {} && \ 
                        Length[answer]\  == \ 1, \ 
                      Print[StringForm[
"Finished last round (``) successfully. 1 root has been found and the set is complete.", 
count]]]\);

If[ 
                      trQ  &&  reg  == {} &&  
                        Length[answer]  ==  0,  
                      Print[
                        StringForm[
"Finished last round \
(``) successfully. There are no roots.", 
count]]];


If[ 
                      trQ  &&  reg  ≠  {}  &&  
                        count  <  maxiter,  
                      Print[StringForm[
"Finished last \
round (``) successfully. `` roots have been found, but there are `` \
unresolved regions.", count,  Length[answer], Length[reg]]]];

If[ 
                      trQ  &&  reg  ≠  {}  &&  
                        count  ==  maxiter,  
                      Print[
"Reached maxiterations \
without resolution."]];

If[count == maxiter,  
                      Message[FindRootsIntervalNewton::maxiters]];
{answer,  reg}  /. 
                      x_Interval  :>   If[formQ,  
                            IntervalForm,  Identity][x]]];

FindRootsIntervalNewton::maxiters  =  "Max iterations reached; possibly \
more roots. Consider setting MergeTolerance to 10^-6, or smaller to merge \
boxes that might share a root near their border. But there is the possibility \
of a multiple root that has not been verified.";

Algorithm 4.5: ValidateRoots

ValidateRoots::usage  =  "ValidateRoots[ff, iters , rlist \
(, opts)] returns True if each point in rlist is provably a root of to the \
accuracy specified in the AccuracyGoal option. In the case of one dimension \
ff can be a single expression, as opposed to a list, as it must be in higher \
dimensions.";

RootQ::usage  =  "RootQ[f, vbles, r] uses Krawczyk \
condition and returns True if r is a zero of function f from R^n -> R^n. If \
it returns False, then that means that the test fails for the specified \
accuracy goal, which does not mean that r is not a root.";

CheckCorrectness::usage  =  "CheckCorrectness is an option \
to ValidateRoots that asks that the ϵ-inflation method be used to \
check correctness of the roots.";

CheckCompleteness::usage  =  "CheckCompleteness is an \
option to ValidateRoots that asks that the Krawczyk method be used to check \
that the number of roots agrees with the number of roots submitted.";

ValidateRoots::notodo =  "There is nothing to check because \
both tests have been turned off.";




Options[RootQ]  =  {AccuracyGoal  ->  6,  
            WorkingPrecision -> $MachinePrecision};

RootQ[ff_, vbles_List , 
          root_?VectorQ,   opts___Rule]  :=  (
f = 
            If[! ListQ[ff],  {ff},  ff]; 
{wp,  
              ag}  =  {WorkingPrecision,  
                  AccuracyGoal}  /.  {opts}  /.  
              Options[RootQ]; 
tol  = 
            10^-ag; 
dim = 
            Length[uvbles]; 
der = 
            Transpose[D[f,  #]  &  /@  
                vbles]; 
KrawczykQ[f,  
            Sequence @@  
              rectToIteratorList[ RealToInterval[root,  2 tol/(Sqrt[dim]) ],  
                vbles], 
    der])


Options[ValidateRoots]  =  {CheckCompleteness -> True,  
            CheckCorrectness  ->  True,  
            WorkingPrecision ->  $MachinePrecision,  TraceQ -> False,  
            AccuracyGoal -> 6,  
            MaxIterations -> 40,   MergeTolerance  ->  0, 
            ShowIntervals -> False,  
            PartialValidation  ->  False};

ValidateRoots[f_, iters__List , 
          rlist_?VectorQ,   opts___Rule]  :=  
        ValidateRoots[f, iters , List  /@  rlist,   opts]  

ValidateRoots[ff_, iters__List , rlist_?MatrixQ,  
            opts___Rule] := Block[{Experimental`$EqualTolerance  =  0},Module[{   wp, reg, grad, int,  vals,  pos, len2,  
                    maxiter, invder, det, count, trQ, t   },
{ag, trQ,  wp,  maxiter,  artQ, mtol,  
                        corQ,  
                        pvQ}  = { AccuracyGoal , TraceQ,  
                            WorkingPrecision,  MaxIterations,  
                            ShowIntervals,  
                            MergeTolerance,   CheckCorrectness,  
                            PartialValidation}  /.  {opts}  /.  
                        Options[ValidateRoots];
f = 
                      If[! ListQ[ff],  {ff},  ff];
reg =  {Interval  /@  Rest  /@  {iters}};
uvbles  =  
                      First  /@  {iters};
tol  =  10^-ag;
dim = Length[uvbles];
rootlist = 
                      Select[rlist, intervalmemberQbox[reg[[1]], #] &];
   
If[! corQ  &&  pvQ,  
                      Message[ValidateRoots::notodo];  
                      Return[Short[
                          HoldForm[
                            ValidateRoots[ff, iters , rlist,  opts]]]]];

h[m_:  ∞] := 
                      Function[Evaluate[uvbles], Evaluate[N[f, m]]];
der = 
                      Transpose[D[f,  #]  &  /@  uvbles];
ans1 = ans2 = True;
  
If[
                      corQ  &&  ! (ans1 = (And @@ 
                                  Map[
KrawczykQ[f,  
                                        Sequence @@  
                                        rectToIteratorList[ 
                                        RealToInterval[#,  
                                        2 tol/(Sqrt[dim]) ],  
                                        uvbles], 
    der] & , rootlist])),  
Message[
                        ValidateRoots::fail1]];
If[! pvQ,
count = 0;   good = 0;
While[reg  !=  {}  &&  
                                count  <  maxiter,count++;
reg = 
                                  subdivide[reg,  wp];
len1 = Length[reg];
If[
                                  dim == 2  &&  artQ, 
                                  Show[Graphics[
                                      reg  /. { Interval[{aa_,  bb_}], 
                                        Interval[{cc_,  
                                        dd_}]} :> {Hue[
                                        Random[]],  
                                        Rectangle[{aa, cc}, {bb  , 
                                        dd}]}], 
PlotRange -> Rest /@  Take[{iters},  2],  AspectRatio -> Automatic,  
                                    Frame -> True]];

  
pos = 
                                  Flatten [ 
                                    Position[
                                      Apply[h[wp], reg, 1],  _?
                                        containszerovector]] ;
reg = 
                                  reg[[
                                    pos ]];    
len2 = 
                                  Length[reg];

  
If[NumericQ[mtol] && Length[reg]  >  1 && Max[diam/@ reg] < mtol, reg  =  
                                      FixedPoint[combine,  reg]];
If[ trQ,  
                                  Print[
                                    StringForm[
"Round \
``. Out of `` regions searched, `` might contain zeros and need to be \
searched. `` validated so far", count, len1, len2,  good]]];
  

goodpos =  
                                  Table[
 
                                    kk = Krawczyk[f, 
                                        Sequence @@ 
                                        rectToIteratorList[reg[[ii]],  
                                        uvbles] ,  
                                        der]; 
If[
                                      kk[[1]], good++];  
                                    If[kk[[1]]  ||  kk[[2]],  
                                      ii], {ii,  Length[reg]}];
reg  =  
                                  reg[[Complement[Range[Length[reg]],  
                                        DeleteCases[goodpos, Null]]]]];
If[ trQ,  
                            Print[StringForm[
"End of \
round ``. Out of `` regions searched, `` might contain zeros and need to be \
searched. `` validated so far", count, len1, len2,  good]]];If[count == maxiter,  
                            Message[FindRootsIntervalNewton::maxiters]];
ans2  =  
                            Length[rootlist] == good];
ans1  &&  ans2]];

\(ValidateRoots::fail1\  = \ "The set of roots did not meet Krawczyk's condition at the given tolerance.";\)

Algorithm 4.6: FindAndValidateRoots     traceQ needs roots I think

FindAndValidateRoots::usage  =  "FindAndValidateRoots[f, \
iters (, opts)] finds the zeros of the system f on the domain indicated by \
iters and proves they are correct. It uses Krawczyk to get seeds, then \
Newton's method, then ϵ-inflation to verify.";

Options[
            FindAndValidateRoots]  =  {WorkingPrecision ->  
              MachinePrecision,  TraceQ -> False,  
            AccuracyGoal -> 6,  
            MaxIterations -> 40,   MergeTolerance  ->  0, 
            ShowIntervals -> False,  
            PartialValidation  ->  False};

FindAndValidateRoots[f_, iters__List ,   opts___Rule]  :=  
          ValidateRoots[f, iters , List  /@  rlist,   opts]  ;

FindAndValidateRoots[ff_, iters__List ,  opts___Rule] := Module[{  f },
{ag, trQ,  wp,  maxiter,  artQ, 
                  mtol}  = { AccuracyGoal , TraceQ,  
                      WorkingPrecision,  MaxIterations,  ShowIntervals,  
                      MergeTolerance}  /.  {opts}  /.  
                  Options[FindAndValidateRoots];
f = 
                If[! ListQ[ff],  {ff},  ff];
reg =  {Interval  /@  Rest  /@  {iters}};
uvbles  =  First  /@  {iters};
tol  =  10^-ag;
dim = Length[uvbles];seeds = {};
  
count = 0;   good = 0;
h[m_:  ∞] := 
                Function[Evaluate[uvbles], Evaluate[N[f, m]]];
der = 
                Transpose[D[f,  #]  &  /@  uvbles];
invder = Inverse[der];

While[reg  !=  {}  &&  count  <  maxiter,count++;
reg = subdivide[reg,  wp];
len1 = Length[reg];
If[dim == 2  &&  artQ, 
                      Show[Graphics[
                          reg  /. { Interval[{aa_,  bb_}], 
                                Interval[{cc_,  
                                    dd_}]} :> 
{Hue[Random[]],  
                                Rectangle[{aa, cc}, {bb  , 
                                    dd}]}], 
PlotRange -> Rest /@  Take[{iters},  2],  AspectRatio -> Automatic,  
                        Frame -> True]];

  
pos = 
                      Flatten [ 
                        Position[
                          Apply[h[wp], reg, 1],  _?containszerovector]] ;
reg = 
                      reg[[pos ]];    
len2 = Length[reg];

  
If[NumericQ[mtol] && Length[reg]  >  1 && Max[Diam/@ reg] < mtol, reg  =  FixedPoint[combine,  reg]];

If[ trQ,  
                      Print[StringForm[
"Round ``. Out of \
`` regions searched, `` might contain zeros and need to be searched. `` \
validated so far", count, len1, len2,  good]]];
  

goodpos =  
                      Table[
 
                        kk = Krawczyk[f, 
                            Sequence @@ 
                              rectToIteratorList[reg[[ii]],  
                                uvbles] ,  
                            der]; 

If[
                          kk[[1]], good++;  
                          AppendTo[seeds,  
                            kk[[3]]]]; 
 
                        If[kk[[1]]  ||  kk[[2]],  ii], {ii,  
                          Length[reg]}];
reg  =  
                      reg[[Complement[Range[Length[reg]],  
                            DeleteCases[goodpos, Null]]]]];

  

rlist  =  (uvbles  /.  
                        FindRoot[
                          f,  
Evaluate[
                            Sequence  @@ (( 
                                    rectToIteratorList[#,  
                                      uvbles]) /.  {s_,  {t_,  
                                        r_}} :>  {s,  
                                      1/2 (t + 
                                        r)}), 
Jacobian ->  der]  ]  &)  /@  seeds;
  
vv = 
                ValidateRoots[f, iters , rlist, 
                  PartialValidation -> True] ;
rlist];

Table 4.1

The following code uses the FindRootsContourMethod code at the top of this file, which uses the contour method for 2-dimensional root-finding. It generates the data for Table 4.1.

Off[FindRoot::lstol];
Do[Print[
{res, 
Sum[
            Length[FindRootsContourMethod[{fx[x, y], fy[x, y]} , {x,  a, 
                  a + 0.25},  {y,  b,  b + 0.25}, 
                ContourResolution -> res,  
                SeedsOnly -> True]],  {a,  -1,  0.75,  
              0.25}, {b,  -1,  0.75,  0.25}], 
Sum[
            Length[FindRootsContourMethod[{fx[x, y], fy[x, y]} , {x,  a, 
                  a + 0.25},  {y,  b,  b +  .25}, 
                ContourResolution -> res]],  {a,  -1,  0.75,  
              0.25}, {b,  -1,  0.75,  0.25}]}],  {res,  20,  400,  
        20}]