辛普森方法的数值积分

Question

我正在尝试通过 simpsons 方法解决此积分问题，并将结果绘制在图中。该图仅从 for 循环中获取 P0= -6 的值。当我设置 I(K,P) 时，它给出了错误：

试图访问 P0(0)；索引必须是正整数或逻辑

我该如何解决？

alpha = 45;     
beta = 185;     
gamma_e = 116;  

% Gain values

G_ei = -18.96;    
G_ee = 18.52;
G_sr = -0.26;
G_rs = 16.92;
G_es = 2.55;
G_re = 4.67;
G_se = 0.73;
G_sn = 2.78;

G_esre = G_es*G_sr*G_re;
G_srs = G_sr*G_rs;
G_ese = G_es*G_se;
G_esn = G_es*G_sn;

t_0 = 0.085;    % corticothalamic loop delay in second
r_e = 0.086;    % Excitatory axon range in metre
f = linspace(-40,40,500);   % f = frequency in Hz
w = 2*pi*f;                 % angular frequency in radian per second
delt_P = 0.5;

L=zeros(1,500);
Q=repmat(L,1);
P=repmat(L,1);

%%%%%%%%%%%%% integration %%%%%%%%%%%%

a = -80*pi;
b = 80*pi;
n=500;

I=repmat(L,1);
P_initial = repmat(L,1);
P_shift = repmat(L,1);
p = repmat(L,1);

for k = 1:length(w)
    for P0 = [6 -6]

        L_initial = @(w1) (1-((1i*w1)/alpha))^(-1)*(1-((1i*w1)/beta))^(-1);                                                                 

        Q_initial = @(w1) (1/(r_e^2))*((1-((1i*w1)/gamma_e))^(2) - (1/(1-G_ei*L_initial(w1)))*....
            (L_initial(w1)*G_ee + (exp(1i*w1*t_0)*(L_initial(w1)^2*G_ese +L_initial(w1)^3*G_esre))/(1-L_initial(w1)^2*G_srs)));                  

        P_initial =  @(w1) (pi/r_e^4)* (abs((L_initial(w1)^2*G_esn)/((1-L_initial(w1)^2*G_srs)*....
            (1-G_ei*L_initial(w1)))))^2 * abs((atan2((imag(Q_initial(w1))),(real(Q_initial(w1)))))/imag(Q_initial(w1)));                   

        G =  150*exp(- (f - P0).^2./(2*(delt_P).^2));   

        P2 = @(w1) G(k) + P_initial(w1);

        L_shift =  @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1);                               

        Q_shift  = @(w1)  (1/(r_e^2))*((1-((1i*(w(k)-w1))/gamma_e))^(2) - (1/(1-G_ei*L_shift(w1)))*...
            (L_shift(w1)*G_ee + (exp(1i*(w(k)-w1)*t_0)*(L_shift(w1)^2*G_ese +L_shift(w1)^3*G_esre))/(1-L_shift(w1)^2*G_srs)));     

        P_shift =  @(w1)  (pi/r_e^4)* (abs((L_shift(w1)^2*G_esn)/((1-L_shift(w1)^2*G_srs)*(1-G_ei*L_shift(w1)))))^2 *....
            abs((atan2((imag(Q_shift(w1))),(real(Q_shift(w1)))))/imag(Q_shift(w1)));                                

        p =  @(w1)  P2(w1)*P_shift(w1);        %  Power spectrum formula for P(w1)*p(w-w1)

        I(k) = simprl(p,a,b,n);

    end
end

figure(1)
plot(f,I,'r--')

figure(2)
plot(f,G,'k')

Answer 1

目前您只使用P0 = -6 的结果，因为您将它们保存在I(k) 中。首先将结果保存为P0 = 6，然后覆盖它并保存另一个。 P0 = 6 的结果既不使用也不保存。如果您按以下方式编写代码，这将得到澄清。

for k = 1:length(w)

    L_shift =  @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1);                               

    Q_shift  = @(w1)  (1/(r_e^2))*((1-((1i*(w(k)-w1))/gamma_e))^(2) - (1/(1-G_ei*L_shift(w1)))*...
        (L_shift(w1)*G_ee + (exp(1i*(w(k)-w1)*t_0)*(L_shift(w1)^2*G_ese +L_shift(w1)^3*G_esre))/(1-L_shift(w1)^2*G_srs)));     

    P_shift =  @(w1)  (pi/r_e^4)* (abs((L_shift(w1)^2*G_esn)/((1-L_shift(w1)^2*G_srs)*(1-G_ei*L_shift(w1)))))^2 *....
        abs((atan2((imag(Q_shift(w1))),(real(Q_shift(w1)))))/imag(Q_shift(w1)));           


    for P0 = [6 -6]
        G =  150*exp(- (f - P0).^2./(2*(delt_P).^2));   
        P2 = @(w1) G(k) + P_initial(w1);
        p =  @(w1)  P2(w1)*P_shift(w1);        
        I(k) = simprl(p,a,b,n);
    end
end

您无法访问I(k,P)，因为I 是向量而不是矩阵。但是，这会让您索引超出矩阵维度。您可以将P0 = -6 的结果保存在一个变量中，将P0 = 6 的结果保存在另一个变量中，因为您的代码中的结果不会相互依赖。