# How can I solve a fourth order system of PDEs?

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.

• "But is it wrong?" is too vague a question. Can you be specific and describe what is going wrong in detail? Commented Jan 13, 2021 at 7:49
• @m_goldberg It do not work, when I run the code, it quit the kernel Commented Jan 13, 2021 at 8:01
• @yunshi Where did you get this system of equations? Commented Jan 13, 2021 at 14:56
• @bbgodfrey the initial condition and pde are from a paper ong wavelength instabilities of square patterns, Physica D 67 (1993) 198 223 Commented Jan 16, 2021 at 8:45
• @yunshi This system of equations describes some temperature pattern in 2D convection flow in square cavity in some approximation. It is good for analytical research, but not for numerical. Nevertheless there are examples of numerical solution with pseudospectral method on a grid 100 by 100 in a paper sited. Commented Jan 16, 2021 at 11:46

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]