1

我正在尝试使用 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

我试图寻找函数句柄数组和一些不同的东西，但到目前为止似乎没有任何问题可以解决。

4

1 回答 1

1

您的代码中的错误是您i在循环中同时用作虚数单位和迭代变量for。有两种方法可以解决此问题：

  1. 将迭代变量更改i为ii。
  2. 更改i为fun，1i以便您明确使用虚数单位 i。

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

  1. 您可以移出fun循环for，k因为fun不依赖于k.
  2. mtimes并且mpower很少使用。使用*and^代替。
于 2016-07-26T23:57:45.273 回答