【发布时间】：2015-06-06 10:38:24
【问题描述】：
我有以下代码，关于牛顿在多变量矢量函数中计算零点的方法：
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
【问题讨论】：
-
Inverse_Matrices 模块是标准的/记录的吗？
-
是的，但只有当我更改主程序中的内容时，解决方案才会有所不同。整个模块也是双精度的。我几乎可以肯定这与双精度表达式在主程序中的编写方式有关。
-
write引起了一些明显不相关的问题，涉及数组越界或使用未初始化的值或类似问题，因此很可能与模块功能有关。尝试打开完全边界检查等。顺便说一句，最好将X上的索引反转为X(3,*)，以便每个向量的组件在内存中对齐。 -
h(x)是g(x)的派生词？ -
你能告诉我们预期的结果是什么，你得到了什么。
标签： fortran double-precision newtons-method