Using NDSolve with EventLocator to solve higher-order differential equations with jumping nonlinearity

how to figure the Fučík spectrum of a fourth-order boundary value problem, related to a model of suspension bridge

Jiří Benedikt

Západočeská univerzita v Plzni, Plzeň, Czech Republic

Introduction

Mathematical model of a suspension bridge

There are many mathematical models of suspension bridge, at different levels of simplification. One of the most simplified models is the fourth-order Dirichlet boundary value problem

(u^((4)))(t)=f(u(t)),          t ∈ [0, 1],
u(
              0) = u' (0) = u(1) = u' (1) = 0,

where u : [0, 1]->R (the solution) stands for the displacement of the bridge deck and f:R->R describes the external force. The boundary conditions are called Dirichlet and mean that the ends of the bridge deck are clamped.

The external force f depends on the displacement u, as it is typical for the suspension bridge - if the bridge deck moves downwards (i.e. u>0 - we consider the downwards direction to be positive), the cable stays pull it back, if it moves upwards (u<0), the cable stays do nothing.

Let us assume that f depends linearly on u, but with different constants for u positive and u negative, i.e. f=(SuperPlus[μu])-(SuperMinus[νu]), where \(TraditionalForm\`\(u\^+\) = max {u, 0}\) is the positive and \(TraditionalForm\`\(u\^-\) = max {\(-u\), 0}\) the negative part of u (hence, \(u = \(u\^+\) - \(u\^-\)\)). This is satisfied, if the cable stays behave like a spring with the modulus μ for u positive and the modulus ν for u negative (ν=0 in our case). Note that such f is called a jumping nonlinearity.

PosPart[a_] := (1/2)*(a + Abs[a]); 
NegPart[a_] := (1/2)*(Abs[a] - a); 
Show[GraphicsArray[
  (Plot[#1[[1]]*PosPart[t] - #1[[2]]*NegPart[t], 
    {t, -1, 1}, AspectRatio -> Automatic, 
    DisplayFunction -> Identity] & ) /@ 
  {{0.2, 0.7}, {1, 0}, {-0.5, -0.5}}]];

The Fučík spectrum

We arrive at the fourth-order Dirichlet boundary value generalized eigenvalue problem

(u^((4)))(t)=(SuperPlus[μu])(t)-(SuperMinus[νu])(t),          t ∈ [0, 1],
u(
              0) = u' (0) = u(1) = u' (1) = 0.

the Fučík (generalized) spectrum is defined as the set of all couples (μ,ν) ∈ (R^2) such that the problem has a nontrivial solution (eigenfunction).

If f is at least asymptotically equal to (SuperPlus[μu])-(SuperMinus[νu]), i.e.

(lim_(s -> +∞)) ((f(s))/s)=μ   and   (lim_(s -> -∞)) ((f(s))/s)=ν,

and (μ, ν) does not belong to the Fučík spectrum, then the topological degree theory can be used to prove that the solution is a priori bounded (otherwise there may exist an unbounded sequence of solutions). Hence it is important to know the Fučík spectrum well. Our aim is to figure the Fučík spectrum by Mathematica.

Shooting method

Initial value problem

Let us consider the fourth-order initial value problem (IVP)

(u^((4)))(t)=(SuperPlus[μu])(t)-(SuperMinus[νu])(t),          t ∈ [0, 1],
u(0) = u' (0) = 0, u'' (0) = 1, u''' (0) = δ,

where μ, ν > 0 and δ∈R. The solution can be proved to be unique, and so we may denote by V=V(μ, ν, δ) the value of the solution at t = 1. Obviously, V is an increasing continuous function of δ, and (lim_(δ -> ±∞))\(V(μ, ν, δ)\)=±∞. Hence there exists exactly one constant (Subscript[δ, 0]) = (Subscript[δ, 0])(μ, ν), such that V(μ, ν, (Subscript[δ, 0])(μ, ν))=0. Since clearly V(μ, ν, δ)>0 for any δ≥0, it must be (Subscript[δ, 0])(μ, ν)<0 for all μ,ν>0. Let D=D(μ, ν) denote the first derivative of the solution of the IVP with δ=(Subscript[δ, 0])(μ, ν). If D(μ, ν)=0, then (μ, ν) belongs to the Fučík spectrum. Moreover, (ν, μ) belongs to the Fučík spectrum, too, and the corresponding eigenfunction is -u.

On the other hand, let us assume that (μ, ν), μ, ν > 0, belongs to the Fučík spectrum and let u be a corresponding eigenfunction. It can be proved that u'' (0) ≠ 0. Let us assume first u'' (0) > 0. Since αu, α > 0, is an eigenfunction, too, we may suppose that u'' (0) = 1. Then, clearly, u is the solution of the IVP with δ=(Subscript[δ, 0])(μ, ν), and so D(μ, ν)=0. If u'' (0) < 0, then -u is an eigenfunction, corresponding to (ν, μ) with (-u)'' (0) > 0. Consequently, D(ν, μ) = 0 in this case.

We showed that (μ, ν), μ, ν > 0, belongs to the Fučík spectrum, if and only if D(μ, ν)=0 or D(ν, μ)=0. Hence it may be figured as a union of the zero contours of D(μ, ν) and D(ν, μ). Our aim now reduces to implementation of D in Mathematica.

Implementation

The first problem: precision

From now on we write (μ^4) and (ν^4) instead of μ and ν in the differential equation. First we implement the function V. The simplest one might be

ValueAtOne[μ_, ν_, δ_] := 
  u[1] /. NDSolve[
    {Derivative[4][u][t] == 
      μ^4*PosPart[u[t]] - ν^4*NegPart[u[t]], 
      u[0] == 0, 
      Derivative[1][u][0] == 0, 
      Derivative[2][u][0] == 1, 
      Derivative[3][u][0] == δ}, 
    u, {t, 0, 1}][[1]];

Now we look for (Subscript[δ, 0])(μ, ν) as a root of V(μ, ν, δ).

p1 = Plot[ValueAtOne[40, 30, t], {t, -45, -35}, 
       DisplayFunction -> Identity]; 
p2 = Plot[ValueAtOne[40, 30, t], 
       {t, -40.37441848174, -40.37441848177}, 
       DisplayFunction -> Identity]; 
Show[GraphicsArray[{p1, p2}]];

We see that if we want to guarantee, e.g., |V(μ, ν, (Subscript[δ, 0])(μ, ν))|<(10^(-6)), then the "machine precision" is not enough. So let us improve the implementation of V.

ValueAtOne[μ_, ν_, δ_, goal_, working_] := 
  u[1] /. NDSolve[
    {Derivative[4][u][t] == 
      SetPrecision[μ^4, working]*PosPart[u[t]] -
      SetPrecision[ν^4, working]*NegPart[u[t]], 
      u[0] == 0, 
      Derivative[1][u][0] == 0, 
      Derivative[2][u][0] == 1, 
      Derivative[3][u][0] == SetPrecision[δ, working]}, 
    u, {t, 0, 1}, PrecisionGoal -> goal, 
    WorkingPrecision -> working][[1]];

p1 = Plot[ValueAtOne[40, 30, t, 20, 40], 
       {t, -45, -35}, 
       DisplayFunction -> Identity]; 
p2 = Plot[ValueAtOne[40, 30, t, 20, 40], 
       {t, -40.37441848174, -40.37441848177}, 
       DisplayFunction -> Identity]; 
Show[GraphicsArray[{p1, p2}]];

Increasing of the working precision caused a significant change of the obtained values. To find an optimal precision we can compare the values for different precisions.

TableForm[
  Table[ValueAtOne[40, 30, -40.37441848177, n, 2*n], 
   {n, 15, 25}], 
  TableHeadings -> {Table[2*n, {n, 15, 25}], 
    {"working", "result"}}]

Obviously, PrecisionGoal 20 and WorkingPrecision 40 is enough if we want the (absolute) precision of the result to be less than (10^(-6)). But what happened when the WorkingPrecision was 48!? Let us plot the solution for this case.

sol = u /. NDSolve[
    {Derivative[4][u][t] == 
      SetPrecision[40^4, 48]*PosPart[u[t]] -
      SetPrecision[30^4, 48]*NegPart[u[t]], 
      u[0] == 0, 
      Derivative[1][u][0] == 0, 
      Derivative[2][u][0] == 1, 
      Derivative[3][u][0] == 
        SetPrecision[-40.37441848177, 48]}, 
    u, {t, 0, 1}, PrecisionGoal -> 24, 
    WorkingPrecision -> 48][[1]]
Plot[sol[t], {t, 0, 1}, 
     PlotRange -> {{0, 1}, {-0.01, 0.01}}];

The problem is that the solution was not computed on the whole interval [0, 1] and extrapolation was used. We can fix the problem using the StiffnessSwitching method.

ValueAtOne[μ_, ν_, δ_, goal_, working_] := 
  u[1] /. NDSolve[
    {Derivative[4][u][t] == 
      SetPrecision[μ^4, working]*PosPart[u[t]] -
      SetPrecision[ν^4, working]*NegPart[u[t]], 
      u[0] == 0, 
      Derivative[1][u][0] == 0, 
      Derivative[2][u][0] == 1, 
      Derivative[3][u][0] == SetPrecision[δ, working]}, 
    u, {t, 0, 1}, PrecisionGoal -> goal, 
    WorkingPrecision -> working, 
    Method -> {StiffnessSwitching}][[1]];

TableForm[
  Table[ValueAtOne[40, 30, -40.37441848177, n, 2*n], 
   {n, 15, 45}], 
  TableHeadings -> {Table[2*n, {n, 15, 45}], 
    {"working", "result"}}]

The equation was always solved on the whole [0, 1]. However, the results are obviously not consistent.

The second problem: jumping nonlinearity

The reason why we got such nonconsistent results is that our f (the right-hand side) is the jumping nonlinearity which is not differentiable at 0. It brings unpredictable errors into the solution as it changes its sign. Consequently, we must stop the integration whenever the solution attains the zero value, and restart the integration. This is easily accomplished with the EventLocator method, new in Mathematica 5.1.

state = First[NDSolve`ProcessEquations[
  {Derivative[4][u][t] == 
    SetPrecision[40^4, 130]*PosPart[u[t]] - 
    SetPrecision[30^4, 130]*NegPart[u[t]], 
    u[0] == 0, 
    Derivative[1][u][0] == 0, 
    Derivative[2][u][0] == 1, 
    Derivative[3][u][0] == 
      SetPrecision[-40.37441848177, 130]}, u, t, 
  Method -> {EventLocator, "Event" -> u[t], 
    Method -> StiffnessSwitching}, 
  WorkingPrecision -> 130, AccuracyGoal -> 30, 
  PrecisionGoal -> 30, MaxSteps -> Infinity]]; 
result = Reap[
  While[Sow[state["CurrentTime"["Forward"]]] < 1, 
    state = NDSolve`ReinitializeVector[state, 
      state["CurrentTime"["Forward"]], 
      state["SolutionVector"["Forward"]]];   
    NDSolve`Iterate[state, 1]; ]; 
  state["SolutionVector"["Forward"]]]; 
SetPrecision[result[[1]][[1]],30]

The used precision appears to be sufficient for any μ, ν ∈ [0, 50]. Using the EventLocator is not appropriate when μ = ν, and so we finally define

ValueAtOne[order_, μ_, ν_, δ_, goal_, working_] := 
  If[Abs[μ - ν] < 10^(-10), 
    SetPrecision[Derivative[order][u][1] /. 
      NDSolve[{Derivative[4][u][t] == 
        SetPrecision[μ^4, working]*PosPart[u[t]] - 
        SetPrecision[ν^4, working]*NegPart[u[t]], 
        u[0] == 0, 
        Derivative[1][u][0] == 0, 
        Derivative[2][u][0] == 1, 
        Derivative[3][u][0] == SetPrecision[δ, working]}, 
        u, {t, 0, 1}, 
        Method -> StiffnessSwitching, 
        WorkingPrecision -> working, 
        AccuracyGoal -> goal, PrecisionGoal -> goal, 
        MaxSteps -> Infinity][[1]][[1]], goal], 
    Module[{state, result}, 
      state = First[NDSolve`ProcessEquations[
        {Derivative[4][u][t] == 
          SetPrecision[μ^4, working]*PosPart[u[t]] - 
          SetPrecision[ν^4, working]*NegPart[u[t]], 
          u[0] == 0, 
          Derivative[1][u][0] == 0, 
          Derivative[2][u][0] == 1, 
          Derivative[3][u][0] == 
            SetPrecision[δ, working]}, 
        u, t, Method -> {EventLocator, "Event" -> u[t], 
          Method -> StiffnessSwitching}, 
        WorkingPrecision -> working, 
        AccuracyGoal -> goal, PrecisionGoal -> goal, 
        MaxSteps -> Infinity]]; 
      result = Reap[While[Sow[state["CurrentTime"[
        "Forward"]]] < 1, 
        state = NDSolve`ReinitializeVector[state, 
          state["CurrentTime"["Forward"]], 
          state["SolutionVector"["Forward"]]]; 
        NDSolve`Iterate[state, 1]; ]; 
        state["SolutionVector"["Forward"]]]; 
      SetPrecision[result[[1]][[order + 1]], goal]]];

ValueAtOne[0, 40, 30, -40.37441848177, 30, 130] // Timing

Having the function V implemented, we can easily compute (Subscript[δ, 0])(μ, ν) and D(μ, ν).

LookForNeg[fction_] := 
  NestWhile[#1*2 & , -1, fction[#1] > 0 & ]; 
InitInt[x_, prec_] := SetPrecision[
  If[x == -1, {-1, 0}, {x, x/2}], prec]; 
Halving[{l_, r_, f_}] := If[f[(l + r)/2] < 0, 
  {(l + r)/2, r, f}, {l, (l + r)/2, f}]; 
MyFindRoot[fction_, goal_, working_] := 
  NestWhile[Halving, Append[
    InitInt[LookForNeg[fction], working], fction], 
    #1[[2]] - #1[[1]] > 10^(-goal) & ][[1]]; 
DerZero[μ_, ν_] := ValueAtOne[1, μ, ν, 
  MyFindRoot[ValueAtOne[0, μ, ν, #1, 30, 130] & , 
  25, 35], 30, 130];

MyFindRoot[ValueAtOne[0, 40, 30, #, 30, 130] & , 25, 35]
DerZero[40, 30]

Finally, the Fučík spectrum may be figured by the ContourPlot function as follows
ContourPlot[DerZero[μ, ν]*DerZero[μ, ν],
  {μ, 0, 50}, {ν, 0, 50}, ContourShading -> None, 
  Contours -> {0}, PlotPoints -> 200]
But it would take months...

Optimization

The ContourPlot method evaluates the function uniformly at all n×n points. We develop an adaptive iterative algorithm for figuring the zero contour of a function. It evaluates the function more densely near already found lines to make them smoother, but still it searches for another components of the contour.

Let us consider the following test function.

f[x_, y_] = Sin[x + Cos[x + y]]*Sin[y + 0.2] + 
  x / 25 - y / 8 - 0.08;
p = ContourPlot[f[x, y], {x, 0, 10}, {y, 0, 10},
  ContourShading -> None, Contours -> {0}, 
  ContourStyle -> Hue[0], PlotPoints -> 65];

To draw the picture, ContourPlot needed to evaluate the function at 65×65 = 4225 points. Our algorithm follows.

size = 10;
proporref = .9;
InsertNew[main_, pos_] := Module[{am, new, numnew},
  am = main;
  new = Complement[pos, am[[2]] /. {x_?NumericQ, y_, f_} -> {x, y}];
  numnew = Length[new];
  am = {am[[1]], Join[am[[2]], new]};
  am = am /. {x_?NumericQ, y_} -> {x, y, f[x, y]};
  am = {am[[1]], Sort[am[[2]], #1[[2]] < #2[[2]] || #1[[2]] == #2[[2]] && #1[[1]] < #2[[1]] & ]};
  {numnew, am}];
Uni[main_] := Module[{am, posnew, numnew},
  am = MapAt[# + 1 & , main, {1, 1}];
  posnew = Flatten[Table[{j * size, i * size}, {i, 0, 1, 2 ^ - am[[1]][[1]]}, {j, 0, 1, 2 ^ - am[[1]][[1]]}], 1];
  {numnew, am} = InsertNew[am, posnew];
  am = MapAt[# - ((2^am[[1]][[1]] + 1) ^ 2 - (am[[1]][[3]] + numnew)) & , am, {1, 4}];
  ReplacePart[am, (2^am[[1]][[1]] + 1) ^ 2, {1, 3}]];
RefDiv[main_, lev_] := Module[{am, posnew, numnew},
  am = main;
  amchoice = Select[am[[2]], Mod[Take[#, 2],size * 2 ^ -lev] == {0, 0} & ];
  posnew = Select[amchoice, Module[{thenext}, thenext = Select[amchoice, Function[a, Take[#, 2] + {size, 0} * 2^ - (lev - 1) ==
    Take[a, 2]], 1]; thenext =!= {} && #[[3]] * thenext[[1]][[3]] <= 0] & ] /. {x_?NumberQ, y_, f_} -> {x + size * 2 ^ -lev, y};
  posnew = Join[posnew, Select[amchoice, Module[{thenext}, thenext = Select[amchoice, Function[a, Take[#, 2] + {0, size} * 2 ^ -(lev - 
    1) == Take[a, 2]], 1];thenext =!= {} && #[[3]] * thenext[[1]][[3]] <= 0] & ] /. {x_?NumberQ, y_, f_} -> {x, y + size * 2 ^ -lev}];
  {numnew, am} = InsertNew[am, posnew];
  MapAt[# + numnew &, am, {1, 4}]];
RefAdd[main_, lev_] := Module[{am, posnew, numnew},
  am = main;
  amchoice = Select[am[[2]], Mod[Take[#, 2],size * 2 ^ -lev] == {0, 0} & ];
  posnew = Flatten[Select[am[[2]], Module[{thenext}, thenext = Select[amchoice, Function[a, Take[#, 2] + {size, 0} * 2 ^ -lev == 
    Take[a, 2]], 1]; thenext =!= {} && #[[3]] * thenext[[1]][[3]] <= 0] & ] /. {x_?NumberQ, y_, f_} -> {{x, y - size * 2 ^ -lev}, 
    {x, y + size * 2 ^ -lev}, {x + size * 2 ^ -lev, y - size * 2 ^ -lev}, {x + size * 2 ^ -lev, y + size * 2 ^ -lev}}, 1];
  posnew = Join[posnew, Flatten[Select[am[[2]], Module[{thenext}, thenext = Select[amchoice, Function[a, Take[#, 2] + {0, size} * 2 ^ 
    -lev == Take[a, 2]], 1]; thenext =!= {} && #[[3]] * thenext[[1]][[3]] <= 0] & ] /. {x_?NumberQ, y_, f_} -> {{x - size * 2 ^ -lev, y}, 
    {x + size * 2 ^ -lev, y}, {x - size * 2 ^ -lev, y + size * 2 ^ -lev}, {x + size * 2 ^ -lev, y + size * 2 ^ -lev}}, 1]];
  posnew = Union[Select[posnew, #[[1]] ≥ 0 && #[[1]] ≤ size && #[[2]] ≥ 0 && #[[2]] ≤ size & ]];
  {numnew, am} = InsertNew[am, posnew];
  MapAt[# + numnew & , am, {1, 4}]];
MyRefine[main_] := Module[{am}, am = main; For[i = 1, i ≤ am[[1]][[2]], i++, am = RefDiv[am, i]; am = RefAdd[am, i]]; am];
Step[main_] := Module[{am, ifrefine, orig}, am = main; ifrefine = am[[1]][[4]] / Plus @@ am[[1]][[{3, 4}]] < proporref; 
  If[ifrefine, am = MapAt[# + 1 & , am, {1, 2}]; orig = am[[1]][[4]]; am = FixedPoint[MyRefine[#] & , am]]; 
  If[!ifrefine || orig == am[[1]][[4]], am=Uni[am]; am = FixedPoint[MyRefine[#] & , am]; am = MapAt[Max[#, am[[1]][[1]]] &, am, {1, 2}]];
  am];
interp[main_, x_, y_] := main[[2]][[Position[#, Min[#]][[1, 1]] & [Norm[{x, y} - Take[#, 2]] & /@ main[[2]]]]][[3]]
depict[main_] := Module[{p1, p2},
  p1 = ContourPlot[interp[main, x, y], {x, 0, size}, {y, 0, size}, PlotPoints -> 2 ^ main[[1]][[2]] + 1, 
    Contours -> {0}, ContourShading -> None, DisplayFunction -> Identity];
  p2 = Graphics[{PointSize[.012], main[[2]] /. {x_?NumberQ, y_, f_} -> Point[{x, y}]}];
  Show[p, p1, p2, DisplayFunction -> Identity, PlotRange -> {{-.1, size + .1}, {-.1, size + .1}}]];

Val = {{{0, 0, 4, 0}, Flatten[Table[{j * size, i * size, f[j * size, i * size]}, {i, 0, 1}, {j, 0, 1}], 1]}};
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
lastpic = Show[depict[Last[Val]], DisplayFunction -> $DisplayFunction];

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;

We get the same result, needing only 19% of the points, that standard ContourPlot needed. Of course, we can continue in the refinement.

Val = Append[Val, Step[Last[Val]]];
Block[{eval = Plus @@ Last[Val][[1]][[{3, 4}]], allc = (2 ^ Last[Val][[1]][[2]] + 1) ^ 2},
  ToString[Last[Val][[1]]] <> ", " <> ToString[eval] <> " / " <> ToString[allc] <> " (" <> ToString[eval / allc * 100.] <> "%)"]
newpic = depict[Last[Val]];
Show[GraphicsArray[{lastpic, newpic}]];
lastpic = newpic;