FORTRAN-正确的数值以简单的精度计算，但错误的解决方案以双倍精度计算

Question

我有以下代码，关于在多变量矢量函数中计算零的牛顿方法：

   program Newton_Method_R_N_R_N_1_Numeric
use Inverse_Matrices
implicit none
integer  :: i, j                            !INDEX
real (kind = 8) :: X(0:501,3)               !VECTORS IN NEWTON METHOD
real (kind = 8) :: A(3,3)                   !JACOBI'S MATRIX FOR EACH VECTOR
integer , parameter :: n = 3                !DIMENSION OF JACOBI'S MATRIZ
real (kind = 8) :: Inv(3,3)                 !INVERSE JACOBI'S MATRIX
real (kind = 8), parameter :: eps = 1D-8    !NEAR-ZERO VALUE
real (kind = 8) :: det                      !DETERMINANT
X(0,:) = (/2.D0,2.D0,2.D0/)
do i = 0, 500
    A = h(X(i,:))
    call LU_Factorization (n, A)
    call Determinant (n, A, det)
    if (abs(det)<eps) then
        exit
    end if
    call Inverse_Matrix (n, A, Inv)
    X(i + 1,:) = X(i,:) - matmul(Inv,g(X(i,:)))
end do
write(*,*) ' A zero of the function is: ', X(i,:)
    read(*,*)
  contains
function g(t)
    implicit none
    real (kind = 8), intent(in) :: t(3)
        real (kind = 8) :: g(3)
        g(1) = t(1)**2 - t(2)**2 + 1.0D0
        g(2) = 2.0D0*t(1)*t(2)
        g(3) = t(3)**2 - 2.0D0
end function g
function g_1(t_1)
    implicit none
    real (kind = 8), intent(in) :: t_1(3)
        real (kind = 8) :: g_1
        g_1 = t_1(1)**2 - t_1(2)**2 + 1.0D0
end function g_1
function g_2(t_2)
    implicit none
    real (kind = 8), intent(in) :: t_2(3)
        real (kind = 8) :: g_2
        g_2 = (2.0D0)*t_2(1)*t_2(2)
end function g_2
function g_3(t_3)
    implicit none
    real (kind = 8), intent(in) :: t_3(3)
        real (kind = 8) :: g_3
        g_3 = t_3(3)**2 - 2.0D0
end function g_3
function h(q)
    implicit none
    real (kind = 8), intent(in) :: q(3)
        real (kind = 8) :: h(3,3)
        real (kind = 8), parameter :: dx = 1D-4 
        real (kind = 8) :: a(3)
        real (kind = 8) :: b(3)
        real (kind = 8) :: c(3) 
        a(1) = dx
        a(2) = 0.D0
        a(3) = 0.D0
        b(1) = 0.D0
        b(2) = dx
        b(3) = 0.D0
        c(1) = 0.D0
        c(2) = 0.D0
        c(3) = dx
        h(1,1) = ((g_1(q + a)) - g_1(q)) / dx
        h(1,2) = ((g_1(q + b)) - g_1(q)) / dx
        h(1,3) = ((g_1(q + c)) - g_1(q)) / dx
        h(2,1) = ((g_2(q + a)) - g_2(q)) / dx
        h(2,2) = ((g_2(q + b)) - g_2(q)) / dx
        h(2,3) = ((g_2(q + c)) - g_2(q)) / dx
        h(3,1) = ((g_3(q + a)) - g_3(q)) / dx
        h(3,2) = ((g_3(q + b)) - g_3(q)) / dx
        h(3,3) = ((g_3(q + c)) - g_3(q)) / dx
end function h
end program

简单精确地显示解决方案。奇怪的是，当我让程序为每一步写向量时，最终的解决方案是正确的，但简单地从do循环中删除write语句，最终的方案就不是一个数字。双精度似乎扰乱了结果，因为在简单精度中它可以正常工作，我找不到可能的错误。也许我在双精度中写错了函数？

如有任何帮助，我们将不胜感激。

---

如果有用的话，这就是模块：

module Matrices_Inversas_2
    implicit none
contains    
    subroutine Factorizacion_LU (n, A)
    implicit none
    integer, intent(in) :: n            !DIMENSIÓN DEL SISTEMA (matriz de coef.)
    real (kind = 8), intent(inout) :: A(n,n)    !MATRIZ DE COEFICIENTES DEL SISTEMA
        integer :: i, j, k              !ÍNDICES
        real (kind = 8) :: V(n), W(n)               !VECTORES AUXILIARES
        !FACTORIZACIÓN LU
        do k = 1, n
            !COLUMNAS
            do i = 1, (k - 1)
                V(i) = A(k,i)
            end do
            do j = k, n
                do i = 1, (k - 1)
                    W(i) = A(i,j)
                end do
                A(k,j) = A(k,j) - dot_product(V,W)
            end do
            !FILAS
            do i = 1, (k - 1)
                W(i) = A(i,k)
            end do
            do i = (k + 1), n
                do j = 1, (k - 1)
                    V(j) = A(i,j)
                end do
                A(i,k) = (A(i,k) - dot_product(V,W)) / A(k,k)
            end do
        end do
    end subroutine Factorizacion_LU
    subroutine Sustitucion_Directa (n, A, B, C, V)
    implicit none
    integer, intent(in) :: n                !DIMENSIÓN DEL SISTEMA  
    real (kind = 8), intent(inout) :: A(n,n)        !MATRICES DEL SISTEMA
    real (kind = 8), intent (in) :: B(n)            !MATRIZ DE TÉRMINOS INDEPENDIENTES
    real (kind = 8), intent(out) :: C(n)            !MATRIZ DE INCÓGNITAS
    real (kind = 8), intent(out) :: V(n)            !VECTOR AUXILIAR
        integer :: i, j                     !ÍNDICES EN LAS MATRICES
        real (kind = 8) :: s                            !SUMATORIO
        do i = 1, n
            V(i) = A(i,i)
            A(i,i) = 1.D0
        end do
        !SUSTITUCIÓN DIRECTA
        C(1) = B(1) / A(1,1)
            do i = 2, n
                s = 0.D0
                do j = 1, (i - 1)
                    s = s + A(i,j)*C(j)
                end do
                C(i) = (B(i) - s) / A(i,i)
            end do
    end subroutine Sustitucion_Directa
    subroutine Sustitucion_Regresiva (n, A, C, X, V)
    implicit none
    integer, intent(in) :: n            !DIMENSIÓN DEL SISTEMA
    real (kind = 8), intent(inout) :: A(n,n)    !MATRIZ DE COEFICIENTES
    real (kind = 8), intent(in) :: C(n)         !MATRIZ DE TÉRMINOS INDEPENDIENTES
    real (kind = 8), intent(out) :: X(n)        !MATRIZ DE INCÓGNITAS
    real (kind = 8), intent(in) :: V(n)         !VECTOR AUXILIAR
        integer :: i,j                  !ÍNDICES EN LAS MATRICES
        real (kind = 8) :: s                        !SUMATORIO EN SUSTITUCIÓN REGRESIVA
        do i = 1, n
            A(i,i) = V(i)
        end do
        !SUSTITCUIÓN REGRESIVA
        X(n) = C(n) / A(n,n)
        do i = (n - 1), 1, -1
            s = 0.D0
            do j = (i + 1), n
                s = s + A(i,j) * X(j)
            end do
            X(i) = (C(i) - s) / A(i,i)
        end do
    end subroutine Sustitucion_Regresiva
    subroutine Matriz_inversa (n, A, Inv)
    implicit none
    integer, intent(in) :: n            !DIMENSIÓN DEL SISTEMA
    real (kind = 8), intent(inout) :: A(n,n)    !MATRIZ DE LA QUE HALLAR LA INVERSA
    real (kind = 8), intent(out) :: Inv(n,n)    !INVERSA
        integer :: i, j                 !ÍNDICES
        real (kind = 8) :: B(n)                     !MATRIZ DE TÉRMINOS INDEPENDIENTES
        real (kind = 8) :: C(n)                     !MATRIZ INTERMEDIA DE INCÓGNITAS
        real (kind = 8) :: X(n)                     !MATRIZ DE INCÓGNITAS
        real (kind = 8) :: V(n)                     !VECTOR AUXILIAR
        do i = 1, n
            B = 0.D0
            B(i) = 1.D0
            call Sustitucion_Directa (n, A, B, C, V)
            call Sustitucion_Regresiva (n, A, C, X, V)
            do j = 1, n
                Inv(j,i) = X(j)
            end do
        end do
    end subroutine Matriz_inversa
    subroutine Determinante (n, A, det)
    integer, intent(in) :: n        !DIMENSIÓN DEL SISTEMA
    real (kind = 8), intent(inout) :: A(n,n)!MATRIZ DE ENTRADA
    real(kind = 8), intent(out) :: det      !VALOR DEL DETERMINANTE
        integer :: i                            !ÍNDICE
        det = 1.D0
        do i = 1, n
            det = det * A(i,i)
        end do
    end subroutine Determinante
    end module

Answer 1

您为我们提供了一个模块Matrices_Inversas_2，其中显然包含西班牙语命名的过程，但您的程序使用了一个名为Inverse_Matrices的模块，因此我们仍然没有一套完全可编译的代码。

您也没有指定您认为"正确显示"的内容。

我对您的模块代码进行了明显的更改，以便使您的过程名称匹配并使用nagfor -kind=byte -u -C=all -C=undefined -gline进行编译。正在运行，我得到

mod.f90:
[NAG Fortran Compiler normal termination]
test.f90:
Warning: test.f90, line 93: Unused local variable J
[NAG Fortran Compiler normal termination, 1 warning]
Loading...
Runtime Error: mod.f90, line 27: Reference to undefined variable V
Program terminated by fatal error
mod.f90, line 27: Error occurred in MATRICES_INVERSAS_2:FACTORIZACION_LU
test.f90, line 18: Called by NEWTON_METHOD_R_N_R_N_1_NUMERIC
Abort (core dumped)

您正在尝试一个由两个大小的n数组（V和W）组成的dot_product，从您的代码中可以看出，只有第一个(k-1)元素被初始化。如果我用dot_product(V(1:(k-1)),W(1:(k-1)))替换你的两个dot_product调用，那么程序运行到完成并显示

A zero of the function is:    0.0000000000000000   1.0000000000000000   1.4142135623730951