Creating arrayfun with two variables

Hello.

I'd like to create cell array with results of a function S for variables q and k, which is defined in a code below.

function z=comparison_sum_integral_trapz
tic
format long

tt=-0.000689609;t=0.242731; muu=0.365908;

[m,NN]=meshgrid(0:100,-500:1:500);

    y1= @(N,q,k) t*q./k.*log((-k.^2+2*k.*q-q.^2+muu+1i*(2*pi*N.*t-(2*m(1,:)+1)*pi*t))./(-k.^2-2*k.*q-...
        q.^2+muu+1i*(2*pi*N.*t-(2*m(1,:)+1)*pi*t)))./(tt*pi+integral(@(a)a.*tanh((a.^2-muu)./(2*t)).*log((2*a.^2+2*a.*q+...
        q.^2-2*muu-1i*2*pi*N*t)./(2*a.^2-2*a.*q+q.^2-2*muu-1i*2*pi*N*t))./q-2,0,10000,'AbsTol',1e-6,'RelTol',1e-3,'ArrayValued',true));
    
      R1=@(q,k) integral(@(N)y1(N,q,k),500,10^6,'AbsTol',1e-6,'RelTol',1e-3,'ArrayValued',true);
      R11=@(q,k) integral(@(N)y1(N,q,k),-10^6,-500,'AbsTol',1e-6,'RelTol',1e-3,'ArrayValued',true);
    
    y2=@(q,k) t*q./k.*log((-k.^2+2*k.*q-q.^2+muu+1i*(2*pi*NN(:,1).*t-(2*m(1,:)+1)*pi*t))./(-k.^2-2*k.*q-...
        q.^2+muu+1i*(2*pi*NN(:,1).*t-(2*m(1,:)+1)*pi*t)))./(tt*pi+integral(@(a)a.*tanh((a.^2-muu)./(2*t)).*log((2*a.^2+2*a.*q+...
        q.^2-2*muu-1i*2*pi*NN(:,1).*t)./(2*a.^2-2*a.*q+q.^2-2*muu-1i*2*pi*NN(:,1).*t))./q-2,0,10000,'AbsTol',1e-6,'RelTol',1e-3,'ArrayValued',true));
    
    R2=@(q,k) sum(y2(q,k));
    
    S=@(q,k) R2(q,k)+R11(q,k)+R1(q,k)-4*sqrt(2)/pi*(1/1000)/(pi^(3/2)*sqrt(t))*q.^2;

q=0.001:1:7;
k=0.001:1:7;

out=arrayfun(S,k,q,'UniformOutput',false)
end

My problem is that I expected to get a cell array 7×7 but instead of I obtained 1×7.

out =

  1×7 cell array

  Columns 1 through 5

    [1×101 double]    [1×101 double]    [1×101 double]    [1×101 double]    [1×101 double]

  Columns 6 through 7

    [1×101 double]    [1×101 double]

Elapsed time is 3.165576 seconds.

I'd like to avoid an expoitation of

meshgrid

due to long time calculations (on the next step I will calculate the cell array out for larger arrays of q and k ).

What I did it wrong?

I will apreciate for any suggestions.

[Rik]

Your S function returns an imaginary array of 101 elements. Run the code below to see what happens for a scalar input. Your code doesn't have any comments, so I don't understand what it is doing, so I can't determine what to fix.

temp=S(q(1),k(1));
x=real(temp);y=imag(temp);
plot(x,y)

[Yuriy Yerin]

I want to get a cell array for each values of k and q. This means that I will have a 7×7 cell array not 1×7 one like in my case. Btw I know that values of the function S are complex as it should be.

Answer 1

[Guillaume]

You seem to be expecting implicit expansion out of arrayfun (the way bsxfun works). arrayfun doesn't do that (and bsxfun can't create cell arrays, so you can't use that either. I'm afraid that you don't have any other choice than using meshgrid or ndgrid:

[k, q] = ndgrid(0.001:7);
out = arrayfun(S, k, q, 'UniformOutput', false)

[Yuriy Yerin]

A lot of thanks for your reponse! Maybe it's possible to use logical operators like

k=0.001:7;
for j=1:7
    
end

to calculate values of S and get cell array out for all values of q?

Before I tried to implement that but I didn't succeded.

[Guillaume]

I have no idea what you're implying. I don't see how logical operators would help in any way.

You could replace arrayfun by a double for loop. That would possibly be marginally faster, although usually the more costly part of arrayfun is the anonymous function call, which you would incur with the loop since S is already an anonymous function.

k = 0.001:7; q = 0.001:7;
out = cell(numel(k), numel(q));
for row = 1:numl(k)
    for col = 1:nume(q)
        out{row, col} = S(k(row), q(col));
    end
end

edit: Note that I've not tried to understand what you're doing with S. The best solution would probably be to modify S so that it can work on array inputs instead of scalars.

[Yuriy Yerin]

Thank you for the comment!

My final goal is to apply trapezoid integration over q for the given number m and given k. Then to make summation over m and finally to plot a result as a function of k.

Earlier I realized above mentioned procedures via a function 'integral' with a parameter arrayvalued true. But for some input parameters (tt, t and mu) during the integration over q I had essential discontinuity of a integrand and the 'integral' works very bed, namely it takes a lot of time. I decided to do the same using trapz function and avoid that problem.