Bifurcation Diagrams Using webMathematica

drawing the zero contour of a function with a minimal number of evaluations

Jiří Benedikt

Katedra matematiky
Západočeská univerzita v Plzni
Univerzitní 22
306 14 Plzeň
Czech Republic

Intruduction

Many types of spectra and bifurcation diagrams of (not only) boundary value problems for ordinary differential equations coincide with a zero contour of a real function of two real variables which can be easily implemented using the shooting method. Evaluation of such function may be very slow, especially in the case of nonlinear higher-order differential equations. Hence, our aim is to draw the zero contour using as few evaluations as possible. We show it on the example of the so-called Fučík spectrum of the fourth-order Dirichlet boundary value problem

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

where (SuperPlus[u])=max{u,0} and (SuperMinus[u])=max{-u,0}, i.e. u=(SuperPlus[u])-(SuperMinus[u]). The Fučík spectrum is defined as the set of all couples (μ,ν) ∈ (ℝ^2) such that (1) has a nontrivial solution. Existence of a stable solution of an initial–boundary value problem that describes the oscillation of a suspension bridge is related to the shape of the Fučík spectrum of (1).

Mathematical Model of a Suspension Bridge and the Fučík Spectrum

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

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

where u:[0,1]->ℝ (the solution) stands for the deflection of the bridge deck and f:ℝ->ℝ describes the external force. The boundary conditions are called Dirichlet and they mean that the ends of the bridge deck are clamped. The solution of (2) bears on the stationary solution of a less simplified initial–boundary value problem.

The external force f depends on the deflection 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]). 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. Examples of a graph of such function f follow.

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

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

(3)   (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 of (3) is defined as the set of all couples (μ,ν) ∈ (ℝ^2) such that (3) has a nontrivial solution (eigenfunction).

Shooting Method

Let us consider the fourth-order initial value problem 

(4)   (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 δ∈ℝ. 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 (4) with δ=(Subscript[δ, 0])(μ,ν). If D(μ,ν)=0, then (μ,ν) belongs to the Fučík spectrum of (3). Moreover, (ν,μ) belongs to the Fučík spectrum of (3), too, and the corresponding eigenfunction is -u.

On the other hand, let us assume that (μ,ν), μ,ν>0, belongs to the Fučík spectrum of (3) 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 (4) 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 of (3), 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.

Implementation

To get a reasonable picture of the Fučík spectrum, we first replace μ and ν by (μ^4) and (ν^4) in (4), respectively. If we define the function V using the NDSolve method with the implicit machine working precision:

ValueAtOne[μ_?NumericQ,ν_?NumericQ,δ_?NumericQ]:=Block[{t},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]]];

we cannot find (Subscript[δ, 0]) with sufficient precision:

Plot[ValueAtOne[40,30,δ],{δ,-40.37441848174,-40.37441848177}]

Increased Working Precision

We have to increase the working precision, but we realize that the solution is then not computed on the whole interval [0,1]:

ValueAtOne[μ_?NumericQ,ν_?NumericQ,δ_?NumericQ,goal_?NumericQ,working_?NumericQ]:=Block[{t},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,MaxSteps->Infinity][[1]]];
TableForm[Prepend[Table[{2*n,ValueAtOne[40,30,-40.37441848174,n,2*n]},{n,20,30}],{"working","result"}]]

StiffnessSwitching Method

We fix the problem using the StiffnessSwitching method.

ValueAtOne[order_?NumericQ,μ_?NumericQ,ν_?NumericQ,δ_?NumericQ,goal_?NumericQ,working_?NumericQ]:=Block[{t},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},PrecisionGoal->goal,WorkingPrecision->working,Method->{StiffnessSwitching}][[1]]];
TableForm[Prepend[Table[{2*n,ValueAtOne[0,40,30,-40.37441848174,n,2*n]},{n,20,30}],{"working","result"}]]

Now, we can define the function D:

DerZero[μ_?NumericQ,ν_?NumericQ]:=ValueAtOne[1,μ,ν,δ/.FindRoot[ValueAtOne[0,μ,ν,δ,40,50]==0,{δ,-(μ+ν+5),0},WorkingPrecision->30,AccuracyGoal->10,PrecisionGoal->Infinity,MaxIterations->1000][[1]],40,50];

The time of evaluation of D at one point may attain a minute (running on Intel Core 2 Duo 1.8 GHz with 1GB RAM):

DerZero[40,30]//Timing

Optimization Algorithm

To draw the zero countour of D in a reasonable time, we develop an algorithm for drawing the zero contour with a minimal number of evaluations.

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,8}],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++,ls=ZLRefDiv[ls,i];ls=ZLRefNei[ls,i]];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]];ls[[2]]=ZLUnEvalCom[ls[[2]],ls[[3]]];
ls];
ZLShow[state_,rules_]:=ListContourPlot[Flatten[Table[{state[[1,2]]+i*2^(-state[[1,8]])(state[[1,3]]-state[[1,2]]),state[[1,4]]+j*2^(-state[[1,8]])(state[[1,5]]-state[[1,4]]),state[[2,j+1,i+1]]},{i,0,2^state[[1,8]]},{j,0,2^state[[1,8]]}],1],Contours->{0},ContourShading->False,rules,PlotRange->{Take[state[[1]],{2,3}],Take[state[[1]],{4,5}]},InterpolationOrder->1]
ZLShow[state_]:=ZLShow[state,{}]
ZLShowSym[state_]:=Show[ZLShow[state],ZLShow[{state[[1]],Transpose[state[[2]]],state[[3]]}]]
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.]<>"%)"]

Demonstration

Since the time of evaluation of D increases with μ and ν, we demontrate the algorithm on drawing the Fučík spectrum of (3) in the square [5,20]×[5,20] only. First, we call the function ZLInit to get an initial state. Then we apply the function ZLStep on the previous state until we are contented with the result. The result is shown using the ZLShow function. The function ZLShowPts shows, moreover, the points where the function was evaluated. The constant 0.9, given in ZLInit, means that the algorithm refines the curves of the zero contour that it has already found only if less than 90% of all evaluations was used for this purpose before. Otherwise it tries to find other components of the zero contour in the next step.

total=Timing[state=ZLInit[DerZero,5.,20.,5.,20.,0.9]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.03]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.015]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.015]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.01]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0.01]

total+=Timing[state=ZLStep[state]][[1]]
ZLInfo[state]
ZLShowPts[state,0]

In the last picture, the curves seem to be sufficiently smooth, so we finally display the Fučík spectrum of (3) using the ZLShowSym function that shows a symmetrization of the zero contour (union of the zero contours of D(μ,ν) and D(ν,μ)).

ZLShowSym[state]

Conclusion

The function ZLInfo informs us that the function was evaluated 2093 times while the classical ContourPlot (from, e.g., Mathematica 5.2.2.0) would need 16641 evaluations to give the same result. Since our algorithm took approximately 4 hours, we can deduce that the classical ContourPlot would take about 1.5 day. The version of ContourPlot contained in Mathematica 6.0.1.0 also uses refinement of the curves (set by the MaxRecursion option). But still it has several disadvantages comparing to our algorithm - the user has to choose the value of PlotPoints and MaxRecursion options in advance. If the computation is too slow and he interrupts it, he gets nothing. On the other hand, if the computation finishes and he realize that the curves are not smooth enough, he starts another computation without a possibility to use the previous result.

Computation on the Web

The author collaborates on a web-project based on webMathematica by means of which any Internet user can apply the algorithm (using an HTML form) on several predefined types of boundary value problems. It is called CML (see http://webmath.zcu.cz) and it is primarily designed for scientific computations. It allows the users to perform predefined long-lasting computations on a cluster of computers located at our university using an HTML form (he obtains an e-mail with a link to the result when the computation is finished). Hence, our aim for the nearest future is to develop a parallel version of the algorithm and to include it in the CML portal.

Acknowledgment

This work has been supported by the Research Plan MSM 4977751301 and grant 1N04078 of the Ministry of Education, Youth and Sports of the Czech Republic.