Code from the Package Chapter4Code

Definition of f(x, y) and other functions related 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;

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[uvbles]; 
            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++\); 

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[
                      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]]]]; 

 (*\ if\ there\ is\ one\ box\ left, \ 
                      try\ to\ use\ rootfinding\ *) 
If[
                      NumericQ[mtol]\  && \ Length[rects]\  > \ 1\  && \ 
                        Max[diam /@ \ rects]\  < \ mtol, \ 
                      rects\  = \ 
                        FixedPoint[combine, \ rects]]; 
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]; 
If[trQ, \ 
                          Print[
StringForm["Finished; no. of current boxes = ``, total rects = ``, ``; total function evals = ``", Length[rects] + Length[smallrects], Flatten[\ rectcount], \ 
                              Plus\  @@ \ Flatten[rectcount], \ 
                              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[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]]];

Sitzung S. 115

IntervalMinimize[f[x, y], {x, -1, 1}, {y, -1, 1}, -3.24, 
      AccuracyGoal -> 102, WorkingPrecision -> 110]