Ask Your Question

Fermi-Dirac integral of half order

asked 2013-09-06 14:20:20 +0200

wlp2 gravatar image

updated 2024-04-14 16:06:34 +0200

FrédéricC gravatar image

I am trying to implement a 1D model for semiconductor pn-junctions. This involves evaluating the Fermi-Dirac integral of half order, which can only be done numerically. I saw that the GSL library has an algorithm to perform this integration, but I can't figure out how to call it in Sage, is this possible? If not, what are some alternatives to evaluate this integral? Thanks in advance.

edit retag flag offensive close merge delete

4 Answers

Sort by » oldest newest most voted

answered 2013-11-07 23:10:40 +0200

Mohan gravatar image

Please try fermi.m from matlab site matlab file exchange. This will give very accurate values(14 digits minimum) for any half-order Fd integral.

edit flag offensive delete link more

answered 2013-09-27 02:32:26 +0200

wlp2 gravatar image

Sorry it took me so long to respond, life circumstances intervened.

First, sorry if the question wasn't clear enough. This is the first time that I have posted here so I am unfamiliar with the level of expertise in the community, which makes it difficult to write a specific, targeted question.

My original question had two parts:

  1. Can I use Sage to run the GSL algorithm to evaluate the Fermi-Dirac integral of half-order? I am only interested in the solution of this specific integral, not the solution of any other Fermi-Dirac integral. Other orders have analytical solutions, but the Fermi-Dirac integral of half order does not; it must be done numerically. Here is one reference asserting the same claim (click on it, press ctrl+f and type in numerical to get to the specific sentence):

  2. If I cannot use the GSL function, what other options are there?

The first has gone unanswered. The second has been soundly answered and I am grateful for your response! Barring an affirmative answer to the first question, I will use the numerical_integral command to evaluate this integral.

I am, however, troubled by your comments in general. I like SAGE because it is a free and open-source alternative to very expensive software packages like MM. Any work that I do in MM I lose as soon as I don't work somewhere that maintains a MM license, which is almost anywhere other than academia. Suggesting that I use MM seems to go against the stated purpose of SAGE and it's creator's original intent. I understand that in certain circumstances it is a more capable tool. That's fine, but SAGE includes at least 3 packages dedicated to numerical computation (octave, numpy, and scilab). I thought I would at least get a suggestion about how to use one of those tools as an alternative to the GSL algorithm before running to MM with my tail between my legs.

Further, your suggestion that MM "knows" how to perform this integral analytically and does so is an egregious error. The Fermi-Dirac integral of negative half order (used in your example) must also be done numerically (quick google search). Even the answer given by MM, a number and not an expression, suggests that it is not employing a symbolic solution, but is indeed switching to a numerical method. It's equally possible that MM is programmed to recognize Fermi-Dirac integrals and respond accordingly. It should be noted that I found no evidence in the MM documentation to suggest that either one of our answers is correct. This is the core problem of the MM black box, and another reason why SAGE is valuable. In this case, the Maxima engine in SAGE is outperforming MM by providing a more accurate and predictable (although useless) outcome. The only thing the MM result confirms is that the numerical_integral method in SAGE provides the same answer. This suggests that MM is correct because it can be compared against ... (more)

edit flag offensive delete link more

answered 2013-09-27 04:15:03 +0200

Mark gravatar image

updated 2013-09-27 06:45:36 +0200

Further, your suggestion that MM "knows" how to perform this integral analytically and does so is an egregious error. The Fermi-Dirac integral of negative half order (used in your example) must also be done numerically (quick google search). Even the answer given by MM, a number and not an expression, suggests that it is not employing a symbolic solution, but is indeed switching to a numerical method

A) unrelated to sage: as I already said in my previous post, I'd be careful making all-out statements. (Even more so based on 'quick' google seraches). So once again, to falsify your claims:

[1] The Analytical Evaluation of the Half-Order Fermi-Dirac Integrals Jerry A. Selvaggi, and Jerry P. Selvaggi The Open Mathematics Journal, 2012, 5, 1-7

[2] Solutions to the Fermi-Dirac Integrals in Semiconductor Physics Using Polylogarithms M. Ulrich, W. Seng, P. Barnes Journal of Computational Electronics 1: 431–434, 2002

and several other papers. You should also realize, that there is a difference between scientific papers, like the preceding, which have undergone peer review, and some unpublished and unrefereed 'notes', which pop up on the internet - like the one you point to @

From the preceeding, I guess it is quite clear, where or to whom to apply vocabulary like 'egregious'. Moreover, it is Ref. [2] which tells you that things work down even to y>-1. Just pls. read the papers before posting.

B) related to sage: MM is not '..programmed to recognize Fermi-Dirac integrals and respond accordingly..'. Basically it uses eqn. (5) of Ref. [2], which is a paper free to be read with results free to be used by anyone - also outside of academia. Obviously MM has simply cared to include this result into its code - while sage has not. As simple as that.

with my tail between my legs

C) related to this forum I like this forum because it tends to be free of such statements. Pls. keep it that way.

edit flag offensive delete link more

answered 2013-09-07 05:53:21 +0200

Mark gravatar image

updated 2013-09-07 06:02:33 +0200

First, when asking questions, it might speed up a potential answer if you would care to give enough detail about your problem ... so, I only guess that this is what you want $$F(y,\mu) = \int_0^\infty \frac{x^{y}}{e^{x+\mu}+1} dx$$

Second, I'd be very careful about www-ing statements like

"...which can only be done numerically..."

since, to my recollection $F(y,\mu)$ is just another way of writing the polylog for all $y>-1$, i.e down to those values of $y$ where the singularity of the integrand remains harmless. Therefore, at least I do not understand your claim, since in 1D the DOS in a semiconductor is $\propto x^{-1/2}$ (which you also did not specify that in your question). So it can be done analytically ...

Third, if still you really want to do this numerically in sage ... just look here and e.g. do this

def F(y,mu):
    ul = 100+mu # increase the 100 if you really need
    return numerical_integral(x^y/(e^(x+mu)+1),0,ul)

(1.072154920510084, 8.286087609477639e-07)

Finally, and only to my very humble understanding: sage is a chain of existing tools, i.e. it is kind of as strong as it's weakest link. Since in plain sage things like Maxima are the link for many of those calculations which 90%-99% of non-mathematician might call 'symbolic', plain sage frequently is of little help if compared to MM or Maple once you come to such 'symbolic' problems.

Applied to your case this means:

sage: var('y,mu')
(y, mu)
sage: assume(y+1>0)
sage: assume(y,'integer')
sage: integral(x^y/(e^(x+mu)+1),x,0,oo)
integrate(x^y/(e^(mu + x) + 1), x, 0, +Infinity)

so, nothing ... or if you try special values

sage: integral(x^(-1/2)/(e^(x+mu)+1),x,0,oo)
integrate(1/((e^(mu + x) + 1)*sqrt(x)), x, 0, +Infinity)

sage: integral(x^1/(e^(x+mu)+1),x,0,oo)
polylog(2, -e^mu) + limit(1/2*x^2 - x*log(e^(mu + x) + 1) - polylog(2, -e^(mu + x)), x, +Infinity, minus)

again, not too helpful. Maybe there are 'tricks' to get further but I don't know.

If however you use MM then, out of the box:

sage: mathematica('Integrate[x^a/ (Exp[x + b] + 1), {x, 0, Infinity}, Assumptions -> {Re[a] > -1}]')
-(a*Gamma[a]*PolyLog[1 + a, -E^(-b)])

so MM 'knows' that your integral can be done analytically and provides the answer, and

sage: mathematica('-(a*Gamma[a]*PolyLog[1 + a, -E^(-b)])/. {a -> -1/2, b -> 0} // N')

which is the numerical value we got previously.

In conclusion: To answer if that implies you should better be using MM, there are more powerful experts around here.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2013-09-06 14:20:20 +0200

Seen: 2,550 times

Last updated: Nov 07 '13