Inverse Fermi Dirac Integral using Octave/Matlab

I need to obtain the inverse of the fermi dirac integral of order -1/2 in my code. This is what I do. I create a file inv_fermi_dirac_mhalf.m

function y=inv_fermi_dirac_mhalf(fval,init_guess)
global fvalfermi;
%get length of vector
Nel=length(fval);

%initialize
y=zeros(Nel,1);

%loop over all data
for nloop=1:Nel
fvalfermi=fval(nloop);
y(nloop)=fzero('myfunfermi',init_guess(nloop));
end
end

The loop enable the function to take vector data. It uses optimization function fzero to obtain the inverse of the fermi dirac integral. You will need to install Optim and Miscellaneous package into your octave before using fzero function. This function calls another function which will find its zero. So in this case myfunfermi.m is

function y=myfunfermi(x)
global fvalfermi;

y=fermi_dirac_mhalf(x)-fvalfermi;

Note that I need to use a global function to pass into myfunfermi.m. This is because the function that fzero can call must be in the form

f(x)=0

And the function fermi_dirac_mhalf.m is from GSL package of Octave.

Advertisements

4 Responses

  1. […] Posted on April 25, 2008 by kurniawano Regarding my previous post on obtaining inverse fermi dirac integral. I found a paper by Antia which uses rational function to solve the inverse integral. Antia, H. M. […]

  2. How do you call this function? What value do you give to fval?

  3. How do you decide the initial guess, init_guess?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: