How can I evaluate a complex definite integral?

Hello;

I would appreciate if anyone can help me with this problem. I am working in a project, and I was asked to evaluate the next definite integral

The constants are:

k= 1.380649e-23

h=6.6226069e-34

T=6000

fg2= 484e12

I tried first to evaluate the integral and then multiply the result with the rest of the equation, but I cannot do it. I am a bit new in Matlab, I tried using Live Script. The result should be 0.2952

Thanks for the help in advance

syms f 

k= 1.380649e-23

k = 1.3806e-23

h=6.6226069e-34 

h = 6.6226e-34

T=6000 

T = 6000

fg2= 484e12 

fg2 = 4.8400e+14

n= (f^2/(exp(h*f/k*T)-1))

n = 

%Integral 

n2=int(n,0,inf)

n2 = 

%Multiply 

n3= ((8.17e-43)*((fg2)/(T^4)))*(n2)

n3 = 

var = vpa(n3)

var = 

[John D'Errico]

When you post a picture of your code, you make it far more difficult to get help, since then someone needs to type in everything you did from scratch. Is there a good reason why you want to make it more difficult, since you are asking for help? It is just as easy, if not easier to paste in text as a picture. But the txt is trivially easy for us to then bring into MATLAB.

If you want help, then make it easy to get help.

[Valentin Augusto Simbron Mckay]

Sorry I forgot , thanks for the advice

[John D'Errico]

That makes it far more easy to offer help. Thanks.

Answer 1

[Paul]

Hi Valentin,

Shouldn't the lower bound on the integral be fg2, not 0? Don't see how the reult will be 0.2952. The integrand is basically zero over the interval 1e8 < f < inf

syms f real
k = sym(1.380649e-23);
h = sym(6.6226069e-34);
T = sym(6000); 
fg2 = sym(484e12); 
n(f) = (f^2/(exp(h*f/k*T)-1));
fplot(n(f),[0 1e8])

So the integral from fg2 to inf will be tiny

n2 = vpaintegral(n(f),f,fg2,inf)

n2 = 

7.55952e-60495963

And the scale factor is much less than 1

vpa(sym(8.17e-43)*fg2/T^4)

ans = 

3.0511419753086418658056274018784e-43

[Valentin Augusto Simbron Mckay]

Yes, you have right, thanks for the help

[Paul]

After reading @David Goodmanson's Answer below, I realize I had a mistake in my code, which I copied from Question. But the code in the Question doesn't represent the actual equation in the Question.

syms f real
k = sym(1.380649e-23);
h = sym(6.6226069e-34);
T = sym(6000); 
fg2 = sym(484e12); 
%n(f) = (f^2/(exp(h*f/k*T)-1));  % should have divided by T !!!!!
n(f) = (f^2/(exp(h*f/k/T)-1));
n2 = 8.17e-43*fg2/T^4*vpaintegral(n(f),f,fg2,inf)

n2 = 

0.31052671103833630524482060960562

Matches @David Goodmanson's result, naturally.

Answer 2

[Walter Roberson]

You will need to switch to numeric operations. The definite integral involves polylog(2) and polylog(3) and the limit of those as f approaches 0 and infinity. When I test out the formula on some larger f values, I seem to encounter complex intermediate results, which I do not get numerically.

format long g
Q = @(v) sym(v);

syms f 

syms T fg2 positive
k = Q(1380649)*sym(10)^(-29)

k = 

h = Q(66226069)*sym(10)^(-41)

h = 

T = sym(6000)

T = 

6000

fg2 = sym(484)*sym(10)^(12)

fg2 = 

484000000000000

h_k = h/k

h_k = 

%assume(0 < h_k & h_k < 1)
n = (f^2/(exp(f*h_k*T)-1))

n = 

%fplot(nfun, [0 2/(h_k*T)])
%Integral 
nfun = matlabFunction(n)

nfun = function_handle with value:
    @(f)f.^2./(exp(f.*2.878040790961352e-7)-1.0)

n2fun = @(F) arrayfun(@(U) integral(nfun, 0, U), F);


F = linspace(0, 1e8, 250);

N = n2fun(F);

Warning: Inf or NaN value encountered.

plot(F, N)

N(end-2:end)

ans = 1×3

1.0e+00 * 1.00847279626234e+20 1.0084727962773e+20 1.00847279629073e+20

Answer 3

[David Goodmanson]

Hi Valentin,

It's not strictly necessary (see Paul's vpa solution) but it is advantageous to scale things properly before doing the integration. Once that is done, the answer is more apparent conceptually. The result below comes out .3105, slightly different from .2952.

result = 8.17e-43 fg2/T^4 Integral{fg2,inf}  f^2/(exp(hf/(kT)) -1) df 

To scale this, let f = fg2*u, df = fg2*du, then

result = 8.17e-43 fg2^4/T^4 Integral{1,inf}  u^2/(exp( [h fg2/(kT)] u ) -1) du

which goes as follows

k = 1.380649e-23;
h = 6.6226069e-34; 
T = 6000; 
fg2 = 484e12; 
C1 = h*fg2/(k*T)
C2 = 8.17e-43*fg2^4/T^4
I = integral(@(u) fun(u,C1), 1, inf)
result = C2*I

function y = fun(u,C1);
y = u.^2./(exp(C1*u)-1);
end

% scaling works
%   C1 = 3.8694
%   C2 = 34.5938
%   I = 0.0088

result = 0.3105