Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB equation with integral

  1. Apr 18, 2013 #1
    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. jcsd
  3. Apr 18, 2013 #2

    kreil

    User Avatar
    Gold Member

    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
    [tex]I =8T \int_0^{\infty} \frac{\sqrt{\epsilon}}{e^{\epsilon - \mu}} d\epsilon = \frac{8T\sqrt{\pi}}{2}e^{\mu}[/tex]
    So plotting I is kind of boring since it is just an exponential.
    Solving for mu tells us that
    [tex]\mu(T) \approx \ln \left( \frac{I(\mu)}{T} \right)[/tex]

    2) Use Leibniz's Integral Rule: http://en.wikipedia.org/wiki/Leibniz_integral_rule
     
    Last edited: Apr 18, 2013
  4. Apr 19, 2013 #3
    that 8T is actually under the exp- like that:

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

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

    and what would the analitical solution be?
     
  5. Apr 19, 2013 #4
    [tex]\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}[/tex]

    Now substitute
    [tex] t := \frac{ε}{8T} [/tex]

    and deploy the integral form of the Γ- function, like kreil did above.
     
  6. Apr 20, 2013 #5
    all right. that seems easy. but what if the integral looked like that? (and I wanted to solve it by matlab?):

    [tex]I=\int_{0}^{inf} dε \frac{\sqrt(ε)}{\exp(\frac{ε-μ}{8T})+1}[/tex]
     
  7. Apr 20, 2013 #6
    What have you yet tried to solve that problem on your own?
     
    Last edited: Apr 20, 2013
  8. Apr 22, 2013 #7
    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
  9. Apr 22, 2013 #8
    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);


    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 [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
  10. Apr 22, 2013 #9
    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
     
  11. Apr 23, 2013 #10
    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 [tex]\zeta(\nu)[/tex]?
     
  12. Apr 23, 2013 #11
  13. Apr 23, 2013 #12
    @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?
     
  14. Apr 23, 2013 #13

    kreil

    User Avatar
    Gold Member

    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
     
  15. Apr 23, 2013 #14

    kreil

    User Avatar
    Gold Member

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: MATLAB equation with integral
  1. Integration in Matlab (Replies: 3)

  2. Integration in Matlab (Replies: 1)

  3. Matlab integral (Replies: 0)

  4. Integration in matlab (Replies: 0)

  5. Integration in matlab (Replies: 3)

Loading...