7
$\begingroup$

I'm trying to solve for a steady state Fick's second law equation with a simple reaction term.

However, when I compare both DSolve and NDSolve they both give very different plots. DSolve Plot NDSolve Plot

Why is that? Is there something wrong with the way I am setting up my code?

Clear["Global`*"]

DiffCo = 1*10^-6; (*Diffusion coefficient*)
k = 0.25;

eqn = (DiffCo*(D[u[x], x, x]) + (k*u[x])) == 0

sol = DSolve[{eqn, DirichletCondition[u[x] == 0, x == 100], 
   DirichletCondition[u[x] == 1, x == 0]}, u, {x, 0, 100}]

sol1 = DSolve[eqn, u, {x, 0, 100}]

Plot[u[x] /. sol, {x, 0, 100}, PlotRange -> Full]

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

Plot[u[x] /. sol2, {x, 0, 100}, PlotRange -> Full]

Any help would be greatly appreciated. Thank you!

$\endgroup$
  • $\begingroup$ Normaly diffusion coefficient is given in m^2/s. Are you sure, that you have to regard the first 100 m for x ? You get something like {{u -> Function[{x}, Cos[500 x] - Cot[50000] Sin[500 x]]}} . I think x should be measurd in mm. $\endgroup$ – Akku14 Jun 5 '18 at 5:09
9
$\begingroup$

My understanding is that when using DirichletCondition, NDSolve will automatically use FEM. So it was using FEM already. The issue seems to be in the grid size used, by default in FEM. since you did not specify one, it used large one and that gave inaccurate result.

k = 1/4;
diffCo = 0.01;
bc = {DirichletCondition[u[x] == 1, x == 0], 
   DirichletCondition[u[x] == 0, x == 100]};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
solNumeric = NDSolve[{eqn == 0, bc},u,{x, 0, 100},Method ->{"FiniteElement"}];

Plot[u[x] /. solNumeric, {x, 0, 100}, PlotRange -> Full]

Mathematica graphics

The smaller the diffCo, the smaller the maxCellMeasure is needed. This manipulate shows this more clearly. When diffCo is large, the default will work OK:

k = 1/4;
diffCo = 1;
bc = {DirichletCondition[u[x] == 1, x == 0], 
   DirichletCondition[u[x] == 0, x == 100]};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
solNumeric = 
  NDSolve[{eqn == 0, bc}, u, {x, 0, 100}, Method -> {"FiniteElement"}];
Plot[u[x] /. solNumeric, {x, 0, 100}, PlotRange -> Full]

Mathematica graphics

Which now matches the DSolve solution.

enter image description here

Added

btw, your ode is of form $\delta u''(x) + k u(x)=0$ where $\delta$ is very small. This calls for boundary layer theory and asymptotic expansion, which is used in these cases, and it gives good result as long as one is close to the point of expansion.

k = 1/4;
diffCo = 1/10^6;
bc = {u[0] == 1, u[100] == 0};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], x];
Plot[Evaluate@solAnalytic, {x, 0, .05}, PlotRange -> Full, 
 ImageSize -> 300]

Mathematica graphics

(*use 50 terms in expansion around 0 *)
sol = AsymptoticDSolveValue[{eqn, bc}, u[x], {x, 0, 50}];
Quiet[Plot[Evaluate@sol, {x, 0, .05}, PlotRange -> Full, ImageSize -> 300]]

Mathematica graphics

Added Manipulate code

Manipulate[
 Module[{k, u, x, bc, eqn, solAnalytic, solNumeric},
  k = 1/4;
  bc = {DirichletCondition[u[x] == 1, x == 0], DirichletCondition[u[x] == 0, x == 100]};
  eqn = diffCo*u''[x] + k*u[x] == 0;
  solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
  solNumeric = NDSolve[{eqn == 0, bc}, u, {x, 0, 100}, 
    Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> maxCellMeasure}];

  Grid[{
    {Plot[solAnalytic, {x, 0, 100}, PlotRange -> Full, ImageSize -> 300]},
    {Plot[Evaluate[u[x] /. solNumeric], {x, 0, 100}, 
      PlotRange -> Full, ImageSize -> 300]}
    }, Frame -> All, Spacings -> {1, 1}, Alignment -> Center]
  ]
 ,

 {{diffCo, 0.001, "diffCo"}, 0.001, 1, 0.01, Appearance -> "Labeled"},

 {{maxCellMeasure, 1, "MaxCellMeasure"}, 0.001, 1, 0.001, Appearance -> "Labeled"},
 TrackedSymbols :> {diffCo, maxCellMeasure}
 ]
$\endgroup$
  • $\begingroup$ Thank you so much! This was really helpful! Question, how did you place the two graphs on top of each other? $\endgroup$ – AhWee Jun 5 '18 at 7:05
  • $\begingroup$ Also, how did you set up your Manipulate plot? When I try, I can't seem to get the graph to change. This is the code I'm using: Manipulate[ Plot[u[x] /. solAnalytic, {x, 0, 100}, PlotRange -> Full], {diffCo, 1*10^-6, 1, Appearance -> "Labeled"}, {MaxCellMeasure, 0, 1, Appearance -> "Labeled"}] $\endgroup$ – AhWee Jun 5 '18 at 7:15
  • $\begingroup$ @J.Yang sorry, forgot to add Manipulate code. it is there now $\endgroup$ – Nasser Jun 5 '18 at 7:22
9
$\begingroup$

Using FEM will certainly give you the same solution as that of DSolve,

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

Plot[u[x] /. sol2, {x, 0, 100}, PlotRange -> Full]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.