# Diffusion equation with concentration dependent diffusivity, NDSolve has problem at C0> 10^20

I am using NDSolve to solve the diffusion equation:

D[Cp[x, t], t] == D[Dc D[Cp[x, t], x], x]


where Dc is a function of Cp: Dc = Dc(Cp). The analytical expression for Dc is given. Here is my program.

Clear[DPi, DP, DIi, DI, DPI, DPVi, DPV, Dtotal, Tc, ni, n, Cp, nn, n1, n2, h]

q = 1.6 10^-19;kB = 1.38 10^-23;

ni[Tc_]:=1.8 10^21 Exp[-((0.66 q)/(kB Tc))];

n[Cp_, Tc_]:= 0.5 (Cp + Sqrt[Cp^2 + 4 ni[Tc]^2]);

h[Cp_, Tc_]:= 1 + n[Cp, Tc]/(2 ni[Tc]) ((n[Cp, Tc]/(2 ni[Tc]))^2 + 1)^-0.5;

DP[Cp_, Tc_]:=4.9 10^-5 Exp[-2.1 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^0.5 h[Cp, Tc];

DI[Cp_, Tc_]:=0.23 Exp[-2.6 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^-1.8 h[Cp, Tc];

DPI[Cp_, Tc_]:= (DP[Cp, Tc] DI[Cp, Tc])/(DP[Cp, Tc] + DI[Cp, Tc]);

DPV[Cp_, Tc_] :=7.6 10^4 Exp[-5.2 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^2 h[Cp, Tc];

Dtotal[Cp_, Tc_] := DPI[Cp, Tc] + DPV[Cp, Tc];

Dc = Dtotal[Cp, 900 + 273] /. Cp -> Cp[x, t];

LogLogPlot[Dtotal[Cp, 850 + 273], {Cp, 10^15, 7 10^20},
PlotRange -> {{10^16, 10^21}, {10^-15, 10^-12}}, Frame -> True]

C0 = 3 10^20;

sol = NDSolve[{D[Cp[x, t], t] == D[ Dc D[Cp[x, t], x], x],
Cp[x, 0] == C0 UnitStep[-x], Cp[0, t] == C0,
Cp[3. 10^-4, t] == 0}, Cp, {x, 0, 3. 10^-4}, {t, 0, 5. 10^3},
Method -> {"MethodOfLines"}, MaxSteps -> 10000, PrecisionGoal -> 2];

LogPlot[Evaluate[Cp[x, 2400] /. sol], {x, 0, 1 10^-4},
PlotRange -> {{0, 1 10^-4}, {10^15, 1 10^21}}, Frame -> True ]


At low concentration from 10^18-10^20, I dont have any problem, but with concentrations higher than 10^20, the solution looks strange.

Could you please help?

-
As a general rule, it's always good to keep the numbers you're working with to values ~unity. Many of the algorithms that Mathematica uses under the hood can fail if the arguments passed to them are too large/small. –  Guillochon Jan 14 at 21:53
Make everything in the input exact and you can do e.g. sol = NDSolve[{D[Cp[x, t], t] == D[Dc D[Cp[x, t], x], x], Cp[x, 0] == C0 UnitStep[-x], Cp[0, t] == C0, Cp[3 10^-4, t] == 0}, Cp, {x, 0, 3 10^-4}, {t, 0, 5 10^3}, Method -> {"MethodOfLines"}, MaxSteps -> 10000, PrecisionGoal -> 2, WorkingPrecision -> 40]; The gap in the plot is from taking logs of negatives and thus getting an imaginary part. This might or might not indicate a problem in the model and/or solution thereof. –  Daniel Lichtblau Jan 15 at 0:50
Thank you, Guillochon and Daniel ! @ Guillochon, if I put concentration to be 1, then I have to rescale diffusivity to be very small, say 10^-30. @ Daniel, right, the solution discontinuity is from negative value in log polt, but why? I tried to increase WorkingPrecision as you suggested, but for conentration higher than 2x10^20, the discontinuity still exist. –  Pham Nam Jan 18 at 21:11
add comment