How do I plot an integral which is repeated over an interval?

2 次查看（过去 30 天）
Mitchell
Mitchell 2024-9-30
评论： Voss 2024-10-1
I am attempting to plot the absporptance (alpha_s_bar) of an uncontaminated spacecraft as a function of the incident light wavelenght. The absorptance is found by integrating two different functions of wavelength over an interval from zero to the wavelength value of interest, on an interval from 0 to 0.7 microns. The code I have attached below runs and returns different values of absorptance as the wavelenght of interest (L) changes, but does not save the values before iterating, so I end up with one scalar value as alpha_s_bar. I am interested in saving the values of alpha_s_bar as a vector so I can plot them against the wavelength.
%Note that I will be using x for the wavelength for ease of typing
% Start by defining the constants for the absorptance for typical
% outgassing products, called alpha_c in problem statement, where
% alpha_c=M10*exp(M11*L)
M10=34.898;
M11=-8.6481;
% Next define the constants for solar absorptivity, alpha_s, where
% aplha_s=M0+M1*L+M2*L^2+M3*L^3
M0=0.28779;
M1=0.06028;
M2=-0.046055;
M3=0.020211;
% Next, define the constants for the solar spectrum (adjusted so units are microns), S, where
% S=2*C1/L^5[exp(C2/L*T)-1]
T=5500;
C1=5.594E-7;
C2=1438.8;
for L=0:0.01:0.7;
fnum=@(x) (M0+M1.*x+M2.*(x.^2)+M3.*(x.^3))./((x.^5).*exp(C2./(x.*T)));
fden=@(x) 1./((x.^5).*exp(C2./(x.*T)));
alpha_bar_s=integral(fnum,0,L)./integral(fden,0,L)
plot(L,alpha_bar_s)
hold on
end

请先登录，再回答此问题。

回答（2 个）

Abhas
Abhas 2024-10-1
Ran in:
Hi Mitchell,
To save the values of "alpha_bar_s" as a vector and plot them against the wavelength, you need to store each computed value in an array during the loop. After the loop, you can plot the results.
Here's the modified MATLAB code to achieve the same:
% Define constants
M10 = 34.898;
M11 = -8.6481;
M0 = 0.28779;
M1 = 0.06028;
M2 = -0.046055;
M3 = 0.020211;
T = 5500;
C1 = 5.594E-7;
C2 = 1438.8;
% Initialize vectors to store results
L_values = 0:0.01:0.7;
alpha_bar_s_values = zeros(size(L_values));
% Loop over wavelength values
for i = 1:length(L_values)
L = L_values(i);
if L == 0
% Handle the case where L is zero
alpha_bar_s_values(i) = 0;
else
% Define the integrand functions
fnum = @(x) (M0 + M1.*x + M2.*(x.^2) + M3.*(x.^3)) ./ ((x.^5) .* exp(C2 ./ (x .* T)));
fden = @(x) 1 ./ ((x.^5) .* exp(C2 ./ (x .* T)));
% Calculate alpha_bar_s for each L
alpha_bar_s_values(i) = integral(fnum, 0.001, L) / integral(fden, 0.001, L);
end
end
% Plot the results
plot(L_values, alpha_bar_s_values, 'b-', 'LineWidth', 2);
xlabel('Wavelength (microns)');
ylabel('Absorptance (\alpha_{s\_bar})');
title('Absorptance as a Function of Wavelength');
grid on;
You can address the situation when "L" is zero according to your needs. If unhandled, this might trigger a warning indicating that the integrand functions "fnum" or "fden" are producing "Inf" or "NaN" values at certain points, likely due to division by zero or overflow. This issue generally arises when the wavelength "L" is zero, as the functions involve division by "x^5", which becomes problematic when "x" is zero.
You may refer to the below MathWorks documentation links to have a better understanding on numerical integration and function handles:
  1. https://www.mathworks.com/help/matlab/ref/integral.html
  2. https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html
Voss
Voss 2024-10-1

Make L a vector, use one element of L on each iteration of the for loop, pre-allocate alpha_bar_s before the loop to be the same size as L, and store one element of alpha_bar_s on each iteration. Also, since the definitions of fnum and fden don't depend on the value of L, they can be set once before the loop.

M10=34.898;
M11=-8.6481;
% Next define the constants for solar absorptivity, alpha_s, where
% aplha_s=M0+M1*L+M2*L^2+M3*L^3
M0=0.28779;
M1=0.06028;
M2=-0.046055;
M3=0.020211;
% Next, define the constants for the solar spectrum (adjusted so units are microns), S, where
% S=2*C1/L^5[exp(C2/L*T)-1]
T=5500;
C1=5.594E-7;
C2=1438.8;
L=0:0.01:0.7; 
fnum=@(x) (M0+M1.*x+M2.*(x.^2)+M3.*(x.^3))./((x.^5).*exp(C2./(x.*T)));
fden=@(x) 1./((x.^5).*exp(C2./(x.*T)));
alpha_bar_s = zeros(size(L));
for ii = 1:numel(L) 
  alpha_bar_s(ii)=integral(fnum,0,L(ii))./integral(fden,0,L(ii));
end
plot(L,alpha_bar_s)

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by