矩阵方程没有及时正确更新

Question

作为一个简单的例子来说明我的观点，我试图从 t=0 开始求解以下方程 f(t+1) = f(t) + f(t)*Tr (f^2)，其中 Tr 是矩阵的轨迹（对角线元素的总和）。下面我提供一个基本代码。我的代码编译没有错误，但没有按照我的意愿更新解决方案。我的预期结果也低于我手工计算的结果（通过矩阵乘法手工检查非常容易）。

在我下面的示例代码中，我有两个存储解决方案的变量，g 用于我实现的 f(t=0)，然后我将 f(t+1) 存储为 f。

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

假设g=f(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

我已经检查过这个结果确实是我想要的，所以我知道 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

这是我意识到我的代码没有正确实现的时候。
当我打印 f(k,l) 时，我希望 t=1 得到 0.224*identity matrix 的结果，现在我明白了。但是对于 t=2 输出是不正确的。所以我的代码在第一次迭代时被正确更新，但之后就没有了。

我正在寻找有关如何正确执行方程式的解决方案，我想获得我期望的结果。

Answer 1

我将回答您似乎遇到的一些问题。第一，痕迹。 3x3 矩阵的轨迹是 A(1,1)+A(2,2)+A(3,3)。第一和第二个索引相同，所以我们使用 one 循环变量。计算 NxN 矩阵的迹 A:

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

我认为你正试图循环 p 和 q 来计算你的轨迹，这是不正确的。在那个总和中，您将添加 A(2,3) 之类的术语，这是错误的。

其次，要计算更新，我建议您将更新后的 f 计算为 fNew，然后您的代码将类似于：

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

这种方法有两个优点。首先，您的代码实际上看起来像您尝试做的数学运算，而且很容易理解。这对于不是那么简单的现实问题非常重要。其次，如果 fNew(i,j) 依赖于 f(i+1,j)，例如，您不会更新到下一个时间级别，而仍需要使用当前时间级别值。

矩阵方程没有及时正确更新

Matrix equation not properly updating in time

fortran

matrix

matrix-multiplication