代码不会在 MATLAB 中产生定积分的值
Code wont produce the value of a definite integral in MATLAB
我的代码在尝试进行积分计算时遇到了问题，但它不会提供强大的功能，P2。
我试过使用匿名函数句柄在 MATLAB 上使用 integral() 函数以及仅使用 int()，但它仍然无法计算。这些值是否太小以至于 MATLAB 无法集成，或者我只是遗漏了一些小东西？
如果能将我推向正确的方向，我们将不胜感激。谢谢！
代码中的问题在标记为 "Power Calculations" 的部分的底部。如果这有所不同，我的积分也会变得非常混乱。
%%%%%%%%%%% Parameters %%%%%%%%%%%%
n0 = 1; %air
n1 = 1.4; %layer 1
n2 = 2.62; %layer 2
n3 = 3.5; %silicon
L0 = 650*10^(-9); %centre wavelength
L1 = 200*10^(-9): 10*10^(-9): 2200*10^(-9); %lambda from 200nm to 2200nm
x = ((pi./2).*(L0./L1)); %layer phase thickness
r01 = ((n0 - n1)./(n0 + n1)); %reflection coefficient 01
r12 = ((n1 - n2)./(n1 + n2)); %reflection coefficient 12
r23 = ((n2 - n3)./(n2 + n3)); %reflection coefficient 23
t01 = ((2.*n0)./(n0 + n1)); %transmission coefficient 01
t12 = ((2.*n1)./(n1 + n2)); %transmission coefficient 12
t23 = ((2.*n2)./(n2 + n3)); %transmission coefficient 23
Q1 = [1 r01; r01 1]; %Matrix Q1
Q2 = [1 r12; r12 1]; %Matrix Q2
Q3 = [1 r23; r23 1]; %Matrix Q3
%%%%%%%%%%%% Graph of L vs R %%%%%%%%%%%
R = zeros(size(x));
for i = 1:length(x)
P = [exp(j.*x(i)) 0; 0 exp(-j.*x(i))]; %General Matrix P
T = ((1./(t01.*t12.*t23)).*(Q1*P*Q2*P*Q3)); %Transmission
T11 = T(1,1); %T11 value
T21 = T(2,1); %T21 value
R(i) = ((abs(T21./T11))^2).*100; %Percent reflectivity
end
plot(L1,R)
title('Percent Reflectance vs. wavelength for 2 Layers')
xlabel('Wavelength (m)')
ylabel('Reflectance (%)')
%%%%%%%%%%% Power Calculation %%%%%%%%%%
syms L; %General lamda
y = ((pi./2).*(L0./L)); %Layer phase thickness with variable Lamda
P1 = [exp(j.*y) 0; 0 exp(-j.*y)]; %Matrix P with variable Lambda
T1 = ((1./(t01.*t12.*t23)).*(Q1*P1*Q2*P1*Q3)); %Transmittivity matrix T1
I = ((6.16^(15))./((L.^(5)).*exp(2484./L) - 1)); %Blackbody Irradiance
Tf11 = T1(1,1); %New T11 section of matrix with variable Lambda
Tf2 = (((abs(1./Tf11))^2).*(n3./n0)); %final transmittivity
P1 = Tf2.*I; %Power before integration
L_initial = 200*10^(-9); %Initial wavelength
L_final = 2200*10^(-9); %Final wavelength
P2 = int(P1, L, L_initial, L_final) %Power production
我重构了你的代码
- 为了更容易阅读
- 提高代码重用
- 提高性能
- 为了更容易理解
为什么要使用这么多不必要的括号？！
无论如何，我在您的代码中发现了一些问题。
你用i作为循环变量，j作为虚数单位。这一次还可以，但勉强如此。以后虚数单位最好用1i或者1j，and/orm或者ii什么的other 而不是 i 或 j 作为循环索引变量。您正在帮助自己和您的同事；这样就不会那么混乱了。
到最后，您连续两次使用了变量名 P1，并且使用了两种不同的方式。虽然它在这里有效，但它令人困惑！我花了一段时间才弄明白为什么一个矩阵生成函数却生成标量...
但到目前为止，您的代码中最大的问题是黑体辐照度计算的数值问题。术语
L⁵ · exp(2484/L) - 1
对于 λ₀ = 200 · 10⁻⁹ m 将需要计算数量
exp(1.242 · 10¹⁰)
不用说，这对计算机来说是相当困难的:)实际上，你的计算有两个问题。首先，取幂肯定超出了 64 位 IEEE-754 双精度的范围，因此会导致 ∞。第二，括号错误；普朗克定律应为
C/L⁵ · 1/(exp(D) - 1)
用C和D常数（涉及普朗克常数、光速和玻尔兹曼常数），你大概已经预先计算好了（我没有检查这些值。我do 知道单位的选择会搞砸这些，所以最好检查一下）。
因此，除了愚蠢的括号错误外，我怀疑主要问题是您只是忘记了将 λ 重新调整为 nm。将黑体方程中的所有内容更改为 nm 并更正这些括号给出代码
I = 6.16^(15) / ( (L*1e+9)^5 * (exp(2484/(L*1e+9)) - 1) );
有了这个，我得到了
的积分的有限值
P2 = 1.052916498836486e-010
但是，再说一次，你最好仔细检查所有内容。
请注意，我使用了 quadgk()，因为它是 R2010a 上可用的更好的之一（我坚持使用它），但您可以轻松地将其替换为 integral()任何比 R2012a 更新的东西。
这是我最终得到的代码：
function pwr = my_fcn()
% Parameters
n0 = 1; % air
n1 = 1.4; % layer 1
n2 = 2.62; % layer 2
n3 = 3.5; % silicon
L0 = 650e-9; % centre wavelength
% Reflection coefficients
r01 = (n0 - n1)/(n0 + n1);
r12 = (n1 - n2)/(n1 + n2);
r23 = (n2 - n3)/(n2 + n3);
% Transmission coefficients
t01 = (2*n0) / (n0 + n1);
t12 = (2*n1) / (n1 + n2);
t23 = (2*n2) / (n2 + n3);
% Quality factors
Q1 = [1 r01; r01 1];
Q2 = [1 r12; r12 1];
Q3 = [1 r23; r23 1];
% Initial & Final wavelengths
L_initial = 200e-9;
L_final = 2200e-9;
% plot reflectivity for selected lambda range
plot_reflectivity(L_initial, L_final, 1000);
% Compute power production
pwr = quadgk(@power_production, L_initial, L_final);
% Helper functions
% ========================================
% Graph of lambda vs reflectivity
function plot_reflectivity(L_initial, L_final, N)
L = linspace(L_initial, L_final, N);
R = zeros(size(L));
for ii = 1:numel(L)
% Transmission
T = transmittivity(L(ii));
% Percent reflectivity
R(ii) = 100 * abs(T(2,1)/T(1,1))^2 ;
end
plot(L, R)
title('Percent Reflectance vs. wavelength for 2 Layers')
xlabel('Wavelength (m)')
ylabel('Reflectance (%)')
end
% Compute transmittivity matrix for a single wavelength
function T = transmittivity(L)
% Layer phase thickness with variable Lamda
y = pi/2 * L0/L;
% Matrix P with variable Lambda
P1 = [exp(+1j*y) 0
0 exp(-1j*y)];
% Transmittivity matrix T1
T = 1/(t01*t12*t23) * Q1*P1*Q2*P1*Q3;
end
% Power for a specific wavelength. Note that this function
% accepts vector-valued wavelengths; needed for quadgk()
function pwr = power_production(L)
pwr = zeros(size(L));
for ii = 1:numel(L)
% Transmittivity matrix
T1 = transmittivity(L(ii));
% Blackbody Irradiance
I = 6.16^(15) / ( (L(ii)*1e+9)^5 * (exp(2484/(L(ii)*1e+9)) - 1) );
% final transmittivity
Tf2 = abs(1/T1(1))^2 * n3/n0;
% Power before integration
pwr(ii) = Tf2 * I;
end
end
end
我的代码在尝试进行积分计算时遇到了问题，但它不会提供强大的功能，P2。
我试过使用匿名函数句柄在 MATLAB 上使用 integral() 函数以及仅使用 int()，但它仍然无法计算。这些值是否太小以至于 MATLAB 无法集成，或者我只是遗漏了一些小东西？
如果能将我推向正确的方向，我们将不胜感激。谢谢！
代码中的问题在标记为 "Power Calculations" 的部分的底部。如果这有所不同，我的积分也会变得非常混乱。
%%%%%%%%%%% Parameters %%%%%%%%%%%%
n0 = 1; %air
n1 = 1.4; %layer 1
n2 = 2.62; %layer 2
n3 = 3.5; %silicon
L0 = 650*10^(-9); %centre wavelength
L1 = 200*10^(-9): 10*10^(-9): 2200*10^(-9); %lambda from 200nm to 2200nm
x = ((pi./2).*(L0./L1)); %layer phase thickness
r01 = ((n0 - n1)./(n0 + n1)); %reflection coefficient 01
r12 = ((n1 - n2)./(n1 + n2)); %reflection coefficient 12
r23 = ((n2 - n3)./(n2 + n3)); %reflection coefficient 23
t01 = ((2.*n0)./(n0 + n1)); %transmission coefficient 01
t12 = ((2.*n1)./(n1 + n2)); %transmission coefficient 12
t23 = ((2.*n2)./(n2 + n3)); %transmission coefficient 23
Q1 = [1 r01; r01 1]; %Matrix Q1
Q2 = [1 r12; r12 1]; %Matrix Q2
Q3 = [1 r23; r23 1]; %Matrix Q3
%%%%%%%%%%%% Graph of L vs R %%%%%%%%%%%
R = zeros(size(x));
for i = 1:length(x)
P = [exp(j.*x(i)) 0; 0 exp(-j.*x(i))]; %General Matrix P
T = ((1./(t01.*t12.*t23)).*(Q1*P*Q2*P*Q3)); %Transmission
T11 = T(1,1); %T11 value
T21 = T(2,1); %T21 value
R(i) = ((abs(T21./T11))^2).*100; %Percent reflectivity
end
plot(L1,R)
title('Percent Reflectance vs. wavelength for 2 Layers')
xlabel('Wavelength (m)')
ylabel('Reflectance (%)')
%%%%%%%%%%% Power Calculation %%%%%%%%%%
syms L; %General lamda
y = ((pi./2).*(L0./L)); %Layer phase thickness with variable Lamda
P1 = [exp(j.*y) 0; 0 exp(-j.*y)]; %Matrix P with variable Lambda
T1 = ((1./(t01.*t12.*t23)).*(Q1*P1*Q2*P1*Q3)); %Transmittivity matrix T1
I = ((6.16^(15))./((L.^(5)).*exp(2484./L) - 1)); %Blackbody Irradiance
Tf11 = T1(1,1); %New T11 section of matrix with variable Lambda
Tf2 = (((abs(1./Tf11))^2).*(n3./n0)); %final transmittivity
P1 = Tf2.*I; %Power before integration
L_initial = 200*10^(-9); %Initial wavelength
L_final = 2200*10^(-9); %Final wavelength
P2 = int(P1, L, L_initial, L_final) %Power production
我重构了你的代码
- 为了更容易阅读
- 提高代码重用
- 提高性能
- 为了更容易理解
为什么要使用这么多不必要的括号？！
无论如何，我在您的代码中发现了一些问题。
你用
i作为循环变量，j作为虚数单位。这一次还可以，但勉强如此。以后虚数单位最好用1i或者1j，and/orm或者ii什么的other 而不是i或j作为循环索引变量。您正在帮助自己和您的同事；这样就不会那么混乱了。到最后，您连续两次使用了变量名
P1，并且使用了两种不同的方式。虽然它在这里有效，但它令人困惑！我花了一段时间才弄明白为什么一个矩阵生成函数却生成标量...但到目前为止，您的代码中最大的问题是黑体辐照度计算的数值问题。术语
L⁵ · exp(2484/L) - 1对于 λ₀ = 200 · 10⁻⁹ m 将需要计算数量
exp(1.242 · 10¹⁰)不用说，这对计算机来说是相当困难的:)实际上，你的计算有两个问题。首先，取幂肯定超出了 64 位 IEEE-754 双精度的范围，因此会导致 ∞。第二，括号错误；普朗克定律应为
C/L⁵ · 1/(exp(D) - 1)用
C和D常数（涉及普朗克常数、光速和玻尔兹曼常数），你大概已经预先计算好了（我没有检查这些值。我do 知道单位的选择会搞砸这些，所以最好检查一下）。
因此，除了愚蠢的括号错误外，我怀疑主要问题是您只是忘记了将 λ 重新调整为 nm。将黑体方程中的所有内容更改为 nm 并更正这些括号给出代码
I = 6.16^(15) / ( (L*1e+9)^5 * (exp(2484/(L*1e+9)) - 1) );
有了这个，我得到了
的积分的有限值P2 = 1.052916498836486e-010
但是，再说一次，你最好仔细检查所有内容。
请注意，我使用了 quadgk()，因为它是 R2010a 上可用的更好的之一（我坚持使用它），但您可以轻松地将其替换为 integral()任何比 R2012a 更新的东西。
这是我最终得到的代码：
function pwr = my_fcn()
% Parameters
n0 = 1; % air
n1 = 1.4; % layer 1
n2 = 2.62; % layer 2
n3 = 3.5; % silicon
L0 = 650e-9; % centre wavelength
% Reflection coefficients
r01 = (n0 - n1)/(n0 + n1);
r12 = (n1 - n2)/(n1 + n2);
r23 = (n2 - n3)/(n2 + n3);
% Transmission coefficients
t01 = (2*n0) / (n0 + n1);
t12 = (2*n1) / (n1 + n2);
t23 = (2*n2) / (n2 + n3);
% Quality factors
Q1 = [1 r01; r01 1];
Q2 = [1 r12; r12 1];
Q3 = [1 r23; r23 1];
% Initial & Final wavelengths
L_initial = 200e-9;
L_final = 2200e-9;
% plot reflectivity for selected lambda range
plot_reflectivity(L_initial, L_final, 1000);
% Compute power production
pwr = quadgk(@power_production, L_initial, L_final);
% Helper functions
% ========================================
% Graph of lambda vs reflectivity
function plot_reflectivity(L_initial, L_final, N)
L = linspace(L_initial, L_final, N);
R = zeros(size(L));
for ii = 1:numel(L)
% Transmission
T = transmittivity(L(ii));
% Percent reflectivity
R(ii) = 100 * abs(T(2,1)/T(1,1))^2 ;
end
plot(L, R)
title('Percent Reflectance vs. wavelength for 2 Layers')
xlabel('Wavelength (m)')
ylabel('Reflectance (%)')
end
% Compute transmittivity matrix for a single wavelength
function T = transmittivity(L)
% Layer phase thickness with variable Lamda
y = pi/2 * L0/L;
% Matrix P with variable Lambda
P1 = [exp(+1j*y) 0
0 exp(-1j*y)];
% Transmittivity matrix T1
T = 1/(t01*t12*t23) * Q1*P1*Q2*P1*Q3;
end
% Power for a specific wavelength. Note that this function
% accepts vector-valued wavelengths; needed for quadgk()
function pwr = power_production(L)
pwr = zeros(size(L));
for ii = 1:numel(L)
% Transmittivity matrix
T1 = transmittivity(L(ii));
% Blackbody Irradiance
I = 6.16^(15) / ( (L(ii)*1e+9)^5 * (exp(2484/(L(ii)*1e+9)) - 1) );
% final transmittivity
Tf2 = abs(1/T1(1))^2 * n3/n0;
% Power before integration
pwr(ii) = Tf2 * I;
end
end
end