Search code examples
matlabfunction-handle

iterating over quadgk with multiple parameters

I am trying to evaluate a numerical integration using quadgk, as I am not expert in matlab, I am having a hard time to get the following code to work. I have matrix g(i,j) where I am evaluating an integral over parameter phi for each element of g. This part of code is working properly but the problem starts when I want to change the size of matrix g,in this case, only the first value is correct, and it is returning zero for all elements of g for higher sizes(k). 

clear;
alpha=2.0;
h=1.0;
lmax=12;
for k=2:2:4
  fun = @(phi,t,s) (exp(-i.*(t-s).*phi).*(exp(-i.*phi)-1)./sqrt((1-h.*exp(i.*phi)).*(1-h.*exp(-i.* phi))))./(2.*pi);
     for i=1:k
       for j=1:k  
        F=@(phi) fun(phi,i,j);
        g(i,j)=real(quadgk(F,0,2.*pi));
       end
     end 
  Y1=mtimes(transpose(g),g)
  Y2=mpower(Y1,1./2.);
  Z1 = 0.5.*(eye(k) - Y2);
  Z2 = 0.5.*(eye(k) + Y2);
  C1 = mpower(Z1, alpha) + mpower(Z2, alpha);
  M1=diag(log(eig(C1)));
  s(k/2)=k;
  ent(k/2)=real(trace(M1))./(1-alpha);
end

here is output for k=2 and 4,

  k =

     2


g =

  -0.636619772367581   0.636619772367581
  -0.212206590789194  -0.636619772367581


Y1 =

   0.450316371743723  -0.270189823046234
  -0.270189823046234   0.810569469138702


k =

     4


g =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0


Y1 =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

I tried to look for the array of function handles, and a couple of different things, but nothing seems to solve the problem so far.

Solution

  • The bug in your code is that you are using i as both an imaginary unit and an iteration variable in a for loop. There are two ways to fix this:

    1. Change the iteration variable i to ii.
    2. Change i in fun to 1i so you are using the imaginary unit i explicitly.

    Incidentally, there are some additional things that you should change and improve in your code:

    1. You can move fun out of the for loop over k since fun does not depend on k.
    2. mtimes and mpower are rarely used. Use * and ^ instead.