¡Descarga pdf de mates y más Apuntes en PDF de Matemáticas solo en Docsity!
5 Ceros de funciones 89
5 Ceros de funciones
Ob jetivos
Describir tres tecnicas numericas iterativas para hallar ceros de funciones (f (x) = 0): meto do de la biseccion, meto do de Newton y meto do de la secante.
Estudiar y comparar los tres meto dos mediante algunos ejemplos numericos.
Explicar las funciones externas FUNCTION de FORTRAN.
5.1 Intro duccion
Muchos problemas pueden mo delarse matematicamente como una ecuacion
f (x) = 0 (5:1)
donde f es una cierta funcion de una variable x. Se trata pues de hallar los valores de x que
satisfacen la ecuacion 5.1. Estos valores se llaman ceros de la funcion f o races de la ecuacion
f (x) = 0, y se denotan p or x. Gra camente, los ceros de una funcion son los puntos de
interseccion de la gra ca y = f (x) con el eje de las x.
Para algunos casos sencillos, la ecuacion 5.1 puede resolverse analticamente. Supongase,
p or ejemplo, que f es un p olinomio de segundo grado, f (x) = ax^2 + bx + c. Entonces, el n umero
de ceros (reales) dep ende del valor del discriminante = b^2 4 ac; para > 0, la funcion f tiene dos ceros x = ( b
p
)= 2 a.
En un problema mas general, si f es una funcion cualquiera, la ecuacion 5.1 no puede
resolverse analticamente. De hecho, ni siquiera se sab e a priori cuantos ceros tiene f : >varios,
uno, ninguno? En estos casos, es necesario utilizar una tecnica numerica iterativa: a partir
de una aproximacion inicial x^0 a un cero x de f , se construye iterativamente una sucesion de
90 Metodos numericos
aproximaciones fxk^ g. El sup erndice k es el contador de iteraciones: en la primera iteracion,
se calcula x^1 ; en la segunda, x^2 , y as sucesivamente. El pro ceso iterativo se detiene cuando,
para un cierto valor de k , el valor xk^ es una aproximacion su cientemente buena a x.
Puede verse pues que para obtener numericamente un cero de f hay que resp onder las tres
preguntas siguientes:
1. >Como se elige la aproximacion inicial x^0?
Como se construye la sucesion fxk^ g de aproximaciones?
3. >Como se decide si xk^ es una aproximacion su cientemente buena a x?
Estas preguntas se resp onden a lo largo del captulo. De momento, y para terminar este apartado de intro duccion, se presentan dos ejemplos de la ecuacion 5.1.
5.1.1 Calculo de races cuadradas
Un ingeniero necesita calcular la raz cuadrada x de un n umero s, x =
p
s, haciendo
unicamente op eraciones aritmeticas elementales (suma, resta, pro ducto y division). Esta es la situacion real en el dise ~no de algunos ordenadores, puesto que solo estas cuatro op eraciones estan incorp oradas a nivel de hardware, y las demas op eraciones deb en hacerse a partir de ellas.
Dado que no se puede calcular directamente la raz cuadrada
p
s, se utiliza una estrategia
alternativa. Elevando la expresion x =
p
s al cuadrado y pasando s a la izquierda de la igualdad,
puede escribirse f (x) = x^2 s = 0 (5:2)
En la ecuacion 5.2 queda claro que el calculo de
p
s equivale a obtener el cero de la funcion
f (x) = x^2 s. En otras palabras, se trata de hallar la interseccion de la gra ca y = f (x) con
el eje de las x (vease la gura 5.1). Como se vera en los apartados siguientes, la ecuacion 5.
puede resolverse de manera iterativa empleando unicamen te las cuatro op eraciones aritmeticas elementales.
5.1.2 Como jugar al billar en una mesa circular
La ultima mo da entre los a cionados al billar es la mesa circular (vease la gura 5.2). Para
los principiantes, el juego consiste simplemente en golp ear la b ola Q con la b ola P despues de
un impacto I en la banda. Los parametros del problema pueden verse en la gura 5.2: la mesa
tiene radio R , la p osicion de las b olas P y Q queda determinada p or las co ordenadas cartesianas
(xP ; yP ) y (xQ ; yQ ), y el punto de impacto I viene de nido p or el angulo .
92 Metodos numericos
Mediante consideraciones geometricas sencillas (que se dejan como ejercicio al lector intere-
sado) puede verse que los valores de que prop orcionan los puntos de impacto I son las races
de la ecuacion
f ( ) =
xP sin yP cos
p
(R cos xP )^2 + (R sin yP )^2
xQ sin yQ cos
p
(R cos xQ )^2 + (R sin yQ )^2
Notese que es la unica incognita de la ecuacion 5.3. Los valores de R , xP , yP , xQ e yQ
son datos del problema. Para resolver la ecuacion 5.3 hay que utilizar una tecnica numerica
iterativa, que construya una sucesion f k^ g de aproximaciones a un cero de la funcion f.
5.2 Meto do de la biseccion
La primera tecnica iterativa para hallar ceros de funciones que se presenta aqu es el metodo de la biseccion. Se ilustrara mediante el calculo de
p
2 a partir de op eraciones aritmeticas
elemenp tales. Se trata, pues, de tomar s = 2 en la ecuacion 5.2. Tal como ya se ha comentado,
2 es el cero de la funcion f (x) = x^2 2 (vease la gura 5.3).
Fig. 5.3 Calculo de
p
2 por el metodo de la biseccion
5 Ceros de funciones 93
El meto do de la biseccion consiste en:
1. Inicializar el contador de iteraciones a cero (k = 0).
2. Elegir una aproximacion inicial x^0 y otro valor a de manera que el intervalo que de nen estos
dos puntos (que sera
x^0 ; a
para x^0 < a y
a; x^0
para x^0 > a) contenga el cero buscado x
y solo ese cero. Para el problema que se esta estudiando es muy sencillo ver que la funcion
f (x) = x^2 2 tiene un unico cero x =
p
2 en el intervalo [1; 2]. Se toma, p or tanto, x^0 = 1
y a = 2. Para un problema general, se utiliza el siguiente control para elegir x^0 y a: si f
es una funcion continua y el intervalo de extremos x^0 y a contiene un cero de f , entonces
f (x^0 ) y f (a) tienen signos distintos (es decir, f (x^0 )f (a) < 0). En la gura 5.3 puede verse
como, efectivamente, f (x^0 ) < 0 y f (a) > 0.
3. Tomar el punto medio del intervalo, xk^ +1^ = (xk^ + a)=2, como siguiente aproximacion a x.
Notese que xk^ +1^ divide al intervalo de extremos xk^ y a en dos nuevos intervalos con la mitad
de longitud. Por este motivo, se habla de meto do de la biseccion.
4. Decidir si xk^ +1^ es una aproximacion su cientemente buena a x. En caso a rmativo, detener
el pro ceso iterativo y tomar x xk^ +1^. En caso negativo, seguir iterando. Para tomar la
decision es necesario emplear alg un criterio de convergencia (ver apartado 5.3).
5. Detectar cual de los dos intervalos obtenidos en el paso 3 contiene x. Puede hacerse de
manera sistematica, sin necesidad de dibujar la gra ca de la funcion, a partir del signo de
f (xk^ +1^ ). Si f (xk^ +1^ ) y f (xk^ ) tienen signos distintos, entonces x esta en el intervalo de
extremos xk^ y xk^ +1^. Si, p or el contrario, el cambio de signo de f se pro duce entre xk^ +1^ y
a, el cero x esta entre xk^ +1^ y a. Por ultimo, puede o currir que f (xk^ +1^ ) = 0; en este caso,
xk^ +1^ es el cero x y se detiene el pro ceso iterativo de biseccion.
- Tomar el intervalo escogido en el paso 5 como nuevo intervalo de traba jo. Para obtener
un algoritmo mas compacto y facilitar su programacion, interesa denotar p or a uno de los
extremos del intervalo (el otro extremo es xk^ +1^ ) durante to do el pro ceso. Para ello se adopta
el siguiente criterio:
si
f (xk^ +1^ )f (xk^ ) < 0 ; a xk
f (xk^ +1^ )f (xk^ ) 0 ; a a
7. Incrementar en 1 el contador de iteraciones (k k + 1) y volver al paso 3.
En resumen, en el meto do de la biseccion se parte de un intervalo inicial que contiene el
cero x , y se va sub dividiendo este intervalo hasta \encerrar" a x en un intervalo tan p eque ~no
como se desee.
Observese que el algoritmo que se acaba de presentar ab orda las tres cuestiones planteadas
en el apartado 5.1: eleccion de la aproximacion inicial x^0 (paso 2), construccion de la sucesion
fxk^ g de aproximaciones a x (pasos 3, 5, 6 y 7) y nalizacion de las iteraciones (paso 4).
5 Ceros de funciones 95
Notese que, con estas tolerancias, se obtienen o cho cifras signi cativas en los resultados. Puede
verse en la tabla 5.1 como el intervalo inicial [1; 2] se va sub dividiendo hasta llegar, despues
de 29 iteraciones, a x^29 = 1 :4142136, que se toma como aproximacion a
p
- Otra entrada del programa es la variable MAXITER, que representa el n umero de iteraciones que se desea realizar como maximo. En caso de alcanzarse este valor, el programa nalizara sin ning un mensa je de convergencia. Habitualmente, al detectarse este fenomeno se deb e ejecutar de nuevo el programa inicializandose el meto do mas cerca de la solucion.
Considerese ahora al problema del billar circular. Tomando como punto de partida el programa 5.1, se escrib e el programa 5.2 (apartado 5.7) que resuelve la ecuacion 5.3 en lugar de la ecuacion 5.2. Se toman los valores R = 1, xP = 0 :6, yP = 0, xQ = 0 :6, yQ = 0
(vease la gura 5.4). Un punto de impacto I viene dado entonces p or = =2. Se puede
capturar esta solucion tomando x^0 = 1 : 5 < =2, a = 1 : 6 > =2. Si se mantienen las tolerancias
tolx = T O Lf = 0 : 5 10 ^8 , el meto do de la biseccion prop orciona los resultados de la tabla 5.2. Efectivamente, en 23 iteraciones se obtiene 23 = 1 : 5707963 =2.
Fig. 5.4 Problema del bil lar para R = 1 , xP = 0 : 6 , yP = 0 , xQ = 0 : 6 e yQ = 0
En cada iteracion del meto do de la biseccion es necesario evaluar f (xk^ ) (paso 5). Para ello,
en el programa FORTRAN se traba ja con una funcion externa (FUNCTION). Se ha visto en el captulo 3 que el FORTRAN disp one de una biblioteca de funciones intrnsecas (trigonometricas, logartmicas, exp onenciales, etc.) ya incorp oradas. Las FUNCTIONs, p or el contrario, son funciones de usuario que se pueden de nir a voluntad para resolver un problema concreto.
Por ejemplo, para el problema del billar circular se ha de nido f ( ) (ver ecuacion 5.3) como
una FUNCTION. En el apartado 5.7 se explica como de nir y traba jar con funciones externas FUNCTION.
Problema 5.1: a) Determinar, p or simple insp eccion visual, cuales son los demas ceros de la
funcion f ( ) (ecuacion 5.3) para los datos de la gura 5.4 (Pista: hay un total
96 Metodos numericos
de cuatro ceros). >Son to dos estos ceros soluciones validas desde un punto de vista fsico?
b) Veri car que los valores obtenidos en el apartado a son realmente ceros de f ( )
utilizando el meto do de la biseccion. Justi car razonadamente los resultados obtenidos.
Tabla 5.2 Obtencion de la solucion = = 2 por el metodo de la biseccion
Iteracion Extremo a Aproximacion x f(x) Error relativo en x ========= ========== ============== ========== =================== 0 1.6000000 1.5000000 3.211D-02 -3.226D- 1 1.6000000 1.5500000 9.440D-03 -1.587D- 2 1.5500000 1.5750000 -1.908D-03 8.000D- 3 1.5750000 1.5625000 3.766D-03 -3.984D- 4 1.5750000 1.5687500 9.290D-04 -1.988D- 5 1.5687500 1.5718750 -4.897D-04 9.950D- 6 1.5718750 1.5703125 2.196D-04 -4.973D- 7 1.5703125 1.5710938 -1.350D-04 2.487D- 8 1.5710938 1.5707031 4.231D-05 -1.243D- 9 1.5707031 1.5708984 -4.635D-05 6.217D- 10 1.5707031 1.5708008 -2.022D-06 3.109D- 11 1.5708008 1.5707520 2.014D-05 -1.554D- 12 1.5708008 1.5707764 9.061D-06 -7.771D- 13 1.5708008 1.5707886 3.519D-06 -3.886D- 14 1.5708008 1.5707947 7.486D-07 -1.943D- 15 1.5707947 1.5707977 -6.368D-07 9.714D- 16 1.5707977 1.5707962 5.592D-08 -4.857D- 17 1.5707962 1.5707970 -2.904D-07 2.429D- 18 1.5707962 1.5707966 -1.173D-07 1.214D- 19 1.5707962 1.5707964 -3.067D-08 6.071D- 20 1.5707964 1.5707963 1.263D-08 -3.036D- 21 1.5707963 1.5707963 -9.020D-09 1.518D- 22 1.5707963 1.5707963 1.803D-09 -7.589D- 23 1.5707963 1.5707963 -3.608D-09 3.795D-
Convergencia en la iteracion 23 Solucion para theta= 1.
5.3 Criterios de convergencia
Se dice que la sucesion fxk^ g converge a x si
lim k!
xk^ = x
98 Metodos numericos
Debido a la aproximacion hecha en la ecuacion 5.7, puede o currir en algunos casos que la
condicion de convergencia dada p or la ecuacion 5.8 se cumpla estando xk^ lejos de x. Para
evitar estos problemas, se complementa el criterio relativo en x (ecuacion 5.8) con un criterio
absoluto en f.
Para ello, basta darse cuenta de que, si fxk^ g converge a x , se veri ca tambien
lim k!
f (xk^ ) = 0
puesto que f (x ) = 0. Esto signi ca que f (xk^ ) es directamente el error absoluto en f. Para
aceptar xk^ como aproximacion nal a x se exige |ademas de la condicion 5.8| que este error
absoluto sea, en valor absoluto, inferior a una tolerancia T O Lf preestablecida:
jf (xk^ )j < T O Lf
5.4 Meto do de Newton
Se ha comprobado en el apartado 5.2 que el meto do de la biseccion es una tecnica robusta
para hallar ceros de funciones: basta que f sea una funcion continua, que el intervalo inicial
de nido p or x^0 y a contenga un cero x y que f tenga signos distintos en los extremos del
intervalo (f (x^0 )f (a) < 0) para garantizar que el meto do va \encerrando" a x en un intervalo
cada vez mas p eque ~no. La longitud del intervalo nal puede controlarse mediante las tolerancias de convergencia. Sin embargo, los dos ejemplos del apartado 5.2 tambien muestran que la biseccion es una tecnica lenta: para las tolerancias exigidas, han sido necesarias entre 20 y 30 iteraciones para alcanzar la convergencia.
Una tecnica mas rapida (aunque no tan robusta, como se vera) es el metodo de Newton. Se hara en primer lugar una deduccion analtica del meto do y luego una deduccion gra ca.
5.4.1 Deduccion analtica del meto do de Newton
Supongase que xk^ es una aproximacion a un cero x de una cierta funcion f. Puesto que
xk^6 = x , resulta f (xk^ ) 6 = 0
Dado que xk^ no es el cero x buscado, se intenta que la siguiente aproximacion xk^ +1^ s lo
sea. Para ello se de ne xk^ +1^ como
xk^ +1^ = xk^ + xk^ +1^ (5:10)
donde xk^ +1^ es la correccion que se hace a xk^ para obtener xk^ +1^. El criterio para calcular esta
correccion xk^ +1^ es precisamente imp oner que xk^ +1^ sea un cero de f , es decir, f (xk^ +1^ ) = 0.
Teniendo en cuenta la ecuacion 5.10, esto se escrib e como
f (xk^ + xk^ +1^ ) = 0 (5:11)
5 Ceros de funciones 99
Esta claro que para una funcion f arbitraria no es p osible desp ejar xk^ +1^ en la ecuacion 5.11.
Por este motivo, se hace un desarrol lo en serie de Taylor de primer orden de f alrededor de
xk^ , y se obtiene
f (xk^ + xk^ +1^ ) f (xk^ ) + f 0 (xk^ )xk^ +1^ (5:12) Es imp ortante resaltar que la ecuacion 5.12 es una aproximacion y no una igualdad, p orque se han despreciado los terminos del desarrollo de Taylor con derivadas de orden sup erior a uno.
Si ahora se sustituye f (xk^ + xk^ +1^ ) en la ecuacion 5.11 p or la aproximacion obtenida en
5.12, resulta
f (xk^ ) + f 0 (xk^ )xk^ +1^ = 0
de donde puede aislarse xk^ +1^ como
xk^ +1^ =
f (xk^ )
f 0 (xk^ )
siempre que f 0 (xk^ ) 6 = 0. Finalmente, reemplazando esta expresion de xk^ +1^ en la ecuacion 5.10 se llega a la expresion del meto do de Newton:
xk^ +1^ = xk^
f (xk^ )
f 0 (xk^ )
La ecuacion 5.14 prop orciona una estrategia para construir la sucesion de aproximaciones fxk^ g
a un cero x. Como ya se ha comentado, para completar el meto do es necesario elegir una
aproximacion inicial x^0 y un criterio de nalizacion de iteraciones. En cuanto a este ultimo
punto, se emplean los mismos criterios de convergencia (relativo en x y absoluto en f ) que para
el meto do de la biseccion (ver apartado 5.3).
5.4.2 Deduccion gra ca del meto do de Newton
El meto do de Newton puede deducirse tambien de manera gra ca, tal y como se muestra en
la gura 5.5. La idea es la siguiente: dada una cierta xk^ , se aproxima la funcion f p or la recta
tangente a la curva y = f (x) en el punto (xk^ ; f (xk^ )). La p endiente de esta recta es justamente
la derivada de f en xk^. A continuacion se toma la interseccion de esta recta con el eje de las x
como siguiente aproximacion xk^ +1^. En la gura 5.5 se puede observar que f (xk^ ), xk^ +1^ y la
p endiente f 0 (xk^ ) estan relacionados seg un
f (xk^ )
xk^ +^
= f 0 (xk^ )
que es equivalente a la ecuacion 5.13 obtenida en la deduccion analtica del meto do.
Se ha escrito un programa en FORTRAN para calcular races cuadradas a partir de op e- raciones elementales mediante el meto do de Newton (programa 5.3, apartado 5.7). En este
caso, es necesario de nir dos funciones externas FUNCTION: una para la funcion f y otra para
su derivada f 0.
5 Ceros de funciones 101
Las tablas 5.1 y 5.3 p onen de mani esto que la convergencia a
p 2 es mucho mas rapida para el meto do de Newton que para el meto do de la biseccion. Esto es debido a que el meto do
de Newton se basa en una estrategia muy \inteligente": a medida que la aproximacion xk^ se
va acercando al cero x de f , la recta tangente se va pareciendo parece cada vez mas a la curva
y = f (x), hasta confundirse con ella (vease la gura 5.6).
Debido a su rapidez, el meto do de Newton es ampliamente utilizado en la practica. Sin embargo, no es tan robusto como el meto do de la biseccion, tal como se ilustra en el problema 5.2.
Problema 5.2: Se desea calcular
p
2 tomando x^0 = 0 como aproximacion inicial. Veri car que:
a) puede hacerse sin di cultades mediante el meto do de la biseccion (con a = 2,
p or ejemplo). b) el meto do de Newton falla. >Por que? Problema 5.3: a) Escribir un programa en FORTRAN que resuelva el problema del billar cir- cular mediante el meto do de Newton. b) Utilizar el programa para hallar las soluciones con los datos de la gura 5.4.
Fig. 5.6 <El metodo de Newton es rapido!
102 Metodos numericos
5.5 Meto do de la secante
Se ha visto en el apartado anterior que la rapidez del meto do de Newton es debida a la utilizacion de la recta tangente en cada punto. Ello obliga a efectuar dos evaluaciones funcionales en cada iteracion (la funcion y su derivada), en lugar de traba jar unicamen te con valores de la funcion. Esto hace que el coste computacional p or iteracion sea mas elevado para el meto do de
Newton que para el meto do de la biseccion. En los ejemplos expuestos aqu el calculo de f y
f 0 no es excesivamente costoso, p ero es bastante com un en problemas reales que la evaluacion
de f y f 0 requiera resolver un problema complejo cada vez.
Por otro lado, puede o currir que para un determinado problema el calculo de derivadas sea muy farragoso o incluso imp osible (p or ejemplo, si no se disp one de una expresion analtica de
f (x)). En estos casos, se puede utilizar el metodo de la secante. La idea basica de este meto do se
ilustra en la gura 5.7: la recta tangente a la curva y = f (x) en el punto (xk^ ; f (xk^ )) se aproxima
mediante la recta secante que pasa p or este punto y el punto (xk^ ^1 ; f (xk^ ^1 )), obtenido en la
iteracion anterior. En otras palabras, la derivada f 0 (xk^ ) no se calcula, sino que se aproxima p or
f 0 (xk^ )
f (xk^ ) f (xk^ ^1 ) xk^ xk^ ^1
La interseccion de esta recta secante con el eje de las x se toma como siguiente aproximacion
xk^ +1^. Este meto do necesita dos aproximaciones iniciales (x^0 y x^1 ) para p o der trazar la primera
recta secante.
El programa 5.4 (ver apartado 5.7) es un programa FORTRAN que halla races cuadradas mediante el meto do de la secante. Al calcular
p
2 con x^0 = 1, x^1 = 1 : 5 (es decir, las dos
primeras aproximaciones del meto do de Newton, vease la tabla 5.3) y tolx = T O Lf = 0 : 5 10 ^8 , se obtienen los resultados de la tabla 5.4. Comparando las tablas 5.1, 5.3 y 5.4 puede verse que la convergencia del meto do de la secante es ligeramente mas lenta que la del meto do de Newton, p ero claramente mas rapida que la del meto do de la biseccion.
Tabla 5.4 Calculo de
p
2 por el metodo de la secante a partir de x^0 = 1 y x^1 = 1 : 5
Iteracion Aproximacion x f(x) Error relativo en x ========= ============== ========== =================== 0 1.0000000 -1.000D+00 -3.333D- 1 1.5000000 2.500D-01 7.143D- 2 1.4000000 -4.000D-02 -9.756D- 3 1.4137931 -1.189D-03 -2.988D- 4 1.4142157 6.007D-06 1.502D- 5 1.4142136 -8.931D-10 -2.233D-
Convergencia en la iteracion 5 La raiz cuadrada de 2.0000000 es 1.
104 Metodos numericos
biseccion, con un decrecimiento mucho mas lento del error.
1,E-
1,E-
1,E-
1,E-
1,E-
1,E-
1,E+
0 5 10 15 20 25 30
Iteración
Error relativo
Bisección Secante Newton
Fig. 5.8 Comparacion de los tres metodos para el calculo de
p 2
Problema 5.5: Representar en gra cas de convergencia los resultados obtenidos en los problemas 5.1, 5.3 y 5.4, corresp ondientes a la resolucion del problema del billar circular p or los meto dos de biseccion, Newton y secante. Comparar razonadamente las distintas curvas.
Problema 5.6: Escribir un programa FORTRAN para obtener el(los) cero(s) de la funcion f (x) = x^2 2 x + 1 mediante los meto dos de la biseccion, de Newton y de la secante (el usuario deb e p o der elegir, al principio de la ejecucion, que meto do desea emplear). Comentar razonadamente los resultados obtenidos. >Presentan los meto dos su comp ortamiento tpico?
5.7 Asp ectos computacionales: las funciones externas FUNCTION en
FORTRAN
La sentencia FUNCTION p ermite al programador de nir sus propias funciones. Estas funciones se llaman externas para distinguirlas de las funciones intrnsecas, ya incorp oradas en el lengua je FORTRAN.
Supongase que en un programa es necesario calcular el area de muchos triangulos distintos
A a partir de su base y su altura (A = (bh)=2). Para no tener que programar esta formula
5 Ceros de funciones 105
rep etidamente a lo largo del programa principal, se de ne como una funcion externa en un subprograma:
function area(base,altura)
area = (base*altura)/2. return end
Cada vez que sea necesario evaluar el area de un triangulo, desde el programa principal puede utilizarse la funcion externa haciendo:
a = area(b,h)
Tal como se puede observar en el ejemplo:
La funcion externa se declara con la sentencia FUNCTION seguida del nombre de la funcion y de sus argumentos entre parentesis.
Puesto que se trata de un subprograma indep endiente del programa principal, es necesario terminar la funcion externa con las sentencias RETURN (para volver al programa principal) y END (para indicar el nal del subprograma). Amb os programas (el principal y el subpro- grama) pueden estar en un mismo archivo o en archivos distintos.
Los nombres de los argumentos de la funcion pueden ser distintos en la llamada desde el programa principal (b, h) y en el subprograma (base, altura). Lo unico que imp orta es el orden que o cupan en la lista de argumentos, su nombre es indiferente.
Deb e existir concordancia de tip o entre la funcion y el resultado que devuelve. Es decir que, p or ejemplo, utilizando declaraciones p or defecto, si el resultado arro jado p or la funcion es un INTEGER4, el nombre de la funcion deb e emp ezar p or I, J, K, L, M o N, mientras que si es un REAL4 la inicial del nombre de la funcion deb e estar comprendida entre A-H y O-Z. Si el resultado es de otro tip o (INTEGER2, REAL8, REAL16, COMPLEX, LOGICAL o CHARACTER), deb e constar explcitamente en la declaracion de la funcion. As, p or ejemplo, para calcular el area de un triangulo en doble precision se puede utilizar la funcion externa real8 function area(base,altura)
implicit real8 (a-h,o-z) area = (basealtura)/2.d return end
El programa 5.1, empleado en el apartado 5.2 para calcular races de ecuaciones del tip o
f (x) = 0 p or el meto do de la biseccion, utiliza una FUNCTION para programar la funcion f.
Como puede observarse, la funcion externa FUNCTION que eval ua la funcion f es del mismo tip o
(REAL*8) que la variable F a la que esta asignada, y en el interior es necesario volver a declarar
5 Ceros de funciones 107
. 'Aproximacion x', 4x, ' f(x) ', 4x, . 'Error relativo en x') write (6,520) 520 format (1x, '=========', 3x, '==========', . 3x, '==============', 4x, '==========', . 4x, '===================')
c___Control para decidir si se puede iniciar la biseccion c___Tienen f(a) y f(x0) signos opuestos? fa = f(a,s) fx = f(x,s) if( (fafx) .gt. 0.d0 ) then write (6,) ' No se cumple la condicion f(a)*f(x0) < 0' stop endif
c___Inicio del proceso iterativo c___del metodo de la biseccion
do k=0, maxiter
pmedio = (a+x)/2.d fmedio = f(pmedio,s)
c___Calculo del error relativo en x rel_x = (x - pmedio)/pmedio
c___Impresion de resultados write (6,600) k,a,x,fx,rel_x 600 format (3x, i3, 7x, f10.7, 5x, f10.7, 6x, 1pd10.3,
. 8x, 1pd10.3)
c___Control de convergencia if ( (abs(fx).lt.tol_f) .and. (abs(rel_x).lt.tol_x) ) then write (6,700) k 700 format (/,1x, 'Convergencia en la iteracion', i3) write (6,800) s,x 800 format (1x,'La raiz cuadrada de ', f10.7, ' es ', f10.7) stop endif
c___Eleccion del nuevo intervalo segun el valor de f(pmedio) if ( (fmedio*fx) .lt. 0.d0 ) then a = x x = pmedio fx = fmedio
108 Metodos numericos
else x = pmedio fx = fmedio endif
enddo
end
c___________________Function f(x,s) real8 function f(x,s) implicit real8 (a-h,o-z) f = x*x - s return end
Prog. 5.1 Calculo de races cuadradas por el metodo de la biseccion
Como se puede ver, la utilizacion de las sentencias FUNCTION p ermite realizar de manera natural una programacion en modulos (programacion estructurada); este asp ecto se comentara en detalle en el captulo 7. En este sentido, es interesante resaltar que si se desea utilizar el programa 5.1 para hallar el cero de otra funcion distinta, solo es necesario mo di car la FUNCTION corresp ondiente (y la entrada de datos necesarios para este caso). Por ejemplo, el programa 5.2 se ha obtenido mo di cando adecuadamente el programa 5.1 presentado anteriormente para hallar las races de la ecuacion 5.3 (juego del billar en una mesa circular) mediante el meto do de la biseccion.
c c Este programa resuelve el problema del billar circular c mediante el metodo de la BISECCION c c___________________________________________________________________ implicit real8 (a-h,o-z) c___Radio R de la mesa de billar y posicion de las dos bolas P y Q write (6,100) read (5,) r write (6,110) read (5,) xp,yp write (6,120) read (5,) xq,yq
100 format (1x,'Radio de la mesa:') 110 format (1x,'Coordenadas de la bola P:') 120 format (1x,'Coordenadas de la bola Q:')
