# NDSolve::femper for system of non-linear PDEs

I am trying to solve a system of PDE in 1 dimension. They are as follows:

My system of equations looks like follows:

$$\frac{\partial \bar{c^{T}_E}}{\partial \bar{t}} = \frac{\partial^2 \bar{c^{T}_E}}{\partial \bar{x}^2} + \alpha \frac{\partial^2}{\partial \bar{x}^2}\Big(\frac{\bar{c^{T}_E} \bar{c_S}}{K_M/c_{S0} + \bar{c_S}}\Big)$$

$$\frac{\partial \bar{c_S}}{\partial \bar{t}} = \beta \frac{\partial^2 \bar{c_S}}{\partial \bar{x}^2} - k_2 \frac{c_{E0}}{c_{S0}}\frac{\bar{c^{T}_E} \bar{c_S}}{K_M/c_{S0} + \bar{c_S}}\frac{d^2}{D_E}$$

The boundary conditions that I want to impose is the following:

$$\frac{\partial c_E}{\partial x} = 0 \hspace{0.5 cm}x = 0; \hspace{0.5cm} \frac{\partial c_E}{\partial x} = 0 \hspace{0.5 cm}x = 1$$

$$c_S(t,x=0) = 1 ; \hspace{0.5cm} c_S(t,x=1) = 0$$

The initial conditions are as follows:

$$c_E(0,x) = 1; \hspace{0.5cm} c_S(0,x=0) = 1 ; c_S(0,x>0) = 0$$

When I try to implement FEM on this system, I am getting the following issue:

NDSolveValue::femper: PDE parsing error of Grad[(Ca Cs)/(2/5+Cs)]. Inconsistent equation dimensions.

Here is the code:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
Needs["NDSolveFEM"]
(*Enzyme-substrate dyamics - hexokinase parameters*)
Clear["Global'*"];
Clear[Derivative];
Ca0 = 10^(-6);(*Initial Protein Concentration*)
d = 180*10^(-6);(*Half-width of channel*)
Da = 7.2*10^(-11);
CB0  = 10*10^(-3);
Km = 4*10^(-3);(*Michealis Menten Constant*)
beta  = 10;(*Subscript[D, s]/Subscript[D, E]*)
gamma = (Ca0/CB0)*(d^2/Da);
k2 = 8000;
alpha = -0.5;
op1 = D[Ca[t, x], t] -
alpha*Inactive[Div][
Inactive[
Grad][(Ca[t, x]*Cs[t, x])/(Km/CB0 + Cs[t, x]), {x}], {x}];
op2 = D[Cs[t, x], t] -
k2*gamma*(Ca[t, x]*Cs[t, x])/(Km/CB0 + Cs[t, x]);
IC = {Ca[0, x] == 1, Cs[0, x] == Exp[-1000 x]};
Subscript[\[CapitalGamma], N1] = NeumannValue[0, x == 0];
Subscript[\[CapitalGamma], N2] = NeumannValue[0, x == 1];
Subscript[\[CapitalGamma], NS1] =
DirichletCondition[Cs[t, x] == 1, x == 0];
Subscript[\[CapitalGamma], NS2] =
DirichletCondition[Cs[t, x] == 0, x == 1];
eqn1 = {op1 ==
Subscript[\[CapitalGamma], N1] + Subscript[\[CapitalGamma], N2]};
sol = NDSolveValue[{op1 ==
Subscript[\[CapitalGamma], N1] + Subscript[\[CapitalGamma], N2],
op2 == 0, IC, Subscript[\[CapitalGamma], NS1],
Subscript[\[CapitalGamma], NS2]}, {Ca, Cs}, {t, 0, 1}, {x, 0, 1}]


What is the reason for this error? Can Mathematica solve this kind of PDEs? My PDE does have both Neumann and Dirichlet BC at same boundary but for different variables.

Edit1:

I simplified the 2nd term of my first equation. Here is what I wrote:

op1 = D[Ca[t, x], t] -
alpha* (Inactive[Grad][Cs[t, x]/(Km/CB0 + Cs[t, x]), {x}] .
Ca[t, x], {x}] + (Cs[t, x]/(Km/CB0 + Cs[t, x])) Inactive[Div][
Inactive[Grad][Cs[t, x]/(Km/CB0 + Cs[t, x]), {x}] +
Ca[t, x] Inactive[Div][
Inactive[Grad][Cs[t, x]/(Km/CB0 + Cs[t, x]), {x}], {x}]);


I now get the following error:

NDSolveValue::femper: PDE parsing error of 0.5 ((Cs Div$$22707)/(2/5+Cs)+Ca Div$$22709+Grad$$22706.Grad$$22708+Grad$$22708.Grad$$22706). Inconsistent equation dimensions.

What is the issue here?

You could try this:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
Needs["NDSolveFEM"]
(*Enzyme-substrate dyamics-hexokinase parameters*)
Clear["Global'*"];
Clear[Derivative];
Ca0 = 10^(-6);(*Initial Protein Concentration*)d =
180*10^(-6);(*Half-width of channel*)Da = 7.2*10^(-11);
CB0 = 10*10^(-3);
Km = 4*10^(-3);(*Michealis Menten Constant*)beta = \
10;(*Subscript[D,s]/Subscript[D,E]*)gamma = (Ca0/CB0)*(d^2/Da);
k2 = 8000;
alpha = -0.5;
op1 = D[Ca[t, x], t] -
Inactive[Div][{{alpha*(Cs[t, x])/(Km/CB0 + Cs[t, x])}} .
op2 = D[Cs[t, x], t] -
k2*gamma*(Ca[t, x]*Cs[t, x])/(Km/CB0 + Cs[t, x]);
IC = {Ca[0, x] == 1, Cs[0, x] == Exp[-1000 x]};
Subscript[\[CapitalGamma], N1] = NeumannValue[0, x == 0];
Subscript[\[CapitalGamma], N2] = NeumannValue[0, x == 1];
Subscript[\[CapitalGamma], NS1] =
DirichletCondition[Cs[t, x] == 1, x == 0];
Subscript[\[CapitalGamma], NS2] =
DirichletCondition[Cs[t, x] == 0, x == 1];
eqn1 = {op1 ==
Subscript[\[CapitalGamma], N1] + Subscript[\[CapitalGamma], N2],
op2 == 0, IC, Subscript[\[CapitalGamma], NS1],
Subscript[\[CapitalGamma], NS2]};
sol = NDSolveValue[eqn1, {Ca, Cs}, {t, 0, 1}, {x, 0, 1}]


Your first equation, the third term was a bit off. Double check that I did not make a mistake transcribing your equations.

• So the equation you have used is not equivalent to mine. In that process, some of the derivative terms will be ignored. Jul 9, 2021 at 12:57
• @Brownian_Motion can you, please, specify how it is not equivalent to your equation? This will help to aid user21 to improve their answer and, in the end, to better help you. Jul 9, 2021 at 13:08
• @user_21 what is the dot operator you have used in the expression? I think the reason I say the equation is not equivalent is that the term Cs[t,x]/(Km/CB0+Cs[t,x]) is now out of the grad term. There should one more term Ca[t,x]*Grad[Cs[t,x]/(Km/CB0+Cs[t,x])] Jul 9, 2021 at 13:12
• @Brownian_Motion, I see. The problem is that inactive version Grad[f[t, x]*g[t, x], {x}] is no implemented. I you could try to do this manually, though, what happens when you do not use Inactive on this part of the equation? Jul 12, 2021 at 6:55