代码不会在MATLAB中产生定积分的值

Question

由于我试图进行积分计算，所以我的代码出现了问题，但它不适合功率P2。

我试过使用匿名函数句柄在MATLAB上使用integral()函数以及仅使用int()，但它仍然不会计算。对于MATLAB来说，这些值是否太小？或者我只是错过了一些小的东西？

任何帮助或建议，将不胜感激，促使我在正确的方向。谢谢！

代码中的问题位于标有“功率计算”部分的底部。如果这有所作为，我的积分也会变得相当混乱。

%%%%%%%%%%% 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 

Answer 1

我已经重构你的代码

1. ，使其更易于阅读
2. ，提高代码的重用
3. 以提高性能
4. ，使其更容易理解

为什么你使用这么多不必要的括号？

无论如何，我在代码中看到了一些问题。

1. 您使用i作为循环变量，和j作为虚数单位。这个例子可以，但几乎没有。将来最好使用1i或1j作为虚数单元，和/或m或ii或其他其他比i或j作为循环索引变量。你在帮助你和你的同事;这种方式不那么容易混淆。

2. 最后，您连续两次使用变量名P1，并以两种不同方式使用。虽然它在这里工作，但它是混淆！花了我一段时间，以解开为什么一个矩阵生产函数产生的标量，而不是...

3. 但到目前为止，您的代码中最大的问题是与黑体辐照度计算的数值问题。对于λ₀  = · 10⁻⁹ 米术语

   L⁵ · exp(2484/L) - 1 
   

   将需要计算量

   exp(1.242 · 10¹⁰) 
   

   其中，不用说，是相当困难的用于计算机:)其实，你的计算问题是双重的。首先，求幂绝对超出64位的IEEE-754双精度范围，因此将导致∞。其次，括号是错误的;普朗克定律应该读

   C/L⁵ · 1/(exp(D) - 1) 
   

   与C和D常数（涉及普朗克常数，光速和玻尔兹曼常数），你已经大概预计算的（我没有检查的价值观。我做知道单位的选择可以搞砸了，所以更好的检查）。

所以，除了愚蠢的括号错误，我怀疑主要的问题是，你只是忘了把λ重新调整到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