Plotting a function in matlab involving an integral

Question

I'm trying to plot a function that contains a definite integral. My code uses all anonymous functions. When I run the file, it gives me an error. My code is below:

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

It seems to me that when the program runs the plot command, it should put each value of rng into sig(y), and that value will be used as the y value in A_integral and Q_integral. However, matlab throws an error when I try to run the program.

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

Any ideas about what I'm doing wrong?

Answer 1

You have

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

all 60001 elements get passed down to

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

which calls

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

(similar for Q). The thing is, integral is an adaptive quadrature method, meaning (roughly) that the amount of x's it will insert into A_integral varies with how A_integral behaves at certain x.

Therefore, the amount of elements in y will generally be different from the elements in x in the call to A_integral. This is why y-x +1i*gamma_t fails.

Considering the complexity of what you're trying to do, I think it is best to re-define all anonymous functions as proper functions, and integrate a few of them into single functions. Look into the documentation of bsxfun to see if that can help (e.g., bsxfun(@minus, y.', x) instead of y-x could perhaps fix a few of these issues), otherwise, vectorize only in x and loop over y.