# Solving a nonlinear PDE with Mathematica10 FEM Solver

I am trying to solve a system of coupled nonlinear PDEs in a rectangular region with the new FEM solver in Mathematica 10. However, I come across an error stating

NDSolveValue::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

Is there any workaround? I find it hard to believe that this task should not be possible. Do I have to adjust the method of the solver?

Here is my complete code:

Ω = Rectangle[{0, 0}, {2, 2}];
pdes = {-Inactive[Laplacian][u[t, x, y], {x, y}] + Derivative[1, 0, 0][u][t, x, y] ==
0.6*u[t, x, y] - u[t, x, y]^3 - v[t, x, y],
-Inactive[Laplacian][v[t, x, y], {x, y}] + Derivative[1, 0, 0][v][t, x, y] ==
1.5*u[t, x, y] - 2*v[t, x, y]};
c = {
(** ic **)
u[0, x, y] == Exp[-5 ((x - 3/2)^2 + (y - 3/2)^2)],
v[0, x, y] == Exp[-5 ((x - 1/2)^2 + (y - 1/2)^2)]
(** blank bc = noflux **)
};
{usol, vsol} = NDSolveValue[{pdes, c},
{u, v},
{x, y} ∈ Ω,
{t, 0, 2 π}
];


If you replace your input with

pdes = {-Derivative[0, 0, 2][u][t, x, y] - Derivative[0, 2, 0][u][t, x, y]
+ Derivative[1, 0, 0][u][t, x, y] == 0.6*u[t, x, y] - u[t, x, y]^3 - v[t, x, y],
-Derivative[0, 0, 2][v][t, x, y] - Derivative[0, 2, 0][v][t, x, y] +
Derivative[1, 0, 0][v][t, x, y] == 1.5*u[t, x, y] - 2*v[t, x, y]};

c = {
(** ic **)

u[0, x, y] == Exp[-5 ((x - 3/2)^2 + (y - 3/2)^2)],
v[0, x, y] == Exp[-5 ((x - 1/2)^2 + (y - 1/2)^2)]
(** blank bc = noflux **)
};
{usol, vsol} =
NDSolveValue[{pdes, c}, {u, v}, {t, 0, 2 π}, {x, 0, 2}, {y, 0,
2}];


NDSolveValue::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution.

NDSolveValue::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable y. Artificial boundary effects may be present in the solution.

NDSolveValue::mxst: Maximum number of 14516 steps reached at the point t == 2.9421825263433456

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 577.1533477585843at t = 2.9421825263433456 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 35 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

But you get something back. However, you'd need to specify boundary conditions to get something decent.

Concerning FEM, no, currently (V10) can not deal with non-linear stuff out of the box, it will in a future version. You can, however, use the low level FEM functions to code one up your self today. Give that a shot. There is an example of how to write a non linear FEM solver.

• That is somewhat disappointing... What does "future version" mean? Version 10.1 or 13.0? What are the low level FEM functions you are referring to? – Oscillon Aug 12 '14 at 19:30
• @Oscillum, I don't know. Writing a non-linear solve that works generally from scratch is not entirely trivial. If you provide a problem set that works with NDSolve on a rectangular region it could be possible to write such a non-linear solve with the low level FEM functions. Have you tried that? – user21 Aug 12 '14 at 19:40
• @Oscillum, added a link to the low level fem functions – user21 Aug 15 '14 at 9:09
• @m_goldberg, is there a tool to do the formatting? – user21 Aug 15 '14 at 9:38
• Is FEM a must? I have had a lot of success with non linear PDEs with NDSolve. Yes, you will have to deal with stiffness but the Method->LSODA is quite the "golden hammer" for stiff non linear problems. – dearN Sep 19 '15 at 2:48

In Version 12.0 you can solve this:

pdes = {-Derivative[0, 0, 2][u][t, x, y] -
Derivative[0, 2, 0][u][t, x, y] +
Derivative[1, 0, 0][u][t, x, y] ==
0.6*u[t, x, y] - u[t, x, y]^3 -
v[t, x, y], -Derivative[0, 0, 2][v][t, x, y] -
Derivative[0, 2, 0][v][t, x, y] +
Derivative[1, 0, 0][v][t, x, y] == 1.5*u[t, x, y] - 2*v[t, x, y]};

c = {(**ic**)
u[0, x, y] == Exp[-5 ((x - 3/2)^2 + (y - 3/2)^2)],
v[0, x, y] == Exp[-5 ((x - 1/2)^2 + (y - 1/2)^2)]
(**blank bc=noflux**)};

{usol, vsol} =
NDSolveValue[{pdes, c}, {u, v}, {t, 0, 2 \[Pi]}, {x, 0, 2}, {y, 0,
2}, Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}];
Plot3D[usol[2 \[Pi], x, y], {x, 0, 2}, {y, 0, 1}]
`