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!

Improve this question
$\endgroup$
1
  • $\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, 2018 at 5:09

2 Answers 2

Reset to default
10
$\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}
 ]
Improve this answer
$\endgroup$
3
  • $\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, 2018 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, 2018 at 7:15
  • $\begingroup$ @J.Yang sorry, forgot to add Manipulate code. it is there now $\endgroup$
    – Nasser
    Jun 5, 2018 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]
Improve this answer
$\endgroup$

Your Answer

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

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