# Using DSolve and NDSolve for 1D steady-state diffusion equation with a reaction term

I am trying to solve for a 1D steady-state diffusion equation (Fick's Second Law/Heat Equation) with a reaction term. I wanted to compare using both DSolve and NDSolve.

However, when running the code, I get multiple errors:

(*L+rUnderoverscript[\[Equilibrium], d1, k1]R*)
k1 = 0.00193
d1 = 0.00700
K1 = d1/k1
NT = 1.7
(*Total Receptors(NT) = r+R*)
R = (NT*u[x])/(K1 + u[x])
r = NT - R
RL = (d1*R) - (k1*u[x]*r)

diffCo = 0.0001; (*Diffusion coefficient*)
bc = {DirichletCondition[u[x] == 1, x == 0], DirichletCondition[u[x] == 0, x == 100]};
eqn = diffCo*u''[x] + RL == 0;

solDSolve = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}]

solNDSolve =
NDSolve[{eqn, bc}, u, {x, 0, 100},
Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> 0.001}]


Am I setting something up incorrectly? Thank you.

• Your code does not return those errors. In particular, notice that you use u[x, t] in the definitions, and you mix that with u[x] in thee equations and boundary conditions. Please amend the code to show exactly what you used. – MarcoB Jun 7 '18 at 19:54
• @MarcoB Sorry about that. I initially tried solving with time dependency before switching to steady state and forgot to adjust those terms. I fixed them and am getting the errors I mentioned. – AhWee Jun 7 '18 at 20:44
• Your problems seem to stem from the choice of the FEM method. NDSolve seems to work fine if you remove the Method specifications. Do you have to / want to use FEM specifically? – MarcoB Jun 7 '18 at 20:58
• @MarcoB even when removing it, I still get the CoefficientArrays and NDSolve error. I don't need it to use FEM specifically. – AhWee Jun 7 '18 at 21:01

The error message CoefficientArrays::poly signifies that FEM only works for linear ODE's, which is a bit strange to me. Or maybe I am not understanding it.

This type of error message has been reported quite a few times but for PDE's not ODE's.

As suggested in the comment you should remove the Method but still you and I (V. 11.2) are getting the same error. I think, it's because the way the boundary conditions are defined (I suspect, NDSolve is calling FEM by default). If we consider the boundary conditions in the more conventional way, everything seems to work just fine,

solNDSolve = NDSolve[{eqn, u[0] == 1, u[100] == 0}, u, {x, 0, 100}]

• Although, when I plot that solution, it's just a linear line when I thought it should be nonlinear. Unless I'm missing something? – AhWee Jun 8 '18 at 4:18
• @AhWee Why you think this ode must have a nonlinear solution? – zhk Jun 8 '18 at 4:19
• Unless I'm misunderstanding I expected the concentration at x=0 to decay/deplete at the rate RL. I tried making the other end of the boundary much larger in the hopes to see at what point did the concentration reach 0. However, when I try this boundary condition, it just moves the 0 point to whichever x, I choose to make it. Also, the y axis along the plot just shows "1." Unless I'm setting up the plot wrong. I put Plot[Evaluate[u[x] /. solNDSolve], {x, 0, 100}, PlotRange -> Full] – AhWee Jun 8 '18 at 4:27
• @AhWee Try this u[x] /. solNDSolve /. x -> 100 you will see that the bc at x=100 is satisfied. – zhk Jun 8 '18 at 4:38
• Plot[Evaluate[u[x] /. solNDSolve /. x -> 100], {x, 0, 100}, PlotRange -> Full]' correct? In this case, I just get a horizontal line at 1. – AhWee Jun 8 '18 at 5:36

Things are very simple. Your differential equation reduces to (u''[x] == 0 .

k1 = 0.00193;
d1 = 0.00700;
K1 = d1/k1;
NT = 1.7;
R = (NT*u[x])/(K1 + u[x]);
r = NT - R;
RL = (d1*R) - (k1*u[x]*r);
diffCo = 0.0001;

bc = {u[0] == 1, u[100] == 0};
eqn = Rationalize[diffCo*u''[x] + RL == 0, 0] // Simplify

(*   (u''[x] == 0   *)

sol = First@DSolve[{eqn, bc}, u, x]

(*   {u -> Function[{x}, (100 - x)/100]}   *)

ndsol = NDSolve[{eqn, bc}, u, {x, 0, 100}]

Plot[u[x] /. ndsol, {x, 0, 100}]


Use Simplify on the equations and "MaxCellMeasure" (as a String) and things work out as expected.

solNDSolve =
NDSolve[{eqn // Simplify, bc}, u, {x, 0, 100},
Method -> {"FiniteElement",
"MeshOptions" -> MaxCellMeasure -> 0.001}]
`