使用多个参数迭代 quadgk

Question

我正在尝试使用 quadgk 评估数值积分，因为我不是 matlab 专家，我很难让以下代码工作。 我有矩阵 g(i,j)，其中我正在评估 g 的每个元素的参数 phi 的积分。这部分代码工作正常，但是当我想更改矩阵 g 的大小时问题就开始了，在这种情况下，只有第一个值是正确的，并且对于更大的尺寸（k），它对于 g 的所有元素都返回零。

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. 将 fun 中的 i 更改为 1i，这样您就可以明确使用虚数单位 i。

顺便说一句，您应该在代码中更改和改进一些额外的东西：

1. 您可以将fun 移出for 循环，覆盖k，因为fun 不依赖于k。
2. mtimes 和 mpower 很少使用。请改用* 和^。