2
$\begingroup$

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: enter image description here

(*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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Jun 7, 2018 at 19:54
  • $\begingroup$ @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. $\endgroup$
    – AhWee
    Commented Jun 7, 2018 at 20:44
  • $\begingroup$ 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? $\endgroup$
    – MarcoB
    Commented Jun 7, 2018 at 20:58
  • $\begingroup$ @MarcoB even when removing it, I still get the CoefficientArrays and NDSolve error. I don't need it to use FEM specifically. $\endgroup$
    – AhWee
    Commented Jun 7, 2018 at 21:01

3 Answers 3

Reset to default
1
$\begingroup$

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}]
Improve this answer
$\endgroup$
6
  • $\begingroup$ Although, when I plot that solution, it's just a linear line when I thought it should be nonlinear. Unless I'm missing something? $\endgroup$
    – AhWee
    Commented Jun 8, 2018 at 4:18
  • $\begingroup$ @AhWee Why you think this ode must have a nonlinear solution? $\endgroup$
    – zhk
    Commented Jun 8, 2018 at 4:19
  • $\begingroup$ 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] $\endgroup$
    – AhWee
    Commented Jun 8, 2018 at 4:27
  • $\begingroup$ @AhWee Try this u[x] /. solNDSolve /. x -> 100 you will see that the bc at x=100 is satisfied. $\endgroup$
    – zhk
    Commented Jun 8, 2018 at 4:38
  • $\begingroup$ `Plot[Evaluate[u[x] /. solNDSolve /. x -> 100], {x, 0, 100}, PlotRange -> Full]' correct? In this case, I just get a horizontal line at 1. $\endgroup$
    – AhWee
    Commented Jun 8, 2018 at 5:36
1
$\begingroup$

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}]

enter image description here

Improve this answer
$\endgroup$
1
$\begingroup$

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}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.