Матричное уравнение, которое должным образом не обновляется во времени

Question

Как простой пример, иллюстрирующий мою точку зрения, я пытаюсь решить следующее уравнение f(t+1) = f(t) + f(t)*Tr (f^2), начинающееся с t = 0, где 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 

Я проверил этот результат действительно то, что я хочу, чтобы это было, так что я знаю, е при Т = 0 является правильно определена.

Теперь я пытаюсь использовать эту матрицу при 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). Первый и второй индексы одинаковы, поэтому мы используем одну переменную цикла. Для того, чтобы вычислить след в 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), например, вы не обновляетесь до следующего уровня времени, пока значения текущего уровня времени еще не используются.