# MATLAB equation with integral

1. Apr 18, 2013

### _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

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

2. Apr 18, 2013

### kreil

1) Use the function arrayfun:

Code (Text):
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

Last edited: Apr 18, 2013
3. Apr 19, 2013

### _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?

4. Apr 19, 2013

### Solkar

\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.

5. Apr 20, 2013

### _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}$$

6. Apr 20, 2013

### Solkar

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

Last edited: Apr 20, 2013
7. Apr 22, 2013

### _Matt87_

ok ok
I started with something like that:

Code (Text):
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:

Code (Text):
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 :(

Last edited: Apr 22, 2013
8. Apr 22, 2013

### _Matt87_

hi all! I've got something!

Code (Text):
function u = proba(T)

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

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

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 [Broken]) where ef is my u0 and x axis should be temperature T.

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

Last edited by a moderator: May 6, 2017
9. Apr 22, 2013

### 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
Afais this
could be useful.

Regards, Solkar

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

10. Apr 23, 2013

### _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)$$?

11. Apr 23, 2013

### Solkar

12. Apr 23, 2013

### 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?

13. Apr 23, 2013

### kreil

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

14. Apr 23, 2013