Codigo Newton Fortran
Codigo Newton Fortran
Codigo Newton Fortran
(x
n
)
, n = 1, 2, . . .
Nuestra implementacin de la subrutina correspondiente requiere que se pase tambin como argumento
no slo la funcion f(x), sino tambin su derivada f
b
m
= a
m
b
k
= a
k
+ b
k+1
x
n
k = m 1, . . . , 0
siendo, entonces
b
0
= f(x
n
),
en tanto que f
(x
n
) es computada haciendo
c
m
= b
m
c
k
= b
k
+ c
k+1
x
n
k = m 1, . . . , 1
siendo, entonces
c
1
= f
(x
n
).
El mtodo de Newton se reduce as a
x
n+1
= x
n
b
0
c
1
El procedimiento resultante se conoce a menudo como mtodo de BirgeVieta.
Nuestra implementacin en la siguiente subrutina pasa los coecientes del polinomio en un arreglo a de
(m+1) componentes, siendo m el grado del polinomio (valor que tambin es requerido como argumento). En
la subrutina, para simplicar el tratamiento de los subndices, el arreglo es declarado con el ndice inferior 0,
no 1, de manera que a(0) = a
0
, a(1) = a
1
, . . . , a(M) = a
m
.
Cdigo 4. Implementacin del mtodo de BirgeVieta
SUBROUTINE birge_vieta(a,m,x0,n,tol,raiz,clave)
! ---------------------------------------------------
! METODO DE BIRGE-VIETA para resolver ECUACIONES
! ALGEBRAICAS: P (x) = 0 donde P es un polinomio de
! grado m de coeficientes reales.
! El mtodo se basa en el mtodo de Newton-Raphson
! implementando el esquema de Horner para la evalua-
! cin del polinomio y su derivada.
! ---------------------------------------------------
! Bloque de declaracin de argumentos
! ---------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: m ! Grado del polinomio
REAL(WP), INTENT(IN) :: a(0:m) ! Vector de m+1 elementos conteniendo
! los coeficientes del polinomio
REAL(WP), INTENT(IN) :: x0 ! Aproximacin inicial a la raz
REAL(WP), INTENT(IN) :: tol ! Tolerancia para el error relativo
INTEGER, INTENT(INOUT) :: n ! Limite de iteraciones/iteraciones realizadas
REAL(WP), INTENT(OUT) :: raiz ! Aproximacin a la raiz
INTEGER, INTENT(OUT) :: clave ! Clave de xito:
! 0 : xito
! >0 : iteraciones excedidas
! ---------------------------------------------------
! Bloque de declaracin de variables locales
! ---------------------------------------------------
INTEGER :: i, j
REAL(WP):: xx0,b,c
! ---------------------------------------------------
! Bloque de procedimiento
! ---------------------------------------------------
clave = 1
xx0 = x0
DO i=1,n
! -------------------------
! Esquema de Horner
! -------------------------
5
b = a(m)
c = a(m)
DO j=m-1,1,-1
b = b
*
xx0+a(j)
c = c
*
xx0+b
ENDDO
b = b
*
xx0+a(0)
! -------------------------
! Mtodo de Newton
! -------------------------
raiz = xx0 - b/c
IF (ABS(raiz-xx0) < tol
*
ABS(raiz)) THEN
clave = 0
n = i
EXIT
END IF
xx0 = raiz
END DO
END SUBROUTINE birge_vieta
5. Mtodo de la secante
El mtodo de la secante procede a partir de dos aproximaciones iniciales obteniendo la aproximacin
x
n+1
como la abscisa del punto de interseccin del eje x con la recta secante que pasa por los puntos
(x
n1
, f(x
n1
)) y (x
n
, f(x
n
)). La frmula de iteracin es entonces
x
n+1
= x
n
f(x
n
)
(x
n
x
n1
)
f(x
n
) f(x
n1
)
, n = 2, 3, . . .
El mtodo de la secante, si bien no converge tan rpido como Newton, tiene la gran ventaja de no requerir
la derivada de la funcin. Eso s, ahora se necesitan dos aproximaciones iniciales para arrancar el mtodo.
Cdigo 5. Implementacin del mtodo de la secante
SUBROUTINE secante(f,x0,x1,n,tol,raiz,clave)
! ---------------------------------------------------
! ALGORITMO DE LA SECANTE para encontrar una solucin
! de f(x)=0, siendo f una funcin continua, dada las
! aproximaciones iniciales x0 y x1.
! ---------------------------------------------------
! Bloque de declaracin de argumentos
! ---------------------------------------------------
IMPLICIT NONE
INTERFACE
FUNCTION f(x) ! Funcin que define la ecuacin
IMPORT :: WP
IMPLICIT NONE
REAL(WP) :: f
REAL(WP), INTENT(IN) :: x
END FUNCTION f
END INTERFACE
REAL(WP), INTENT(IN) :: x0,x1 ! Aproximaciones iniciales a la raz
INTEGER, INTENT(INOUT) :: n ! Lmite mximo de iteraciones/ iteraciones
realizadas
REAL(WP), INTENT(IN) :: tol ! Tolerancia para el error relativo
REAL(WP), INTENT(OUT) :: raiz ! Aproximacin a la raiz
INTEGER, INTENT(OUT) :: clave ! Clave de xito:
! 0 : xito
! >0 : iteraciones excedidas
! ---------------------------------------------------
! Bloque de declaracin de variables locales
! ---------------------------------------------------
6
INTEGER :: i
REAL(WP):: xx0, xx1, fx0, fx1
! ---------------------------------------------------
! Bloque de procesamiento
! ---------------------------------------------------
clave = 1
xx0 = x0
xx1 = x1
fx0 = f(x0)
fx1 = f(x1)
DO i= 2,n
raiz = xx1 - fx1
*
((xx1-xx0)/(fx1-fx0))
IF (ABS(raiz-xx1) < tol
*
ABS(raiz)) THEN
clave = 0
n = i
EXIT
ENDIF
xx0 = xx1
fx0 = fx1
xx1 = raiz
fx1 = f(raiz)
END DO
END SUBROUTINE secante
6. Iteracin de punto jo
El mtodo de punto jo requiere que se reescriba la ecuacin f(x) = 0 en la forma
x = (x)
y entonces, a partir de una aproximacin inicial x
0
, se obtiene la sucesin de aproximaciones x
1
, x
2
, ... segn
x
n+1
= (x
n
), n = 1, 2, . . .
En la siguiente implementacin, es importante recordar que la funcin que es pasada por argumento es
ahora (x) y no f(x), funcin que debe ser implementada por el usuario como una FUNCTION.
Cdigo 6. Implementacin del mtodo de punto jo
SUBROUTINE punto_fijo(f,x0,n,tol,raiz,clave)
! ---------------------------------------------------
! ALGORITMO DE PUNTO FIJO o DE APROXIMACIONES SUCESIVAS
! para encontrar una solucin de x=f(x) dada una
! aproximacin inicial x0.
! ---------------------------------------------------
! Bloque de declaracin de argumentos
! ---------------------------------------------------
IMPLICIT NONE
INTERFACE
FUNCTION f(x) ! Funcin que define la ecuacin
IMPORT :: WP
IMPLICIT NONE
REAL(WP) :: f
REAL(WP), INTENT(IN) :: x
END FUNCTION f
END INTERFACE
REAL(WP), INTENT(IN) :: x0 ! Aproximacin inicial a la raz
INTEGER, INTENT(INOUT) :: n ! Limite de iteraciones/ iteraciones
realizadas
REAL(WP), INTENT(IN) :: tol ! Tolerancia para el error relativo
REAL(WP), INTENT(OUT) :: raiz ! Aproximacin a la raiz
INTEGER, INTENT(OUT) :: clave ! Clave de xito:
! 0 : xito
7
! >0 : iteraciones excedidas
! ---------------------------------------------------
! Bloque de declaracin de variables locales
! ---------------------------------------------------
INTEGER :: i
REAL(WP) :: xx0
! ---------------------------------------------------
! Bloque de procesamiento
! ---------------------------------------------------
clave = 1
xx0 = x0
DO i=1,n
raiz = f(xx0)
IF (ABS(raiz-xx0) < tol
*
ABS(raiz)) THEN
clave = 0
n = i
EXIT
ENDIF
xx0 = raiz
END DO
END SUBROUTINE punto_fijo
7. Ejemplo
Como ejemplo del uso de nuestro mdulo roots consideremos el problema de determinar la raz de la
ecuacin
cos x x = 0.
Al gracar las curvas y = x e y = cos x, vemos en la g. 1, que la intereseccin de las mismas ocurre dentro
del intervalo [0.6, 0.8]. As pues, la raz buscada se encuentra en el interior de este intervalo y con l podemos
proceder a calcular la raz con el mtodo de biseccin asumiendo una tolerancia para el error absoluto de
1
2
10
6
y un nmero mximo de, digamos, 100 iteraciones. El siguiente programa implementa lo requerido.
Cdigo 7. Determinacin de la raz de la ec. cos x x = 0
PROGRAM ejemplo
USE precision, WP => DP
USE roots, ONLY: biseccion
IMPLICIT NONE
INTERFACE
FUNCTION f(x)
IMPORT :: WP
IMPLICIT NONE
REAL(WP) :: f
REAL(WP), INTENT(IN) :: x
END FUNCTION f
END INTERFACE
REAL(WP) :: a, b, tol, raiz
INTEGER :: n, clave
! Datos de iniciales
a = 0.6_WP
b = 0.8_WP
tol = 0.5E-6_WP
n = 100
! Determinar la raz
CALL biseccion(f,a,b,n,tol,raiz,clave)
IF (clave == 0) THEN
WRITE(
*
,
*
) 'Raz = ', raiz
WRITE(
*
,
*
) 'Nmero de iteraciones realizadas =', n
ELSE
WRITE(
*
,
*
) 'Error =', clave
ENDIF
END PROGRAM ejemplo
8
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4
x
cos(x)
Figura 1. Determinacin de una aproximacin inicial de la raz de la funcin
cos x x.
FUNCTION f(x)
! Funcin que define la ecuacin
USE precision, WP => DP
IMPLICIT NONE
REAL(WP) :: f
REAL(WP), INTENT(IN) :: x
f = cos(x) - x
RETURN
END FUNCTION f
La compilacin del mismo procede en la lnea de comandos como sigue:
$ gfortran -Wall -o ejemplo roots.f90 ejemplo.f90
Su ejecucin arroja entonces:
$ ./ejemplo
Raz = 0.73908500671386723
Nmero de iteraciones realizadas = 19
As, pues, la raz buscada, con seis decimales correctos, es 0.739085.
9