# Help in advection-dispersion equation using NDSolveValue

Dz = 0.0000738739; as = 2.21622*10^-6;
Dz*D[u[z, t], z, z] - as*D[u[z, t], z] == D[u[z, t], t]


The boundary and initial conditions are:

u[z, 0] == 0.10,
u[0, t] == 0.30, , (D[u[z, t], z] /. z -> Infinity) == 0


As i wanted to solve numerically, i changed "infinity" to 10^5 and tried to solve the pde numerically with the following code:

 pde = {Dz*D[u[z, t], z, z] - as*D[u[z, t], z] - D[u[z, t], t] == 0,
u[z, 0] == 0.1,
u[0, t] ==
0.30 + (0.10 - 0.30)*Exp[-100000*t], (D[u[z, t], z] /.
z -> 10^5) == 0}

solution = NDSolveValue[pde, u, {z, 0, 100000}, {t, 0, 100000}]

Plot[solution[0.6, t], {t, 0, 100000}, PlotRange -> All]


The Exp[-100000*t] term was used to not give inconsisting boundary condition values. The problem is that when i plot the solution, the results are totally wrong comparing to the reality and also to the analitic solution of this pde. I've tried to use the NeumannValue and DirichletCondition functions, i get better results, close to the analitical solution, but still very different from the analitic solution. And i don't know why the results change when i use the NeumannValue and DirichletCondition functions, in my mind they should be the same, as i wrote the code correctly. Here is the analitic solution of the proposed PDE.

z1[z_, t_, as_, Dz_] := (z + as*t)/(2*Sqrt[Dz*t])

z2[z_, t_, as_, Dz_] := (z - as*t)/(2*Sqrt[Dz*t])

A[z_, t_, as_, Dz_] :=
1/2 (Erfc[z1[z, t, as, Dz]] + Exp[((as*z)/Dz)]*Erfc[z2[z, t, as, Dz]])

u[z_, t_, u0_, ui_, as_, Dz_] := ui + (u0 - ui)*A[z, t, as, Dz]

Plot[{u[0.6, t, 0.30, 0.10, as, Dz]}, {t, 0.01, 1000000}, PlotRange -> All]


Can anybody, please, help me to find the solution to this problem? Thanks!

• Hi, I would help if you add Dz and as to the code and delete the code that that has syntax errors in it. That will make it more likely that someone will that the time to look at your issue. Use the eidit button under your post for that. – user21 Mar 18 '19 at 6:45
• Hello. Thanks for your comment. The Dz value is literal on the first equation, but, on the third, i substituted its value. And in the text i also explain and put its numerical value. The code is exactly as i posted and it runs perfectly. Thanks. – shewlong Mar 18 '19 at 18:11
• ……The analytic solution is incorrect. Just try Clear[as, Dz]; func = {z, t} \[Function] u[z, t, u0, ui, as, Dz]; error = Dz*D[u[z, t], z, z] - as*D[u[z, t], z] - D[u[z, t], t] /. u -> func // FullSimplify Plot3D[Block[{u0 = 1, ui = 3, as = 1, Dz = 1}, error] // Evaluate, {t, 0, 1}, {z, 0, 1}] – xzczd Mar 19 '19 at 6:57
• Even if it is incorrent, is there any way to explain why the ndsolvevalue gives inconsisted values comparing to what is expected? – shewlong Mar 19 '19 at 13:31
• Let me be clearer. I found the analytic solution in a scientific paper, and it may be wrong. But that's not exactly the point, it was just to illustrate the expected process. The numerical solution should, when plotted, be able to show what we expected: different points converging to a limit of 0.30. If i change "as" and "dz", it doesn't change the plotted results. If i change z, it doesn't change neither. I don't know what is happening to this solution. It is going too fast for every z point. I've tried everything that was in my abilities. I appreciate the time looking at it. – shewlong Mar 19 '19 at 19:39

First of all, the given analytic solution is incorrect:

Clear[as, Dz];
func = {z, t} \[Function] u[z, t, u0, ui, as, Dz];
error = Dz*D[u[z, t], z, z] - as*D[u[z, t], z] - D[u[z, t], t] /. u -> func;
Plot3D[
Block[{u0 = 1, ui = 3, as = 1, Dz = 1}, error] // Evaluate, {t, 0, 1}, {z, 0, 1},
PlotPoints -> 50] Then the problem is related to this one. The numeric solution given by NDSolve is inaccurate, because the domain of z isn't properly chosen. Notice the values of as and Dz are very small, this suggest the variation of the solution will only happen in a small area near z == 0, and the default spatial grid (only 25 points in z direction) is not dense enough to capture the variation. So the simplest solution is to chose a smaller number for approximation of $$\infty$$:

Dz = 0.0000738739; as = 2.21622*10^-6;
tend = 1/100; zend = 1/200;
pde = {Dz D[u[z, t], z, z] - as D[u[z, t], z] - D[u[z, t], t] == 0, u[z, 0] == 0.1,
u[0, t] == 0.30 + (0.10 - 0.30) Exp[-100000 t], (D[u[z, t], z] /. z -> zend) == 0};

solution = NDSolveValue[pde, u, {z, 0, zend}, {t, 0, tend}];

Plot3D[solution[z, t], {z, 0, zend}, {t, 0, tend}, PlotRange -> All] Hello @xzczd thank you for the time you spent in this problem. I have some points to take into consideration.

Fortunately, yesterday i was struggling with this issue and i came up to a solution very similar to yours, because of the discussion we previously had. I realized that the domain was too big and, as the initial conditions are constant, Mathematica thinks that the error is negligible. So i increase the refinement of the grid until it stabilizes. Here is the code:

u0 = 0.30;
ui = 0.10;
ContDireito = 1000;
Dz = 0.00007387387387387386;
as = 2.2162162162162158*^-6;

EqNumerica = {Dz*Derivative[2, 0][u][z, t] -
as*Derivative[1, 0][u][z, t] - Derivative[0, 1][u][z, t] == 0,
u[z, 0] == If[z == 0, u0, ui],
u[0, t] == u0, (D[u[z, t], z] /. z -> ContDireito) == 0}

SolutionNew =
NDSolveValue[EqNumerica, u, {z, 0, 1}, {t, 0, 10000000},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 2000}}]


After that, i plotted SolutionNew: Of course, this is not the point that deserves an answer. What i would like to discuss is about the analytical solution. First of all, the right-Boundary condition is in infinite, so we need to make sure that the domain is large enough comparing to the points we analyze, that's why i put the domain in 1000. After seeing te convergence (which we could calculate, but it was easier for me to just keep trying until the plotting remains constant) i compared it to the analytic solution, which, according to your code, is obviously wrong. Here is the comparison: So, what's the point? I couldn't see the proof of this analytic solution, but i suppose it uses some approximations in its deduction, so that for small values of as and Dz it works perfectly. The values of as and Dz that you used were a little bit unrealistic for the studied case. I appreciate If you could do any further collaboration about that. Thanks for the valuable discussion.

• Interesting. Sadly I can't figure out what approximation has been used here. The only thing I know is, the exact analytic solution is the inverse Laplace transform of (E^(1/2 (as/Dz - Sqrt[as^2 + 4 Dz s]/Dz) z) (u0 - ui))/s + ui/s, which unfortunately cannot be calculated with InverseLaplaceTransform. – xzczd Mar 20 '19 at 14:05
• Hello @xzczd if you are interested in the solution of this equation, i have some interesting papers which derive its analytic solution, but not with the exactly right-side boundary condition. Those are: pubs.usgs.gov/pp/0411a/report.pdf, journals.aps.org/pr/abstract/10.1103/PhysRev.23.412. Both are beautiful papers that present a very similar solution, but not with the exact right-side boundary condition. Thanks for the debate! – shewlong Mar 21 '19 at 3:45