在matlab中繪製一個涉及積分的函數

Question

我想繪製一個包含定積分的函數。我的代碼使用所有匿名函數。當我運行該文件時，它給了我一個錯誤。我的代碼如下：

%%% List of Parameters %%% 
    gamma_sp = 1; 
    cap_gamma = 15; 
    gamma_ph = 0; 
    omega_0 = -750; 
    d_omega_0 = 400; 
    omega_inh = 100; 
    d_omega_inh = 1000; 

    %%% Formulae %%% 
    gamma_t = gamma_sp/2 + cap_gamma/2 + gamma_ph; 
    G = @(x) exp(-(x-omega_inh).^2./(2*d_omega_inh.^2))./(sqrt(2*pi)*d_omega_inh); 
    F = @(x) exp(-(x-omega_0).^2./(2*d_omega_0.^2))./(sqrt(2*pi)*d_omega_0); 
    A_integral = @(x,y) G(x)./(y - x + 1i*gamma_t); 
    Q_integral = @(x,y) F(x)./(y - x + 1i*gamma_t); 
    A = @(y) integral(@(x)A_integral(x,y),-1000,1000); 
    Q = @(y) integral(@(x)Q_integral(x,y),-3000,0); 

    P1 = @(y) -1./(1i.*(gamma_sp + cap_gamma)).*(1./(y + 2.*1i.*gamma_t)*(A(y)-conj(A(0)))-1./y.*(A(y)-A(0))+cap_gamma./gamma_sp.*Q(y).*(A(0)-conj(A(0)))); 

    P2 = @(y) conj(P1(y)); 
    P = @(y) P1(y) - P2(y); 
    sig = @(y) abs(P(y)).^2; 

    rng = -2000:0.05:1000; 
    plot(rng,sig(rng)) 

在我看來，該程序運行繪圖命令時，應該把RNG到SIG（Y）的每一個值，而該值將被用作A_integral y值和Q_積分。但是，當我嘗試運行該程序時，matlab會引發錯誤。

Error using - 
Matrix dimensions must agree. 

Error in @(x,y)G(x)./(y-x+1i*gamma_t) 

Error in @(x)A_integral(x,y) 

Error in integralCalc/iterateScalarValued (line 314) 
      fx = FUN(t); 

Error in integralCalc/vadapt (line 133) 
     [q,errbnd] = iterateScalarValued(u,tinterval,pathlen); 

Error in integralCalc (line 76) 
    [q,errbnd] = vadapt(@AtoBInvTransform,interval); 

Error in integral (line 89) 
Q = integralCalc(fun,a,b,opstruct); 

Error in @(y)integral(@(x)A_integral(x,y),-1000,1000) 

Error in 
@(y)-1./(1i.*(gamma_sp+cap_gamma)).*(1./(y+2.*1i.*gamma_t)*(A(y)-conj(A(0)))-1. /y.*(A(y)-A(0))+cap_gamma./gamma_sp.*Q(y).*(A(0)-conj(A(0)))) 

Error in @(y)P1(y)-P2(y) 

Error in @(y)abs(P(y)).^2 

Error in fwm_spec_diff_paper_eqn (line 26) 
plot(rng,sig(rng)) 

有關我在做什麼錯的任何想法？

Answer 1

你有

>> rng = -2000:0.05:1000; 
>> numel(rng) 
ans = 
    60001 

所有60001個元素，會傳下來的

A = @(y) integral(@(x)A_integral(x,y),-1000,1000); 

這就要求

A_integral = @(x,y) G(x)./(y - x + 1i*gamma_t); 

（類似於爲Q）。問題是，integral是一種自適應正交方法，意思是（大致）它將插入到A_integral中的x的數量隨A_integral在某些x處的行爲方式而變化。

因此，y中元素的數量通常會與A_integral調用中的x中的元素不同。這就是爲什麼y-x +1i*gamma_t失敗。

考慮到你想要做的事情的複雜性，我認爲最好將所有匿名函數重新定義爲適當的函數，並將其中的一些函數集成到單個函數中。查看bsxfun的文檔以查看是否有幫助（例如，bsxfun(@minus, y.', x)而不是y-x可能可以解決這些問題中的一些問題），否則，僅在x中進行矢量化並且在y上循環。

Answer 2

謝謝Rody，這對我有意義。我一直試圖使用像mathematica一樣的matlab，我忘記了matlab如何做。我稍微修改了一下代碼，併產生了正確的結果。積分評估非常粗略，但應該很容易解決這個問題。我已在下面發佈我的修改後的代碼。

%%% List of Parameters %%% 
gamma_sp = 1; 
cap_gamma = 15; 
gamma_ph = 0; 
omega_0 = -750; 
d_omega_0 = 400; 
omega_inh = 100; 
d_omega_inh = 1000; 

%%% Formulae %%% 
gamma_t = gamma_sp/2 + cap_gamma/2 + gamma_ph; 
G = @(x) exp(-(x-omega_inh).^2./(2*d_omega_inh.^2))./(sqrt(2*pi)*d_omega_inh); 
F = @(x) exp(-(x-omega_0).^2./(2*d_omega_0.^2))./(sqrt(2*pi)*d_omega_0); 
A_integral = @(x,y) G(x)./(y - x + 1i*gamma_t); 
Q_integral = @(x,y) F(x)./(y - x + 1i*gamma_t); 

w = -2000:0.05:1000; 
sigplot = zeros(size(w)); 
P1plot = zeros(size(w)); 
P2plot = zeros(size(w)); 
Pplot = zeros(size(w)); 
aInt_range = -1000:0.1:1200; 
qInt_range = -2000:0.1:100; 
A_0 = sum(A_integral(aInt_range,0).*0.1); 

for k=1:size(w,2) 

P1plot(k) = -1./(1i*(gamma_sp + cap_gamma)).*(1./(w(k)+2.*1i.*gamma_t).*(sum(A_integral(aInt_range,w(k)).*0.1)-conj(A_0))-1./w(k).*(sum(A_integral(aInt_range,w(k)).*0.1)-A_0)+cap_gamma./gamma_sp.*sum(Q_integral(qInt_range,w(k)).*0.1).*(A_0-conj(A_0))); 
P2plot(k) = conj(P1plot(k)); 
Pplot(k) = P1plot(k) - P2plot(k); 
sigplot(k) = abs(Pplot(k)).^2; 

end 

plot(w,sigplot)