How can I solve an equation with an integral in MATLAB?

[_Matt87_]

hi people,
I've just recently started using MATLAB (last week) and 've already got plenty of problems. Could anyone help me with this for example [??]:

1)

3=∫(0,inf) (sqrt(ε)/exp[(ε-μ)/(8*T)] dε

so this is the equation. I would like to solve in a way that μ=... , and then plot the μ(T) with T=0:1:10000;

2)

f=∫(-inf,inf) d[ 1/{exp[(ε-μ)/(8*T)]+1} ]/dT dε

where μ is a functions from the previous. and there is a partial derivative with respect to T.
also plot it as f(T).

 

[kreil]

1) Use the function arrayfun:

I = @(u) arrayfun(@(mu) integral(@(e) sqrt(e).*exp(-e+mu),0,inf),u);

This solves the integral for whatever values are specified in u. Analytically, we know that
I =8T \int_0^{\infty} \frac{\sqrt{\epsilon}}{e^{\epsilon - \mu}} d\epsilon = \frac{8T\sqrt{\pi}}{2}e^{\mu}
So plotting I is kind of boring since it is just an exponential.
Solving for mu tells us that
\mu(T) \approx \ln \left( \frac{I(\mu)}{T} \right)

2) Use Leibniz's Integral Rule: http://en.wikipedia.org/wiki/Leibniz_integral_rule

 

[_Matt87_]

that 8T is actually under the exp- like that:

I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})}

does it make a big difference? or I just add it to the code?

and what would the analitical solution be?

 

[Solkar]

that 8T is actually under the exp- like that:

I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})}

does it make a big difference? or I just add it to the code?

and what would the analitical solution be?

\begin{align}\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})}&= \int_{0}^{inf} dε \sqrt(ε) \exp(-\frac{ε-μ}{8T})\\ &= \int_{0}^{inf} ε^{1/2}\, \exp(\frac{-ε}{8T}) \exp(\frac{μ}{8T})\, dε\\&= \exp(\frac{μ}{8T}) \int_{0}^{inf} ε^{1/2}\, \exp(\frac{-ε}{8T})\, dε \end{align}

Now substitute
t := \frac{ε}{8T}

and deploy the integral form of the Γ- function, like kreil did above.

 

[_Matt87_]

all right. that seems easy. but what if the integral looked like that? (and I wanted to solve it by matlab?):

I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})+1}

 

[Solkar]

[...]but what if the integral looked like that? (and I wanted to solve it by matlab?):
I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})+1}

What have you yet tried to solve that problem on your own?

 

[_Matt87_]

ok ok
I started with something like that:

function u= chempot

k=1.380648813e-23;
u0=7.370012199e-19;
I=@(e,T,u) sqrt(e)/(exp((e-u)/(k*T))+1);
    II=integral(I,0,inf);

for T=1:1:10000
    x0=[-10*u0];
    u=fsolve(@(u) 2*u0^(3/2)/3-II,x0);

    plot(T,u)
    hold on
end

doesn't work. plenty errors, don't know what's going on.
is that 'integral' function actually integrate over e?/////
on the second thought it maybe should look more like this:

function u = proba(T)
k=1.380648813e-23;
u0=7.370012199e-19;

I=@(e,T,u) sqrt(e)/(exp((e-u)/(k*T))+1);
    II=integral(I,0,inf);
    
    u=fsolve(@(u) 2*u0^(3/2)/3-II,u0);

with temperature being set from outside of the function. It still doesn't work though :(

 

[_Matt87_]

hi all! I've got something!

function u = proba(T)

k=1.380648813e-23;
u0=7.370012199e-19;

options=optimset('TolFun',1.0e-16,'TolX',1.0e-16);


u=fsolve(@(u)(2*u0^(3/2)/3-quad(@(e)(sqrt(e)./(exp((e-u)./(k.*T))+1)),0,1000000,1e-20)),u0,options);
plot(T,u);
hold on
end

only that it doesn't work :) but it's very close.
I open the function in a loop for T=1:1:10000

but! the numbers seem to be too small for 1) quad 2)fsolve.
so that quad (or 'integral') gives 0. and whole fsolve give just the constant value of u0.

Could anyone help me now with it? as We're so close to the solution ?

It should be giving a curve like the one in the picture below for Three dimensions (I found it here http://www.lcst-cn.org/Solid%20State%20Physics/Ch63.html ) where ef is my u0 and x-axis should be temperature T.

http://www.lcst-cn.org/Solid%20State%20Physics/Ch63.files/image010.gif

 

[Solkar]

OK, that looks like work.

But before I consult Matlab, Maple, Maxima etc vX.Y I tend to consult Brain v1.0 and a good integral table

I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})+1}

Afais this

\int_0^\infty \frac{x^{\nu-1}}{e^{\mu x} + 1}\, dx = \frac{1}{\mu^\nu} \left(1 - 2^{1-\nu}\right) \Gamma(\nu) \zeta(\nu)\qquad \left[\Re(\mu) > 0, \Re(\nu)> 0 \right]

could be useful.

Regards, Solkar


[JZ07] Jeffrey, A. & Zwillinger, D. Table of Integrals, Series, and Products. Elsevier Science, 2007

 

[_Matt87_]

Thank you.
Only that I am in the same time trying to learn Matlab. And as I've got some other functions similar to mentioned, which I need to solve numerically, I would need someone just to tell me how to write a proper function for that. Or at least tell me why the one I put before do not work?

as for the Table of Integrals. Could you tell me what is \zeta(\nu)?

 

[Solkar]

Could you tell me what is \zeta(\nu)?

http://en.wikipedia.org/wiki/Riemann_zeta_function

 

[Solkar]

@kreil

The
exp(-µ/kT) (=: C)
part in the denominator is the main burden.

What do you thing about simply replacing this with i and then picking Im() of the result after plain complex-valued numerical integration?

What could hurt?

 

[kreil]

all right. that seems easy. but what if the integral looked like that? (and I wanted to solve it by matlab?):

I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})+1}

I should have guessed this is what the integral actually was. It's the Fermi-Dirac distribution.

I recommend using this m-file I found on the exchange and editing it to meet your needs.
http://www.mathworks.com/matlabcentral/fileexchange/13616-fermi

 

[kreil]

For more information on how to solve this,

https://docs.google.com/viewer?a=v&...I-rOHe&sig=AHIEtbSEX_yjg8KCNGBvHxS0z1i9UYiNYQ