How can I solve a fourth order system of PDEs?

Question

My equations are:

$A_t=r A-\mu |A|^2A-\nu |B|^2A+(d/dx+\frac{e}{2i}(\partial^2/\partial_{xx}+\partial^2/\partial_{yy}))^2A$ $B_t=r B-\mu |B|^2B-\nu |A|^2B+(d/dy+\frac{e}{2i}(\partial^2/\partial_{xx}+\partial^2/\partial_{yy}))^2B$

where $r=40.0$, $\mu=2,\nu=1,e=0.1,i=\sqrt{-1}$

The initial conditions are:

$A_0=\sqrt{\frac{r-q^2(1+eq/2)^2}{\mu+\nu}}\exp(iqx)$

$B_0=\sqrt{\frac{r-q^2(1+eq/2)^2}{\mu+\nu}}\exp(iqy)$

The period condition is used here.

eq1 = 
  D[A[x, y, t], t] == 
    r A[x, y, t] - μ Abs[A[x, y, t]]^2 A[x, y, t] - 
    ν Abs[B[x, y, t]]^2 A[x, y, t] + D[A[x, y, t], {x, 2}] + 
    ϵ/I (D[A[x, y, t], {x, 3}] + D[D[A[x, y, t], {x, 1}], {y, 2}]) + 
    (ϵ/(2 I))^2 (D[A[x, y, t], {x, 4}] + D[A[x, y, t], {y, 4}]);

 eq2 = 
   D[B[x, y, t], t] == 
     r B[x, y, t] - μ Abs[B[x, y, t]]^2 B[x, y, t] - 
     ν Abs[A[x, y, t]]^2 B[x, y, t] + D[B[x, y, t], {y, 2}] + 
     ϵ/I (D[B[x, y, t], {y, 3}] + D[D[B[x, y, t], {y, 1}], {x, 2}]) + 
     (ϵ/(2 I))^2 (D[B[x, y, t], {x, 4}] + D[B[x, y, t], {y, 4}]);

q = -2; ϵ = 0.1; r = 40.0; μ = 2.0; ν = 1.0; L = 20;
R0 = Sqrt[(r - q^2 (1 + ϵ q/2)^2)/(μ + ν)];
As = R0 Exp[I q x]; Bs = R0 Exp[I q y]; 

sol = 
  NDSolve[
    {eq1, eq2, 
     A[x, y, 0] == As, B[x, y, 0] == Bs, 
     A[x, -L/2, t] == A[x, L/2, t], A[-L/2, y, t] == A[L/2, y, t], 
     B[x, -L/2, t] == B[x, L/2, t], B[-L/2, y, t] == B[L/2, y, t]}, 
    {A, B}, {x, -L/2, L/2}, {y, -L/2, L/2}, {t, 0, 10}, 
    MaxSteps -> 10]

The code causes a kernel crash and I do not know how to deal with it.

Answer 1

To demonstrate numerical solution we need to put $L=2 n \pi, n=1, 2, ..., 15, ...$ (for example in a paper $n=15$). This numerical solution with $n=1$ demonstrates zigzag instabilities

eq1 = D[A[x, y, t], t] == 
   r A[x, y, t] - \[Mu] Abs[A[x, y, t]]^2 A[x, y, 
      t] - \[Nu] Abs[B[x, y, t]]^2 A[x, y, t] + 
    D[A[x, y, t], {x, 2}] + \[Epsilon]/
      I (D[A[x, y, t], {x, 3}] + 
       D[D[A[x, y, t], {x, 1}], {y, 2}]) + (\[Epsilon]/(2 I))^2 (D[
        A[x, y, t], {x, 4}] + D[A[x, y, t], {y, 4}]);

eq2 = D[B[x, y, t], t] == 
   r B[x, y, t] - \[Mu] Abs[B[x, y, t]]^2 B[x, y, 
      t] - \[Nu] Abs[A[x, y, t]]^2 B[x, y, t] + 
    D[B[x, y, t], {y, 2}] + \[Epsilon]/
      I (D[B[x, y, t], {y, 3}] + 
       D[D[B[x, y, t], {y, 1}], {x, 2}]) + (\[Epsilon]/(2 I))^2 (D[
        B[x, y, t], {x, 4}] + D[B[x, y, t], {y, 4}]);

q = -2; \[Epsilon] = 0.1; r = 40.0; \[Mu] = 2.0; \[Nu] = 1.0; L = 2 Pi;
R0 = Sqrt[(r - q^2 (1 + \[Epsilon] q/2)^2)/(\[Mu] + \[Nu])];
As = R0 Exp[I q x]; Bs = R0 Exp[I q y];


{aS, bS} = 
 NDSolveValue[{eq1, eq2, A[x, y, 0] == As, B[x, y, 0] == Bs, 
   A[x, -L/2, t] == A[x, L/2, t], A[-L/2, y, t] == A[L/2, y, t], 
   B[x, -L/2, t] == B[x, L/2, t], B[-L/2, y, t] == B[L/2, y, t]}, {A, 
   B}, {x, -L/2, L/2}, {y, -L/2, L/2}, {t, 0, 1}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 40, "MaxPoints" -> 80, 
      "DifferenceOrder" -> "Pseudospectral"}}]

ContourPlot[
 Re[aS[x, y, 1] + bS[x, y, 1]], {x, -L/2, L/2}, {y, -L/2, L/2}, 
 ColorFunction -> Hue]