I have the following code, regarding Newton's method for calculating zeros in multivariable vectorial functions:
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
In simple precision the solutions are correctly displayed. The curious thing is that when I make the program write the vectors for each step, the final solution comes up correctly, but simply by eliminating the write statement from the do loop, the final solution comes up as not a number. Double precision seems to be messing with the results, as in simple precision it works properly, and I cannot find what the error might be. Perhaps I am writing the functions wrong in double precision?
Any help would be greatly appreciated.
This is the module, if helpful:
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
