8
$\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 2 days ago
  • $\begingroup$ @Alex Trounev Yes. $\endgroup$ Commented 2 days ago
  • $\begingroup$ @AlexeiBoulbitch Actually we can solve this problem using linear FEM without this small addition. $\endgroup$ Commented 2 days ago
  • $\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 2 days ago

2 Answers 2

Reset to default
7
$\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$
7
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.