Итерирование по 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 

здесь выводится при к = 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, чтобы вы использовали мнимую единицу i явно.

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

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