How do I plot a integration result?

14 views (last 30 days)
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')

Accepted Answer

Star Strider
Star Strider on 1 Feb 2022
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.
.
  2 Comments
Md. Golam Zakaria
Md. Golam Zakaria on 1 Feb 2022
Would you please redo this using for loop. It will be benificial for me as I may need this soon.
Star Strider
Star Strider on 1 Feb 2022
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.
.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!