行列の方程式が正しく更新されない

Question

簡単な例として、私は次の式を解くことを試みています。はt = 0から始まり、Trは行列（対角要素の和）のトレースです。以下は、私は基本的なコードを提供します。私のコードはエラーなしでコンパイルされますが、私が望むようにソリューションを更新していません。私の予想される結果は、私が手で計算したものでもあります（行列の乗算によって手作業でチェックするのはとても簡単です）。

下のサンプルコードでは、ソリューションを格納する変数が2つあります。gはf(t=0)です。次にf(t+1)をfとして保存します。

complex,dimension(3,3) :: f,g 
integer :: k,l,m,p,q 

（T = 0）が、私はこの結果をチェックした

do l=1,3 !matrix index loops 
    do k=1,3 !matrix index loops 

     if (k == l) then 

      g(k,l) = cmplx(0.2,0) 

     else if (k /= l) then 

      g(k,l) = cmplx(0,0) 

     end if 


    end do 
end do 

以下のように定義されたグラム= Fを想定し、私はそれになりたいものを実際にあるので、私はt = 0でfは知っています正しく定義されています。

今私はt = 0でこの行列を使用しようとしており、式f(t+1) = f(t)+f(t)*Tr(f^2)によって支配されている行列を常に見つけようとしていますが、これは私が望むコードを正しく実装していない場所です。この結果を印刷

do m=1,3 !loop for 3 time iterations 

     do p=1,3 !loops for dummy indices for matrix trace 
     do q=1,3 

       g(1,1) = g(1,1) + g(1,1)*g(p,q)*g(p,q) !compute trace here 
       f(1,1) = g(1,1) 

       !f(2,2) = g(2,2) + g(2,2)*g(p,q)*g(p,q) 
       !f(3,3) = g(3,3) + g(3,3)*g(p,q)*g(p,q) 

       !assume all other matrix elements are zero except diagonal 


    end do 
    end do 

end do 

は、私は自分のコードが正しく実装されていない実現したときにこれがある

print*, "calculated f where m=", m 
do k=1,3 
    print*, (f(k,l), l=1,3) 
end do 

によって行われます。
を印刷すると、t = 1の結果が0.224*identity matrixとなることが予想されます。しかし、t = 2の場合、出力は正しくない。だから、私のコードは最初の反復で正しく更新されていますが、その後は更新されません。

私は期待している結果を得るために、方程式を正しく実装する方法についてのソリューションを探しています。

Answer 1

問題が発生しているように思われるいくつかのことにお答えします。まず、トレース。 3×3行列のトレースはA(1,1)+A(2,2)+A(3,3)です。 1番目と2番目のインデックスは同じですので、とのループ変数を使用します。 NxN行列Aのトレース計算するには：

trace = 0. 
do i=1,N 
    trace = trace + A(i,i) 
enddo 

を私はあなたが間違っている、あなたのトレースを計算するためにpとqにわたってループ処理をしようとしていると思います。その合計で、A(2,3)のような条件で追加しますが間違っています。

第二には、更新を計算するために、私はあなたがfNewに更新fを計算お勧めして、あなたのコードのようなものになります。

do m=1,3 ! time 
    ! -- Compute f^2 (with loops not shown) 
    f2 = ... 

    ! -- Compute trace of f2 (with loop not shown) 
    trace = ... 

    ! -- Compute new f 
    do j=1,3 
     do i=1,3 
     fNew(i,j) = f(i,j) + trace*f(i,j) 
     enddo 
    enddo 

    ! -- Now update f, perhaps recording fNew-f for some residual 
    ! -- The LHS and RHS are both arrays of dimension (3,3), 
    ! -- so fortran will automatically perform an array operation 
    f = fNew 
enddo 

をこの方法では、2つの利点があります。最初に、あなたのコードは実際にあなたがしようとしている数学のように見え、従うのは簡単です。これは非常に簡単ではない現実的なproblesmにとって重要な非常にです。次に、fNew(i,j)がf(i+1,j)に依存している場合、現在の時間レベル値を使用する必要がある間は、次の時間レベルに更新していません。