通過辛普森方法數值積分

Question

我想通過辛普森方法來解決這種積分，並將結果繪製在一個圖中。該圖僅取自for循環的P0 = -6的值。當我設置我（K，P）時會給出錯誤：

Attempted to access P0(0); index must be a positive integer or logical

我該如何解決？

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

目前只爲你將它們保存在I(k)使用P0 = -6結果。首先你將結果保存爲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保存在另一個變量中，因爲您的代碼中的結果不會相互依賴。