# Solving a system of temporal non-linear (reaction-diffusion) PDEs over a region using Neumann conditions

I am trying to solve a system of PDEs with non-linear terms:

$\frac{\partial a(x,y,z,t)}{\partial t}=\color{red}{-\text{$\tau_2 $} a(x,y,z,t) h(x,y,z,t)}+\text{$\tau_1 $} d(x,y,z,t) \\\frac{\partial b(x,y,z,t)}{\partial t}=\color{red}{-\text{$\tau_2 $} b(x,y,z,t) i(x,y,z,t)}+\text{$\tau_1$} e(x,y,z,t) \\\frac{\partial c(x,y,z,t)}{\partial t}=\color{red}{-\text{$\tau_2 $} c(x,y,z,t) g(x,y,z,t)}+\text{$\tau_1$} f(x,y,z,t) \\\frac{\partial d(x,y,z,t)}{\partial t}=\color{red}{\text{$\tau_2 $} a(x,y,z,t) h(x,y,z,t)}-\text{$\tau_1 $} d(x,y,z,t) \\\frac{\partial e(x,y,z,t)}{\partial t}=\color{red}{\text{$\tau_2 $} b(x,y,z,t) i(x,y,z,t)}-\text{$\tau_1 $} e(x,y,z,t) \\\frac{\partial f(x,y,z,t)}{\partial t}=\color{red}{\text{$\tau_2 $} c(x,y,z,t) g(x,y,z,t)}-\text{$\tau_1 $} f(x,y,z,t) \\\frac{\partial g(x,y,z,t)}{\partial t}=\color{blue}{\mathscr{D} \nabla _{\{x,y,z\}}^{}g(x,y,z,t)}+\text{$\tau_3$} a(x,y,z,t)-\frac{g(x,y,z,t)}{\text{$\tau $4 }}+\text{$\tau_1 $} f(x,y,z,t) \\\frac{\partial h(x,y,z,t)}{\partial t}=\color{blue}{\mathscr{D} \nabla _{\{x,y,z\}}^{}h(x,y,z,t)}+\text{$\tau_3$} b(x,y,z,t)-\frac{h(x,y,z,t)}{\text{$\tau $4 }}+\text{$\tau_1 $} d(x,y,z,t) \\\frac{\partial i(x,y,z,t)}{\partial t}=\color{blue}{\mathscr{D} \nabla _{\{x,y,z\}}^{}i(x,y,z,t)}+\text{$\tau_3 $} c(x,y,z,t)-\frac{i(x,y,z,t)}{\text{$\tau $4 }}+\text{$\tau_1 $} e(x,y,z,t)$

with non-linear terms in $\color{red}{red}$ and spatial terms in $\color{blue}{blue}$

i.e.

pdes = {
Derivative[0, 0, 0, 1][a][x, y, z, t] == 0.05*d[x, y, z, t] - 0.05*a[x, y, z, t]*h[x, y, z, t],
Derivative[0, 0, 0, 1][b][x, y, z, t] == 0.05*e[x, y, z, t] - 0.05*b[x, y, z, t]*i[x, y, z, t],
Derivative[0, 0, 0, 1][c][x, y, z, t] == 0.05*f[x, y, z, t] - 0.05*c[x, y, z, t]*g[x, y, z, t],
Derivative[0, 0, 0, 1][d][x, y, z, t] == -0.05*d[x, y, z, t] + 0.05*a[x, y, z, t]*h[x, y, z, t],
Derivative[0, 0, 0, 1][e][x, y, z, t] == -0.05*e[x, y, z, t] + 0.05*b[x, y, z, t]*i[x, y, z, t],
Derivative[0, 0, 0, 1][f][x, y, z, t] == -0.05*f[x, y, z, t] + 0.05*c[x, y, z, t]*g[x, y, z, t],
Derivative[0, 0, 0, 1][g][x, y, z, t] == 100*a[x, y, z, t] + 0.05*f[x, y, z, t] +
0.05*(Derivative[0, 0, 2, 0][g][x, y, z, t] + Derivative[0, 2, 0, 0][g][x, y, z, t] + Derivative[2, 0, 0, 0][g][x, y, z, t]),
Derivative[0, 0, 0, 1][h][x, y, z, t] == 100*b[x, y, z, t] + 0.05*d[x, y, z, t] +
0.05*(Derivative[0, 0, 2, 0][h][x, y, z, t] + Derivative[0, 2, 0, 0][h][x, y, z, t] + Derivative[2, 0, 0, 0][h][x, y, z, t]),
Derivative[0, 0, 0, 1][i][x, y, z, t] == 100*c[x, y, z, t] + 0.05*e[x, y, z, t]  +
0.05*(Derivative[0, 0, 2, 0][i][x, y, z, t] + Derivative[0, 2, 0, 0][i][x, y, z, t] + Derivative[2, 0, 0, 0][i][x, y, z, t])
};


with the following intitial conditions:

initcs = {
a[x, y, z, 0] == (Sqrt[40/Pi])/
E^(40*((0.5 + x)^2 + y^2 + z^2)),
b[x, y, z, 0] == (Sqrt[40/Pi])/E^(40*(x^2 + y^2 + z^2)),
c[x, y, z, 0] == (Sqrt[40/Pi])/E^(40*((-0.5 + x)^2 + y^2 + z^2)),
d[x, y, z, 0] == 0, e[x, y, z, 0] == 0, f[x, y, z, 0] == 0,
g[x, y, z, 0] == 0, h[x, y, z, 0] == 0, i[x, y, z, 0] == 0
};


if I solve this in a cubic region I DO get an answer (although it tells me that the step-size might be too large):

sol = NDSolve[
Flatten[{pdes, initcs}], {a, b, c, d, e, f, g, h, i}, {x, -1,
1}, {y, -1, 1}, {z, -1, 1}, {t, 0, 1}]


to plot:

    Export["disks.gif",
ListDensityPlot3D /@
Transpose[sol[[1, 9, 2]]["ValuesOnGrid"], {2, 3, 4, 1}]] However, I want to solve it in a specific region (a complex curved region). Lets take a cuboid region as an example since it should give the exact same solution:

 sol2 = NDSolve[
Flatten[{pdes, initcs}], {a, b, c, d, e, f, g, h,
i}, {x, y, z} \[Element] Cuboid[{-1, -1, -1}, {1, 1, 1}], {t, 0,
1}]


this gives me an error, even though it is the exact same problem

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


Why does the second method not work when the second one does? How can I solve my problem?

Edit: I have been suggested to look at the amazing answer of user21 to solving the naiver-stokes equation. This seems like the right way to start, but this solves the steady state instead of the here required time-resolved solution.

After linearizion (see chapter 4 and 5) I come to:

alfabet = {a, b, c, d, e, f, g, h, i};
coords = {x, y, z};
rulefunct = # -> #[x, y, z] & /@ alfabet;
alfabet2 = alfabet /. rulefunct;
F = { #} & /@ -{a*h - τ1 d, b*i - τ1 e,
c*g - τ1 f, -a*h + τ1 d, -b*i + τ1 e, -c*
g + τ1 f, τ3*g, τ3 h, τ3 i} /.
rulefunct /. {τ1 -> 1, τ3 -> 1};
A = Table[-D[F[[α]], alfabet2[[β]]], {α, 9}, {β,9}];
σ = -Normal[SparseArray[Table[{i, i, j, j} -> - d, {i, 7, 9}, {j, 1, 3}] // Flatten[# , 1] &]];
Γ =  Join[ ConstantArray[0, {6, 3}],Table[-d D[alfabet2[[α]],coords[[β]]], {α, 7,9}, {β, 3}]];
τ = IdentityMatrix


to implement in:

nr = ToNumericalRegion[Ball[]];
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {alfabet,
coords}];
sd = NDSolveSolutionData["Space" -> nr];
nlPdeCoeff = InitializePDECoefficients[vd, sd, "LoadCoefficients" ->(*F*)F,
"ReactionCoefficients" ->(*a*)A,
"DampingCoefficients" -> IdentityMatrix,
"DiffusionCoefficients" -> σ]


I do not yet see a way to give the right initial conditions(and initialize a 4D region?) such that the coeficients can indeed be scalar given I require the temporal solution.

• The error message is quite clear, FEM still doesn't support Nonlinear coefficients.
– zhk
Jul 10, 2017 at 11:40
• In sol1 NDSolve doesn't use FEM but by default uses MOL that's why it solve the system with some warnings.
– zhk
Jul 10, 2017 at 12:03
• "Or is the method of lines unavailable for arbitrary shapes?" No, it's "TensorProductGrid" method that is unavailable for arbitrary shapes. "TensorProductGrid" and "FiniteElement" can both be used for spatial discretization of "MethodOfLines", the former can handle nonlinear coefficient but can't handle irregular domain, while the latter can handle irregular domain but cannot handle nonlinear coefficient (at least now). Here is a related post. Jul 22, 2017 at 5:46
• If you want to solve nonlinear PDE in irregular domain, then some relatively low level programming is needed. There already exist several examples in this site. You can start from this post. Jul 22, 2017 at 5:48
• BTW, if low-level FEM programming is too hard, there exists another possible method, which is not that accurate but maybe acceptable: please check this answer. Jul 27, 2017 at 12:04

Needs["NDSolveFEM"] 