Edit - Mathematica Stack Exchange

You are not logged in. Your edit will be placed in a queue until it is peer reviewed.

We welcome edits that make the post easier to understand and more valuable for readers. Because community members review edits, please try to make the post substantially better than how you found it, for example, by fixing grammar or adding additional resources and hyperlinks.

Required fields*

Rev

Title*

Body*

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][1]][1]
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][2]][2]

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][3]][3]

**Any idea?**


  [1]: https://i.sstatic.net/cwvPttKg.png
  [2]: https://i.sstatic.net/lGmbFxA9.png
  [3]: https://i.sstatic.net/XIvkjmjc.png