# Parsing error [bug?] with systems of nonlinear PDEs, Mathematica 12.0.0

Bug introduced in 12.0 and persisting through 12.1.1 - Fixed in Version: 12.2

I find the following behavior using 12.0.0. (EDIT: I asked a friend to try it in 12.1.1 and he finds the same result.)

Here is a system of coupled nonlinear PDEs that NDSolve cannot parse (don't worry about the system itself, it is just the simplest example I've found that produces this behavior):

c = {{1, 0}, {0, v[x, y]}};
alpha = {0, -u[x, y]};
NDSolveValue[{-Inactive[Div][c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
-Inactive[Div][Inactive[Times][alpha, v[x, y]], {x, y}] ==0}, {u[x, y], v[x, y]},
Element[{x, y}, Disk[]]]


The output is

NDSolveValue::femper: PDE parsing error of Div[{{1,0},{0,v}}.Grad$9730]. Inconsistent equation dimensions. It seems like extremely basic functionality to be a bug. Maybe I am missing something simple. I tried including copious calls to Inactive so the parser would have no trouble identifying the coefficients, which might otherwise be ambiguous with nonlinear systems. Note that parsing is no problem with a similar linear problem: c = {{1, 0}, {0, 1}}; alpha = {0, -1};  There is also no trouble parsing when reducing the dependent variables to one dimension ({u} instead of {u,v}), but keeping it nonlinear (e.g. c1 = {{1, 0}, {0, u[x, y]}}). So this trouble seems to be due to the combination of being coupled and nonlinear. Also note there is an obvious workaround, which is to go the "FEM programming" route and just specify the pde coefficients via InitializePDECoefficients. But still... what is up here? • I am looking at this right now, and it smells like a bug. I am not at the bottom of this yet. I think you best bet right now would be to use InitializePDECoefficients. I'll keep digging. Also, the behavior is the same in 12.1 Commented Aug 19, 2020 at 7:29 • Thanks, @user21 Commented Aug 19, 2020 at 8:03 • Do you by any chance have a solution and/or boundary conditions for this. I'd like to test a potential fix. Commented Aug 21, 2020 at 6:44 • Good question. Actually, I am working on a research project involving a more complicated system, and the nature/existence of the solutions is an open question. So this would not be useful for you, I think. I am working on a simple limit, which I can forward to you when ready. Commented Aug 23, 2020 at 6:35 • I have committed a fix for this. This is going to be available in V12.2. Thanks again for reporting. Commented Nov 3, 2020 at 6:27 ## 2 Answers This is a bug and has been fixed in version 12.2 c = {{1, 0}, {0, v[x, y]}}; alpha = {0, -u[x, y]}; NDSolveValue[{-Inactive[Div][ c . Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0, -Inactive[Div][Inactive[Times][alpha, v[x, y]], {x, y}] == 0}, {u[x, y], v[x, y]}, Element[{x, y}, Disk[]]]  You'll get an expected warning about missing boundary conditions but other than that it returns a solution. I have the same bug with NDSolve. For Mathematica 11, it works. But the same notebook is run in Mathematica 12.0, it produces: NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.. In[2]:=$Version

Out[2]= "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)"
`
• I doubt that this is the same issue: Version 11 did not have a nonlinear FEM solver and the error message tag is not the same. You must have solved some other problem in version 11. Commented Nov 4, 2020 at 7:05