2
$\begingroup$

I am something of a novice with Mathematica. I have written the following code;

Clear[band, en, w, fermi, k, T, S, a]
k = 86*10^-6;
T = 4000;
band[en_, w_] := 10000 Exp[-en^2/2 w^2];
fermi[en_, ef_] := 1/(Exp[(en - ef)/(k T)] + 1);
Plot[band[en, w] /. w -> 1, {en, -4, 4}]
Plot[fermi[en, ef] /. ef -> 0, {en, -4, 4}]
Plot[{band[en, w] /. w -> 1, (band[en, w]*fermi[en, ef]) /. {ef -> 0, w -> 1}}, {en,-4,4}]
S[ef_] := -k NIntegrate[(band[en, w]*fermi[en, ef]*Log[fermi[en, ef]]) /. w -> 1, {en, -Infinity, Infinity}]
Plot[S[ef], {ef, -4, 4}, PlotRange -> All]
a[ef_] := D[S[ef], ef]
Plot[a[ef], {ef, -4, 4}]

It works very nicely until I try to plot a, which is simply the derivative of S wrt ef - even though I can plot S vs ef a few lines above. Plotting a gives error messages along the lines of;

General::ivar: "-3.83657 is not a valid variable."

This is clearly something to do with NIntegrate evaluating too early, or not being a function of ef in the way I would expect it to be, or similar. But I have no idea how to fix it.

I would be very grateful for any suggestions anyone could give. I think this issue is showing up my fundamental lack of understanding how Mathematica works, and so replies along the lines of "You should read and understand such-and-such article..." would also be very welcome!

Many thanks,

Nick

$\endgroup$
  • $\begingroup$ Many thanks belisarius, that's a useful package to know about. $\endgroup$ – NLambert Dec 5 '13 at 16:31
2
$\begingroup$

Use a numerical derivative:

Clear[band, en, w, fermi, k, T, S, a]
<< NumericalCalculus`
k = 86*10^-6;
T = 4000;
band[en_, w_] := 10000 Exp[-en^2/2 w^2];
fermi[en_, ef_] := 1/(Exp[(en - ef)/(k T)] + 1);
S[ef_?NumericQ] := -k NIntegrate[(band[en, w] fermi[en, ef] Log[fermi[en, ef]]) /. w -> 1, 
                   {en, -Infinity, Infinity}]

Plot[{S[x], ND[S[ef], ef, x]}, {x, -4, 4}]  

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.