How do I plot a integration result?

Hello everyone

I need to plot the result of the integral as I tried in the code below (figure 2). Obviously I am doing something wrong or missing something; I am not that experienced. The Result of the integral returns a single value. How can I plot the result of the integral elementwise. It would be a great help.

Thank You.

clc
clear all

Fs= 2.16*10^-6*pi;            % Geometrical Factor
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun

E=0:0.0001: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(1)
plot(E,bs),xlabel('Photon Energy'),ylabel('Spectral Flux Density')
figure(2)
plot(E,Ps),xlabel('Photon Energy'),ylabel('Power  Density')

Answer 1

[Star Strider]

The ‘Ps’ integration produceds onlly one value, the value of the integral between the limits of integration.

To get a vector for it, either use a for loop or arrayfun (which is a for loop in disguise) as I did here.

% clc
% clear all

Fs= 2.16*10^-6*pi;            % Geometrical Factor
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun

% E=0:0.0001:4;
E = linspace(0, 6, 250);

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=arrayfun(@(ul)integral(fun,0,ul), E);

Warning: Inf or NaN value encountered.

figure(1)
plot(E,bs),xlabel('Photon Energy'),ylabel('Spectral Flux Density')

figure(2)
plot(E,Ps),xlabel('Photon Energy'),ylabel('Power  Density')

I changed the end value of ‘E’ just to see what the curve did, and discovered that it became asymptotic at about 6.. Change it back if necessary. Also, 250 elements of it are enough to produce a smooth curve.

.

[Md. Golam Zakaria]

Would you please redo this using for loop. It will be benificial for me as I may need this soon.

[Star Strider]

Sure!

% clc
% clear all

Fs= 2.16*10^-6*pi;            % Geometrical Factor
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun

% E=0:0.0001:4;
E = linspace(0, 6, 250);

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)));

for k = 1:numel(E)
    Ps(k) = integral(fun,0,E(k));
end

Warning: Inf or NaN value encountered.

figure(1)
plot(E,bs),xlabel('Photon Energy'),ylabel('Spectral Flux Density')

figure(2)
plot(E,Ps),xlabel('Photon Energy'),ylabel('Power  Density')

The same result.

I am not certain whether arrayfun or the for loop is faster, since I did not time them here.

.