Итерация по quadgk с несколькими параметрами

Question

Я пытаюсь оценить числовую интеграцию, используя quadgk, так как я не эксперт в Matlab, мне трудно заставить следующий код работать. У меня есть матрица g(i,j), где я оцениваю интеграл по параметру phi для каждого элемента g. Эта часть кода работает правильно, но проблема начинается, когда я хочу изменить размер матрицы g, в этом случае верным является только первое значение, и оно возвращает ноль для всех элементов g для более высоких размеров (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

здесь вывод для k=2 и 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

Я попытался найти массив дескрипторов функций и пару разных вещей, но пока что ничего не решает проблему.

Answer 1

Ошибка в вашем коде в том, что вы используете i как мнимая единица и переменная итерации в for петля. Есть два способа исправить это:

1. Изменить переменную итерации i в ii,
2. + Изменить i в fun в 1i так что вы используете мнимую единицу я явно.

Кстати, есть некоторые дополнительные вещи, которые вы должны изменить и улучшить в своем коде:

1. Вы можете двигаться fun вне for зацикливаться k поскольку fun не зависит от k,
2. mtimes а также mpower редко используются. использование * а также ^ вместо.