2

私は Fortran の初心者で、4x4 行列を解くためにガウス消去法コードを書く必要があります。以下のコードは間違った結果を返し、問題をデバッグできませんでした。助けていただければ幸いです。

      common /grid/ A(100,100), NEQ, C(100), X(100)

      open(10, file="NEQ.txt", status='unknown')
      read(10,*) NEQ
      close (10)

      open(12, file="C1.txt", status='unknown')
        do i=1,NEQ
        read(12,*) C(i)
        enddo
      close (12)

      open(11, file="A1.txt", status='unknown')
        do i=1,NEQ
        read(11,*) (A(i,k), k=1,NEQ)
        enddo
      close (11)

      call SOL

      open(13, file="X.txt", status='unknown')
        do i=1,NEQ
        write(13,*) X(i)
        enddo
      close (13)

      stop
      end

      subroutine SOL
      common /grid/ A(100,100), NEQ, C(100), X(100)


c     Forward Reduction Phase:

      do 10 K=2,NEQ
      do 10 I=K,NEQ
      R=A(I,K-1)/A(K-1,K-1)

      C(I)=C(I)-R*C(K-1)

      do 10 J=K-1,NEQ
10    A(I,J)=A(I,J)-R*A(K-1,J)

c     Back Substitution Phase:

      X(NEQ)=C(NEQ)/A(NEQ,NEQ)
      do 30 K=NEQ-1,1,-1
      X(K)=C(K)
      do 20 J=K+1,NEQ

20    X(K)=X(K)-A(K,J)*X(J)
30    X(K)=X(K)/A(K,K)

      return
      end

私の場合、NEQ はテキスト ファイルから 4 として読み取られます。私の A1.txt は次のとおりです。

18, -6, -6, 0
-6, 12, 0, -6
-6, 0, 12, -6
0, -6, -6, 18

C1.txt は次のとおりです。

60
0
20
0

結果の X 行列は次のようになります。

8.3333330    
6.6666665    
8.3333321    
4.9999995

それ以外の：

13.13
14.17
15.83
15.00
4

1 に答える 1

4

Vladimir F がコメントしたように、何か新しいコードを書かなければならない場合、少なくとも Fortran 90 を使用することは非常に役に立ちます。Fortran 77 よりもはるかに優れたオプションを提供して、コードを構造化および整理 (および読み取り可能) に保ちます。implicit noneエラーを減らすための非常に貴重なステートメントでもあることを付け加えたいと思います。

そうは言っても、Fortran 90 でアルゴリズムがどのように見えるかの例を次に示します (関数として記述しました)。

function gaussian_elimination(A, C) result(X)
  implicit none
  real, intent(inout) :: C(:), A(size(C), size(C))
  real :: X(size(C))

  real    :: R(size(C))
  integer :: i, j, neq

  neq = size(C)

  ! Forward reduction, only two loops since reduction is now row by row
  do i = 1, neq
    R = A(:,i)/A(i,i)

    do j = i+1, neq
      A(j,:) = A(j,:) - R(j)*A(i,:)
      C(j) = C(j) - R(j)*C(i)
    enddo
  enddo

  ! Back substitution, only one loop
  do i = neq, 1, -1
    x(i) = (C(i) - sum(A(i, i+1:) * x(i+1:))) / A(i,i)
  enddo
end function gaussian_elimination

これは、なぜあなたの結果X = [8.33, 6.67, 8.33, 5.00]が間違っていると思うのか疑問に思うだけです? これは正しい解であり、行列 A を掛けることで確認できます:matmul(A, X)は (ほぼ) に等しいはずCです。

于 2012-12-08T17:52:25.600 に答える