fortran - 使用高斯方法求解线性系统返回错误结果

Question

谢谢你们，你们所有的建议都是相关的，但我放弃了 Fortran 并将我的代码翻译成 Python。在那里，只有 1 或 2 个小错误，我的代码运行良好

我编写了一个程序来使用高斯方法求解线性系统。我写了所有的算法，前向消除和后向替换，我做了很多其他的子程序，我不知道出了什么问题，我不知道我的代码是否有问题，或者 Fortran 编程是否有问题，因为我是这门语言的新手。我将把我的代码和我应该找到解决方案的线性系统放在下面

PROGRAM metodo_Gauss
    IMPLICIT NONE

    REAL :: det_a_piv
    INTEGER :: n, i, j
    REAL, DIMENSION(:,:), ALLOCATABLE :: a, a_piv 
    INTEGER, DIMENSION(:), ALLOCATABLE :: p 
    REAL, DIMENSION(:), ALLOCATABLE :: b, x 

    PRINT*, "Entre com a dimensão n do sistema a ser resolvido"
    READ*, n

    ! allocate memory
    ALLOCATE(a(n, n))
    ALLOCATE(a_piv(n, n))
    ALLOCATE(p(n))
    ALLOCATE(b(n))
    ALLOCATE(x(n))

    CALL matriz_a(n, a)
    CALL vetor_b(n, b)

    a_piv(1:n, 1:n) = a(1:n, 1:n)

    DO i = 1, n
        x(i) = 0
    END DO

    CALL eliminacao(n, a, a_piv, p)

    det_a_piv = (-1) ** n
    DO j = 1, n
        det_a_piv = det_a_piv * a_piv(j, j)
    END DO

    IF (det_a_piv == 0) THEN

        PRINT*, "O sistema linear é indeterminado"

    ELSE IF (abs(det_a_piv) <= 1) THEN

        PRINT*, "O sistema linear é mal-condicionado"
    ELSE

        CALL substituicao(n, a_piv, p, b, x)
        PRINT*, "A solução do sistema é:"
        PRINT*, x
    END IF

END PROGRAM metodo_Gauss

SUBROUTINE matriz_a(n, a)
    IMPLICIT NONE

        INTEGER, INTENT(in) :: n
        REAL, DIMENSION(n,n), INTENT(inout) :: a


        INTEGER :: i, j !Indícios usados em loops para percorrer os arrays

    PRINT*, "Por favor digite os valores do elementos da matriz sistema linear seguindo pela ordem das linhas até o final:"

        DO i = 1, n
            DO j = 1, n
                READ*, a(i,j)
            END DO
        END DO

END SUBROUTINE matriz_a

SUBROUTINE vetor_b(n, b)
    IMPLICIT NONE

    INTEGER, INTENT(in) :: n
    REAL, DIMENSION(n), INTENT(inout) :: b

    INTEGER :: i 

    PRINT*, "Por favor entre com os elementos do vetor b:"

    DO i = 1, n
        READ*, b(i)
    END DO

END SUBROUTINE vetor_b

SUBROUTINE eliminacao(n, a, a_piv, p)
    IMPLICIT NONE


    INTEGER, INTENT(in) :: n
    REAL, DIMENSION(n, n), INTENT(in) :: a
    REAL, DIMENSION(n, n), INTENT(out) :: a_piv
    INTEGER, DIMENSION(n), INTENT(out) :: p


    INTEGER :: i, j, local, dim 
    REAL :: mult 

    DO i = 1, (n - 1)

        dim = n - 1

        CALL local_pivo(dim, a(i:n, i), local)

        a_piv(i, i:n) = a(local, i:n)
        a_piv(local, i:n) = a(i, i:n)

        p(i) = local

        DO j = (i + 1), n
                mult = (-1) * (a_piv(j,i) / a_piv(local,i))
                a_piv(j,i) = mult
                a_piv(j, j:n) = a_piv(j, j:n) + mult * a_piv(i, j:n)
        END DO

    END DO

END SUBROUTINE eliminacao

SUBROUTINE local_pivo(n, a, local)
    IMPLICIT NONE

    INTEGER, INTENT(in) :: n
    REAL, DIMENSION(n), INTENT(in) :: a
    INTEGER, INTENT(inout) :: local

    INTEGER :: i 

    local = 1


    DO i = 2, n
        IF ((ABS(a(i))) > ABS(a(local))) THEN
            local = i
        END IF
    END DO

END SUBROUTINE local_pivo

SUBROUTINE substituicao(n, a_piv, p, b, x)
    IMPLICIT NONE


    INTEGER, INTENT(in) :: n
    REAL, DIMENSION(n, n), INTENT(in) :: a_piv
    REAL, DIMENSION(n), INTENT(out) :: b, x
    INTEGER, DIMENSION(n), INTENT(in) :: p


    INTEGER :: i, j, k, l, pivo 
    REAL :: aux 


    DO i = 1, (n - 1)


        pivo = p(i)


        IF (pivo /= i) THEN
            aux = b(i)
            b(i) = b(pivo)
            b(pivo) = aux
        END IF

        DO j = (i + 1), n
            b(j) = a_piv(j, i) * b(j) + b(i)
        END DO

    END DO

    DO k = n, 1, -1
        IF (k == n) THEN
            x(n) = b(n) / a_piv(n, n)
        ELSE
            x(k) = (b(k) + a_piv(k, n) * x(n)) / a_piv(k, k)
            DO l = n, k, -1
                x(l) = x(l) + (a_piv(k, l) * x(l)) / a_piv(k, k)
            END DO
        END IF
    END DO



END SUBROUTINE substituicao

这是我要解决的系统

我的输入是：

4

4
3
2
2
2
1
1
2
2
2
2
4
6
1
1
2

5
8
3
1

我的输出是：

-40.5000000      -40.2500000      -3.75000000      -37.5000000 

但输出应该是：

6.500000
-44.000000
72.000000
-16.500000

标签： fortranlinear-algebranumerical-methodsfortran90