I am having trouble doing integration for the following code.

I am having the following warning.

Warning: Infinite or Not-a-Number value encountered.

> In integralCalc/iterateScalarValued (line 349)

In integralCalc/vadapt (line 132)

In integralCalc (line 83)

In integral (line 88)

In Q2 (line 17)

clc
clear all

Fs= 2.16*10^-5*pi;         % Geometrical Factor
h= 6.626*10^-34;           % Plancks Constant
c= 3*10^8;                 % Speed of light
K = 1.38*10^-23;           % Boltzmanns Constant
Ts=5760;                   % Temparature of the sun

E=0:0.5:4;

bs=((2*Fs)/((h^3)*(c^2))).*(E.^2./(exp(E./((K*Ts)-1))));

fun= @(E)(E.*((2*Fs)/((h^3)*(c^2))).*(E.^2./(exp(E./((K*Ts)-1)))));
Ps = integral (fun,0,Inf) ;

figure(2)
plot(E,Ps)

Answer 1

[Star Strider]

The warning is likely the result of the second term, that becomes Inf at higher values.

Fs= 2.16*10^-5*pi;         % Geometrical Factor
h= 6.626*10^-34;           % Plancks Constant
c= 3*10^8;                 % Speed of light
K = 1.38*10^-23;           % Boltzmanns Constant
Ts=5760;                   % Temparature of the sun

fun= @(E)(E.*((2*Fs)/((h^3)*(c^2))).*(E.^2./(exp(E./((K*Ts)-1)))));
Result1 = fun(0)

Result1 = 0

Result2 = fun(realmax)

Result2 = Inf

Check = @(E) E.^2/exp(E./((K*Ts)-1))

Check = function_handle with value:
    @(E)E.^2/exp(E./((K*Ts)-1))

Check(realmax)

ans = Inf

I am not certain that there is any solution for that, other than restricting the upper integration limit to something smaller.

.

Answer 2

[Walter Roberson]

Q = @(v) sym(v);
Fs = Q(2.16)*Q(10^-5)*Q(pi)         % Geometrical Factor

Fs = 

h = Q(6.626)*Q(10^-34)           % Plancks Constant

h = 

c = Q(3)*Q(10^8)                 % Speed of light

c = 

300000000

K = Q(1.38)*Q(10^-23)           % Boltzmanns Constant

K = 

Ts = Q(5760)                   % Temparature of the sun

Ts = 

5760

E = Q(0):Q(0.5):Q(4);
bs = ((2*Fs)/((h^3)*(c^2))).*(E.^2./(exp(E./((K*Ts)-1))));

fun = @(E)(E.*((2*Fs)/((h^3)*(c^2))).*(E.^2./(exp(E./((K*Ts)-1)))));
syms x b real
assume(b >= 0)
Ps = int(fun(x), x, 0, b)

Ps = 

limit(Ps, b, inf)

ans = 

∞