An efficient way to figure spectra and bifurcation diagrams of boundary value problems

Jiří Benedikt

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

The author is supported by grants 1N04078 and MSM 4977751301 of the Ministery of Education, Youth and Sports of the Czech Republic.

In the Wolfram Technology Conference 2005, we showed that the so-called Fučík spectrum of the fourth-order Dirichlet boundary value problem

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

which bears on stability of suspension bridges, can be figured as the zero-contour of a function of two variables. In the same way, spectra and bifurcation diagrams of many boundary value problems can be depicted. As it showed that evaluation of this function at one point may take even minutes (in the case of the fourth-order problem), it is not appropriate to use the built-in ContourPlot function which needs values at all n×n points.

We present an algorithm which needs values at much fewer points — the dependence of the number of evalution points on n is linear instead of the quadratic one (under several assumptions on the zero-contour). Our algorithm is iterative — the user may apply it on the previous state until he is contented with the result. The algorithm chooses evaluation points in order to both refine the curves and search for another components of the zero-contour, with user-predefined proportion.

The standard 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 close to already found lines to make them smoother, but still it searches for another components of the zero contour.

Let us consider the following test function.

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

To draw the picture, ContourPlot needed to evaluate the function at 257×257 = 66049 points. Implementation of our algorithm follows (below, the algorithm is demonstrated on an example of the above function f).

ZLInit[f_, xmin_?NumericQ, xmax_?NumericQ, ymin_?NumericQ, 
          ymax_?NumericQ, proporref_?NumericQ] := 
        Module[{state = {{f, xmin, xmax, ymin, ymax, proporref, 0, 0, 0, 0}, 
                Table[0, {2}, {2}], Table[0, {2}, {2}]}}, 
          ZLEval[state, Flatten[Table[{i, j}, {i, 1, 2}, {j, 1, 2}], 1], 
            True]];
ZLAddLev[state_] := 
        Module[{ls = 
              state}, ls[[1, 8]]++; 
ls[[2]] = 
            ls[[2]] /. 
              vec_?VectorQ :> 
                Drop[Flatten[
                    Transpose[{vec, 
                        Table[0, {2^(ls[[1, 8]] - 1) + 1}]}], 
                    1], -1]; 
ls[[2]] = 
            Drop[Flatten[
                Transpose[{ls[[2]], 
                    Table[0, {2^(ls[[1, 8]] - 1) + 
                          1}, {2^ls[[1, 8]] + 1}]}], 
                1], -1]; 
ls[[3]] = 
            ls[[3]] /. 
              vec_?VectorQ :> 
                Drop[Flatten[
                    Transpose[{vec, 
                        Table[0, {2^(ls[[1, 8]] - 1) + 1}]}], 
                    1], -1]; 
ls[[3]] = 
            Drop[Flatten[
                Transpose[{ls[[3]], 
                    Table[0, {2^(ls[[1, 8]] - 1) + 
                          1}, {2^ls[[1, 8]] + 1}]}], 
                1], -1]; 
ls];
ZLRefDiv[state_, lev_?IntegerQ] := 
        Module[{dif, indeces}, 
dif = 
            2^(state[[1, 8]] - lev); 
indeces = 
            Flatten[Table[
                Select[Flatten[
                    Table[{i*dif + 1, j*dif + 1}, {j, dir, 2^lev - dir, 
                        2}, {i, 1 - dir, 2^lev - (1 - dir), 2}], 1], 
                  state[[3, #[[2]] + 
                              If[dir == 0, 
                                0, {-1, 0, 1}*dif], #[[1]] + 
                              If[dir == 0, {-1, 0, 1}*dif, 
                                0]]] == {1, 0, 1} && 
                      Times @@ state[[2, #[[2]] + 
                                If[dir == 0, 
                                  0, {-1, 1}*dif], #[[1]] + 
                                If[dir == 0, {-1, 1}*dif, 0]]] < 
                        0 &], {dir, 0, 1}], 1]; 
ZLEval[
            state, indeces, False]];
ZLRefNei[state_, lev_?IntegerQ] := 
        Module[{dif, pos, indeces}, 
dif = 
            2^(state[[1, 8]] - lev); 
pos = 
            Flatten[{{#, Map[Reverse, #, {2}]} &[
                    Flatten[
                      Table[{{1*(-1)^j, 0}, {1*(-1)^j, 
                            1*(-1)^i}}, {i, 0, 1}, {j, 0, 1}], 1]], 
                  Table[{{(-1)^j, 0}, {0, (-1)^i}}, {i, 0, 
                      1}, {j, 0, 1}]}, 2]*dif; 
indeces = 
            Select[Flatten[
                Table[{i*dif + 1, j*dif + 1}, {j, 0, 2^lev}, {i, 0, 2^lev}], 
                1], Function[index, 
                state[[3, index[[2]], index[[1]]]] == 0 && 
                  Or @@ ((#[[1, 2]] == 
                                1 && #[[2, 2]] == 
                                1 && #[[1, 1]]*#[[2, 1]] < 
                                0) & /@ (Table[
                                state[[{2, 
                                      3}, #[[i, 2]], #[[i, 
                                        1]]]], {i, 1, 2}] & /@ 
                            Select[
                              Table[{#, #} &[{index[[1]], 
                                      index[[2]]}] + pos[[i]], {i, 1, 
                                  Length[pos]}], 
                              Min[#] ≥ 1 && 
                                  Max[#] ≤ 
                                    2^state[[1, 8]] + 
                                      1 &]))]]; 
ZLEval[
            state, indeces, False]];
ZLUni[state_] := 
        Module[{ls = state, lev, dif, 
            indeces}, 
lev = ++ls[[1, 
                  7]]; 
dif = 
            2^(ls[[1, 8]] - lev); 
indeces = 
            Select[Flatten[
                Table[{i*dif + 1, j*dif + 1}, {dir, 0, 1}, {j, dir, 
                    2^lev - dir, 2}, {i, 1 - dir, 2^lev - (1 - dir), 
                    2 - dir}], 2], 
              ls[[3, #[[2]], #[[1]]]] == 
                  0 &]; 
ZLEval[ls, indeces, 
            True]];
ZLEval[state_, ind_, uni_] := 
        Module[{ls = 
              state}, ({ls[[2, #[[2]], #[[1]]]], 
                      ls[[3, #[[2]], #[[1]]]]} = {ls[[1, 
                            1]][
                        ls[[1, 
                              2]] + (ls[[1, 3]] - 
                                ls[[1, 2]])*(#[[1]] - 1)/
                              2^ls[[1, 8]], 
                        ls[[1, 
                              4]] + (ls[[1, 5]] - 
                                ls[[1, 4]])*(#[[2]] - 1)/
                              2^ls[[1, 8]]], 1};) & /@ 
            ind; 
If[uni, 
            ls[[1, 9]] = (2^ls[[1, 7]] + 1)^2; 
            If[ls[[1, 7]] > 0, 
              ls[[1, 
                    10]] -= ((2^ls[[1, 7]] + 
                            1)^2 - (2^(ls[[1, 7]] - 1) + 
                            1)^2) - Length[ind]], 
            ls[[1, 10]] += 
              Length[ind]]; 
  
ls];
\(\(ZLRefine[state_] := 
        Module[{ls = state}, 
          For[i = state[\([1, 7]\)] + 1, 
            i ≤ 
              state[\([1, 8]\)], \(i++\),  (*\(Print\)\([\)\("RefDiv " <> 
                  ToString[
                    i] <> " " <> \(\(ToString\)\([\)\(\(Timing\)\([\)\)\)\)*) ls = ZLRefDiv[ls, 
                i] (*\(\(\(]\) \([\([1]\)]\)\)\(]\)\)\(]\)*) ;  (*\(Print\)\([\)\("RefNei " <> 
                  ToString[
                    i] <> " " <> \(\(ToString\)\([\)\(\(Timing\)\([\)\)\)\)*) ls = ZLRefNei[ls, i] (*\(\(\(]\) \([\([1]\)]\)\)\(]\)\)\(]\)*) ]; 
          ls];\)\)
ZLRefineRep[state_] := 
        Module[{ls = state, pts}, pts = Plus @@ ls[[1, {9, 10}]]; 
          ls = ZLRefine[ls]; 
          While[pts < Plus @@ ls[[1, {9, 10}]], 
            pts = Plus @@ ls[[1, {9, 10}]]; ls = ZLRefine[ls]]; 
          ls];
ZLUnEval[values_, bools_] := 
        Module[{ls = values, maxi, val, fi}, 
maxi = 
            Length[ls[[1]]]; 
For[i = 1, 
            i ≤ maxi, i++, 
            val = Null; 
For[j = 1, 
              j ≤ 
                maxi, j++, If[bools[[i, j]] == 1, 
                  If[val === Null, fi = j]; val = ls[[i, j]], 
                  ls[[i, j]] = val];]; 
If[
              val =!= Null, val = ls[[i, fi]]; 
              For[j = fi - 1, 
                j ≥ 
                  1, j--, ls[[i, j]] = 
                    val;];];]; 
For[j = 1, 
            j ≤ maxi, j++, val = Null; 
            For[i = 1, 
              i ≤ 
                maxi, i++, If[ls[[i, j]] =!= Null, 
                  If[val === Null, fi = i]; 
                  val = ls[[i, j]], ls[[i, j]] = 
                      val;];]; 
val = ls[[fi, j]]; 
            For[i = fi - 1, 
              i ≥ 
                1, i--, ls[[i, j]] = 
                  val;];]; 
ls];
ZLUnEvalCom = 
        Compile[{{values, _Real, 2}, {bools, _Integer, 2}}, 
          ZLUnEval[values, bools]];
\(\(ZLStep[state_] := 
        Module[{ls = state, ifrefine, pts}, 
          pts = Plus @@ ls[\([1, {9, 10}]\)]; 
          ifrefine = ls[\([1, 10]\)]/pts < state[\([1, 6]\)]; 
          If[ifrefine, ls = ZLAddLev[ls]; 
            ls = ZLRefineRep[ls]]; 
If[\(! ifrefine\) || 
              ls[\([1, 10]\)] == 0, ls = ZLUni[ls]; 
            ls = ZLRefineRep[
                ls]]; 
 (*\(Print\)\([\)\("UnEval " <> \(\(ToString\)\([\)\(\(Timing\)\([\)\)\)\)*) ls[\([2]\)] = 
            ZLUnEvalCom[ls[\([2]\)], 
              ls[\([3]\)]] (*\(\(\(]\) \([\([1]\)]\)\)\(]\)\)\(]\)*) ; 
ls];\)\)
ZLShow[state_, rules_] := 
        Show[ContourGraphics[state[[2]], Contours -> {0}, 
            ContourShading -> False, 
            MeshRange -> {Take[state[[1]], {2, 3}], 
                Take[state[[1]], {4, 5}]}], rules, 
          PlotRange -> {Take[state[[1]], {2, 3}], 
              Take[state[[1]], {4, 5}], {-1, 
                1}}];
ZLShow[state_] := ZLShow[state, {}];
ZLShowPts[state_, ptsize_, rules_] := 
        Show[ZLShow[state, DisplayFunction -> Identity], 
          Graphics[{PointSize[
                ptsize], Point[{state[[1, 
                            2]] + #[[2]]/(2^
                                state[[1, 8]])*(state[[1, 3]] - 
                              state[[1, 2]]), 
                      state[[1, 
                            4]] + #[[1]]/(2^
                                state[[1, 8]])*(state[[1, 5]] - 
                              state[[1, 4]])}] & /@ 
                Select[Flatten[
                    Table[{i, j}, {i, 0, 2^state[[1, 8]]}, {j, 0, 
                        2^state[[1, 8]]}], 1], 
                  state[[3, #[[1]] + 1, #[[2]] + 1]] == 
                      1 &]}], rules, 
          PlotRange -> ({{1.02, -0.02}, {-0.02, 
                        1.02}} . # & /@ {Take[state[[1]], {2, 3}], 
                  Take[state[[1]], {4, 5}]}), 
          DisplayFunction -> $DisplayFunction];
ZLShowPts[state_, ptsize_] := 
        ZLShowPts[state, ptsize, {}];
ZLInfo[state_] := 
      Block[{evals = Plus @@ state[[1, {9, 10}]], 
          allc = (2^state[[1, 8]] + 1)^2}, 
        ToString[state[[1, {7, 8, 9, 10}]]] <> ", " <> 
          ToString[evals] <> " / " <> ToString[allc] <> " (" <> 
          ToString[evals/allc*100. ] <> "%)"]

state = 
        ZLInit[f, 0. , 10. , 0. , 10. , 0.8];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.03, 
              DisplayFunction -> Identity]}]];

The function ZLInit returns the first iteration. The parameters are the name of the function, the zero contour of which we want to figure (f), the domain ([0, 10]×[0, 10]), and a proportion between refinement of already found lines and searching for another components of the zero countour (the value 0.8 means that if the algorithm used more than 80% of the evaluations for the refinement, then it tries to find another component in the next step). The function ZLInit starts with evaluation of the function at the four corners of the given domain (see the right picture). Since all the four values are positive in this case, we have no line for the present (as we see in the left picture — by red colour the zero contour obtained by the starndard ContourPlot is shown).

The function ZLInfo provides us with information about the computation. The first number means how many times the algorithm searched for the components of the zero contour (except the first "try" of ZLInit), the second one how many times it refined the lines, the third one is the number of evaluations used to find components and the fourth one is the number of evaluations used to refine the already found lines. Then ZLInfo prints how effective the algorithm is at this moment — 4 / 4 means that we needed four evaluations and standard ContourPlot would need four evaluation to produce a similar picture (of course, we did not save any evaluation in this first step).

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.03, 
              DisplayFunction -> Identity]}]];

In the second step, returned by ZLStep, the algorithm refines the grid so that the neighbours (by neighbours we mean two points in the same row (column) and neighbouring columns (rows) of the grid) have half the distance (we have 3×3 = 9 points now). First it evaluates the function at those of the new five points which are between two points which were neighbours in the previous grid (this is satisfied by four points in the middle of each side of the domain — they are between the corners) and at which the function has values with oposite signs (this is satisfied at none of the four points since the function is positive at all the corners). The reason is that if the function is continous then its zero contour must intersect the abscissa between two points where the values have oposite signs. Adding a point between them helps to specify the position of the zero contour.

Since the algorithm found no point where to evaluate the function, it refines the grid where the function is evaluated uniformly (like the standard ContourPlot does), so it evaluates the function at all the new five points. ZLInfo now tells that the grid for the uniform evaluation was refined once, and the lines are not refined on a finer grid (the second number is one, too). In fact, we have not found any line yet, since the function is also positive at all the new five points.

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.03, 
              DisplayFunction -> Identity]}]];

Similarly as in the previous step, the algorithm refines the grid, now adding new 16 new points (hence we have 5×5 = 25 points). 12 of them are between neighbours from the previous grid, but none of these 12 couples has oposite signs. So the algorithm evaluates the function on the whole grid again and it founds some lines because the value of the function is negative at several points.

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.03, 
              DisplayFunction -> Identity]}]];

After evaluation at the points between the neighbours from the previous grid (explained above), the algoritm continues by evaluation at points close to neighbours in the new (9×9 = 81) grid having values with oposite signs. For example, if the values at the blue points in the following picture (where the function has been already evaluated) have oposite signs, the algorithm evaluates the function at the green points (if it was not evaluated there before, of course). The reason is that (assuming continuity of the function) the zero contour goes between the blue points and it must leave the squares above and below intersecting at least one another side of each square. This rule is applied recursively until no more point to evaluate the function at is found. Note that ZLInfo now states that the "uniform grid" (on which the function is evaluated everywhere) was refined twice, while the grid for further refinement of the found lines was refined once more. Consequently, the algorithm begins saving some evaluations (31 out of 81 at this moment).

Show[p, 
        Graphics[{RGBColor[0, 0, 1], PointSize[0.03], 
            Point[{2/8*10, 6/8*10}], Point[{3/8*10, 6/8*10}], 
            RGBColor[0, 1, 0], Point[{2/8*10, 5/8*10}], 
            Point[{3/8*10, 5/8*10}], Point[{2/8*10, 7/8*10}], 
            Point[{3/8*10, 7/8*10}]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.02, 
              DisplayFunction -> Identity]}]];

In the previous step, 25 of the total 50 evaluation points were used for refinement. It is less than 80%, so the refinement continues. Now we have 17×17 = 289 points, but the function was evaluated at 119 points only, which is about 41% of points which standard ContourPlot would need. It is a substantial saving, but it is compensated by an obvious disadvantage — the algorithm may omit some components of the zero contour. On the other hand, if we continue in the iteration process, they will be found.

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.02, 
              DisplayFunction -> Identity]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.02, 
              DisplayFunction -> Identity]}]];

After the previous step, the algorithm used 242 out of 267 points for the refinement of the lines, which is more than 80%, so it now refines the grid for uniform evaluation to the third level (9×9), and it found one of the two remaining components. To find the small one, we need several more iterations.

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.01, 
              DisplayFunction -> Identity]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.01, 
              DisplayFunction -> Identity]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.01, 
              DisplayFunction -> Identity]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.01, 
              DisplayFunction -> Identity]}]];

state = ZLStep[state];
ZLInfo[state]
Show[
        GraphicsArray[{Show[p, 
              ZLShow[state, DisplayFunction -> Identity], 
              DisplayFunction -> Identity], 
            ZLShowPts[state, 0.01, 
              DisplayFunction -> Identity]}]];

Now we obtain exactly the same picture as that produced by the standard ContourPlot, needing only 5.8% of the evaluations. This example show that if the user of this algorithm does not know the shape of the zero contour and he cannot foreclose existence of such small components, he should choose the last parameter given to ZLInit closer to zero. If he is sure that all the components will be found on a rough grid, he may choose a number closer to one (may be even equal to one) in order not to waste time finding new components. Of course, we can now continue in the refinement...

state = ZLStep[state];
ZLInfo[state]
ZLShow[state];