How can I evaluate a complex definite integral?

조회 수: 1(최근 30일)
Valentin Augusto Simbron Mckay 2022년 11월 18일
편집: David Goodmanson 2022년 11월 29일
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 =
댓글 수: 3표시숨기기 이전 댓글 수: 2
John D'Errico 2022년 11월 19일
That makes it far more easy to offer help. Thanks.

댓글을 달려면 로그인하십시오.

채택된 답변

Paul 2022년 11월 19일
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
댓글 수: 2표시숨기기 이전 댓글 수: 1
Paul 2022년 11월 28일
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.

댓글을 달려면 로그인하십시오.

추가 답변(2개)

Walter Roberson 2022년 11월 19일
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
댓글 수: 0표시숨기기 이전 댓글 수: -1

댓글을 달려면 로그인하십시오.

David Goodmanson 2022년 11월 28일
편집: David Goodmanson 2022년 11월 29일
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

댓글을 달려면 로그인하십시오.

범주

Find more on Common Operations in Help Center and File Exchange

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by