# Solving the system of nonlinear PDE

I have never used Mathematica before. I am trying to solve the system of PDE

$$\begin{array}{l} \frac{\partial u}{\partial t} = D_u \nabla^2 u - \beta \chi_u \nabla ( u \nabla \pi_u) + u (\pi_u - \kappa (u+v)) \\ \frac{\partial v}{\partial t} = D_v \nabla^2 v - \beta \chi_v \nabla ( v \nabla \pi_v) + v (\pi_v - \kappa (u+v)) \end{array}$$

where

$$\begin{array}{l} \pi_u = a \frac{u}{u + v} + b \frac{v}{u + v} \\ \pi_v = c \frac{u}{u + v} + d \frac{v}{u + v} \end{array}$$

are the function with respect to $$u$$ and $$v$$.

And I am giving the initial values : $$u$$~Normal(44.444, 0.0001) and $$v$$~Normal(22.222, 0.0001)

And this is the Mathamatica Code I tried

{beta, chiu, chiv, kappa, Du, Dv} = {10, 0.3, 2.4, 0.001, 0.733, 0.733}
{a, b, c, d} = {-0.1, 0.4, 0, 0.2}
piu[u_, v_] := {a*(u/(u + v)) + b*(v/(u + v))}
piv[u_, v_] := {c*(u/(u + v)) + d*(v/(u + v))}
ph = NDSolveValue[{\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[t, x, y]$$\) ==
Du*Inactive[Laplacian][u[t, x, y], {x, y}] -
beta*chiu*
piu[u[t, x, y], v[t, x, y]], {x, y}]), {x, y}] +
u[t, x, y]*(piu[u[t, x, y], v[t, x, y]] -
kappa*(u[t, x, y] + v[t, x, y])), \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$v[t, x, y]$$\) ==
Dv*Inactive[Laplacian][v[t, x, y], {x, y}] -
beta*chiv*
piv[u[t, x, y], v[t, x, y]], {x, y}]), {x, y}] +
v[t, x, y]*(piv[u[t, x, y], v[t, x, y]] -
kappa*(u[t, x, y] + v[t, x, y])),
u[t, x, -1] == u[t, x, 1] == u[t, -1, y] == u[t, 1, y] ==
v[t, x, -1] == v[t, x, 1] == v[t, -1, y] == v[t, 1, y] == 0,
u[0, x, y] == RandomVariate[NormalDistribution[22.222, 0.0001]],
v[0, x, y] ==
RandomVariate[NormalDistribution[44.444, 0.0001]]}, {u, v}, {t, 0,
10}, {x, -1, 1}, {y, -1, 1}]


But the error message says that : Objects of unequal length in ...

• Expression \[Del]\[Del]u[...] gives vectorgradient. Perhaps your pde needs \[Del]\[CenterDot]\[Del]u  Commented Mar 6 at 8:52
• Grad[a Grad[b]] is undefined, probably by $\nabla(a \nabla(b))$ = Div[a Grad[b,{x,y}],{x,y}] is physically and algebraically correct. Commented Mar 6 at 9:35
• @UlrichNeumann Can we solve this problem using NDSolve or not? Commented Mar 7 at 14:13
• @AlexTrounev Ride hand side of the pde's is only scalar if we have the form Grad[...]\[CenterDot] (...Grad[]) ( dot product ) . For this case NDSolve should be able to solve. Still waiting for a reply to my first comment from OP Commented Mar 7 at 14:31
• @random487510 Did you check the correct form of your pde's? By the way the Initial conditions (!=0) and the boundary conditions (==0) are inconsistent! Commented Mar 9 at 7:35

It seems that we can't solve this system of PDEs using NDSolve. Nevertheless we can solve it as system of ODEs using hand made method of lines code as follows

noise = 0.01; conu0 = 22.222; conv0 = 44.444;

{beta, chiu, chiv, kappa, du, dv} = {9.8, 0.3, 2.4, 0.001, 0.733,
0.733};
{a, b, c, d} = {-0.1, 0.4, 0, 0.2};
piu[u_, v_] := a*(u/(u + v)) + b*(v/(u + v));
piv[u_, v_] := c*(u/(u + v)) + d*(v/(u + v));

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{tmax, domain, pts, difforder} = {10, {{-1, 1}, {-1, 1}}, {21, 21}, 4};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, difforder]; {uvar0,
vvar0} = Table[ConstantArray[0, {nx ny, tmax}], {2}]; boundaries =
Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1];
uvar = Table[uu[i][t], {i, nx ny}]; vvar = Table[vv[i][t], {i, nx ny}];

eqnu = du (dx2 + dy2) . uvar +
uvar (piu[uvar, vvar] - kappa (uvar + vvar)) -
beta chiu (dx . (uvar dx . piu[uvar, vvar]) +
dy . (uvar dy . piu[uvar, vvar])); eqnv =
dv (dx2 + dy2) . vvar +
vvar (piv[uvar, vvar] - kappa (uvar + vvar)) -
beta chiv (dx . (vvar dx . piv[uvar, vvar]) +
dy . (vvar dy . piv[uvar, vvar]));

nxy = Complement[Range[1, nx ny], boundaries]; iniu =
Table[uu[i][0] == conu0 + noise*RandomReal[{-1, 1}], {i,
nx ny}]; iniv =
Table[vv[i][0] == conv0 + noise RandomReal[{-1, 1}], {i, nx ny}];


This system of ODEs we can solve with NDSolve. Before do this we should note that system is unstable at $$\beta>9.75$$, while author asked to solve it at $$\beta =10$$. As compromise we solve it at $$\beta=9.8$$, but solution lost stability at t = 0.8939960639052812

sol = NDSolveValue[
Join[Table[uu[i]'[t] == eqnu[[i]], {i, nxy}],
Table[vv[i]'[t] == eqnv[[i]], {i, nxy}],
Table[uu[i]'[t] == 0, {i, boundaries}],
Table[vv[i]'[t] == 0, {i, boundaries}], iniu, iniv],
Join[uvar, vvar], {t, 0, 1}, Method -> "ImplicitRungeKutta",
AccuracyGoal -> 5, PrecisionGoal -> 4, MaxSteps -> 10^6];


Note that we start from random initial conditions that look like

{u, v} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[sol /. t -> .0, Length[sgrid]]];

{ContourPlot[u[x, y], {x, -1, 1}, {y, -1, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
PlotRange -> All],
ContourPlot[v[x, y], {x, -1, 1}, {y, -1, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v", PlotRange -> All]}


Final state looks quite coherent

{u, v} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[sol /. t -> .89, Length[sgrid]]];

{ContourPlot[u[x, y], {x, -1, 1}, {y, -1, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
PlotRange -> All],
ContourPlot[v[x, y], {x, -1, 1}, {y, -1, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v", PlotRange -> All]}


Update 1. We can compare solution computed with code above at beta=2 to solution computed with code proposed by Ulrich Neumann. First note that in his code there is a typo in equation for v, and after clean up his code looks as follows (we generate mesh same as sgrid above)

SeedRandom[1234];
<< NDSolveFEM
t0 = 1;
reg = Rectangle[{-1, -1}, {1, 1}]; mesh =
ToElementMesh[reg, "MaxCellMeasure" -> 1/110,
n = Length[
mesh["Coordinates"]]; noise = 0.01; conu0 = 22.222; conv0 = 44.444;
u0 = ElementMeshInterpolation[{mesh},
conu0 + noise*(RandomReal[{-1, 1}, n])]; v0 =
ElementMeshInterpolation[{mesh},
conv0 + noise*(RandomReal[{-1, 1}, n])];
{beta, chiu, chiv, kappa, Du, Dv} = {2, 0.3, 2.4, 0.001, 0.733, 0.733};
{a, b, c, d} = {-0.1, 0.4, 0, 0.2};
piu[u_, v_] := a*(u/(u + v)) + b*(v/(u + v));
piv[u_, v_] := c*(u/(u + v)) + d*(v/(u + v));

solI = NestList[
NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] ==
Inactive[
Div][(Du -
beta chiu (a -
b) (#[[1]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] + #[[2]][t, x,
y])^2)) Inactive[Grad][u[t, x, y], {x, y}], {x,
y}] + Inactive[
Div][(beta chiu (a -
b) (#[[1]][t, x, y] #[[1]][t, x,
y])/((#[[1]][t, x, y] + #[[2]][t, x, y])^2)) Inactive[
Grad][v[t, x, y], {x, y}], {x, y}] +
u[t, x, y]*(piu[#[[1]][t, x, y], #[[2]][t, x, y]] -
kappa*(#[[1]][t, x, y] + #[[2]][t, x, y])),
Derivative[1, 0, 0][v][t, x, y] ==
Inactive[
Div][(Dv +
beta chiv (c -
d) (#[[1]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] + #[[2]][t, x,
y])^2)) Inactive[Grad][v[t, x, y], {x, y}], {x,
y}] - Inactive[
Div][(beta chiv (c -
d) (#[[2]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] + #[[2]][t, x, y])^2)) Inactive[
Grad][u[t, x, y], {x, y}], {x, y}] +
v[t, x, y]*(piv[#[[1]][t, x, y], #[[2]][t, x, y]] -
kappa*(#[[1]][t, x, y] + #[[2]][t, x, y])),
DirichletCondition[u[t, x, y] == u0[x, y], True],
DirichletCondition[v[t, x, y] == v0[x, y], True],
u[0, x, y] == u0[x, y], v[0, x, y] == v0[x, y]}, {u, v}, {t, 0,
1}, Element[{x, y}, mesh],
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement"}}] &, {u0[x, y] &,
v0[x, y] &}, 3];


Visualization

GraphicsRow[{ContourPlot[solI[[-1, 1]][1, x, y],
Element[{x, y}, region], ColorFunction -> "Rainbow",
Contours -> 20, PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "u", PlotRange -> All],
ContourPlot[solI[[-1, 2]][1, x, y], Element[{x, y}, region],
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v", PlotRange -> All]}]


In my code we put beta=2 and use same pseudorandom generator state SeedRandom[1234], as result we have

Note that pictures look very similar, but we also can compare detailed pictures as

py = {Plot[solI[[-1, 1]][1, 0, y], {y, -1, 1},
PlotStyle -> {Red, Dashed}],
Plot[solI[[-1, 2]][1, 0, y], {y, -1, 1},
PlotStyle -> {Red, Dashed}]};
py1={Plot[u[0, y], {y, -1, 1}], Plot[v[0, y], {y, -1, 1}]};
{Show[py1[[1]],py[[1]]],Show[py1[[1]],py[[1]]]}


Here solid lines computed with my code and dashed lines computed with Ulrich code. The agreement is fine. Therefore we pass test at beta=2.
Test for beta=9 gives last stably solution computed with Ulrich code

Same solution computed with my code

Since solutions look like identical we can conclude that Ulrich code and my code passes test at beta=9 as well. For beta=9.8 Ulrich code runs forever without answer, while with my code we have picture shown above.

• Nice solution (+1) Which pde's did you use? The "Div[Grad[...]]" form? Commented Mar 11 at 7:45
• Yes, we convert "inactive" form in matrix form. Maybe it could be better to use an activated form. Commented Mar 11 at 10:30
• Thanks. I am still trying to find a FEM solution, perhaps iterative , but couldn't find a correct maybe inactive form Commented Mar 11 at 10:45
• Meanwhile I got a nonlinear FEM-solution(see my second answer) , which unfortunately differs form your solution. Have to check the pde-system... Commented Mar 11 at 12:31
• My latest insight: If you remove , Method -> "ImplicitRungeKutta", AccuracyGoal -> 5, PrecisionGoal -> 4, MaxSteps -> 10^6 in sol you get complete stationary solution in the time range 0<t<10 ! Commented Mar 13 at 13:22

To long for a comment, as a reply to @AlexTrounev.

First change the definitions (no curly brakets!)

piu[u_, v_] := a*(u/(u + v)) + b*(v/(u + v))
piv[u_, v_] := c*(u/(u + v)) + d*(v/(u + v))


If we change the pde to Grad[...Grad[]] to Div[...Grad[...]] , which should be confirmed by the OP(!), we get system equations

eqns = {Derivative[1, 0, 0][u][t, x, y] ==
Du* Laplacian [u[t, x, y], {x, y}] -
beta*chiu*
Inactive[Div] [
piu[u[t, x, y], v[t, x, y]], {x, y}] , {x, y}] +
u[t, x, y]*(piu[u[t, x, y], v[t, x, y]] -
kappa*(u[t, x, y] + v[t, x, y])),
Derivative[1, 0, 0][v][t, x, y] ==
Dv* Laplacian [v[t, x, y], {x, y}] -
beta*chiv*
Inactive[Div][
v[t, x, y]* Inactive[ Grad] [
piv[u[t, x, y], v[t, x, y]], {x, y}] , {x, y}] +
v[t, x, y]*(piv[u[t, x, y], v[t, x, y]] -
kappa*(u[t, x, y] + v[t, x, y])) ,
u[t, x, -1] == u[t, x, 1] == u[t, -1, y] == u[t, 1, y] ==
v[t, x, -1] == v[t, x, 1] == v[t, -1, y] == v[t, 1, y] == 0,
u[0, x, y] == RandomVariate[NormalDistribution[22.222, 0.0001]],
v[0, x, y] ==
RandomVariate[NormalDistribution[44.444, 0.0001]]}


Still there are nonlinearities on the right handside, which have to be linearized to make NDSolve FEM solver work.

I tried examplary Inactive[ Grad [piv[u[t, x, y], v[t, x, y]], {x, y}] ] but did not succeed...

If I apply this last modification NDSolve gives following messages

"Grad::argtu: Grad called with 1 argument; 2 or 3 arguments are expected."

and

"NDSolve::femper: PDE parsing error of Grad[ -((0.1 u)/(u+v))+(0.4 v)/(u+v) ]. Inconsistent equation dimensions."

Any ideas?

Hoping my transformation are correct I actually get these pdes:

eqns = {Derivative[1, 0, 0][u][t, x, y] ==
Du* Laplacian [u[t, x, y], {x, y}] -
beta*chiu*

Div [   (
Inactivate[((a - b) v[t, x, y] u[t, x,
y])/(u[t, x, y] + v[t, x, y])^2] Grad   [
u[t, x, y], {x, y}] -
Inactivate[((a - b) u[t, x, y] u[t, x,
y])/(u[t, x, y] + v[t, x, y])^2] Grad   [
v[t, x, y], {x, y}]  ) , {x, y}] +
u[t, x, y]*
Inactivate[(a*u[t, x, y]/(u[t, x, y] + v[t, x, y]) +
b*v[t, x, y]/(u[t, x, y] + v[t, x, y]) -
kappa*(u[t, x, y] + v[t, x, y]))],
Derivative[1, 0, 0][v][t, x, y] ==
Dv* Laplacian [v[t, x, y], {x, y}] -
beta*chiv*

Div  [(Inactivate[((c - d) v[t, x, y] v[t, x,
y])/(u[t, x, y] + v[t, x, y])^2] Grad   [
u[t, x, y], {x, y}] -
Inactivate[((c - d) u[t, x, y] v[t, x,
y])/(u[t, x, y] + v[t, x, y])^2] Grad   [
v[t, x, y], {x, y}]  ), {x, y}] +
v[t, x, y]*
Inactivate[(c*u[t, x, y]/(u[t, x, y] + v[t, x, y]) +
d*v[t, x, y]/(u[t, x, y] + v[t, x, y]) -
kappa*(u[t, x, y] + v[t, x, y])) ],
DirichletCondition[u[t, x, y] == 0, True],
DirichletCondition[v[t, x, y] == 0, True] ,
u[0, x, y] == RandomVariate[NormalDistribution[22.222, 0.0001]],
v[0, x, y] ==
RandomVariate[NormalDistribution[44.444, 0.0001]]}


But solution fails

NDSolve[eqns, {u, v}, {t, 0, 10}, {x, -1, 1}, {y, -1, 1},
Method -> {"MethodOfLines", "TemporalVariable" -> t,"SpatialDiscretization" -> {"FiniteElement"}}]


error message is unclear to me "NDSolve::underdet: There are more dependent variables, {u[t,x,y],v[t,x,y],u[t,x,y],v[t,x,y]}, than equations, so the system is underdetermined."

My handmade transformation give this form of the two pde's. Actually I don't know how to inactivate parts of these equations to get standard form accepted by the FEM solver...

• So, FEM failed. Maybe we can use "MethodOfLines"? Commented Mar 8 at 0:44
• @AlexTrounev Yess, MethodOfLines and FiniteElement. But pde has to be prepared to avoid parsing error Commented Mar 8 at 7:58
• Try Activate[eqns], then you will see real problem of FEM . Commented Mar 8 at 12:44
• @AlexTrounev What exactly do you mean? Mathematica v12.2 shows two nonlinear pde's without error message Commented Mar 8 at 14:07
• Use activation in NDSolve :) Commented Mar 8 at 14:13

modified typo corrected, thanks @AlexTrounev

With manual transformation using

it's possible to find an iterative FEM solution:

region=Rectangle[{-1,-1},{1,1}]
{beta, chiu, chiv, kappa, Du, Dv} = {2, 0.3, 2.4, 0.001, 0.733,0.733}
{a, b, c, d} = {-0.1, 0.4, 0, 0.2}

u0 = RandomVariate[NormalDistribution[22.222, 0.0001]];
v0 = RandomVariate[NormalDistribution[44.444, 0.0001]];

solI = NestList[
NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] ==
Inactive[
Div][(Du -
beta chiu (a - b) (#[[1]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] +  #[[2]][t, x,
y])^2) ) Inactive[Grad][u[t, x, y], {x, y}], {x, y}] +
Inactive[
Div][( beta chiu (a - b) (#[[1]][t, x, y] #[[1]][t, x,
y])/((#[[1]][t, x, y] +  #[[2]][t, x, y])^2) ) Inactive[
Grad][v[t, x, y], {x, y}], {x, y}] +
u[t, x, y]*(piu[#[[1]][t, x, y], #[[2]][t, x, y]] -
kappa*(#[[1]][t, x, y] + #[[2]][t, x, y])),
Derivative[1, 0, 0][v][t, x, y] ==
Inactive[
Div][(Dv +
beta chiv (c - d) (#[[1]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] + #[[2]][t, x, y])^2) ) Inactive[
Grad][v[t, x, y], {x, y}], {x, y}] -
Inactive[
Div][( beta chiv (c - d) (#[[2]][t, x, y] #[[2]][t, x,
y])/((#[[1]][t, x, y] +   #[[2]][t, x,
y])^2) ) Inactive[Grad][u[t, x, y], {x, y}], {x, y}] +
v[t, x,
y]*(piv[#[[1]][t, x, y], #[[2]][t, x, y]] -
kappa*(#[[1]][t, x, y] + #[[2]][t, x, y])),
DirichletCondition[u[t, x, y] == u0, True],
DirichletCondition[v[t, x, y] == v0
, True] , u[0, x, y] == u0, v[0, x, y] == v0
}, {u, v}, {t, 0, 10}, Element[{x, y}, region],
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement"}}] &, {u0 &, v0 &},
3];

GraphicsRow[{ContourPlot[solI[[-1, 1]][10, x, y],
Element[{x, y}, region], ColorFunction -> "Rainbow",
Contours -> 20, PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
PlotRange -> All] ,
ContourPlot[solI[[-1, 2]][10, x, y], Element[{x, y}, region],
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "v",
PlotRange -> All] }]


It would be very helpful to find the inactive form of the nonlinear equations, because we could avoid the slow iterative solution!

Any ideas?

• Your code is not complete. Du, Dv are not defined. Commented Mar 11 at 15:01
• Sorry, I added the definitions! Commented Mar 11 at 15:15
• Thank you (+1). It looks like your code working as an averaging procedure. With my code we can generate same picture for beta=2. Commented Mar 11 at 19:13
• Please check, maybe there is a typo in your code in equation for v. Commented Mar 11 at 22:53
• @AlexTrounev Where exactly do you suspect the error? Commented Mar 12 at 9:58

Sorry, first time for me to give three answers to one very interesting question

If we ara looking for a stationary solution of the two pde's we can give a fast nonlinear Finite Element solution as follows. Here a assume random start conditions InitialSeeding and fix the boundary accordingly.

mesh and interpolation for InitialSeeding:

Needs["NDSolveFEM"]
region = ToElementMesh[Rectangle[{-1, -1}, {1, 1}]];
pts = region["Coordinates"];
ui = Map[RandomVariate[NormalDistribution[22.222, 0.01]] &, pts];
ipu = ElementMeshInterpolation[region, ui](* InitialSeeding u[x,y]*)
vi = Map[ RandomVariate[NormalDistribution[44.444, 0.01]] &, pts];
ipv = ElementMeshInterpolation[region, vi](* InitialSeeding v[x,y]*)


parameters:

{beta, chiu, chiv, kappa, Du, Dv} = {10  , 0.3, 2.4, 0.001, 0.733,0.733}
{a, b, c, d} = {-0.1, 0.4, 0, 0.2}
piu[u_, v_] := a*(u/(u + v)) + b*(v/(u + v))
piv[u_, v_] := c*(u/(u + v)) + d*(v/(u + v))

c11 = Du - beta chiu (a - b) (u[x, y] v[x, y])/(u[x, y] + v[x, y])^2;
c12 = beta chiu (a - b) (u[x, y] u[x, y])/(u[x, y] + v[x, y])^2;
c21 = -beta chiv (c - d) (v[x, y] v[x, y])/(u[x, y] + v[x, y])^2;
c22 = Dv + beta chiv (c - d) (u[x, y] v[x, y])/(u[x, y] + v[x, y])^2;


nonlinear FEM solver

solS = NDSolveValue[{0 ==
Inactive[Div][(- {{c11, 0}, {0, c11}}) .
Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[Div][(- {{c12, 0}, {0, c12}}) .
Inactive[Grad][v[x, y], {x, y}], {x, y}] +
u[x, y]*(piu[u[x, y], v[x, y]] - kappa*(u[x, y] + v [x, y])),
0 == Inactive[Div][ (- {{c21, 0}, {0, c21}}) .
Inactive[Grad][v[x, y], {x, y}], {x, y}] -
Inactive[Div][(- {{c22, 0}, {0, c22}}) .
Inactive[Grad][u[x, y], {x, y}], {x, y}] +
v[x, y]*(piv[u[x, y], v[x, y]] - kappa*(u[x, y] + v[x, y])),
DirichletCondition[u[x, y] == ipu[x, y], True],
DirichletCondition[v[x, y] == ipv[x, y]
, True]
}, {u, v} , Element[{x, y}, region] ,
"InitialSeeding" -> {u[x, y] == ipu[x, y], v[x, y] == ipv[x, y]} ,
Method -> "PDEDiscretization" -> {"FiniteElement"}] ;


solution plot:

GraphicsRow[{
Plot3D[solS[[ 1]][ x, y], Element[{x, y}, region],
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> Automatic, MaxRecursion -> 2, PlotPoints -> 50,
PlotLabel -> "u[x,y] (beta=" <> ToString[beta] <> ")",
PlotRange -> All] ,
Plot3D[solS[[ 2]][ x, y], Element[{x, y}, region],
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> Automatic, MaxRecursion -> 2, PlotPoints -> 50,
PlotLabel -> "v[x,y] (beta=" <> ToString[beta] <> ")",
PlotRange -> All] }]


• Thank you for NDSolve solution (+1), but this is the final stage only. How to force NDSolve` to solve this system on $0\le t \le 1$? Commented Mar 13 at 16:44
• I tried the transient case too, but got error message "LinearSolve::exopt1: The option setting Method -> Multifrontal cannot be used with arbitrary-precision or exact arguments." Commented Mar 13 at 18:36