2016-07-26 2 views
1

私はmatlabの専門家ではないので、quadgkを使用して数値積分を評価しようとしていますが、次のコードを実行するのに苦労しています。 私は行列g(i、j)を持っています。ここでは、gの各要素に対して、パラメータphiの積分を評価しています。この部分のコードは正常に動作していますが、行列gのサイズを変更したいときに問題が発生します。この場合、最初の値のみが正しいので、より大きいサイズ(k)のgのすべての要素についてゼロが返されます。ここ複数のパラメータでquadgkを反復する

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 
     [email protected](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 

ための出力である私は、関数ハンドルの配列、および異なる物事のカップルを探してみましたが、何もそれほど問題を解決するようです遠い

答えて

1

コード内のバグは、forループ内で、iを虚数単位と繰り返し変数として使用していることです。これを修正する方法は2つあります。

  1. iiiに変更します。
  2. ifun1iに変更して、虚数単位iを明示的に使用しています。あなたがkに依存しないfunので、k以上forループの外funを移動することができます

  3. ところで、あなたのコードに変更し、改善すべきいくつかの追加のものがあります。
  4. mtimesおよびmpowerはめったに使用されない。代わりに*^を使用してください。
+0

こんにちはedwinksl、ありがとう、答えは、今は完全に動作している!しかし、私のコードがインデックスkの繰り返しを使ってもkの個々の値に対して正しい結果を出していたので、私は混乱しています。だから私は{i}を使って問題が発生したとは思わなかったのです。 – setareh

関連する問題