PDF de programación - Subrutinas en Fortran 95 para la resolución de ecuaciones no lineales de una variable

Subrutinas en Fortran 95 para la resolución
de ecuaciones no lineales de una variable

Pablo Santamaría
v0.3.1 (Mayo 2014)

1.

Introducción

En general, las raíces de una ecuación no lineal f(x) = 0 no pueden ser obtenidas por fórmulas explícitas
cerradas, con lo que no es posible obtenerlas en forma exacta. De este modo, para resolver la ecuación nos
vemos obligados a obtener soluciones aproximadas a través de algún método numérico.

Estos métodos son iterativos, esto es, a partir de una o más aproximaciones iniciales para la raíz, generan
una sucesión de aproximaciones x0, x1, x2, . . . que esperamos convergan al valor de la raíz α buscada. El
proceso iterativo se continúa hasta que la aproximación se encuentra próxima a la raíz dentro de una
tolerancia
> 0 preestablecida. Como la raíz no es conocida, dicha proximidad, medida por el error absoluto
|xn+1 − α|, no puede ser computada. Sin un conocimiento adicional de la función f(x) o su raíz, el mejor
criterio para detener las iteraciones consiste en proceder hasta que la desigualdad

|xn+1 − xn|

|xn+1|

<


se satisfaga, dado que esta condición estima en cada paso el error relativo. Ahora bien, puede ocurrir en ciertas
circunstancias que la desigualdad anterior nunca se satisfaga, ya sea por que la sucesión de aproximaciones
diverge o bien que la tolerancia escogida no es razonable. En tal caso el método iterativo no se detiene
nunca. Para evitar este problema consideramos además un número máximo de iteraciones a realizarse. Si
este número es excedido entonces el problema debe ser analizado con más cuidado.

¿Cómo se escoge un valor correcto para las aproximaciones iniciales requeridas por los métodos? No existe
una respuesta general para esta cuestión. Para el método de bisección es suficiente conocer un intervalo que
contenga la raíz, pero para el método de Newton, por ejemplo, la aproximación tiene que estar suficientemente
cercana a la raíz para que funcione 1. En cualquier caso primeras aproximaciones iniciales para las raíces
pueden ser obtenidas graficando la función f(x).

En las siguientes secciones presentamos implementaciones de los métodos numéricos usuales como
subrutinas Fortran. Con el fin de proporcionar subrutinas de propósito general, las mismas tienen entre sus
argumentos a la función f involucrada, la cual pude ser entonces implementada por el usuario como un
subprograma FUNCTION con el nombre que quiera. Otros argumentos que requieren estas subrutinas son los
valores para las aproximaciones iniciales que necesite el método, una tolerancia para la aproximación final
de la raíz y un número máximo de iteraciones. La raíz aproximada es devuelta en otro de los argumentos.
Dado que el método puede fallar utilizamos también una variable entera como clave de error para advertir al
programa principal. Por convención tomaremos que si dicha clave es igual a cero, entonces el método funcionó
correctamente y el valor devuelto es la raíz aproximada dentro de la tolerancia preescrita. En cambio, si la
clave de error es distinta de cero, entonces ocurrió un error. La naturaleza del error dependerá del método,
pero un error común a todos ellos es que el número máximo de iteraciones fue alcanzado.

Con el fin de aprovechar la capacidad de Fortran 95 de detectar errores de tipo en los argumentos
al llamar a las subrutinas, debemos hacer explícita la interface de las mismas. La forma más simple y
poderosa de efectuar ésto consiste en agrupar las mismas en un módulo, al que denominaremos roots. En
nuestra implementación todas las cantidades reales serán de la clase de doble precisión, la cual, para máxima
flexibilidad, está definida en forma paramétrica utilizando el modulo intrínseco iso_fortran_env.

1De hecho, efectuando algunas iteraciones del método de bisección podemos obtener una buena aproximación para iniciar el

método de Newton.

1

Código 1. Implementación del módulo roots

MODULE roots

USE, intrinsic :: iso_fortran_env, ONLY: WP => REAL64
IMPLICIT NONE

CONTAINS

SUBROUTINE biseccion(f,a,b,n,tol,raiz,clave)

...

END SUBROUTINE biseccion

SUBROUTINE newton(f,df,x0,n,tol,raiz,clave)

...

END SUBROUTINE newton

SUBROUTINE birge_vieta(a,m,x0,n,tol,raiz,clave)

...

END SUBROUTINE birge_vieta

SUBROUTINE secante(f,x0,x1,n,tol,raiz,clave)

...

END SUBROUTINE secante

SUBROUTINE punto_fijo(f,x0,n,tol,raiz,clave)

...

END SUBROUTINE punto_fijo

END MODULE roots

El código correspondiente a cada subrutina, que debe insertarse donde se encuentran los puntos suspensivos,

será discutido por separado en las siguientes secciones.

2. Método de bisección

El método de bisección comienza con un intervalo [a, b ] que contiene a la raíz. Entonces se computa el
punto medio x0 = (a + b)/2 del mismo y se determina en cual de los dos subintervalos [a, x0] o [x0, b ] se
encuentra la raíz analizando el cambio de signo de f(x) en los extremos. El procedimiento se vuelve a repetir
con el nuevo intervalo así determinado. Es claro que la raíz es acotada en cada paso por el intervalo así
generado y que una estimación del error cometido en aproximar la raíz por el punto medio de dicho intervalo
es igual a la mitad de la longitud del mismo. Esta estimación es utilizada, en la siguiente implementación del
método, como criterio de paro para la sucesión de aproximaciones.

Código 2. Implementación del método de bisección

SUBROUTINE biseccion(f,a,b,n,tol,raiz,clave)

! ---------------------------------------------------
! METODO DE BISECCION para encontrar una solución
! de f(x)=0 dada la función continua f en el intervalo
! [a,b] donde f(a) y f(b) tienen signos opuestos.
! ---------------------------------------------------
! Bloque de declaración de argumentos
! ---------------------------------------------------
INTERFACE

REAL(WP) FUNCTION f(x)

! Función que define la ecuación

IMPORT :: WP
IMPLICIT NONE
REAL(WP), INTENT(IN) :: x

END FUNCTION f

END INTERFACE
REAL(WP), INTENT(IN)
REAL(WP), INTENT(IN)

:: a
:: b

! Extremo izquierdo del intervalo inicial
! Extremo derecho del intervalo inicial

2

INTENT(INOUT) :: n

INTEGER,
REAL(WP), INTENT(IN)
REAL(WP), INTENT(OUT)
INTEGER,
INTENT(OUT)

! Límite de iteraciones/iteraciones realizadas
! Tolerancia para el error absoluto
! Aproximación a la raiz

:: tol
:: raiz
:: clave ! Clave de éxito:

!
!
!

0 : éxito

>0 : iteraciones excedidas
<0 : no se puede proceder (f de igual signo

en a y b)

:: i

! ---------------------------------------------------
! Bloque de declaración de variables locales
! ---------------------------------------------------
INTEGER
REAL(WP) :: xl, xr, error
! ---------------------------------------------------
! Bloque de procesamiento
! ---------------------------------------------------
clave = 1
= a
xl
xr
= b
IF (SIGN(1.0_WP,f(xl))*SIGN(1.0_WP,f(xr)) > 0.0_WP) THEN

clave = -1
RETURN

ENDIF
DO i=1,n

error = (xr-xl)*0.5_WP
raiz
IF (error < tol) THEN

= xl + error

clave = 0
n = i
EXIT

ENDIF
IF ( SIGN(1.0_WP,f(xl))* SIGN(1.0_WP,f(raiz)) > 0.0_WP) THEN

xl = raiz

ELSE

xr = raiz

ENDIF

ENDDO
! ---------------------------------------------------

END SUBROUTINE biseccion

3. Método de Newton–Raphson

El método de Newton comienza con una aproximación inicial x0 dada, a partir de la cual se genera la
sucesión de aproximaciones x1, x2, . . ., siendo xn+1 la abscisa del punto de intersección del eje x con la recta
tangente a f(x) que pasa por el punto (xn, f(xn)). Esto conduce a la fórmula de iteración

xn+1 = xn − f(xn)

f0(xn) , n = 1, 2, . . .

Nuestra implementación de la subrutina correspondiente requiere que se pase también como argumento
no sólo la funcion f(x), sino también su derivada f0(x), la cual debe ser implementada por el usuario, al
igual que f(x), como una FUNCTION.

3

Código 3. Implementación del método de Newton–Raphson

SUBROUTINE newton(f,df,x0,n,tol,raiz,clave)

! ---------------------------------------------------
! Metodo DE NEWTON-RAPHSON para encontrar una
! solución de f(x)=0 dada la función derivable
! f y una aproximación inicial x0.
! ---------------------------------------------------
! Bloque de declaración de argumentos
! ---------------------------------------------------
INTERFACE

REAL(WP) FUNCTION f(x)

! Función que define la ecuación

IMPORT :: WP
IMPLICIT NONE
REAL(WP), INTENT(IN) :: x

END FUNCTION f
REAL(WP) FUNCTION df(x)

IMPORT :: WP
IMPLICIT NONE
REAL(WP), INTENT(IN) :: x

END FUNCTION df

! Derivada de la función
! que define a la ecuación

INTENT(INOUT) :: n

:: x0

END INTERFACE
REAL(WP), INTENT(IN)
INTEGER,
REAL(WP), INTENT(IN)
REAL(WP), INTENT(OUT)
INTEGER,
INTENT(OUT)

! Aproximación inicial a la raíz
! Límite de iteraciones/iteraciones realizadas
! Tolerancia para el error relativo
! Aproximación a la raíz

:: tol
:: raiz
:: clave ! Clave de éxito:

!
!

0 : éxito

>0 : iteraciones excedidas

::

i

! ---------------------------------------------------
! Declaración de variables locales
! ---------------------------------------------------
INTEGER
REAL(WP) :: xx0
! ---------------------------------------------------
! Bloque de procesamiento
! ---------------------------------------------------
clave = 1
xx0
DO i=1,n

= x0

= xx0 - f(xx0)/df(xx0)

raiz
IF (ABS(raiz-xx0) < tol*ABS(raiz) ) THEN

clave = 0
n
= i
EXIT

ENDIF
xx0 = raiz

END DO
! ---------------------------------------------------

END SUBROUTINE newton

4. Método de Newton para ecuaciones algebraicas

En el caso particular en que f(x) es un polinomio, el método de Newton puede ser eficientemente
implementado si la evaluación de f(xn) (y su derivada) es realizada por el método iterativo de Horner. En
efecto, supongamos que f(x) es un polinomio de grado m:

la evaluación de f(xn) por la regla de Horner procede computando

f(x) = a0 + a1x + a2x2 + ··· + amxm,

(

bm = am
bk = ak + bk+1xn

k = m − 1, . . . , 0

4

siendo, entonces

en tanto que f0(xn) es computada haciendo
cm = bm
ck = bk + ck+1xn

(

k = m − 1, . . . , 1

b0 = f(xn),

siendo, entonces

c1 = f0(xn).

El método de Newton