Fortran'da Temel Gauss Eliminasyonu ile Beklenmeyen Davranış

Question

Birçok yaygın olarak kullanılabilen rutinlerden birini, özellikle de bu işlevi gauss (bilinen bir sonuçla birlikte örnek verileriyle birlikte) kullanarak matris formundaki bir lineer denklem sistemini çözmeye çalışıyorum ama Kullandığım eleme stratejisinden bağımsız olarak sürekli olarak yanlış sonuçlar elde ediyorum ve bu noktada nereye döneceğimi bilmiyorum. Benim kodum:

MODULE gaussMod 
    CONTAINS 
    function gauss(a,b) result(x) 
     implicit none 
     real*8 :: b(:), a(size(b), size(b)) 
     real*8 :: x(size(b)) 

     real*8 :: r(size(b)) 
     integer i,j, neq 

     neq = size(b) 

     do i =1, neq 
      r = a(:,i)/a(i,i) 

      do j = i+1, neq 
       a(j,:) = a(j,:) - r(j)*a(i,:) 
       b(j) = b(j) - r(j)*b(j) 
      enddo 
     enddo 

     do i= neq, 1, -1 
      x(i) = (b(i) - sum(a(i, i+1:) * x(i+1:)))/a(i,i) 
     enddo 
    END function 
END MODULE 

SUBROUTINE outFile(n,x) 
    IMPLICIT none 
    INTEGER n 
    INTEGER i 
    REAL*8, DIMENSION(n) :: x 


    OPEN(UNIT=20,FILE="solution.csv",action="write",status="replace") 

    DO i=1, n 
     WRITE (20,"(1(f0.30,',',:))") x(i) 
    END DO 
    CLOSE(20) 
END SUBROUTINE 

PROGRAM elimtest 
     USE gaussMod 
     IMPLICIT none 
     INTEGER, PARAMETER :: coeff_kind = selected_real_kind(p=30, r=99) 
     INTEGER n 
     REAL*8, DIMENSION(:,:), ALLOCATABLE :: a 
     REAL*8, DIMENSION(:), ALLOCATABLE :: b 
     REAL*8, DIMENSION(:), ALLOCATABLE :: x 
     n = 4 

     ALLOCATE(a(n,n)) 
     ALLOCATE(b(n)) 
     ALLOCATE(x(n)) 

      a(1,1) = 18. 
      a(1,2) = -6. 
      a(1,3) = -6. 
      a(1,4) = 0. 
      a(2,1) = -6. 
      a(2,2) = 12. 
      a(2,3) = 0. 
      a(2,4) = -6. 
      a(3,1) = -6. 
      a(3,2) = 0. 
      a(3,3) = 12. 
      a(3,4) = -6. 
      a(4,1) = 0. 
      a(4,2) = -6. 
      a(4,3) = -6. 
      a(4,4) = 18. 

      b(1) = 60. 
      b(2) = 0. 
      b(3) = 20. 
      b(4) = 0. 


     x = gauss(a,b) 
     CALL outFile(n,x) 


END PROGRAM 

Herhangi bir yardım çok takdir edilecektir!

Answer 1

Bu aşağıdaki satırı

b(j) = b(j) - r(j)*b(j) 

benziyor kodunuzu doğru sonucu verir Bu değişiklik ile

b(j) = b(j) - r(j)*b(i) 

bir yazım hatası olduğunu (muhtemelen!): * I

x(:) = 8.3333333333333339 6.6666666666666670 8.3333333333333339 5.0000000000000000