计算定积分会得到一个非常复杂的数值表达式而不是数值

Question

我正在尝试计算可见光谱中太阳辐射的百分比。
我的代码是：

h=6.626e-34;  %planck constant
c=3e8;        %speed of light
k=1.38066e-23; %boltzman constant
T1=5777;       %temperature in surface of sun  

syms f

J1 =(2*pi*h/(c^2))*((f^3)/(exp((h*f)/(k*T1))-1)); %planck radiation law

hold on;
ezplot(J1,[0,1.5e15])
a=int(J1,f,0,inf)  %total energy radiated
b=int(J1,f,4e14,8e14) %energy between 400-800Thz (visible radiation spectrum)
Vp = (b/a)*100  %percentage of visible/total radiation

虽然情节完全符合预期，也就是说我没有完全搞砸，但积分的结果是这样的：

Vp = -1888568826285205004703258735621345426367059580820393216707788800000000000*((73396718487075910602267519716779133887030184268951416015625*log(1 - exp(23642358029674224853515625/7114894705749824515342336)))/205953985278202888163058116890711980917418253877248 - (73396718487075910602267519716779133887030184268951416015625*log(1 - exp(23642358029674224853515625/3557447352874912257671168)))/25744248159775361020382264611338997614677281734656 - (2390487322005187985890576650155251369405251302897930145263671875*polylog(2, exp(23642358029674224853515625/3557447352874912257671168)))/1857466834100924357302864708366291649120175241251782656 + (25952248522181378144831874533777514511468878135769288445192954296875*polylog(3, exp(23642358029674224853515625/3557447352874912257671168)))/67008813354583054015095330308295397033299967000474002915328 - (110058495767576691259417256590823904271799518678985866399923893187011219092083*polylog(4, exp(23642358029674224853515625/3557447352874912257671168)))/1888568826285205004703258735621345426367059580820393216707788800000000 + (2390487322005187985890576650155251369405251302897930145263671875*polylog(2, exp(23642358029674224853515625/7114894705749824515342336)))/7429867336403697429211458833465166596480700965007130624 - (25952248522181378144831874533777514511468878135769288445192954296875*polylog(3, exp(23642358029674224853515625/7114894705749824515342336)))/134017626709166108030190660616590794066599934000948005830656 + (110058495767576691259417256590823904271799518678985866399923893187011219092083*polylog(4, exp(23642358029674224853515625/7114894705749824515342336)))/1888568826285205004703258735621345426367059580820393216707788800000000 + 101409666798338314597227594049400067888200283050537109375/22835963083295358096932575511191922182123945984))/(12228721751952965695490806287869322696866613186553985155547099243001246565787*pi^4)

它只是一个数值表达式（不包含任何常量），但我期待（并且我正在努力寻找）一个单一的值。
任何想法如何克服这个？提前致谢

double(Vp)

returns 44.3072 这正是我一直在寻找的

vpa(Vp)

returns 44.307238264260285485868531074049 + 1.5008384323384489567473242679822e-35*i

norm(vpa(Vp))

给出与 double(Vp) 相同的结果

然而，当我添加几行代码时：

d=int(J1,f,0,4e14);  %infrared energy
e=int(J1,f,8e14,inf); %ultraviolet energy
Ir = (d/a)*100;  %percentage of infrared radiation
Uv = (e/a)*100;  %percentage of ultraviolet radiation

double(Ir) 给出了这个错误：

Error using mupadmex
Error in MuPAD command: DOUBLE cannot convert the input expression into a  double array.

If the input expression contains a symbolic variable, use the VPA function  instead.

Error in sym/double (line 514)
    Xstr = mupadmex('symobj::double', S.s, 0);

Error in symplanckradlaw (line 21)
Infrared=double(Ir)

而 vpa(Ir) 甚至 norm(vpa(Ir)) 被赋予了复杂的复杂数值表达式，就像这样

(abs(7.3340024900713941224341547032249e-56*limit((652255981594442743874468745505068648842285814001516463259648*f^2*polylog(2, exp((3873563939581825*f)/466281739436020499437475332096)))/15004497594028668398955870330625 - (608270107310811468411217054651916403272252179121058584901915330886243977872955782952124416*f*polylog(3, exp((3873563939581825*f)/466281739436020499437475332096)))/58120880811771703455030054556291506586990890625 - f^4/4 + (466281739436020499437475332096*f^3*log(1 - exp((3873563939581825*f)/466281739436020499437475332096)))/3873563939581825 + (283625243683820020974472587101002022201492492463545426687503953676410023626686775978867575742611811818507908158310055936*polylog(4, exp((3873563939581825*f)/466281739436020499437475332096)))/225134948049212098682315198853176286979563186266469146812890625, f == 0, Right) - 146.24549829953781858190522266202 - 2.501397387230748261245540446637e-36*i)^2)^(1/2)

Answer 1

您可以将其转换为双倍使用，嗯：double:

Vp_double = double(Vp)

如果要选择精度也可以用vpa:

Vp_vpa = vpa(Vp)

你得到那个非常长的表达式的原因是因为 MATLAB 设法找到了 "closed form definite integral"。 IE。 MATLAB 设法找到了一个可以准确计算积分且没有舍入误差的表达式。这并不总是可行的，你会得到一个错误：

Warning: Explicit integral could not be found.

如果发生这种情况，那么您应该尝试 this answer.

中给出的方法

Answer 2

出于某种原因，Matlab 能够计算该积分的极限

a=int(J1,f,0,inf)

它在计算此积分中的极限为 0 时遇到问题

d=int(J1,f,0,4e14);

并且在这个积分中无限大

e=int(J1,f,8e14,inf);

我解决这个问题的方法是用一个非常低的值 (1e-45) 代替 0，用一个非常大的值 (1e22) 代替无限。出于我的目的，我得到了很好的结果，但我仍然觉得 matlab 可以在一种情况下评估限制但不能在另一种情况下评估相同的限制真的很奇怪。

计算定积分会得到一个非常复杂的数值表达式而不是数值

Calculating a definite integral results in a really complex numerical expression instead of a numerical value

matlab

symbolic-math