Initialization

Button[Style["Equations initialization",White,Bold,24],
Clear[x];

sz=43;

Ac[expr_]:=
DynamicModule[{c1=0},
               EventHandler[     
               Dynamic[ 
                                   If[c1==0, 
                                      Style[expr,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"]//ExpressionCell,
                                                Dynamic@If[Clock[1,1.0,4]<.5,
                                                Style[expr,RGBColor[0.8,0,0],Plain,sz,Italic,FontFamily->"TimesNewRoman"]//ExpressionCell,
                               Style[expr,Black,Plain,sz,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell
                                       ]
                             ],                                                                              
{"MouseDown":>(c1=c1/.{0->1,1->0})}
                                      ]
                 ];

AcBold[expr_]:=
DynamicModule[{c1=0},
               EventHandler[     
               Dynamic[ 
                                   If[c1==0, 
                                      Style[expr,Black, Bold,sz,Italic,FontFamily->"TimesNewRoman"]//ExpressionCell,
                                                Dynamic@If[Clock[1,1.0,4]<.5,
                                                Style[expr,RGBColor[0.8,0,0],Bold,sz,Italic,FontFamily->"TimesNewRoman"]//ExpressionCell,
                               Style[expr,Black,Bold,sz,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell
                                       ]
                             ],                                                                              
{"MouseDown":>(c1=c1/.{0->1,1->0})}
                                      ]
                 ];

AcSubsBold[expr1_,expr2_]:=
DynamicModule[{c1=0},
               EventHandler[     
               Dynamic[ 
                                   If[c1==0, 
                                        Subscript[ Style[expr1,Black, Bold,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-10,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                                Dynamic@If[Clock[1,1.0,4]<.5,
                                                  Subscript[ Style[expr1,RGBColor[0.8,0,0], Bold,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,RGBColor[0.8,0,0], Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                 Subscript[ Style[expr1,Black, Bold,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]]//ExpressionCell
                                       ]
                             ],                                                                              
{"MouseDown":>(c1=c1/.{0->1,1->0})}
                                      ]
                 ];

AcSubs[expr1_,expr2_]:=
DynamicModule[{c1=0},
               EventHandler[     
               Dynamic[ 
                                   If[c1==0, 
                                        Subscript[ Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                                Dynamic@If[Clock[1,1.0,4]<.5,
                                                  Subscript[ Style[expr1,RGBColor[0.8,0,0], Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,RGBColor[0.8,0,0], Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                 Subscript[ Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-10,Italic,FontFamily->"TimesNewRoman"]]]//ExpressionCell
                                       ]
                             ],                                                                              
{"MouseDown":>(c1=c1/.{0->1,1->0})}
                                      ]
                 ];

AcSuper[expr1_,expr2_]:=
DynamicModule[{c1=0},
               EventHandler[     
               Dynamic[ 
                                   If[c1==0, 
                                        Superscript[ Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                                Dynamic@If[Clock[1,1.0,4]<.5,
                                                  Superscript[ Style[expr1,RGBColor[0.8,0,0], Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,RGBColor[0.8,0,0], Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]//ExpressionCell,
                                 Superscript[ Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"], Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]]//ExpressionCell
                                       ]
                             ],                                                                              
{"MouseDown":>(c1=c1/.{0->1,1->0})}
                                      ]
                 ];



Pl[expr_]:=ExpressionCell[Style[expr,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"]];

PlBold[expr_]:=ExpressionCell[Style[expr,Black, Bold,sz,Italic,FontFamily->"TimesNewRoman"]];

PlSuper[expr1_,expr2_]:=ExpressionCell[Superscript[Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"],Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]];

PlSubs[expr1_,expr2_]:=ExpressionCell[Subscript[Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"],Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]];

PlLarge[expr_]:=ExpressionCell[Style[expr,Black, FontSlant->"ExtraLight",2.5sz,FontFamily->"TimesNewRoman"]];

PlSubSuper[expr1_,expr2_,expr3_]:=ExpressionCell[Subsuperscript[Style[expr1,Black, Plain,sz,Italic,FontFamily->"TimesNewRoman"],Style[expr2,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"],Style[expr3,Black, Plain,sz-15,Italic,FontFamily->"TimesNewRoman"]]],

Background ->RGBColor[1.,0.51,0], 

BaseStyle -> {"GenericButton", Bold},ImageSize->{400,60}]

Button[Style["Evaluate Initialization Cells",White,Bold,24],FrontEndExecute[FrontEndToken["EvaluateInitialization"]],Background ->RGBColor[1.,0.51,0], 
BaseStyle -> {"GenericButton", Bold},ImageSize->{400,60}]

Numerical solution of nonlinear PDEs exhibiting soft bifurcations

Alexei Boubitch, 
Dr. habil., 
Researcher, 
IEE Luxembourg

Join the Conversation #WolframTechConf

ABSTRACT

This talk communicates a numerical pseudo-dynamic approach to solve nonlinear stationary equations exhibiting bifurcations. The solution is achieved by passing from the stationary partial differential equation - to a pseudo-time dependent one. The latter is constructed in such a way that the desired nontrivial solution of the stationary equation represents its fixed point. The numeric solution of the stationary equation which is being looked for is then obtained as the solution of the pseudo-time dependent equation at a high enough value of the pseudo-time.

Scope:

1. I will first tell a few words recalling the theory of soft bifurcations

2. Then formulate the idea of the approach...

3. ...and Illustrate its use by a simple example.

4. Eventual discussion of peculiarities of hard bifurcations.

Introduction: The bifurcation theory revisited

A nonlinear PDE in a most general form

Let us study an equation (or a system of equations):

FlipView[{
Column[{
Row[{"F"//AcBold,"["//Pl,"u"//AcBold,"(x),"//Pl,"λ"//Ac,"]=0"//Pl}],
Spacer[5],
Row[{"u"//AcBold,"(x)(Subscript[|, ∂Ω])=0"//Pl}]
}],

Column[{
Row[{"L"//AcBold,"["//Pl,"u"//AcBold,"(x),"//Pl,"λ"//Ac,"]+"//Pl,"N[u(x)]"//AcBold,"=0"//Pl}],
Spacer[5],
Row[{"u"//AcBold,"(x)(Subscript[|, ∂Ω])=0"//Pl}]
}]

}]

Such that

Row[{"F"//PlBold,"(0, λ)≡0"//Pl}]

A bifurcation

At λ>(Subscript[α, 0])

Only a trivial solution u(x)≡0 of the equation F[u(x),λ]=0 exists

At λ≤(Subscript[α, 0])

One or several nontrivial solutions u(x)≠0 of the equation F[u(x),λ]=0 emerge.

A supercritical (soft) and a subcritical (hard) bifurcation: just to recall

A soft is a bifurcation taking place continuously in the bifurcation point. A hard bifurcation takes place jumpwise.

Row[{
Plot[If[λ<1,(Sqrt[1-λ]),0],{λ,0,2},Exclusions->None,PlotTheme->"Classic",PlotRange->{0,1.1},PlotStyle->{Blue,Thickness[0.005]},AxesLabel->{Style["λ",36,Black],Row[{Style["||",36,Black],Style["u",36,Black,Italic],Style["||",36,Black]}]},Ticks->{{{1,Style["(Subscript[α, 0])",30]}},None},AxesStyle->Arrowheads[0.05],ImageSize->550,
Epilog->Inset[Style["Soft",30],Scaled[{0.5,0.8}]]
],
Spacer[60],
Show[{
Plot[If[λ<1,(Sqrt[1.2-λ]),0],{λ,0,2},PlotTheme->"Classic",PlotRange->{0,1.2},PlotStyle->{Blue,Thickness[0.005]},AxesLabel->{Style["λ",36,Black],Row[{Style["||",36,Black],Style["u",36,Black,Italic],Style["||",36,Black]}]},Ticks->{{{1,Style["(Subscript[α, 0])",30]}},None},AxesStyle->Arrowheads[0.05],ImageSize->550,
Epilog->Inset[Style["Hard",30],Scaled[{0.5,0.8}]]],
Graphics[{Blue,Dashed,Line[{Scaled[{0.502,0}],Scaled[{0.502,0.358}]}]}]
}]
}]

Here I only address soft bifurcations. Hard ones can be briefly discussed in the end.

A problem often arising during solving such equations numerically

A direct numerical solution of the above equation with a nonlinear solver often represents a problem: the solver only returns the trivial solution, u(x)≡0.

Possible explanation: the solver starts from the boundary, where, the BC, (Subscript[u, s])(Subscript[|, ∂Ω])=0,  already belongs to the trivial solution. The solver, thus, “sits” on the true solution, u=0, and is unable to “stand up” from it.

A pseudo-dynamic approach

Instead of the stationary equation F(u,λ)=0 let us study the behavior of the pseudo-time dependent equation.

Column[{
Row[{(1/(Pl["∂t"]))Row[{"∂"//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t)"//Pl}],"="//Pl,AcBold["F"],"[ "//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t), λ]"//Pl,"    ⟺   "//Pl,AcSubsBold["u(x,t)→u","s"],"(x),   t→∞"//Pl}],
Spacer[5],
Row[{"u"//AcBold,"(x)(Subscript[|, ∂Ω])=0"//Pl}],
Spacer[5],
Row[{AcBold["u"],"("//Pl,AcBold["x"],",0)="//Pl,AcSubsBold[u,0],"("//Pl,AcBold["x"],")≢0"//Pl}]
}]

Example

A 1D Ginzburg-Landau equation

Column[{
Row[{(1/(PlSuper["∂x","2"]))Row[{AcSuper["∂","2"],"u(x)"//Ac}],"=2["//Pl,"λ"//Ac,"+"//Pl,Ac["U(x)"],"]"," u(x)"//Ac,"+ "//Pl,AcSuper["u","3"],"(x)"//Pl}],
Spacer[5],
Row[{Ac["u(±∞)=0"]}]
}]

Plot[-2((Cosh[x])^(-2)),{x,-10,10},PlotTheme->"Classic",PlotRange->{-2.1,0.5},PlotStyle->{Red,Thickness[0.005]},AxesLabel -> {Style["x", 36, Italic,Black], Style["U(x)=-(1/((cosh^2) (x)))", 36, Italic,Black]},AxesStyle->Arrowheads[0.05],ImageSize->650]

What the theory of bifurcations gives?

The bifurcation point

(Subscript[α, 0])=(1/2)

The branching off solution

u(x)=({{{, ({{0,, λ>1/2}, {±(((Sqrt[3])(((1/2-λ))^(1/2)))/(cosh(x))),, λ≤1/2}})}})

Row[{
Plot[(1/(Cosh[x])),{x,-10,10},PlotStyle->{Red,Thickness[0.005]},PlotRange->{-0.05,1.2},Ticks->None,AxesLabel -> {Style["x", 36, Italic,Black,FontFamily->"TimesNewRoman"], Style["u(x) @λ<0.5", 36,Black,Italic,FontFamily->"TimesNewRoman"]},AxesStyle->Arrowheads[0.05],ImageSize->550],
Spacer[60],
Plot[If[λ<0.5,(Sqrt[0.5-λ]),0],{λ,0.3,0.6},Exclusions->None,PlotTheme->"Classic",PlotRange->{0,0.5},PlotStyle->{Blue,Thickness[0.005]},AxesLabel->{Style["λ",36,Black],Row[{Style["||",36,Black],Style["u",36,Black,Italic],Style["||",36,Black]}]},Ticks->{{{0.5,Style["0.5",30]}},None},AxesStyle->Arrowheads[0.05],ImageSize->550],

}]

Pseudo-time dependent 1D Ginzburg-Landau equation...

Column[{
Row[{(1/(Ac["∂t"]))Row[{"∂u(x,t)"//Ac}],"="//Pl,(1/(AcSuper["∂x","2"]))Row[{AcSuper["∂","2"],"u(x,t)"//Ac}],"-2[λ - "//Pl,(("1"//Pl)/(Row[{AcSuper["cosh","2"],"(x)"//Pl}])),"]"//Pl," u(x,t)- "//Ac,AcSuper["u","3"],"(x,t)"//Pl}],
Spacer[5],
Row[{Ac["u("],Ac["±L"]," t)=0"//Pl}],
Spacer[5],
Row[{Ac["u(x,0)"],"="//Pl,AcSubs["u","0"],"(x)≠0"//Pl}]
}]

... and its numerical solution

eq[λ_]:=(∂_t) u[x, t]==(∂_({x,2})) u[x, t]- 2(λ  - (1/((Cosh[x])^2)))*u[x, t] - ((u[x, t])^3);

s[λ_,tmax_,L_,ic_]:=NDSolve[{eq[λ], u[-L, t] == 0, u[L, t] == 0, ic}, u[x, t], {x, -L,L}, {t, 0,tmax}, Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MinPoints"->2 L}}];

Evolving in time of the solution: λ>(Subscript[α, 0])≈0.5

Row[{AcSubs["u","0"],"(x)=0.2 exp(-"//Pl,PlSuper["x","2"],")"//Pl}]

Key code element

Evaluate[u[x, t]/.s[0.7,50,10, u[x,0]==0.2*Exp[-(x^2)]]/.t->t1]

u=u(x,t) at  λ=0.7

ic1=u[x,0]==0.2*Exp[-(x^2)];

Animate[
Plot[

Evaluate[u[x, t]/.s[0.7,50,10,ic1]/.t->t1], 

{x, -10, 10},
PlotTheme->"Classic", PlotRange->{{-10.2,11},{0,0.25}},PlotStyle->{Blue,Thick},AxesLabel -> {Style["x", 36,Black, Italic], Style["u(x,t)",36,Black, Italic]},TicksStyle->Directive[18,Black],AxesStyle->Arrowheads[0.03],ImageSize->600],{{t1,0,Dynamic@Panel[Row[{Style["(Subscript[t, 1])=",Blue,30,Italic],Style[NumberForm[t1,{3,2}],Blue,30,Italic]}],ImageSize->{180,60}]},0,15,0.01},AnimationRepetitions->1,AnimationRate->1.5,AnimationRunning->False,ContentSize->620]

Evolving in time of the solution: λ<(Subscript[α, 0])≈0.5

Row[{AcSubs["u","0"],"(x)=0.2 exp(-"//Pl,PlSuper["x","2"],")"//Pl}]

Key code element

Evaluate[u[x, t]/.s[0.3,50,10, u[x,0]==0.2*Exp[-(x^2)]]/.t->t1]

u=u(x,t) at  λ=0.3

ic1=u[x,0]==0.2*Exp[-(x^2)];

Animate[
Plot[{

Evaluate[u[x, t]/.s[0.3,50,10,ic1]/.t->t1],Evaluate[u[x, t]/.s[0.3,250,10,ic1]/.t->20]}, 

{x, -10, 10},
PlotTheme->"Classic", PlotRange->{{-10.2,11},{0,0.9}},PlotStyle->{{Blue,Thick}, {Dashed,Red}},AxesLabel -> {Style["x", 36,Black, Italic], Style["u(x,t)",36,Black, Italic]},TicksStyle->Directive[18,Black],AxesStyle->Arrowheads[0.03],ImageSize->600],{{t1,0,Dynamic@Panel[Row[{Style["(Subscript[t, 1])=",Blue,30,Italic],Style[NumberForm[t1,{3,2}],Blue,30,Italic]}],ImageSize->{180,60}]},0,10,0.01},AnimationRepetitions->1,AnimationRate->1.5,AnimationRunning->False,ContentSize->620]

... with an asymmetric initial condition: λ<0.5

Row[{AcSubs["u","0"],"(x)=0.3×"//Pl,PlSuper["10","-4"]," (x+10)"//Pl,PlSuper["(10-x)","3"]}]

Key code element

Evaluate[u[x, t]/.s[0.3,50,10, u[x,0]==0.3*(10^(-4))*(x+10)(((10-x))^3)]/.t->t1]

u=u(x,t) at a fixed λ=0.3

ic3=u[x,0]==0.3*(10^(-4))*(x+10)(((10-x))^3);

Animate[
Plot[

{Evaluate[u[x, t]/.s[0.3,50,10,ic3]/.t->t1],Evaluate[u[x, t]/.s[0.3,250,10,ic1]/.t->25]},

 {x, -10, 10}, 
  PlotTheme->"Classic",PlotStyle->{{Blue,Thick}, {Dashed,Red}},AxesLabel -> {Style["x", 36, Italic,Black], Style["u(x,t)", 36, Italic,Black]},AxesStyle->Arrowheads[0.03],PlotRange->{{-10.2,11.2},{0,0.9}},ImageSize->600],{{t1,0,Dynamic@Panel[Row[{Style["t=",Blue,30,Italic],Style[NumberForm[t1,{3,2}],Blue,30,Italic]}],ImageSize->{180,60}]},0,7,0.01},AnimationRepetitions->1,AnimationRate->1.5,AnimationRunning->False,ContentSize->620]

The dependence on λ

Key code element

Evaluate[u[x, t]/.s[λ,150,10,u[x,0]==0.2*Exp[-(x^2)]]/.t->150]

u=u(x,λ)

Manipulate[
col=If[λ<0.5,Red,Blue];
Plot[

Evaluate[u[x, t]/.s[λ,150,10,ic1]/.t->150], 

{x, -10, 10}, ImageSize->600,
  PlotTheme->"Classic",PlotRange->{-0.01,0.9},PlotStyle->{col,Thick}, AxesLabel -> {Style["x", 36, Italic], Style["u(x)",col, 36, Italic]},AxesStyle->Arrowheads[0.03],PlotPoints->ControlActive[9,25]],
{{λ,0.7,Dynamic@Panel[Row[{Style["λ=",30,Italic,col],Style[NumberForm[λ,{3,2}],col,30,Italic]}],ImageSize->{120,60}]},0.3,0.7},SaveDefinitions->True,ContentSize->620]

Conclusion

The solution emerges as λ comes close to 0.5 and then gradually increases

Solution norm and the convergence control

The proximity of the solution to the fixed point can be estimated using the Hilbert norm, ||u||≡η(t,λ):

Row[{AcSuper["η","2"],"(t,λ)="//Pl,,(("1"//Pl)/(Pl["2L"])),"(∫_(-L)^L)"//Pl,PlSuper["u","2"],"(x,t)\[DifferentialD]x"//Pl}]

Key code element

Table[
{tmax,(((NIntegrate[((u[x, t])^2)/.s[0.48,tmax,10,ic1]/.t->tmax,{x,-10,10}]/20))^(1/2))},
{tmax,10,1000,10}]

Outcome

lstEtaA=Table[{tmax,(((NIntegrate[((u[x, t])^2)/.s[0.48,tmax,10,ic1]/.t->tmax,{x,-10,10}]/20))^(1/2))},{tmax,10,1000,10}]/.{x_,{y_}}->{x,y};
ListLogLogPlot[lstEtaA, PlotStyle->{Blue,PointSize[0.015]},PlotTheme->"Classic",TicksStyle->{{Black,16},{Black,16}},TicksStyle->{16,Black},AxesLabel -> {Style["(Subscript[t, max])", 36, Italic,Black], Style["η((Subscript[t, max]),0.48)",36, Italic,Black]},AxesStyle->Arrowheads[0.03],ImageSize->650]

In search of the bifurcation point

According to the general theory one expects the behavior:

η(λ)={({{0,                          λ>(Subscript[α, 0])}, {(Subscript[η, 0]) ((((Subscript[α, 0])-λ))^(1/2)),  λ≤(Subscript[α, 0])}})

FlipView[{
Row[{AcSubs["t","max"],"=",((AcSubs["t","max0"])/(Ac["|λ-(Subscript[α, 0])|"]))}],
Row[{AcSubs["t","max"],"=",((AcSubs["t","max0"])/("[(λ-(Subscript[α, 0]))^2+ϵ]^(1/2)"//Pl))}]
}]

Key code elements

ϵ=(10^(-6));   tmax[λ_]:=(8/((((((0.5-λ))^2)+ϵ))^(1/2)));

lstEta=Table[
{λ,(((NIntegrate[((u[x, t])^2)/.s[λ, tmax[λ], L , ic1]/.t->tmax[λ],{x,-L,L}]/(2*L)))^(1/2))}, 
{λ,0.45,0.52,0.001}]/.{x_,{y_}}->{x,y};

FindFit[lstEta, If[λ>=α0, 0, η0*((((((α0-λ))^2)+ϵ))^(1/4))],{η0, {α0,0.49}},λ]

Outcome

Clear[λ];
ϵ=(10^(-6));
tmax[λ_]:=(8/((((((0.5-λ))^2)+ϵ))^(1/2)));
L=10;

lstEta=Table[{λ,(((NIntegrate[((u[x, t])^2)/.s[λ,tmax[λ],L,ic1]/.t->tmax[λ],{x,-L,L}]/(2*L)))^(1/2))},{λ,0.45,0.52,0.001}]/.{x_,{y_}}->{x,y};

fit=FindFit[lstEta, If[λ>=α0,0,η0*((((((α0-λ))^2)+ϵ))^(1/4))],{η0,{α0,0.49}},λ]

Show[{
ListPlot[lstEta, PlotStyle->{Red,PointSize[0.015]},PlotTheme->"Classic", PlotRange->{0,0.14},AxesLabel -> {Style["λ", 34, Italic,Black], Style["η((Subscript[t, max]), λ)",34, Italic,Black]},AxesStyle->Arrowheads[0.03],
PlotLegends->Placed[PointLegend[{Style["   Numeric solution",26,Italic,"TimesNewRoman",Black]}],Scaled[{0.77,0.77}]]
],

Plot[If[λ>α0,0,η0*(((α0-λ))^0.5)]/.fit,{λ,0.46,0.54},PlotStyle->{Blue,Thick},PlotTheme->"Classic",PlotLegends->Placed[LineLegend[{Row[{Style[" Fitting: (Subscript[α, 0])≈",24,Italic,Black],Style[NumberForm[fit[[2,2]],{5,4}],24,Italic,Black]}]}],Scaled[{0.75,0.6}]]]
},ImageSize->650]

Summary

I shared an approach enabling one to solve nonlinear stationary partial differential equations numerically.

The approach is based on solving a pseudo-time dependent equation instead of the stationary one, the initial condition being essentially different from zero. The latter enables the solver to avoid "adhering" to the trivial solution and helps it to pass to the nontrivial one.

The penalty is, however, an increased simulation time.

Take home:

How to operate?

Step 1:  Go from  F(u,λ)=0  to   ∂u/∂t=F(u,λ)

Step 2: Solve it with the standard routine NDSolve using the MethodOfLines

Step 3: Use an appropriate norm to control convergence, and find a good (Subscript[t, max])

Step 4: Keep in mind the that in the vicinity of bifurcation point (if any) the critical slowing down occurs.

Discussion

What sign in front of F(u,λ) to choose: plus or minus?

One can take either

Row[{(1/(Pl["∂t"]))Row[{"∂"//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t)"//Pl}],"="//Pl,AcBold["F"],"[ "//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t), λ]"//Pl}]

or

Row[{(1/(Pl["∂t"]))Row[{"∂"//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t)"//Pl}],"="//Pl,Ac["-"],AcBold["F"],"[ "//Pl,AcBold["u"],"("//Pl,AcBold["x"],", t), λ]"//Pl}]

In one of these two cases the solution converges

Row[{AcSubsBold["u(x,t)→u","s"],"(x),   t→∞"//Pl}]

in the other it diverges

Row[{AcBold["u"],"("//Pl,AcBold["x"],",t)→∞,"//Pl,"   t→∞"//Pl}]

Bya trial and error one choses a right sign.

Nonzero boundary conditions

The approach can be applied also in the case of nonzero boundary conditions.They represent no obstacle.

Dimensionality

Application of the above approach is not limited by the space dimensionality.

Eventual further discussion

If an “energy” functional exists...

Some stationary PDE admit a generating functional, E{u,λ}, (“the energy”) such that the PDE brings it a minimum:

FlipView[{
Row[{((Ac["δ E{u, λ}"])/(Ac["δ u"]))," = 0"//Pl}],
Row[{((Ac["δ E{u, λ}"])/(Ac["δ u"]))," = 0"//Pl,"   ⟺   "//Pl, AcBold["F"],"("//Pl,AcBold["u"],",λ)=0"//Pl}]
}]

The energy functional -> into a function

Substitution u=(Subscript[ξ, 0])(Subscript[ψ, 0])(x)+(Subscript[ξ, 1])(Subscript[ψ, 1])(x)+... 
transforms the energy functional into a function of (Subscript[ξ, 0]), (Subscript[ξ, 1]),...

Row[{"E{u, λ}"//Ac,"="//Pl,Ac["E"],"("//Pl,AcSubs["ξ","0"],","//Pl, AcSubs["ξ","1"],",..., λ)"//Pl}]

Illustration

thermo[x_,label_String,labelSize_]:=Show[{
Graphics3D[{EdgeForm[Gray],FaceForm[{White,Blue,Opacity[0.05]}],Cuboid[{-0.25,0.2,-0.6},{0.25,0.23,0.55}]},ViewPoint->{0.05,-1.5,0.},Boxed->False,ImageSize->{200,400}],
Graphics3D[{Text[Style["(Subscript[λ, c])",16],{-0.15,0,0}],
Text[Style[label,labelSize,Bold],{0,0,0.45}]},ViewPoint->{0,-1.5,0},Boxed->False,ImageSize->{200,400}],
Graphics3D[{LightBlue, Specularity[White,20],Opacity[0.5],Sphere[{0,0,-0.4},0.11]}],
Graphics3D[{Red, Specularity[White,20],Sphere[{0,0,-0.4},0.1]}],
Graphics3D[{LightBlue,Opacity[0.5],Tube[{{0,0,-0.4},{0,0,0.35}},0.05]}],
Graphics3D[{Red,Specularity[White,20],Cylinder[{{0,0,-0.45},{0,0,x}},0.04]}],

Graphics3D[{Thickness[0.015],Line[{{-0.07,0,#},{-0.09,0,#}}]&/@Table[i,{i,-0.3,0.3,0.1}]}]

}];



Manipulate[

Row[{
thermo[-p,"λ",20],


Animate[

Show[{
Plot[(((-p))/2)*(x^2)+((x^4)/4),{x,-1.2,1.2},AxesLabel->{Style["ξ",30,Black],Style["E(ξ,λ)",30,Italic,Black]},PlotRange->{-0.03,0.07}, PlotStyle->{Blue, Thick},Ticks->None, PlotPoints->30,ImageSize->500,AxesStyle->Arrowheads[0.04],PlotPoints->ControlActive[3,10]],
Graphics[{Red,PointSize[0.015],Point[{y,-(p/2)*(y^2)+((y^4)/4)}]}]
   }],{y,FindRoot[-p*(y^2)/2+(y^4)/4==0.05,{y,1}][[1,2]],FindMinimum[-(p/2)*(y^2)+((y^4)/4),{y,0.5}][[2,1,2]]},AnimationRunning->False,AnimationRepetitions->1,Paneled->False]


}],


{{p,-0.3},-0.3,0.3},SaveDefinitions->True]

Soft/hard bifurcation and the energy

In hard bifurcations, the energy may have several minimums in which the equation F(u,λ)=0 holds.

thermo[x_,label_String,labelSize_]:=Show[{
Graphics3D[{EdgeForm[Gray],FaceForm[{White,Blue,Opacity[0.05]}],Cuboid[{-0.25,0.2,-0.6},{0.25,0.23,0.55}]},ViewPoint->{0.05,-1.5,0.},Boxed->False,ImageSize->{200,400}],
Graphics3D[{Text[Style["(Subscript[λ, c])",16],{-0.15,0,0}],
Text[Style[label,labelSize,Bold],{0,0,0.45}]},ViewPoint->{0,-1.5,0},Boxed->False,ImageSize->{200,400}],
Graphics3D[{LightBlue, Specularity[White,20],Opacity[0.5],Sphere[{0,0,-0.4},0.11]}],
Graphics3D[{Red, Specularity[White,20],Sphere[{0,0,-0.4},0.1]}],
Graphics3D[{LightBlue,Opacity[0.5],Tube[{{0,0,-0.4},{0,0,0.35}},0.05]}],
Graphics3D[{Red,Specularity[White,20],Cylinder[{{0,0,-0.45},{0,0,x}},0.04]}],

Graphics3D[{Thickness[0.015],Line[{{-0.07,0,#},{-0.09,0,#}}]&/@Table[i,{i,-0.3,0.3,0.1}]}]

}];



Manipulate[

Row[{
thermo[λ,"λ",20],


Row[{

Show[{
Plot[(((λ))/2)*(x^2)+((x^4)/4),{x,-1.2,1.2},AxesLabel->{Style["ξ",30,Black],Style["E(ξ,λ)",30,Italic,Black]},PlotRange->{-0.03,0.07}, PlotStyle->{Blue, Thick},Ticks->None, PlotPoints->30,ImageSize->500,AxesStyle->Arrowheads[0.04],PlotPoints->ControlActive[3,10]],
Graphics[{Text[Style["Soft",36],Scaled[{0.9,0.95}]]}]
}],

Spacer[50],
Show[{
Plot[(((λ+0.35))/2)*(x^2)-((x^4)/4)+((x^6)/6),{x,-1.2,1.2},AxesLabel->{Style["ξ",30,Black],Style["E(ξ,λ)",30,Italic,Black]},PlotRange->{-0.03,0.07}, PlotStyle->{Blue, Thick},Ticks->None, PlotPoints->30,ImageSize->500,AxesStyle->Arrowheads[0.04],PlotPoints->ControlActive[3,10]],
Graphics[{Text[Style["Hard",36],Scaled[{0.9,0.95}]]}]
}]

}]


}],


{{λ,0.3},-0.3,0.3},SaveDefinitions->True]

Show[{
Plot[(1/2)0.1945*(u^2)-(1/4)*(u^4)+(1/6)*(u^6),{u,-1.,1},PlotTheme->"Classic",Ticks->None,PlotStyle->{Blue,Thickness[0.005]},PlotRange->{0,0.015},AxesLabel->{Style["ξ",36,Black],Style["E",36,Black,Italic]},AxesStyle->Arrowheads[0.05]],
Plot[(1/2)0.1945*(u^2)-(1/4)*(u^4)+(1/6)*(u^6),{u,0.65,0.8577},PlotStyle->Red],
Plot[(1/2)0.1945*(u^2)-(1/4)*(u^4)+(1/6)*(u^6),{u,0,0.4},PlotStyle->RGBColor[0,0.5,0]],

Graphics[{Red,Arrowheads[0.039],Arrow[{{0.80988,0.0032308},{0.8577,0.0026}}],PointSize[0.015],Point[{0.65,0.009}],RGBColor[0,0.5,0], Arrowheads[0.039],Arrow[{{0.00949,0.000066},{0,0}}],Point[{0.4,0.0098}]}]

},ImagePadding->45,ImageSize->450]

Where does the point run?

Whether the solutions corresponding to different minimums belong to the complete solution set of the problem is decided basing on the problem origin. If these solutions are parasitic, the workaround is to try to find initial conditions belonging to the basin of attraction of the desired minimum, that is, as close to the expected solution as possible.