9
$\begingroup$

I am trying to solve a system of 3 nonlinear partial differential equations:

 eq1A = D[p1[y1, y2], {y1, 2}] + D[p1[y1, y2], {y2, 2}] - a*p1[y1, y2] + (y2 p3[y1, y2])/(y1^2 + y2^2 + 10^-6) - p1[y1, y2]^3 == 0; eq2A = D[p2[y1, y2], {y1, 2}] + D[p2[y1, y2], {y2, 2}] - a*p2[y1, y2] - (y1 p3[y1, y2])/(y1^2 + y2^2 + 10^-6) - p2[y1, y2]^3 == 0; eq3A = D[p3[y1, y2], {y1, 2}] + D[p3[y1, y2], {y2, 2}] - a*p3[y1, y2] + (y2 p1[y1, y2])/(y1^2 + y2^2 + 10^-6) - ( y1 p2[y1, y2])/(y1^2 + y2^2 + 10^-6) - p3[y1, y2]^3 == 0;

Here p1, p2 and p3 are the dependent variables, y1 and y2 are the independent ones, a is a parameter, a real number.

I whant to solve it within a disk:

Needs["NDSolve`FEM`"]; bm = ToBoundaryMesh[Disk[{0, 0}, 10]]; mesh1 = ToElementMesh[bm];

Zero Dirichlet conditions are imposed at the boundary:

dC1 = DirichletCondition[p1[y1, y2] == 0, {y1, y2} \[Element] bm]; dC2 = DirichletCondition[p2[y1, y2] == 0, {y1, y2} \[Element] bm]; dC3 = DirichletCondition[p3[y1, y2] == 0, {y1, y2} \[Element] bm];

I expect to see formation of a nontrivial solution p1(y1,y2)!=0, p2(y1,y2)!=0 and p3(y1,y2)!=0 at some small value of the parameter 0<a<1. I do not know, at what exactly a the nontrivial solution should appear, but a=0.1 or a=0.01 seem to be a plausible case.

The plain application of NDSolve

 a = 0.01; ndslA = NDSolve[{eq1A, eq2A, eq3A, dC1, dC2, dC3}, {p1, p2, p3}, {y1, y2} \[Element] mesh1][[1]]

brings nothing (as I expect by experience with such equations):

 Plot3D[{p1[y1, y2] /. ndslA, p3[y1, y2] /. ndslA}, {y1, y2} \[Element] mesh1, PlotStyle -> {Blue, {Orange, Opacity[0.5]}}, PlotRange -> All]

enter image description here Indeed, we see here only the trivial solution.

Further, I make a trial with the Initial Seeding:

ndslB = NDSolveValue[{eq1A, eq2A, eq3A, dC1, dC2, dC3, iC1}, {p1, p2, p3}, {y1, y2} \[Element] mesh1, InitialSeeding -> {p1[y1, y2] == 3, p2[y1, y2] == 2, p3[y1, y2] == 1}]

It returns the following warning: "NDSolveValue::femnlisl: An initial seed should be specified for each dependent variable (6).The number of initial seeds specified is 3".

I do not understand it, since I fixed the initial seeding for three dependent variables, and there are no more dependent ones.

By the way, the Help contains no examples of the initial seeining for a system of several equations. A good idea would be to add one, since the message looks as if I apply an incorrect syntax.

My next trial was to use the so-called, relaxation method. Within this method one instead of the original stationary system of equations formulates a time-dependent one. Its solution at large t must converge to the desired solution of the stationary equation. Here it is:

eq1B = D[p1[t, y1, y2], t] == D[p1[t, y1, y2], {y1, 2}] + D[p1[t, y1, y2], {y2, 2}] - a*p1[t, y1, y2] + ( y2 p3[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - p1[t, y1, y2]^3; eq2B = D[p2[t, y1, y2], t] == D[p2[t, y1, y2], {y1, 2}] + D[p2[t, y1, y2], {y2, 2}] - a*p2[t, y1, y2] - ( y1 p3[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - p2[t, y1, y2]^3; eq3B = D[p3[t, y1, y2], t] == D[p3[t, y1, y2], {y1, 2}] + D[p3[t, y1, y2], {y2, 2}] - a*p3[t, y1, y2] + ( y2 p1[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - (y1 p2[t, y1, y2])/( y1^2 + y2^2 + 10^-6) - p3[t, y1, y2]^3;

The Dirichlet conditions are the same as the ones used previously, but in this case, I also need the initial condition:

iC1 = p1[0, y1, y2] == p2[0, y1, y2] == p3[0, y1, y2] == (100 - (y1^2 + y2^2))*Exp[-(y1^2 + y2^2)/4];

There is an illusion that this works:

a = 0.1; ndslB = NDSolveValue[{eq1B, eq2B, eq3B, dC1, dC2, dC3, iC1}, {p1, p2, p3}, {t, 0, 50}, {y1, y2} \[Element] mesh1]

since it returns a solution. However, looking at this solution we see that it even does not satisfy the boundary conditions:

Plot3D[{ndslB[t, y1, y2][[1]] /. t -> 50, ndsl[t, y1, y2][[3]] /. t -> 50}, {y1, y2} \[Element] mesh1, PlotStyle -> {Blue, Orange, Red}, PlotRange -> All]

enter image description here

My final attempt was to apply the MethodOfLines together with the relaxation method:

ndslC= NDSolve[{eq1B, eq2B, eq3B, bCA1, bCA2, bCA3, ic1, ic2, ic3}, {p1, p2, p3}, {t, 0, 300}, {y1, -20, 20}, {y2, -20, 20}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 200}}]

with the boundary and initial conditions as follows:

bCA1 = p1[t, 20, y2] == p1[t, -20, y2] == p1[t, y1, -20] == p1[t, y1, 20] == 0; bCA2 = p2[t, 20, y2] == p2[t, -20, y2] == p2[t, y1, 20] == p2[t, y1, -20] == 0; bCA3 = p3[t, 20, y2] == p3[t, -20, y2] == p3[t, y1, 20] == p3[t, y1, -20] == 0; ic1 = p1[0, y1, y2] == (y1 + 20)*(y1 - 20)*(y2 + 20)*(y2 - 20)* Exp[-y1^2 - y2^2]; ic2 = p2[0, y1, y2] == (y1 + 20)*(y1 - 20)*(y2 + 20)*(y2 - 20)* Exp[-y1^2 - y2^2]; ic3 = p3[0, y1, y2] == (y1 + 20)*(y1 - 20)*(y2 + 20)*(y2 - 20)* Exp[-y1^2 - y2^2];

also with the result contradicting the boundary and initial conditions:

Plot3D[{sol[300, y1, y2][[1]], sol[300, y1, y2][[3]]}, {y1, -20, 20}, {y2, -20, 20}, PlotStyle -> {Blue, Orange}]

enter image description here

Any idea?

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ There's a (currently) undefined iC1 in ndslB, which looks like the culprit :) . Taking it away, the femnlisl warning doesn't show up, and a non-trivial solution is obtained: i.sstatic.net/WvceWhwX.png $\endgroup$ Commented Oct 20 at 1:35
  • 1
    $\begingroup$ It seems that the small addition in the denominator (10^-6) is done in order to eliminate the singularity. $\endgroup$ Commented Oct 20 at 4:02
  • $\begingroup$ @Alex Trounev Yes. $\endgroup$ Commented Oct 20 at 10:35
  • $\begingroup$ @AlexeiBoulbitch Actually we can solve this problem using linear FEM without this small addition. $\endgroup$ Commented Oct 20 at 12:37
  • $\begingroup$ @Alex Trounev: Thank you for the hint. I added this cut-off term to avoid an occasional infinity in the numeric process. The equation I treated is essentially nonlinear even in its minimal form I gave above. However, the statement of a linear problem is also of acute interest, scince I am looking first of all for the point of the first bifurcation of the nonlinear system. This point is given by the smallest discrete eigenvalue of the linear system. Therefore, any hint (with an example) of such a problem is of a high interest. $\endgroup$ Commented Oct 20 at 13:41

3 Answers 3

Reset to default
8
$\begingroup$

Using linear FEM and the false transient method (see here and here) we have

Needs["NDSolve`FEM`"]; L = 10; bm = ToBoundaryMesh[Disk[{0, 0}, L]]; mesh = ToElementMesh[bm, MeshRefinementFunction -> Function[{vertices, area}, area > 0.005 (0.1 + 10 Norm[Mean[vertices]])]]; mesh["Wireframe"]

Figure 1

a = 0.1; dt = 0.1; tmax = 50; P1[0][x_, y_] := 3 Exp[-x^2 - y^2]; P2[0][x_, y_] := 2 Exp[-x^2 - y^2]; P3[0][x_, y_] := Exp[-x^2 - y^2]; Do[eq1A = D[p1[y1, y2], {y1, 2}] + D[p1[y1, y2], {y2, 2}] - a*p1[y1, y2] + (y2 p3[y1, y2])/(y1^2 + y2^2) - P1[i][y1, y2]^2 p1[y1, y2] == (p1[y1, y2] - P1[i][y1, y2])/dt; eq2A = D[p2[y1, y2], {y1, 2}] + D[p2[y1, y2], {y2, 2}] - a*p2[y1, y2] - (y1 p3[y1, y2])/(y1^2 + y2^2) - P2[i][y1, y2]^2 p2[y1, y2] == (p2[y1, y2] - P2[i][y1, y2])/dt; eq3A = D[p3[y1, y2], {y1, 2}] + D[p3[y1, y2], {y2, 2}] - a*p3[y1, y2] + (y2 p1[y1, y2])/(y1^2 + y2^2) - (y1 p2[y1, y2])/(y1^2 + y2^2) - P3[i][y1, y2]^2 p3[y1, y2] == (p3[y1, y2] - P3[i][y1, y2])/dt; dC1 = DirichletCondition[p1[y1, y2] == 0, {y1, y2} \[Element] bm]; dC2 = DirichletCondition[p2[y1, y2] == 0, {y1, y2} \[Element] bm]; dC3 = DirichletCondition[p3[y1, y2] == 0, {y1, y2} \[Element] bm]; {P1[i + 1], P2[i + 1], P3[i + 1]} = NDSolveValue[{eq1A, eq2A, eq3A, dC1, dC2, dC3}, {p1, p2, p3}, {y1, y2} \[Element] mesh];, {i, 0, tmax}]

Visualization

Table[Plot3D[{P1[n][y1, y2], P3[n][y1, y2]}, {y1, y2} \[Element] mesh, PlotStyle -> {Blue, {Orange, Opacity[0.5]}}, PlotRange -> All, PlotTheme -> "Scientific"], {n, tmax - 1, tmax + 1}]

Figure 2

Asymptotic behavior at $x^2+y^2>>1 $

Table[Plot[{P1[tmax][r Cos[phi], r Sin[phi]], P3[tmax][r Cos[phi], r Sin[phi]]}, {r, 0, L}, PlotStyle -> {Blue, Red}, PlotRange -> All, PlotTheme -> "Scientific"], {phi, 0, 2 Pi, Pi/10}]

Figure 3

Update 1. The same solution we can evaluate without mesh generation (just to avoid boundary mesh usage) as follows

Needs["NDSolve`FEM`"]; L = 10; reg = Disk[{0, 0}, L]; bm = y1^2 + y2^2 == L^2; a = 0.1; dt = 0.1; tmax = 50; P1[0][x_, y_] := 3 Exp[-x^2 - y^2]; P2[0][x_, y_] := 2 Exp[-x^2 - y^2]; P3[0][x_, y_] := Exp[-x^2 - y^2]; Do[eq1A = D[p1[y1, y2], {y1, 2}] + D[p1[y1, y2], {y2, 2}] - a*p1[y1, y2] + (y2 p3[y1, y2])/(y1^2 + y2^2) - P1[i][y1, y2]^2 p1[y1, y2] == (p1[y1, y2] - P1[i][y1, y2])/dt; eq2A = D[p2[y1, y2], {y1, 2}] + D[p2[y1, y2], {y2, 2}] - a*p2[y1, y2] - (y1 p3[y1, y2])/(y1^2 + y2^2) - P2[i][y1, y2]^2 p2[y1, y2] == (p2[y1, y2] - P2[i][y1, y2])/dt; eq3A = D[p3[y1, y2], {y1, 2}] + D[p3[y1, y2], {y2, 2}] - a*p3[y1, y2] + (y2 p1[y1, y2])/(y1^2 + y2^2) - (y1 p2[y1, y2])/(y1^2 + y2^2) - P3[i][y1, y2]^2 p3[y1, y2] == (p3[y1, y2] - P3[i][y1, y2])/dt; dC1 = DirichletCondition[p1[y1, y2] == 0, bm]; dC2 = DirichletCondition[p2[y1, y2] == 0, bm]; dC3 = DirichletCondition[p3[y1, y2] == 0, bm]; {P1[i + 1], P2[i + 1], P3[i + 1]} = NDSolveValue[{eq1A, eq2A, eq3A, dC1, dC2, dC3}, {p1, p2, p3}, {y1, y2} \[Element] reg];, {i, 0, tmax}]

Visualization

Table[Plot3D[{P1[n][y1, y2], P3[n][y1, y2]}, {y1, y2} \[Element] reg, PlotStyle -> {Blue, {Orange, Opacity[0.5]}}, PlotRange -> All, PlotTheme -> "Scientific"], {n, tmax - 1, tmax + 1}]

Figure 4

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you, Alex. I will adopt this approach in my work. Would you happen to have a link to a clear and accessible text where this method is described in more detail? I asked an AI about it and found out that it's possible to adapt deltaT depending on the convergence of the process. Is it feasible to implement this adaptation in the scheme you provided? Once again, many thanks. Best regards, Alexei Boulbitch. $\endgroup$ Commented Oct 23 at 8:55
  • $\begingroup$ @AlexeiBoulbitch At first place this method was introduced in the paper G. D. Mallinson and G. de Vahl Davis, 'The method of the false transient for the solution of coupled partial differential equations', J. Comp. Phys., 12, 435 (1973). Then it has been implemented in the paper Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264. $\endgroup$ Commented Oct 23 at 11:38
8
$\begingroup$

The syntax you use for the initial seeding works just fine for me. What version do you have? Perhaps some stray variables were in the way?

I made a few other corrections as well:

Needs["NDSolve`FEM`"]; \[CapitalOmega] = Disk[{0, 0}, 10]; mesh1 = ToElementMesh[\[CapitalOmega]];

We use markers to specify the boundary - using the boundary mesh is a recipe for disaster. There is not any example of such usage in the documentation.

mesh1["BoundaryElementMarkerUnion"] (* {1} *) dC1 = DirichletCondition[p1[y1, y2] == 0, ElementMarker == 1]; dC2 = DirichletCondition[p2[y1, y2] == 0, ElementMarker == 1]; dC3 = DirichletCondition[p3[y1, y2] == 0, ElementMarker == 1]; a = 0.01; ndslA = NDSolve[{eq1A, eq2A, eq3A, dC1, dC2, dC3}, {p1, p2, p3}, {y1, y2} \[Element] mesh1, InitialSeeding -> {p1[y1, y2] == 3, p2[y1, y2] == 2, p3[y1, y2] == 1}][[1]] Plot3D[{p1[y1, y2] /. ndslA, p3[y1, y2] /. ndslA}, {y1, y2} \[Element] mesh1, PlotStyle -> {Blue, {Orange, Opacity[0.5]}}, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you. Using the boundary mesh indeed led to an incorrect process. However, once I replaced it with a reference to the boundary marker, my relaxation procedure began to function properly. Thank you once again. I remain unconvinced by the approach involving initial seeding, as it results in a solution that approaches the boundary at a finite angle—something clearly visible in the plot included in your response. $\endgroup$ Commented Oct 23 at 8:49
  • $\begingroup$ Continuation: Based on my experience with simpler equations of this type, as well as from analytical asymptotic analysis, I am confident that the true solution should asymptotically approach zero at the boundary. In other words, the nonzero distribution should be localized near the coordinate origin. At some distance it decay dramatically, and exponentially approaches zero as the distance further increases. At a distance of around 10 to 20 units, the solution should become practically indistinguishable from zero. $\endgroup$ Commented Oct 23 at 8:50
  • $\begingroup$ Continuation 2: I am not familiar with the numerical procedure used in the initial seeding method, but I suspect that, when applied to my equations, it does not yield a solution sufficiently close to the one I am aiming for. $\endgroup$ Commented Oct 23 at 8:51
1
$\begingroup$

I present my version of the solution, which combines ideas from @user21 and @Alex Trounev.

We begin with the mesh: it is coarse at the outer boundary—where the solution is expected to be close to zero—and finer near the centre:

Needs["NDSolve`FEM`"]; R = 40.; bM = ToBoundaryMesh[Disk[{0, 0}, R]]; hCenter=0.5; (* ~0.1…1 near r=0*) hBoundary=30; (* ~5…10 near r=R*) h[r_] := hCenter + (hBoundary - hCenter) (r/R)^1.5; Amax[r_] :=Sqrt[3]*(h[r]^2)/4; mrf = Function[{verts, area}, (*refine if the current area is above the radiusdependent target*)area > Amax[Norm[Mean[verts]]]]; mesh2 = ToElementMesh[bM, MeshRefinementFunction  mrf, MaxCellMeasure->(Sqrt[3]/4)*hBoundary^2, "MeshOrder"->1]; mesh2["Wireframe"]

enter image description here

The following system corresponds to the time-dependent system 𝐴 in the original question. Here, a is a numerical parameter taking real values 0<a<0.5.

eq1[a_] := D[p1[t, y1, y2], t] == D[p1[t, y1, y2], {y1, 2}] + D[p1[t, y1, y2], {y2, 2}] - a*p1[t, y1, y2] + (y2 p3[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - p1[t, y1, y2]^3; eq2[a_] := D[p2[t, y1, y2], t] == D[p2[t, y1, y2], {y1, 2}] + D[p2[t, y1, y2], {y2, 2}] - a*p2[t, y1, y2] - (y1 p3[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - p2[t, y1, y2]^3; eq3[a_] := D[p3[t, y1, y2], t] == D[p3[t, y1, y2], {y1, 2}] + D[p3[t, y1, y2], {y2, 2}] - a*p3[t, y1, y2] + (y2 p1[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - ( y1 p2[t, y1, y2])/(y1^2 + y2^2 + 10^-6) - p3[t, y1, y2]^3; dC1 = DirichletCondition[p1[t, y1, y2] == 0, ElementMarker == 1]; dC2 = DirichletCondition[p2[t, y1, y2] == 0, ElementMarker == 1]; dC3 = DirichletCondition[p3[t, y1, y2] == 0, ElementMarker == 1]; iC1A = p1[0, y1, y2] == 0.01*Exp[-(y1^2 + y2^2)/4]; iC2A = p2[0, y1, y2] == 0.02*Exp[-(y1^2 + y2^2)/4]; iC3A = p3[0, y1, y2] == 0.03 Exp[-(y1^2 + y2^2)/4];

The Dirichlet conditions were corrected following @user21's advice. I apply a relaxation method related to the false-transient approach recommended by @Alex Trounev. Unlike the linearised false-transient variant, this relaxation procedure is shorter in practice. On the same mesh it ran about 1.5 times faster in my tests, which motivated its use. More details appear in my article Pseudo-Dynamic Approach to the Numerical Solution of Nonlinear Stationary Partial Differential Equations, Mathematica Journal (2018).

The numerical solver is:

ndsl[a_, tmax_] := NDSolveValue[{eq1[a], eq2[a], eq3[a], dC1, dC2, dC3, iC1A, iC2A, iC3A}, {p1, p2, p3}, {t, 0, tmax}, {y1, y2} \[Element] mesh2];

To monitor convergence I compute the Hilbert norm of 𝑝3 as a function of simulation time:

lst1 = ParallelTable[{tmax, NIntegrate[(ndsl[0.05][[3]][tmax, y1, y2])^2, {y1, y2} \[Element] mesh2]^(1/2)}, {tmax, 0, 200, 10}];

The plot:

ListPlot[lst1, PlotStyle -> Red, ImageSize -> 300, AxesLabel -> {Style["\!\(\*SubscriptBox[\(t\), \(max\)]\)", 16, Italic, Black, FontFamily -> "Times"], Style["||\!\(\*SubscriptBox[\(p\), \(3\)]\)||", 16, Italic, Black, FontFamily -> "Times"]}]

enter image description here

indicates that convergence is satisfactory around tmax=100. To locate the bifurcation point — where computation time typically diverges — I performed a coarse sweep and found the transition near a≈0.22–0.24. Consequently I set tmax=5/Abs[a-0.222]for the refined sweep:

ndsl2[a_] := NDSolveValue[{eq1[a], eq2[a], eq3[a], dC1, dC2, dC3, iC1A, iC2A, iC3A}, {p1, p2, p3}, {t, 0, 5/ Abs[a - 0.222]}, {y1, y2} \[Element] mesh2];

I made a list of values of the parameter a:

aLst = Range[0.05, 0.28, 0.005];

and evaluated the solution over these values:

solA2 = ParallelMap[ndsl2[#] &, aLst];

and computed the Hilbert norm, ||p3|| for each a:

lstHilbertNorm = ParallelTable[{aLst[[m]], NIntegrate[(solA2[[m, 1]][5/Abs[aLst[[m]] - 0.222], y1, y2])^2, {y1, y2} \[Element] mesh2]^(1/2)}, {m, 1, Length[aLst]}];

The resulting plot

enter image description here

shows that the bifurcation lies between a=0.235and a=0.240.

Finally, I display the solution. For example, let us look at p1[y1,y2] at some a value:

enter image description here

The solution is localised near the origin and decays exponentially toward the outer boundary.

I would like to thank @user21 and @Alex Tournev for the discussion that helped me so much. Thank you!

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.