Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 0 down vote favorite

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?

share|improve this question
1  
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
1  
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

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.