# DSolve and NDSolve plots are different

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.  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!

• 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 Sin[500 x]]}}  . I think x should be measurd in mm. – Akku14 Jun 5 '18 at 5:09

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] 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] Which now matches the DSolve solution. 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 == 1, u == 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] (*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]] 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}
]

• Thank you so much! This was really helpful! Question, how did you place the two graphs on top of each other? – AhWee Jun 5 '18 at 7:05
• 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"}] – AhWee Jun 5 '18 at 7:15
• @J.Yang sorry, forgot to add Manipulate code. it is there now – Nasser Jun 5 '18 at 7:22

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