Вычисление определенного интеграла приводит к действительно сложному численному выражению вместо численного значения

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) 

возвращает 44.3072, который является именно то, что я искал в то время как

vpa(Vp) 

возвращает 44.307238264260285485868531074049 + 1.5008384323384489567473242679822e-35 * я, который преодолен путем

norm(vpa(Vp)) 

давая тот же результат, что и двойной (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 

двойной (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) и даже норма (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

По какому-то причине Matlab времени был в состоянии оценить пределы в этом интеграл

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

это были проблемы, оценивающих предел 0 в этом интегральном

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

и предел к бесконечному в этом интегральных

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

Я работал вокруг этого, поставив вместо 0 очень низкое значение (1е-45) и вместо бесконечного очень большого (1e22). Я получил очень хорошие результаты для своих целей, но мне все еще кажется странным, что Matlab может оценить пределы в одном случае, но не может оценить те же пределы в другом случае.

Answer 2

Вы можете преобразовать его в двойной с использованием, а также: double:

Vp_double = double(Vp) 

Вы также можете использовать vpa, если вы хотите выбрать точность:

Vp_vpa = vpa(Vp) 

---

Причина, почему вы получите, что очень долго выражение, потому что MATLAB удалось найти «закрытую форму определенный интеграл». То есть MATLAB удалось найти выражение, которое точно вычисляет интеграл без ошибок округления. Это не всегда возможно, и вы получите сообщение об ошибке:

Warning: Explicit integral could not be found.

Если это произойдет, то вы должны попробовать подход, заданный в this answer.