Creepy numerical artefacts of a closed-form solution

Question

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[[1]]$, $\alpha /.s[[2]]$ and $\alpha /.s[[3]]$. I plotted each of them as a function of temperature to see which one is real and positive. It seemed, that $\alpha /.s[[2]]$ is the candidate, so I tried to plot (with some realistic constants):

$$ Plot[Re[(k T Log[\alpha /. s[[2]]])/. \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[[2]]])/. {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[[2]]])/. \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[[2]]])/. {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?

Thank you in advance.

Answer 1

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[[2]] 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[[2]] /. parametervalues}, 
  {t, 10^-12, 15/1000, 1/10000}
];

ListPlot[plotdata, PlotRange -> Full, AspectRatio -> 1]

Answer 2

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[[2]]]) /. {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.001`78, 0.04`78, 0.001`78}
  ]
(*
  {-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.]