Creepy numerical artefacts of a closed-form solution

I'm wishing you a nice day and here is my problem:

I tried to obtain relation between chemical potential and temperature in a semiconductor and I arrived at this equation:

$$2 e^{-\Delta/kT} \left( \frac{m_c^* kT}{2 \pi \hbar^2} \right)^{3/2} \alpha = 2 \left( \frac{m_v^* kT}{2 \pi \hbar^2} \right)^{3/2} \frac{1}{\alpha} + \frac{N_D}{1+2e^{-\varepsilon_D/kT} \alpha}$$ where: $$e^{\mu/kT} \equiv \alpha$$

So basically, Mathematica's task would be to solve this equation with respect to $\alpha$, for a given set of numerical parameters $\Delta, k, m_c^*, m_v^*, \hbar$ and $\varepsilon_D$ find positive real solution and plot a function $\mu = kT \log \alpha$ as a function of $T > 0$. As it seems, the equation is of the third order, so it is solvable in a closed-form. I obtained those solutions:

$$s = Solve[2 e^{-\Delta/kT} \left( \frac{mc kT}{2 \pi \hbar^2} \right)^{3/2} \alpha == 2 \left( \frac{mv kT}{2 \pi \hbar^2} \right)^{3/2} \frac{1}{\alpha} + \frac{N_D}{1+2e^{-\varepsilon d/kT} \alpha},\alpha]$$ Mathematica code:

eq = 2 E^(-delta/(k T))((mc k T)/(2 Pi hbar^2))^(3/2) alpha == 2 ((mc k T)/(2 Pi hbar^2))^(3/2) (1/alpha) + Nd/(1 + 2 E^(-epsilond/(k T)) alpha);
s = Solve[eq, alpha]


I got three solutions, $\alpha /.s[]$, $\alpha /.s[]$ and $\alpha /.s[]$. I plotted each of them as a function of temperature to see which one is real and positive. It seemed, that $\alpha /.s[]$ is the candidate, so I tried to plot (with some realistic constants):

$$Plot[Re[(k T Log[\alpha /. s[]])/. \left\{\Delta -> \frac{114}{100}, k -> 1, mc -> \frac{2}{10} \frac{910938291}{100000000} 10^{-31}, mv -> \frac{2}{10} \frac{910938291}{100000000} 10^{-31}, Nd -> 10^{21}, \hbar-> \frac{105457173}{100000000} 10^{-34}, \varepsilon d -> \frac{109}{100} \right\}],{T,0,1/10}]$$ ($T, \Delta, \varepsilon d$ and $\mu$ itself are in electronvolts, so mathematica does not have to deal with too large or too small numbers)

Code:

Plot[Re[(k T Log[alpha /. s[]])/. {delta -> 114/100, k -> 1,
mc -> (2/10)(910938291/100000000) 10^-31,
mv -> (2/10)(910938291/100000000) 10^-31,
Nd -> 10^21, hbar-> 105457173/100000000 10^-34, epsilond ->
109/100}],{T,0,1/10}]


So at first it replaces all constants with their numerical values except the temperature and then it tries to plot that. Even the result has to be real, I had to add Re[] to overcome issue with non-zero imaginary part and the result is quite strange. It looks promising, I think the function does what chemical potential would do, it goes from $\varepsilon d$ at $T = 0$ down to $\Delta/2$ at higher temperatures, but there are strange spikes all over the region for $T < 0.03$. It is very strange, because the solution is obtained analytically via Solve, and the function has no reason to behave like that (those are not a high oscillation spikes, it's rather a numerical artefact).

Then I tried to run ContourPlot on the first equation and I got the same spiky result.

Next thing I tried was to create a list of values of $\mu$ for $0 < T < 0.04$ to see where exactly does the oscillation occur:

$$Table[Re[(k T Log[\alpha /. s[]])/. \left\{\Delta -> \frac{114}{100}, k -> 1, mc -> \frac{2}{10} \frac{910938291}{100000000} 10^{-31}, mv -> \frac{2}{10} \frac{910938291}{100000000} 10^{-31}, Nd -> 10^{21}, \hbar-> \frac{105457173}{100000000} 10^{-34}, \varepsilon d -> \frac{109}{100}, T -> i \right\}],{i,0.001,0.04,0.001}]$$

Code:

Table[Re[(k T Log[alpha /. s[]])/. {delta -> 114/100, k -> 1,
mc -> (2/10)(910938291/100000000) 10^-31,
mv -> (2/10)(910938291/100000000) 10^-31,
Nd -> 10^21, hbar-> 105457173/100000000 10^-34, epsilond ->
109/100, T -> i }],{i,0.001,0.04,0.001}]


...and I'm left with:

{1.07225, 1.01546, 1.03827, 0.937135, 1.00205, 0.984445, 0.950138, 0.927387, 0.931658, 0.910465, 0.892761, 0.87479, 0.856831, 0.809198, 0.820963, 0.771908, 0.755787, 0.767006, 0.711472, 0.724181, 0.706065, 0.687321, 0.626185, 0.651721, 0.633072, 0.614729, 0.597576, 0.569833, 0.560556, 0.572135, 0.569195, 0.570001, 0.569994, 0.570001, 0.569991, 0.569992, 0.570002, 0.57, 0.57, 0.57}

You can clearly see the oscillations, which is strange, and in fact, creating list with finer step and then ListPlotting this with respect to the same list with $T$ in same range and step results into a similar spiky mess I obtained before. It does not matter how fine I choose the step, always more spikes show up, which in my opinion means, that this is a numerical issue, not a feature of a closed-form solution.

My question is: how should I attack this problem (which is actually solvable in closed-form, but I won't resist a numerical solution)? How can I get a smooth plottable solution for $\mu$ which I can use further both analytically and/or numerically?

• Could you provide your equations and parameter values as code blocks, rather than in mathematical notation? These formatting instructions may be helpful. I'd like to give your problem a look, but I don't want to re-type your equations by hand... – MarcoB May 4 '15 at 22:06
• This will do? s = Solve[([Alpha]*((kmcT)/[HBar]^2)^(3/2))/(Sqrt*E^([CapitalDelta]/(kT))*Pi^(3/2)) == Nd/(1 + (2*[Alpha])/E^([Epsilon]d/(kT))) + ((kmvT)/[HBar]^2)^(3/2)/(Sqrt*Pi^(3/2)*[Alpha]), [Alpha]] and plot: Plot[Re[kTLog[[Alpha] /. s[]] /. {[CapitalDelta] -> 114/100, k -> 1, mc -> ((2/10)*(910938291/100000000))/10^31, mv -> ((2/10)*(910938291/100000000))/10^31, Nd -> 10^21, [HBar] -> 105457173/100000000/10^34, [Epsilon]d -> 109/100}], {T, 0, 0.05}] is this right? when i paste it back to mathematica, it is good... – user16320 May 4 '15 at 23:32
• Yes that will do for a start, thank you. I will probably get rid of the special symbols though. I think I am on to something. Take a look at the answer I just posted and let me know. – MarcoB May 5 '15 at 0:43
• You can format inline code and code blocks by selecting it and clicking the {} button above the edit window. (Inline code is placed between backticks; code blocks are indented four spaces.) The edit window help button ? is also useful for learning how to format your questions and answers. -- The code in your comment is missing the backslashes to the Greek letter names. Please add the code to your question by editing it. – Michael E2 May 5 '15 at 1:48
• People here generally like users to post code as Mathematica code instead of TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. – Michael E2 May 5 '15 at 1:51

I think your problem originates within the Plot function.

Plot has some amazing internal magic, especially relating to adaptive sampling over the plot range, but this magic goes astray sometimes. In those cases, I attempt to construct a table of values of the function to be plotted, then use ListPlot on this table to obtain what I couldn't get out of Plot.

Let me first rewrite your equation and parameter values in a more friendly fashion, avoiding special symbols or uppercase variables. Avoiding special symbols simply makes for simpler manipulation and transfer. Uppercase variable names or user functions are deprecated in Mathematica, since they could conflict with built-in functions.

equation =
2 Exp[-delta/(k t)]((mc k t)/(2 Pi hbar^2))^(3/2) alpha ==
2 ((mv k t) / (2 Pi hbar^2))^(3/2) 1/alpha +
nD / (1 + 2 Exp[-(epsilond / (k t))] alpha);

parametervalues={
delta->114/100,
k->1,
mc->910938291/50000000000000000000000000000000000000,
mv->910938291/50000000000000000000000000000000000000,
nD->10^21,
hbar->21/200000000000000000000000000000000000,
epsilond->109/100
}


Not as pretty as your version, but much more machine readable! As you can see, I kept all numbers at arbitrary precision.

We can then solve the equation for alpha:

analyticalsolutions = Solve[equation, alpha];


We get three solutions. The second one analyticalsolutions[] is the one with physical significance, as you mentioned as well. I will explicitly construct a table of values $(t,\ kt\ log\ \alpha)$ at arbitrary precision and store it in plotdata, then plot this list of data points using ListPlot.

Two caveats here: 1) you can get arbitrarily close to $t=0$, but the alpha function is not defined there, so trying to evaluate starting from $t=0$ will give rise to all sorts of errors; 2) use PlotRange->Full to make sure that your plotting function is not cutting out portions of your plot.

plotdata = Table[
{t, k t Log[alpha] /. analyticalsolutions[] /. parametervalues},
{t, 10^-12, 15/1000, 1/10000}
];

ListPlot[plotdata, PlotRange -> Full, AspectRatio -> 1] • Thank you very much. I tried to make a table of values as well, but I ended with a mess of some sort-it produced the right plot, but also a lot incorrect points around the curve, even when I tried to rationalize the values of numerical parameters. I see that you rationalized the values too, and your graph is correct, from what I know about real mu(T) dependance. Can you tell me, what can be the reason behind those redundant incorrect points I got from the ListPlot? I think it is related to Michael's answer, but you didn't define precision in your parameters, although the plot is correct. – user16320 May 5 '15 at 12:14

One way to check the numerics is to use high-precision arbitrary-precision numbers. At precisions under 62 digits for i, the expression does not evaluate to a number at all values of i in the OP's table. One has to use 78 digits to get results that are all accurate to machine precision, but if converted first to MachinePrecision, the differences at the tail end of the table would come out to be zero due to their relative magnitude. The values and differences below were computed at an initial precision of 78 digits and converted to machine precision. All of the differences are negative, showing that the solution do not really oscillate as they appear to do at machine precision.

OPsol2 = Re[(k T Log[alpha /. s[]]) /. {delta -> 114/100, k -> 1,
mc -> (2/10) (910938291/100000000) 10^-31, mv -> (2/10) (910938291/100000000) 10^-31,
Nd -> 10^21, hbar -> 105457173/100000000 10^-34, epsilond -> 109/100, T -> i}];
N @ Differences @ Table[OPsol2,
{i, 0.00178, 0.0478, 0.00178}
]
(*
{-0.0698337, -0.0706186, -0.0711283, -0.0715073, -0.0718093,
-0.0720605, -0.0698803, -0.0054057, -2.00043*10^-6, -3.35825*10^-9,
-1.78541*10^-11, -2.24966*10^-13, -5.49761*10^-15, -2.2587*10^-16,
-1.40545*10^-17, -1.22506*10^-18, -1.40944*10^-19, -2.04389*10^-20,
-3.60257*10^-21, -7.49735*10^-22, -1.79951*10^-22, -4.88669*10^-23,
-1.47777*10^-23, -4.91123*10^-24, -1.77394*10^-24, -6.89868*10^-25,
-2.8654*10^-25, -1.26242*10^-25, -5.86458*10^-26, -2.85787*10^-26,
-1.45433*10^-26, -7.69811*10^-27, -4.22369*10^-27, -2.39471*10^-27,
-1.3992*10^-27, -8.40464*10^-28, -5.17876*10^-28, -3.26703*10^-28,
-2.10641*10^-28}
*)


An alternative approach would be to make a table of exact value and use N[_, prec], with a specific number for prec, to calculate the values.

N[Differences@
Table[OPsol2,
{i, 1/1000, 4/100, 1/1000}],
6]
(*  *looks* like the above *)


Applying these approaches to Plot:

Plot[OPsol2 /. i -> SetPrecision[i0, 78], {i0, 0.001, 0.04}]
Plot[N[OPsol2 /. i -> SetPrecision[i0, Infinity], 6], {i0, 0.001, 0.04}]


The second method is noticeably slower. Or one could apply ListPlot to the basic tables above (without the Differences).

[Updated: The solution s` is from the OP's code.]

• Thank you, When I wanted to work further with the function mu(T) I had to specify precision because results were inaccurate. This helped too. – user16320 May 5 '15 at 14:11
• @user16320 You're welcome. – Michael E2 May 5 '15 at 14:12