• Author / Uploaded
• Anonymous S9FeQUSIsQ

Citation preview

´ UNIVERSIDAD AUTONOMA DE CHIAPAS

Facultad de Ciencias en F´ısica y Matem´ aticas

Apuntes del curso de ´ ´ METODOS NUMERICOS

realizados por el

Dr. Roberto Arceo Reyes

Agosto - Diciembre 2017

Tuxtla Guti´ errez, Chiapas, Mexico

iii

Objetivo El an´ alisis num´erico trata de la obtenci´on, descripci´on y an´alisis de algoritmos para el estudio y soluci´ on de problemas matem´aticos. El desarrollo continuo de las m´aquinas computadoras y su cada vez m´as f´acil accesibilidad, aumenta a igual velocidad como la importancia de los m´etodos num´ericos en la soluci´on de problemas en las ciencias y la ingenier´ıa. Gran parte de los egresados de una licenciatura en Ciencias F´ısico-Matem´aticas tratan con problemas que requieren del uso de estos m´etodos. El presente curso pretende dar un panorama amplio de la gama de problemas matem´aticos que se pueden resolver usando los m´etodos num´ericos y se obtienen los algoritmos correspondientes. El curso proporciona material que debe ser conocido por todo cient´ıfico o ingeniero. El curso requiere de los conocimientos de la l´ınea de c´alculo, el manejo de un lenguaje de programaci´ on y el conocimiento de conceptos de ´algebra lineal como espacios vectoriales, normas y matrices. Al t´ermino del curso el estudiante dominar´a los puntos esenciales del mismo, a saber: estudio y clasificaci´on de errores, soluci´on de ecuaciones no lineales, aproximaci´ on e interpolaci´ on, derivaci´on e integraci´on num´erica, soluci´on num´erica de ecuaciones diferenciales ordinarias y soluci´on de sistemas de ecuaciones lineales y habr´ a realizado programas de computadora de los principales m´etodos estudiados.

Prop´ osito Conocer´ a la importancia del an´alisis num´erico por su aplicaci´on a problemas cl´asicos: soluci´ on de ecuaciones, aproximaci´on de funciones, soluci´on num´erica de ecuaciones diferenciales ordinarias y soluci´ on de sistemas de ecuaciones lineales; y los conceptos de algoritmo, convergencia y estabilidad de un algoritmo. Calcular los errores absoluto y relativo de una aproximaci´ on; · Apreciar la utilidad del valor relativo; · Entender los conceptos de cifras significativas; programaci´on del error en operaciones aritm´eticas; programaci´ on del error en la evaluaci´on de funciones y condicionamiento de

iv

un algoritmo. Conocer y diferenciar los m´etodos de bisecci´on, secante, falsa posici´on e iteraci´on de punto fijo; · Comprender los conceptos de aceleraci´on de la convergencia; convergencia lineal y cuadr´ atica; · Aplicar los m´etodos estudiados al c´alculo de ra´ıces de polinomios. · Aplicar el m´etodo de Newton a sistemas de ecuaciones no lineales. Entender los conceptos b´ asicos de aproximaci´on de funciones y ser´a capaz de aproximarlas usando los diferentes criterios estudiados: a saber, m´ınimos cuadrados por polinomios; polinomios ortogonales; funciones splines; e interpolaci´on polinomial. Conocer y aplicar los m´etodos de derivaci´on y sus errores de redondeo y truncamiento; y los de integraci´ on, sus f´ ormulas de error y la aplicaci´on de ´estas para la integraci´ on compuesta. Dominar los m´etodos de Euler, Taylor; Runge-Kutta, los basados en integraci´on num´erica, Adams-Bashfort para solucionar ecuaciones diferenciales ordinarias; · Comprender los conceptos de estabilidad num´erica. Conocer y ser capaz de aplicar los m´etodos directos para la soluci´on de sistemas de ecuaciones lineales; las estrategias de pivoteo; y los m´etodos iterativos. As´ı tambi´en los m´etodos correspondientes para el c´alculo de valores propios.

TABLA DE CONTENIDOS LISTA DE TABLAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LISTA DE FIGURAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. INTRODUCCION

ix 1

1.1.

Problemas cl´ asicos del Anal´ısis Num´erico . . . . . . . . . . . . . . . . . . .

1.2.

Descripci´ on de un Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2. Estudio general del error en un proceso num´ erico

1

13

2.0.1. Exactitud y precisi´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.

Reglas de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3. Resoluci´ on de ecuaciones no lineales

21

3.1.

La serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2.

El m´etodo de bisecci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.3.

El m´etodo de la regla falsa . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.4.

Iteracci´ on de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.5.

El m´etodo de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.6.

El m´etodo de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.7.

Ra´ıces m´ ultiples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4. Aproximaci´ on e interpolaci´ on 4.1.

41

M´ nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 v

vi

5.

4.2.

Interpolaci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.3.

F´ ormulas de integraci´on de Newton-Cotes . . . . . . . . . . . . . . . . . . . 50

4.4.

Regla del trapecio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.5.

Regla de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Soluci´ on num´ erica de ecuaciones diferenciales ordinarias: problema a valores iniciales

55

5.1.

La serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.2.

M´etodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.3.

M´etodo de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.4.

Eliminaci´ on Gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

6. Ejercicios resultos BIBLIOGRAFIA

71 xi

Lista de Tablas

Tabla 1.1 Cinco pasos necesarios para producir y dar soporte a programas de alta ....... 12 Tabla 2.1: Ilustraci´ on de los d´ıgitos retenidos y descartados de un n´ umero con cinco .... 18

vii

Lista de Figuras

Figure 1.1 Representaci´ on de la fuerza que actua sobre un ........ 2 Figure 1.2 Soluci´ on anal´ıtica al problema del paracaidista ......... 6 Figure 1.3 Uso de una diferencia finita para aproximar la .......... 7 Figure 3.1 Esquema gr´ afico del m´etodo de la regla falsa ............ 25 Figure 3.2 Esquema gr´ afico del m´etodo Newton-Raphson .......... 28 Figure 3.3 Esquema gr´ afico del m´etodo de la secante ................. 31 Figure 4.1 El residuo en la regresi´on lineal representa ................ 36 Figure 4.2 Las ´ areas sombreadas muestran tri´angulo .................. 37 Figure 4.3 Diferencia entre f´ ormulas de integraci´on .................... 41 Figure 4.4 Esquema gr´ afico de la regla trapezoidal ..................... 43 Figure 5.1 Esquema gr´ afico del m´etodo de un paso ..................... 46 Figure 5.2 M´etodo de Euler (o m´etodo de Euler-Cauchy.............. 49 Figure 5.3 Soluci´ on del m´etodo de Euler ...................................... 51

ix

Bibliograf´ıa [1] Atkinson, K. (1989). An Introduction to Numerical Analysis. Wiley, 2nd edition. [2] Hamming, R.W. (1987). Numerical Methods for Scientists and Engineers. Dover Publications, 2nd edition. [3] Burden, R.L., Faires, J. D. (2002). Anlisis Numrico. Mxico: International Thomson Editores, S.A. De C.V. [4] Henrici, P. (1980). Elementos de Anlisis Numrico. Mxico: Editorial Trillas. [5] Carnahan, B., Luther, H. A., Wilkes, J.O. (1969). Applied Numerical Methods. Wiley, 1st edition. [6] Kincaid, D.R., Cheney, E.W. (2001). Numerical Analysis: Mathematics of Scientific Computing. Brooks Cole, 3rd edition. [7] Gibbs, W.R. (2001). Computation in Modern Physics. USA: World Scientific Publishing Company, 2nd editition.

xi

Cap´ıtulo 1

Introducci´ on El presente trabajo es un compilaci´on de los apuntes dados en el curso de M´etodos N´ umericos impartido en el semestre Agosto - Diciembre del 2017 en la Licenciatura en F´ısica de la Facultad de Ciencias en F´ısica y Matem´aticas de la Universidad Aut´onoma de Chiapas.

1.1.

Problemas cl´ asicos del Anal´ısis Num´ erico

Recordando la segunda ley de Newton,

F = ma

(1.1)

donde F es la fuerza neta que act´ ua sobre el cuerpo, m es la masa del objeto, y a es su aceleraci´ on.

La ecuaci´ on (1.1) tiene varias caracter´ısticas habituales de los modelos matem´aticos del mundo f´ısico. 1. Describe un sistema o proceso natural en t´erminos matem´ aticos. 2. Representa una idealizaci´on y una simplificaci´on de la realidad. 1

2

3. Conduce a resultados predecibles, y en consecuencia puede emplearse para prop´ositos de predicci´ on. La aceleraci´ on de la ecuaci´ on (1.1) puede despejarse,

a = F/m

(1.2)

esta ecuaci´ on puede usarse para calcular la velocidad final de un cuerpo en ca´ıda libre cerca de la superficie terrestre. El cuerpo en descenso ser´ a un paracaidista.

Figura 1.1: Representaci´ on de las fuerzas que act´ uan sobre un paraca´ıdista en descenso. FD es la fuerza hacia abajo debido a la atracci´on de la gravedad. FU es la fuerza hacia arriba debido a la resistencia del aire. Podemos crear un modelo matem´atico al expresar la aceleraci´on como la raz´on de cambio on (1.1) para dar de la velocidad con respecto al tiempo ( dv dt ) y sustituir en la ecuaci´

m· donde v es la velocidad.

dv =F dt

(1.3)

3

Para un cuerpo que cae dentro del per´ımetro de la Tierra, la fuerza total est´a compuesta por dos fuerzas contrarias. La atracci´on hacia abajo debida a la gravedad FD y la fuerza hacia arriba debida a la resistencia del aire FU .

F = FD + FU

(1.4)

Si a la fuerza hacia abajo se les asiga un signo positivo, se puede usar la segunda ley para formular la fuerza debida a la gravedad como

FD = mg

(1.5)

donde g es la constante de gravitac´on, o la aceleraci´on debida a la gravedad, con un valor de 9.8 m/s2 (mks) o 980 cm/s2 (cgs). La resistencia del aire puede formularse de varias maneras, una es suponer que es linealmente proporcional a la velocidad,

FU = −cv

(1.6)

donde c es una constante llamada el coeficiente de arrastre (en gramos por segundo o kg/s). Combinando las ecuaciones (1.3) a (1.6),

m

dv = mg − cv dt

dv c =g− v dt m

dv = dt, c (g − m v)

(1.7)

(1.8)

4

haciendo u=g−

du = −

dv = − m c du

c v m

c dv m

y

−

m du c u

= dt

c du = − dt u m Integrando !

du c =− u m

!

dt

(en reposo v = 0 en t = 0)

ln(u)|v0 = −

ln(g −

ln(g −

ln(

c ""t t" m 0

c c ""v v)" = − t m 0 m

c c v) − ln(g) = − t m m

c v g−m c )=− t g m

ln(1 −

c c v) = − t mg m

5

(1 −

c c v) = e− m t mg

c c v = (1 − e( m )t ) mg

−→ v(t) =

c mg [1 − e− m t ](1.9) c

Ejemplo: Soluci´ on anal´ıtica al problema del paracaidista que cae. Un paracaidista con una masa de 68, 100 gramos salta de un aeroplano. Apl´ıquese la ecuaci´ on (1.9) para calcular la velocidad antes de abrir el paracaidas. El coeficiente de arrastre c es aproximadamente igual a 12,500 g/s.

Soluci´ on: Al sustituir los valores en la ecuaci´on (1.9) se tiene,

v(t) =

(980cm/s2 )(68, 100g) 12500g/s [1 − exp(− t)] 12, 500g/s 68, 100g

v(t) = (5339.04cm/s)[1 − e−0.18355t/s ] Al dar valores a t se obtiene, t(s)

v(cm/s)

0

0

2

1640.5

4

2776.9

6

3564.2

8

4109.5

6

10

4487.3

12

4749.0

∞

5339.0

Figura 1.2: Soluci´on anal´ıtica al problema del paracaidista.

A la ecuaci´ on (1.9) se le llama una soluci´on anal´ıtica exacta. Si no es posible solucionar la ecuaci´ on se tiene que usar una soluci´on num´erica que se aproxime a la soluci´on exacta, y emplearemos un m´etodo num´erico. La segunda ley de Newton en la aceleraci´on se puede aproximar,

dv ∆v v(ti+1 ) − v(ti ) ≈ = dt ∆t ti+1 − ti

(1.10)

donde ∆v y ∆t son diferencias en la velocidad y el tiempo calculadas sobre intervalos

7

Figura 1.3: Uso de una diferencia finita para aproximar la primera derivada de V con respecto a t.

finitos, v(ti ) es la velocidad en el tiempo inicial ti , y v(ti+1 ) es la velocidad en un tiempo m´as tarde ti+1 . La ecuaci´ on (1.10) es una diferencia finita dividida en el tiempo ti . Al sustituir en la ecuaci´ on (1.8) se da,

c v(ti+1 ) − v(ti ) = g − v(ti ) ti+1 − ti m

(1.11)

al ordenar esta ecuaci´ on d´ a,

v(ti+1 ) = v(ti ) + [g −

c v(ti )](ti+1 − ti ) m

(1.12)

(nuevo valor de v)= (valor anterior de v)+ (valor estimado de la pendendiente) x (incremento del tiempo)

8

Ejemplo 2: Soluci´ on num´erica al problema del paracaidista que cae. Al hacer lo mismo pero con un incremento de 2 s. En t1 = 0 la v(ti )=0 y se puede estimar v(ti+1 ) en ti+1 =2 s.

i = 1 : t1 = 0

y

v(t2 ) = v(t1 ) + [980 −

v(t1 ) = 0

12500 v(t1 )(t2 − t1 ) 68100

t2 = 2s

v(t2 ) = (980)(2 − 0) = 1960

cm s

i=2:

v(t3 ) = v(t2 ) + [980 −

12500 v(t2 )](t3 − t2 ) 65100

t3 = 4s

v(4) = 1960 + [980 −

12500 (1960)](4 − 2) 68100

v(4) = 3200.5

i=3:

cm s

9

v(t4 ) = V (t3 ) + [980 −

12500 v(t3 )](t4 − t3 ) 68100

t4 = 6s

v(6) = 3200.5 + [980 −

12500 3200.5](6 − 4) 68100

v(6) = 3985.6

cm s

v(8) = 4482.5

v(10) = 4796.9

v(12) = 4995.9

∞ = 5339.0 t(s)

v(cm/s)

0

0

2

1960.0

4

3200.5

6

3985.6

8

4482.5

10

4796.9

10

12 ∞

4995.9 5339.0

1.2.

Descripci´ on de un Algoritmo

1.2.1.

Dise˜ no de algoritmos

Un algoritmo es una secuencia l´ogica de pasos necesarios para ejecutar una tarea espec´ıfica tal como la soluci´ on de un problema. Los buenos algoritmos tienen ciertas caracter´ısticas. Siempre deben terminar despu´es de una cantidad finita de pasos y deben ser lo m´ as general posible para tratar cualquier caso particular.

1.2.2.

Composici´ on de un programa

Despu´es de confeccionar un algoritmo, el paso siguiente es expresarlo como una secuenca de declaraciones de programaci´on llamado c´odigo. Lo realizaremos en FORTRAN.

Programa de computadora en F ORT RAN para el problema de la suma simple. c

Version en FORTRAN A=25 B=15 C=A+B WRITE(6,1) C

1

FORMAT(’ ’,F10.3) STOP END ↑ ↑

↑

Columna 7 Columna 5

Columna

1

11

1.2.3.

FORTRAN

Es la contrucci´ on de f´ ormula translation (traducci´on de f´ormulas), y se desarroll´o en la dec´ada de 1950. Debido a que fue expresamente dise˜ nado para c´alculos, ha sido el lenguaje m´as usado en la ingenier´ıa y la ciencia.

1.2.4.

Rastreo y prueba

Despu´es de escribir el c´ odigo del programa, se debe probar para buscar los errores, a los que se les llama bugs. Al proceso de localizar y corregir los errores se les conoce como rastreo.

1.2.5.

Documentaci´ on

Despu´es del rastreo y probado del programa, ´este se debe documentar. La documentaci´ on es la inclusi´ on de comentarios que le permiten al usario implementar el programa m´as f´ acilmente.

1.2.6.

Almacenamiento y mantenimiento

Los pasos finales del programa son el almacenamiento y mantenimiento del mismo. El mantenimiento involucra acondicionar el programa e incluso hacerle cambios que lo hagan accesible a problemas reales. El almacenamiento se refiere a la manera en que los programas se guardan para uso posterior.

Compilaci´ on en fortran: gf ortran

−o

nombre : del : ejecutable

codigo : f uente.f or

12

Tabla 1.1: Cinco pasos necesarios para producir y dar soporte a programas de alta calidad.

Cap´ıtulo 2

Estudio general del error en un proceso num´ erico 2.0.1.

Exactitud y precisi´ on

Los errores asociados con los c´alculos y medidas se pueden caracterizar observando su precisi´ on y exactitud. La precisi´on se refiere a:

1. El n´ umero de cifras significativas en una cantidad. La extensi´on en las lecturas repetidas. 2. La exactitud se refiere a la aproximaci´on de un n´ umero o de una medida al valor verdadero que se supone representa.

Los errores n´ umericos se originan con el uso de aproximaciones. Estos incluyen errores de truncamiento y los errores de redondeo, el primer caso se generan por representar aproximadamente un procedimiento matem´atico exacto y el segundo por representar n´ umeros exactos.

Valor verdadero = Valor aproximado + error 13

14

Error num´erico igual a la diferencia del valor verdadero y el valor aproximado:

Ev = valor verdadero - valor aproximado

Ev : valor exacto del error

Una manera de medir las magnitudes de las cantidades que se est´an evaluando es normalizar el error respecto al valor verdadero, como en:

Error : relativo : f uncional =

error valor : verdadero

El error relativo tambi´en se puede multiplicar por el 100 % para expresarlo como:

!v =

error : verdadero × 100 % valor : verdadero

donde !v denota el error relativo porcentual (real).

Normalizar el error es una alternativa, usando la mejor estimaci´on posible del valor verdadero, esto es, a la aproximaci´on misma, como:

!a =

error : aproximado × 100 % valor : aproximado

donde el sub´ındice a significa que el error est´a normalizado a un valor aproximado.

15

!a : error relativo porcentual aproximado.

Ciertos m´etodos num´ericos usan un esquema iterativo para calcular resultados. Se hace una aproximaci´ on anterior, esto se repite varias veces, en tales casos, el error a menudo se calcula como la diferencia entre la aproximaci´on previa y la actual. Por lo tanto, el error relativo porcentual est´ a dado por:

!a =

aproximaci´ on : actual − aproximaci´ on : previa × 100 % aproximaci´ on : actual

a menudo no importa el signo del error, si no que su valor absoluto sea menor a una tolerancia prefijada !s . | !a |< !s !s =(0.5×102−n ) % si se cumple este u ´ltimo criterio el resultado es correcto en al menos n cifras significativas.

Ejemplo: Estimaci´ on del error para m´ etodos iterativos A menudo se pueden representar las funciones mediante una serie infinita. La funci´ on exponencial se puede representar como,

ex = 1 + x +

x2 x3 x4 + + + ... 2 3 4

(Expansi´ on en series de Maclaurin) Empezando con el primer t´ermino, ex = 1, y agregando m´as t´erminos, est´ımese el valor de e0.5 . Despu´es que se agregue cada t´ermino, calc´ ulesen los errores relativos porcentuales real y aproximado. El valor real de e0.5 = 1.648721271. Agregar t´erminos hasta que el valor absoluto del error aproximado ea sea menor al criterio preestablecido, es que contempla tres cifras significativas.

Soluci´ on:

16

!s =(0.5×102−3 ) % = 0.005 % (primer estimaci´ on)

ex ≈ 1 (segunda estimaci´ on)

ex ≈ 1 + x para x = 0.5

e0.5 ≈ 1 + 0.5 = 1.5 error relativo porcentual real,

!v =

1.648721271 − 1.5 × 100 % = 9.02 % 1.648721271

!a =

1.5 − 1 × 100 % = 33.3 % 1.5

!a no es menor al valor prefijado, !s ; as´ı que se contin´ ua agregando m´as t´erminos y calculando los errores (!a < !s ). (tercera estimaci´ on)

ex ≈ 1 + 0.5 + (

!v =

0.512 = 1.625 2

1.648721271 − 1.6525 × 100 % = 1.438 % 1.648721271

17

!a =

2.1.

1.625 − 1.5 × 100 % = 7.69 % 1.625

Reglas de redondeo

1. En el redondeo, se conservan las cifras significativas y el resto se descarta (ver figura). El u ´ltimo d´ıgito que se conserva se aumenta en uno si el primer d´ıgito descartado es mayor de 5. De otra manera se deja igual. Si el primer d´ıgito descartado es 5 o 5 seguido de ceros, entonces el u ´ltimo d´ıgito retenido se incrementa en 1, s´olo si es impar. 2. En la suma y en las restas, el redondeo se lleva a cabo de forma tal que el u ´ltimo d´ıgito retenido en la respuesta corresponda al u ´ltimo d´ıgito m´as significativo de los n´ umeros que est´ an sumando o restando. N´otese que un d´ıgito en la columna de las cent´esimas es m´ as significativo que uno de la columna de las mil´esimas. 3. Para la multiplicaci´ on y para la divisi´on el redondedo es tal que la cantidad de cifras significativas del resultado es igual al n´ umero m´as peque˜ no de cifras significativas que contiene la cantidad en la operaci´on. 4. Para combinaciones de las operaciones aritm´eticas, existen dos casos generales. Se puede sumar o restar el resultado de las multiplicaciones o de las divisiones.

(Multiplicaci´ on o divisi´on)±(Multiplicaci´on o divisi´on)

o tambi´en se pueden multiplicar o dividir los resultados de las sumas y las restas:

(Suma o resta)×÷(Suma o resta)

18

En ambos casos, se ejecutan las operaciones entre par´entesis y el resultado antes de proceder con otra operaci´on, en vez de redondear u ´nicamente el resultado final.

Tabla 2.1: Ilustraci´ on de los d´ıgitos retenidos y descartados de un n´ umero con cinco cifras significativas.

19

Ejemplo: Ilustraciones de las reglas de redondeo 1. Errores de redondeo 5.6723 −→ 5.67

3 cifras significativas

10.406 −→ 10.41

4 cifras significativas

7.3500 −→ 7.4

2 cifras significativas

88.21650 −→ 88.216

5 cifras significativas

1.25001 −→ 1.3

2 cifras significativas

2. Sumas y restas 2.2 − 1.768 = 0.432 −→ 0.4 3. Multiplicaci´ on y divisi´ on a)

0.0642 × 4.8 = 0.30816 −→ 0.4

b)

945 0.3185

= 2967.032967 −→ 2970

4. Combinaciones

Cap´ıtulo 3

Resoluci´ on de ecuaciones no lineales 3.1.

La serie de Taylor

La serie de Taylor da una formulaci´on para predecir el valor de la funci´on en xi+1 en t´erminos de la funci´ on y de sus derivadas en una vecindad al punto xi . Se construir´ a t´ermino a t´ermino, el primer t´ermino de la serie es:

f (xi+1 ) ≈ f (xi )

(aproximaci´ on de

orden cero)

(3.1)

La ecuaci´ on (3.1) da una buena aproximaci´on si la funci´on que se va a aproximar es una constante. La aproximaci´ on a primer orden se obtiene sumando otro t´ermino al anterior para obtener:

f (xi+1 ) ≈ f (xi ) + f " (xi )(xi+1 − xi )

(3.2)

f " (xi ) es la pendiente multiplicada por la distancia entre xi+1 y xi .

Agregando a la serie un t´ermino de 2do orden para obtener algo sobre la curvatura de la 21

22

funci´ on si es que la tiene:

f (xi+1 ) ≈ f (xi ) + f " (xi )(xi+1 − xi ) +

f "" (xi ) (xi+1 − xi )2 2!

(3.3)

agregando m´ as t´erminos para desarrollar la expansi´on completa de la serie de Taylor:

f (xi+1 ) ≈ f (xi )+f " (xi )(xi+1 −xi )+

f "" (xi ) f """ (xi ) f n (xi ) (xi+1 −xi )2 + (xi+1 −xi )3 +...+ (xi+1 −xi )n +Rn 2! 3! n! (3.4)

se incluye un t´ermino residual para considerar todos los t´erminos desde (n + 1) hasta el infinito: f (n+1) (ξ) (xi+1 − xi )n+1 (n + 1)!

Rn =

(3.5)

donde el sub´ındice n indica que el residuo es de la aproximaci´on a n-´esimo orden y ξ es un valor cualquiera de x que se encuentra en xi y xi+1 .

Rn da una estimaci´ on exacta del error.

h = xi+1 − xi y la ecuaci´ on (3.4) toma la forma:

f (xi+1 ) = f (xi ) + f " (xi )h +

f "" (xi ) 2 f """ (xi ) 3 f (n) (xi ) n h + h + ... + h + Rn 2! 3! n!

en donde

Rn =

f n+1 (ξ) n+1 h (n + 1)!

23

Ejemplo: Usense t´erminos en la serie de Taylor de cero a cuarto orden para aproximar la funci´ on:

f (x) = −0.1x4 − 0.15x3 − 0.5x2 − 0.25x + 1.2 desde el punto xi = 0 y con h = 1. Esto es, predecir el valor de la funci´on en xi+1 = 1.

3.2.

El m´ etodo de bisecci´ on

Si una funci´ on f (x) es real y continua en el intervalo de xl a xu y f (xl ) y f (xu ) tienen signos opuestos, esto es,

f (xl )f (xu ) < 0 entonces hay, al menos una ra´ız real entre xl y xu . El m´etodo de bisecci´ on, conocido tambi´en como de corte binario, de partici´on en dos intervalos iguales o m´etodo de Bolzano, es un m´etodo de b´ usqueda incremental donde el intervalo se divide siempre en dos. Si la funci´on cambia de signo sobre un intervalo, se eval´ ua el valor de la funci´ on en el punto medio. La posici´on de la ra´ız se determina sit´ uandola en el punto medio del subintervalo dentro del cual ocurre un cambio de signo. El proceso se repite hasta obtener una mejor aproximaci´on.

Pasos: 1. Esc´ ojanse los valores iniciales xl y xu de forma tal que la funci´on cambie de signo sobre el intervalo. Esto se puede verificar asegur´andose de que f (xl )f (xu ) < 0. 2. La primera aproximaci´ on a la ra´ız xl se determina como: xr =

xl +xu 2

3. Real´ıcense las siguientes evaluaciones y determ´ınense en que subintervalo cae la ra´ız: a)

Si f (xl )f (xr ) < 0, entonces la ra´ız se encuentra dentro del primer subintervalo.

24

Por lo tanto, resu´elvanse xu = xr y contin´ uese en el paso 4. b)

Si f (xl )f (xr ) > 0, entonces la ra´ız se encuentra dentro del segundo subintervalo.

Por lo tanto, resu´elvanse xl = xr y contin´ uese en el paso 4. c)

Si f (xl )f (xr ) = 0, entonces la ra´ız es igual a xr y se terminan los c´alculos.

4. Calc´ ulese una nueva aproximaci´on a la ra´ız mediante:

xr =

xl + xu 2

5. Dec´ıdase si la nueva aproximaci´on es tan exacta como se desee. Si es as´ı, entonces los c´ alculos terminan, de otra manera, regr´esese al paso 3. Ejemplo (bisecci´ on): Usese el m´etodo de la bisecci´ on para determinar la ra´ız de f (x) = e−x − x.

xl = 0

xr =

y

xu = 1

xl + xu 0+1 1 = = = 0.5 2 2 2

y la cuarta iteraci´ on es:

f (0.5)f (0.625) = −0.010 < 0 la ra´ız est´ a entre 0.5 y 0.625: xu = 0.625

xr =

0.5 + 0.625 = 0.5625 y 2

| !v |= 0.819 %

El m´etodo se puede repetir para obtener mejores estimaciones.

25

El error relativo aproximado !a

| !a |=|

xnueva − xanterior r r | ×100 % xnueva r

xnueva : es la ra´ız de la iteraci´on actual. r

xanterior : es el valor de la ra´ız de la iteraci´on anterior. r

cuando | !a | es menor que un valor previamente fijado, que define el criterio de paro, !s , el programa se detiene.

Programaci´ on del m´ etodo de bisecci´ on (biseccion.f or). implicit real*8 (a-h,o-z) F(X)=EXP(-X)-X OPEN(5,FILE=’biseccion.dat’) READ(5,*)XL,XU,ES,IM 1 FORMAT(3F10.0,I5)

AA=F(XL)*F(XU)

IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM

AA=F(XL)+f(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR

26

XN=(XL+XU)/2 IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 230 XR=XN 240 CONTINUE OPEN(6,FILE=’biseccion.out’) WRITE(6,2) 2 FORMAT(’ ’,’NO SE ENCONTRO LA RAIZ’) WRITE(6,3)XR,EA 3 FORMAT(’ ’,2F10.5) GOTO 310 280 WRITE(6,4)XN,EA,NI 4 FORMAT(’ ’,2F10.5,I5) GOTO 310 300 WRITE(6,5)XR 5 FORMAT(’ ’,’LA RAIZ EXACTA ES =’,F10.5) 310 STOP END

27

3.3.

M´ etodo de la regla falsa

Un defecto del m´etodo de bisecci´on es que al dividir el intervalo xl a xu en mitades iguales, no se toma en consideraci´on la magnitud de f (xl ) y de f (xu ). Por ejemplo, si f (xl ) est´a mucho m´ as cerca de cero que f (xu ), es l´ogico que la ra´ız se encuentra m´as cerca de f (xl ) que de f (xu ). Este m´etodo alternativo aprovecha la idea de unir los puntos con una l´ınea recta. La intersecci´ on de esta l´ınea con el eje x proporciona una mejor estimaci´on de la ra´ız. El reemplazamiento de la curva por una l´ınea recta da una posici´ on f alsa de la ra´ız, de aqu´ı el nombre de m´etodo de la regla falsa o en lat´ın, Regla falsi. Tamb´ıen se le conoce como m´etodo de interpolaci´on lineal. Con el uso de tri´ angulos semejantes, la intersecci´on de la l´ınea recta y el eje x se puede calcular de la siguiente manera: f (xl ) f (xu ) = xr − xl xr − xu Resolvi´endose para xr xr = xu −

f (xu )(xl − xu ) f (xl ) − f (xu )

Est´a es la f´ ormula de la regla falsa. Programaci´ on del m´ etodo de la regla falsa (rf alsa.f or). implicit real*8 (a-h,o-z) F(X)=EXP(-X)-X OPEN(5,FILE=’rfalsa.dat’) READ(5,*)XL,XU,ES,IM 1 FORMAT(3F10.0,I5)

AA=F(XL)*F(XU)

IF(AA.GE.0.0) GOTO 310

28

Figura 3.1: Esquema gr´ afico del m´etodo de la regla falsa. La f´ormula se deriva de los tri´angulos semejantes (´ areas sombreadas).

XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU))

29

IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 230 XR=XN 240 CONTINUE OPEN(6,FILE=’rfalsa.out’) WRITE(6,2) 2 FORMAT(’ ’,’NO SE ENCONTRO LA RAIZ’) WRITE(6,3)XR,EA 3 FORMAT(’ ’,2F10.5) GOTO 310 280 WRITE(6,4)XN,EA,NI 4 FORMAT(’ ’,2F10.5,I5) GOTO 310 300 WRITE(6,5)XR 5 FORMAT(’ ’,’LA RAIZ EXACTA ES =’,F10.5) 310 STOP END

30

Los m´etodos abiertos se basan en f´ormulas que requieren de un solo valor de x o un par de ellos pero no necesariamente encierran la ra´ız, algunas veces divergen o se alejan de la ra´ız a medida que crece el n´ umero de iteraciones.

3.4.

Iteraci´ on de punto fijo

Se emplea una f´ ormula que se puede desarrollar para la iteraci´on de punto fijo o real usando la ecuaci´ on f (x) = 0 de tal forma que x quede del lado requerido de la ecuaci´ on:

x = g(x)

Dada una aproximaci´ on inicial de la ra´ız, xi , se puede obtener una nueva aproximaci´ on xi+1 , expresada por la f´ ormula iterativa:

xi+1 = g(xi )

Se puede estimar el error: | !a |=|

xi+1 − xi | ×100 % xi+1

Ejemplo: Localizar la ra´ız de f (x) = e−x − x. Programaci´ on del m´ etodo de Iteraci´ on de punto fijo (iteracion.f or). implicit real*8(a-h,o-z) F(X)=EXP(-X)-X OPEN(5,FILE=’iteracion.dat’) READ(5,*) XR,ES,IM PRINT *, XR,ES,IM 1 FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=F(XR)

31

IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 170 XR=XN 180 CONTINUE OPEN(6,FILE=’iteracion.out’) WRITE(6,2) 2 FORMAT(’ ’,’NO SE ENCONTRO LA RAIZ’) NI=NI-1 210 WRITE(6,3) XN,EA,NI 3 FORMAT(’ ’,2F15.5,I5) STOP END

32

3.5.

M´ etodo de Newton-Raphson

La f´ ormula de Newton-Raphson es una de las m´as ampliamente usadas. Si el valor inicial de la ra´ız es xi , entonces se puede extender una tangente desde el punto [xi , f (xi )]. El punto donde esta tangente cruza al eje x representa una aproximaci´on mejorada a la ra´ız. La primera derivada en x es equivalente a la pendiente:

f " (xi ) =

f (xi ) − 0 xi − xi+1

ordenando se obtiene:

xi+1 = xi −

f (xi ) f " (xi )

Conocida como la f´ ormula de Newton-Raphson.

3.5.1.

Criterios de paro y estimaci´ on de errores

La derivaci´ on del m´etodo con la serie de Taylor proporciona un conocimiento te´ orico relacionado con la velocidad de convergencia expresado como: Ei+1 = 0(Ei )2 . De esta forma, el error debe ser casi proporcional al cuadrado del error anterior.

La serie de Taylor se puede representar como:

f (xi+1 ) = f (xi ) + f " (xi )(xi+1 − xi ) + ξ ! [xi , xi+1 ]

truncado a la primera derivada

f "" (ξ) (xi+1 − xi )2 2

33

Figura 3.2: Esquema gr´afico del m´etodo de Newton-Raphson.

f (xi+1 ) ≈ f (xi ) + f " (xi )(xi+1 − xi ) en x intersectado f (xi+1 ) = 0,

y

0 ≈ f (xi ) + f " (xi )(xi+1 − xi )

xi+1 = xi − id´entica a la ecuaci´ on de Newton-Raphson.

f (xi ) f " (xi )

(∗)

34

xi+1 = xr en donde xr es el valor exacto de la ra´ız f (xr ) = 0 y

f (xr ) = f (xi ) + f " (xi )(xr − xi ) +

f "" (ξ) (xτ − xi )2 2

(∗∗)

restando (*) y (**), se obtiene:

0 = f " (xi )(xr − xi+1 ) +

f "" (ξ) (xr − xi )2 2

(∗ ∗ ∗)

El error es igual a la diferencia entre xi+1 y el valor real, xr como en:

Ev,i+1 = xr − xi+1

y la ecuaci´ on (***) se expresa como:

0 = f " (xi )Ev.i+1 +

f "" (ξ) 2 Ev,i 2

y si hay convergencia, xi y ξ se aproximan a la ra´ız xr y podemos expresar:

Ev,i+1 ≈ −

f "" (xr ) 2 E 2f " (xr ) v,i

Programaci´ on del m´ etodo Newton-Raphson (nraphson.f or). implicit real*8(a-h,o-z) F(X)=EXP(-X)-X FD(X)=-EXP(-X)-1 OPEN(5,FILE=’nraphson.dat’) READ(5,*) XR,ES,IM 1 FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR)

35

IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 170 XR=XN 180 CONTINUE OPEN(6,FILE=’nraphson.out’) WRITE(6,2) 2 FORMAT(’ ’,’NO SE ENCONTRO LA RAIZ’) NI=NI-1 210 WRITE(6,3) XN,EA,NI 3 FORMAT(’ ’,2F12.5,I5) STOP END

36

3.6.

M´ etodo de la secante

La derivada se puede aproximar mediante una diferencia dividida, como

f " (xi ) ≈

f (xi−1 ) − f (xi ) xi−1 − xi

Se puede sustituir en la ecuaci´on de Newton-Raphson obteniendo la ecuaci´on iterativa:

xi+1 = xi −

xi+1 = xi −

f (xi ) f (xi−1 )−f (xi ) xi−1 −xi

f (xi ) f " (xi )

= xi −

f (xi )(xi−1 − xi ) f (xi−1 ) − f (xi )

F´ormula para el m´etodo de la secante:

xi+1 = xi −

f (xi )(xi−1 − xi ) f (xi−1 ) − f (xi )

No requiere que f (x) cambie de signo entre estos valores, a este m´etodo no se le clasifica como aquellos que usan intervalos. ´ Ejemplo: Usese el m´etodo de la secante para calcular la ra´ız de f (x) = e−x −x. Empi´ecese con los valores iniciales de x−1 = 0 y x0 = 1.0. Programaci´ on del m´ etodo de la secante (secante.f or). implicit real*8(a-h,o-z) F(X)=EXP(-X)-X OPEN(5,FILE=’secante.dat’) READ(5,*) XR,XI,IM 1 FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170

37

Figura 3.3: Esquema gr´afico del m´etodo de la secante. EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 170 XR=XN 180 CONTINUE OPEN(6,FILE=’secante.out’) WRITE(6,2) 2 FORMAT(’ ’,’NO SE ENCONTRO LA RAIZ’) NI=NI-1 210 WRITE(6,3) XN,F(XR),F(XI) 3 FORMAT(’ ’,3F12.5) STOP END

38

3.7.

Ra´ıces M´ ultiples

Una ra´ız multiple corresponde a un punto donde una funci´on es tangencial al eje x. Por ejemplo, dos ra´ıces repetidas resultan de:

f (x) = (x − 3)(x − 1)(x − 1)

f (x) = x3 − 5x2 + 7x − 3

(esta

(multiplicando t´ erminos)

f unci´ on tiene

una

ra´ız

doble)

La ra´ıces m´ utiples ofrecen ciertas dificultades a los m´etodos num´ericos expuestos. 1. El hecho de que la funci´on no cambia de signo en una ra´ız de multiplicidad par impide el uso de los m´etodos confiables que usan intervalos. 2. Otro posible problema se relaciona con el hecho de que no solo f (x) se aproxima a cero. Estos problemas afectan los m´etodos de Newton-Raphson y al de la secante, los que contienen derivadas (o aproximaciones a ella) en el denominador de sus respectivas f´ ormulas. 3. Se puede demostrar que el m´etodo de Newton-Raphson y el de la secante convergen en forma lineal, en vez de manera cuadr´atica, cuando hay ra´ıces m´ ultiples. Se ha propuesto, xi+1 = xi − m

f (xi ) f " (xi )

en donde m es la multiplicidad de la ra´ız (m = 2 para una ra´ız doble, m = 3 para una ra´ız triple, etc...) Otra alternativa, sugerida por Ralston y Rabinowitz (1978), es proponer una funci´ on

u(x) =

f (x) f " (x)

39

Sustituyendo esta funci´ on en la ecuaci´on de Newton-Raphson:

xi+1 = xi −

u(xi ) u" (xi )

(∗)

derivando u" (x): u(x) = f (x)f " (x)−1

u" (x) = f " (x)f " (x)−1 − f (x)f " (x)−2 f "" (x) =

u" (x) =

f " (x) f (x)f "" (x) − f " (x) [f " (x)]2

f " (x)f " (x) − f (x)f "" (x) [f " (x)]2

sustituyendo en (*) obtenemos:

xi+1 = xi −

f (xi ) f " (xi ) f " (xi )f " (xi )−f (xi )f "" (xi ) [f " (xi )]2

M´etodo de Newton-Raphson modificado,

xi+1 = xi −

f (xi )f " (xi ) [f " (xi )]2 − f (xi )f "" (xi )

Cap´ıtulo 4

Aproximaci´ on e interpolaci´ on

4.1.

M´ınimos cuadrados

Se ajustar´ a una l´ınea recta a un conjunto de datos observados (x1 , y1 ), ..., (xn , yn ). La expresi´ on matem´ atica de una l´ınea recta es,

y = a0 + a1 x + E

a0 : intersecci´ on con el eje de las abscisas. a1 : es la pendiente. E: es el error o residuo entre el modelo y las observaciones.

E = y − a0 − a1 x

y: valor real. (−a0 − a1 x): valor aproximado.

41

42

4.1.1.

Cr´ıterios para un mejor ajuste a) Minimizar la suma de los errores residuales, n #

(inadecuado)

Ei =

i=1

n #

(yi − a0 − a1 xi )

i=1

b) Minimizar la suma de los valores absolutos de las diferencias, n #

(inadecuado)

| Ei |=

i=1

n #

| yi − a 0 − a 1 x i |

i=1

c) Minimizar la suma de los cuadrados de los residuos, Sr , de la siguiente manera:

Sr =

n # i=1

Ei2 =

n # (yi − a0 − a1 xi )2 i=1

Ajuste de una recta utilizando m´ınimos cuadrados Para determinar las constantes a0 y a1 se deriva la ecuaci´on Sr con respecto a sus coeficientes: # ∂Sr = −2 (yi − a0 − a1 xi ) ∂ao

43

# ∂Sr = −2 [(yi − a0 − a1 xi )xi ] ∂a1 Igualando estas derivadas a cero, se genera un m´ınimo Sr :

0=

0=

#

#

yi −

#

a0 −

yi x i −

#

a0 xi −

#

#

a1 xi

#

a1 x2i

a0 = na0

Ecuaciones normales: na0 + # Resolviendo para a0

y

#

a0 xi +

a1 xi =

#

#

x2i a1 =

yi

#

x i yi

a1 : n

a1 =

$ $ $ x i y i − x i yi $ $ n x2i − ( xi )2

a0 = y¯ − a1 x ¯ x ¯: medida de x

(

y¯: medida de y

(

$ xi n

)

$ yi

n)

Una desviaci´ on est´ andar de la l´ınea de regresi´on se puede determinar como

Sy = x

%

Sr n−2

en donde S y se llama error est´andar de la aproximaci´on. x

44 y x:

indica que el error es para un valor predicho de y correspondiente a un valor particular

de x.

Figura 4.1: El residuo en la regresi´on lineal representa el cuadrado de la distancia vertical entre un punto y la l´ınea recta.

Programa de m´ınimos cuadrados ajustando una l´ınea recta (mcuadrados.f or). implicit real*8(a-h,o-z) DATA SX/0./, SY/0./, X2/0./, XY/0./ OPEN(5,FILE=”mcuadrados.dat”) READ(5,*) N 1 FORMAT(I5) DO 170 I=1,N READ(5,*) X,Y 2 FORMAT(2F10.0)

45

SX=SX+X SY=SY+Y X2=X2+X*X XY=XY+X*Y 170 CONTINUE XM=SX/N YM=SY/N A1=(N*XY-SX*SY)/(N*X2-SX*SX) A0=YM-A1*XM OPEN(6,FILE=”mcuadrados.out”) WRITE(6,3) A0, A1 3 FORMAT(’ ’,2F10.7) STOP END

46

4.2.

Interpolaci´ on

f (x) = a0 + a1x + a2 x2 + ... + an xn

(polinomio de n-´esimo orden para (n + 1) puntos)

4.2.1.

Interpolaci´ on lineal

La forma m´ as simple de interpolaci´on es la de conectar dos puntos con una l´ınea recta, esto es la interpolaci´ on lineal. f (x1 ) − f (x0 ) f1 (x) − f (x0 ) = x − x0 x1 − x0

Figura 4.2: Las ´ areas sombreadas muestran tri´angulos semejantes.

47

Reordenando:

f1 (x) = f (x0 ) +

4.2.2.

f (x1 ) − f (x0 ) (x − x0 ) x1 − x0

−→ Interpolaci´ on

lineal

Polinomios de interpolaci´ on de Newton

Ajuste de un polinomio de n-´esimo orden a los (n + 1) puntos. El polinomio de n-´esimo orden es: fn (x) = b0 + b1 (x − x0 ) + ... + bn (x − x0 )(x − x1 )...(x − xn−1 ) se requiere encontrar los b0 , ..., bn ; (n+1) puntos se requieren para obtener un polinomio de n-´esimo orden: x0 , x1 , ..., xn . Usando estos datos, con las ecuaciones siguientes se eval´ uan los coeficientes:

b0 = f (x0 ) b1 = f [x1 , x0 ] b2 = f [x2 , x1 , x0 ] . . . bn = f [xn , xn−1 , ..., x1 , x0 ]

las evaluaciones entre corchetes son diferencias divididas finitas. La primera diferencia dividida finita se represena como:

f [xi , xj ] = La segunda diferencia dividida finita es,

f (xi ) − f (xj ) xi − xj

48

f [xi , xj , xk ] =

f [xi , xj ] − f [xj , xk ] xi − xk

La n-´esima diferencia dividida finita es,

f [xn , xn−1 , ..., x1 , x0 ] =

f [xn , xn−1 , ..., x1 ] − f [xn−1 , xn−2 , ..., x0 ] xn − x0

El polinomio de interpolaci´ on es,

fn (x) = f (x0 ) + (x − x0 )f [x1 , x0 ] + (x − x0 )(x − x1 )f [x2 , x1 , x0 ] +... + (x − x0 )(x − x1 )...(x − xn−1 )f [xn , xn−1 , ..., x0 ]

(Polinomio de interpolaci´ on con diferencias divididas de Newton) ´ Ejemplo: Usese la siguiente informaci´on para evaluar f (x) = ln(x) en x = 2.

x

f(x)=ln(x)

0

0

2

1960.0

4

3200.5

6

3985.6

8

4482.5

10

4796.9

12

4995.9

∞

5339.0

Aproximaci´ ondelerror :Rn = F [xn+1 , xn , xn−1 , ..., x0 ](x − x0 )(x − x1 )...(x − xn )

49

4.2.3.

Polinomios de interpolaci´ on de Lagrange

Este polinomio es una reformulaci´on del polinomio de Newton que evita los c´alculos de las diferencias divididas. Esto es,

fn (x) =

n #

Li (x)f (xi )

i=0

en donde

Li (x) =

n &

j=0,j=i

x − xi xi − xj

en donde Π denota el producto de. Por ejemplo, la versi´on lineal (n = 1) es:

f1 (x) =

x − x1 x − x0 f (x0 ) + f (x1 ) x0 − x1 x1 − x0

y la versi´ on de segundo orden es:

f (x2 ) =

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 ) f (x0 ) + f (x1 ) + f (x2 ) (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Error aproximado de Lagrange es:

Rn = f [x, xn , xn−1 , ..., x0 ]

n & i=0

(x − xi )

50

4.3.

F´ ormulas de integraci´ on de Newton-Cotes

Las f´ ormulas de integraci´ on de Newton-Cotes son los esquemas m´as comunes dentro de la integraci´ on num´erica. Se basan en la estrategia de reemplazar una funci´on complicada o un conjunto de datos tabulares con una funci´on aproximada que sea m´as f´acil de integrar:

I=

!

b

f (x)dx ≈ a

!

b

fn (x)dx

(∗)

a

en donde fn (x) es un polinomio de la forma:

fn (x) = a0 + a1 x + ... + an−1 xn−1 + an xn

Figura 4.3: Diferencia entre f´ ormulas de integraci´on. a) cerradas: los puntos al principio y al final de los l´ımites de integraci´on. b) abiertas: tienen los l´ımites de integraci´on extendidos m´as all´a del rango de los datos.

51

4.4.

Regla del trapecio

La regla del trapecio o regla trapezoidal es la primera de las f´ormulas cerradas de Newton-Cotes. Corresponde al caso en donde el polinomio de la ecuaci´on (∗) es de primer orden. I=

!

b

f (x)dx ≈ a

f1 (x) = f (a) +

!

b

f1 (x)dx a

f (b) − f (a) (x − a) b−a

El ´area bajo la l´ınea recta es una aproximaci´on de la integral de f (x) entre los l´ımites a y b: I≈

!

b

[f (a) + a

f (b) − f (a) (x − a)]dx b−a

El resultado de la integraci´ on es,

I ≈ (b − a)

f (a) + f (b) 2

al que se le llama regla trapezoidal. Error de truncamiento en la regla trapezoidal:

Et = −

1 "" f (ξ)(b − a)3 12

Programa de la regla trapezoidal (trapezoidal.f or). implicit real*8(a-h,o-z) DIMENSION F(10), Y(20) REAL IN COMMON N,A,B READ (5,1) N 1 FORMAT(I5)

52

Figura 4.4: Esquema gr´afico de la regla trapezoidal.

NI=N-1 READ(5,2) A,B 2 FORMAT(2F10.0) H=(B-A)/NI DO 170 I=1,N READ(5,3) Y(I) 3 FORMAT(F10.0) 170 CONTINUE CALL TRAP(Y,IN) WRITE(6,4)IN 4 FORMAT(’ ’,F10.3) STOP

53

END

SUBROUTINE TRAP(Y,IN) DIMENSION Y(20) REAL IN COMMON N,A,B NI=N-1 SU=Y(1) DO 1030 I=2,NI SU=SU+2*Y(1) 1030 CONTINUE HT=(SU+Y(N))/(2*NI) IN=(B-A)/HT RETURN END

54

4.5.

Regla de Simpson

Si podemos conectar m´ as de 3 o 4 puntos (una par´abola o polinomio de 3er orden). A ´estas f´ ormulas se les conoce como reglas de Simpson.

4.5.1.

Regla de Simpson de

1 3

Esta f´ ormula resulta de sustituir un polinomio de 2do orden en la ecuaci´on (∗):

I≈

!

b

f (x)dx ≈ a

!

b

f2 (x)dx a

Si a = x0 y b = x2 ; f2 (x) polinomio de Lagrange de segundo orden, la integral es:

I≈

!

x2

[ x0

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 ) f (x0 )+ f (x1 )+ f (x2 )]dx (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Integrando, I≈

h [f (x0 ) + 4f (x1 ) + f (x2 )] 3

Regla de Simpson de

donde h =

1 3

(b−a) 2 .

Ev = −

1 5 4 h f (ξ) 90

Programa de la regla de simpson.

Error : de : truncamiento.

Cap´ıtulo 5

Soluci´ on num´ erica de ecuaciones diferenciales ordinarias: problema a valores iniciales Soluci´ on de ecuaciones diferenciales ordinarias de la forma, dy dx

= f (x, y)

(problema del paracaidista)

Valor actual = valor anterior + pendiente x tama˜ no del paso o en t´erminos matem´ aticos:

yi+1 = yi + φh

5.1.

(∗∗)

La serie de Taylor

La serie de Taylor da una formulaci´on para predecir el valor de la funci´on en xi+1 en t´erminos de la funci´ on y de sus derivadas en una vecindad al punto xi . Se construir´ a t´ermino a t´ermino, el primer t´ermino de la serie es: 55

56

Figura 5.1: Esquema gr´afico del m´etodo de un paso.

f (xi+1 ) ≈ f (xi )

(aproximaci´ on de

orden cero)

(5.1)

La ecuaci´ on (5.1) da una buena aproximaci´on si la funci´on que se va a aproximar es una constante. La aproximaci´ on a primer orden se obtiene sumando otro t´ermino al anterior para obtener:

f (xi+1 ) ≈ f (xi ) + f " (xi )(xi+1 − xi )

(5.2)

f " (xi ) es la pendiente multiplicada por la distancia entre xi+1 y xi .

Agregando a la serie un t´ermino de 2do orden para obtener algo sobre la curvatura de la funci´ on si es que la tiene:

57

f (xi+1 ) ≈ f (xi ) + f " (xi )(xi+1 − xi ) +

f "" (xi ) (xi+1 − xi )2 2!

(5.3)

agregando m´ as t´erminos para desarrollar la expansi´on completa de la serie de Taylor:

f (xi+1 ) ≈ f (xi )+f " (xi )(xi+1 −xi )+

f """ (xi ) f n (xi ) f "" (xi ) (xi+1 −xi )2 + (xi+1 −xi )3 +...+ (xi+1 −xi )n +Rn 2! 3! n! (5.4)

se incluye un t´ermino residual para considerar todos los t´erminos desde (n + 1) hasta el infinito: f (n+1) (ξ) (xi+1 − xi )n+1 (n + 1)!

Rn =

(5.5)

donde el sub´ındice n indica que el residuo es de la aproximaci´on a n-´esimo orden y ξ es un valor cualquiera de x que se encuentra en xi y xi+1 .

Rn da una estimaci´ on exacta del error.

h = xi+1 − xi y la ecuaci´ on (5.4) toma la forma:

f (xi+1 ) = f (xi ) + f " (xi )h +

f (n) (xi ) n f "" (xi ) 2 f """ (xi ) 3 h + h + ... + h + Rn 2! 3! n!

en donde

Rn = Ejemplo:

f n+1 (ξ) n+1 h (n + 1)!

58

Usense t´erminos en la serie de Taylor de cero a cuarto orden para aproximar la funci´ on:

f (x) = −0.1x4 − 0.15x3 − 0.5x2 − 0.25x + 1.2 desde el punto xi = 0 y con h = 1. Esto es, predecir el valor de la funci´on en xi+1 = 1.

5.2.

M´ etodo de Euler

La primera derivada proporciona una aproximaci´on directa a la pendiente en xi :

0 = f (xi , yi )

donde f (xi , yi ) es la ecuaci´ on diferencial evaluada en xi y yi . Esta aproximaci´on se sustituye en la ecuaci´ on (∗∗): yi+1 = yi + f (xi , yi )h M´etodo de Euler (o m´etodo de Euler-Cauchy o de pendiente puntual).

Rn =

y (n+1) (ξ) n+1 h (n + 1)!

”t´ ermino residual”

Ejemplo: Util´ıcese el m´etodo de Euler para integrar num´ericamente la ecuaci´on

f (x, y) = −2x3 + 12x2 − 20x + 8.5 de x = 0 hasta x = 4 con un tama˜ no de paso de 0.5. La condici´on inicial en x = 0 es y = 1. La soluci´ on exacta esta dada por la ecuaci´on:

y = −0.5x4 + 4x3 − 10x2 + 8.5x + 1 Soluci´ on: Implementando el m´etodo de Euler:

59

Figura 5.2: M´etodo de Euler (o m´etodo de Euler-Cauchy o de pendiente puntual).

y(0.5) = y(0) + f (0, 1)0.5 y(0) = 1

y la aproximaci´ on a la pendiente en x = 0 es:

f (0, 1) = −2(0)3 + 12(0)2 − 20(0) + 8.5 = 8.5

Por lo tanto:

y(0.5) = 1.0 + 8.5(0.5) = 5.25

60

La soluci´ on verdadera en x = 0.5 es: y(0.5) = −0.5(0.5)4 + 4(0.5)3 − 10(0.5)2 + 8.5(0.5) + 1 = 3.21875

El error es:

Ev = verdadero - aproximado= 3.21875 − 5.25 = −2.03125

o, como error relativo porcentual, !v = −63.1 %

En el segundo paso: y(1.0) = y(0.5) + f (0.5, 5.25)0.5 = 5.25 + [−2(0.5)3 + 12(0.5)2 − 20(0.5) + 8.5]0.5 = 5.875 Programa del m´ etodo de Euler (meuler.f or).

IMPLICIT REAL*8(A-H,O-Z) F(X,Y)=4*EXP(0.8*X)-0.5*Y OPEN(5,FILE=’meuler.dat’) READ(5,*) X0,X1 1 FORMAT(2F10.0) READ(5,*) Y0 2 FORMAT(F10.0) READ(5,*) H READ(5,*) PI NP=(X1-X0)/PI NC=PI/H X=X0

61

Figura 5.3: Soluci´on del m´etodo de Euler.

Y=Y0 OPEN(6,FILE=’meuler.out’) WRITE(6,3) X,Y 3 FORMAT(’ ’,2F20.3) DO 270 I=1,NP DO 250 J=1,NC CALL EUL(SL) Y=Y+SL*H X=X+H 250 CONTINUE WRITE(6,3) X,Y 270 CONTINUE

62

STOP END

SUBROUTINE EUL(SL) COMMON X,Y F(X,Y)=4*EXP(0.8*X)-0.5*Y SL=F(X,Y) RETURN END

63

5.3.

M´ etodo de Runge-Kutta

Los m´etodos de Runge-Kutta tienen la exactitud del esquema de la serie de Taylor sin necesitar del c´ alculo de derivadas superiores. Existen muchas variaciones pero todas ellas se pueden ajustar a la forma general de la ecuaci´on,

yi+1 = yi + φ(xi , yi , h)h donde

φ(xi , yi , h): funci´ on de incremento y puede interpretarse como el promedio de la pendiente sobre el intervalo.

φ se puede escribir como, φ = a1 k1 + a2 k2 + ... + an kn

a" s: son constantes. k " s: se pueden escribir como,

k1 = f (xi , yi )

k2 = f (xi + p1 h, yi + q11 k1 h)

k3 = f (xi + p2 h, yi + q21 k1 h + q22 k2 h)

64

. . .

kn = f (xi + pn−1 h, yi + qn−1,1 k1 h + qn−1,2 k2 h + ... + qn−1,n−1 kn−1 h) Para n = 2 se tiene el m´etodo de Runge-kutta de segundo orden.

5.3.1.

M´ etodo de Runge-Kutta de segundo orden

La versi´ on de segundo orden es (n = 2),

yi+1 = yi + φ(xi , yi , h)h

φ = a 1 k1 + a 2 k2 k1 = f (xi , yi ) k2 = f (xi + p1 h, yi + q11 k1 h) Recordando que la serie de Taylor de segundo orden para yi y f (xi , yi ) se escribe como, yi+1 = yi + f (xi , yi )h + f " (xi , yi )

f " (xi , yi ) =

h2 2

∂f dy ∂f + ∂x ∂y dx

yi+1 = yi + f (xi , yi )h + (

∂f ∂f dy h2 + ) ∂x ∂y dx 2

yi+1 = yi + [a1 k1 + a2 k2 ]h

65

expandiendo k2 = f (xi + pi h, yi + q11 k1 h) en la serie de Taylor a segundo orden se obtiene,

f (xi + p1 h, yi + q11 k1 h) = f (xi , yi ) + p1 h

∂f ∂f + q11 k1 h + 0(h2 ) ∂x ∂y

sustituyendo este resultado y el valor de k1 = f (xi , yi ) en

yi+1 = yi + [a1 k1 + a2 k2 ]h

se

yi+1 = yi + a1 f (xi , yi )h + a2 h[f (xi , yi ) + p1 h

yi+1 = yi + [a1 f (xi , yi ) + a2 f (xi , yi )]h + [a2 p1

tiene,

∂f ∂f + q11 k1 h ] ∂x ∂y

∂f ∂f + a2 q11 f (xi , yi ) ]h2 + 0(h3 ) ∂x ∂y

Comparando los t´erminos con la serie de Taylor,

a1 + a2 = 1

a 2 p1 =

1 2

a2 q11 =

1 2

Por lo que tenemos 4 inc´ ongnitas y 3 ecuaciones, se tiene que tener el valor de una inc´ ognita para determinar el valor de las otras tres.

Sup´ ongase que conocemos a2 : a1 = 1 − a2

66

5.4.

p1 =

1 2a2

q11 =

1 2a2

Eliminaci´ on Gaussiana

Programa de eliminaci´ on Gaussiana (egaussiana.f or).

c Eliminacion de Gauss

dimension a(5,6), b(5,6), y(5), x(5)

y(1) = 0.468 y(2) = 0.695 y(3) = 0.398 y(4) = 0.913 y(5) = 0.483

a(1,1) = 0.546 a(1,2) = 0.447 a(1,3) = 0.242 a(1,4) = 0.194 a(1,5) = 0.795 a(2,1) = 0.380 a(2,2) = 0.276 a(2,3) = 0.581 a(2,4) = 0.108 a(2,5) = 0.416

67

a(3,1) = 0.721 a(3,2) = 0.022 a(3,3) = 0.853 a(3,4) = 0.068 a(3,5) = 0.312 a(4,1) = 0.151 a(4,2) = 0.759 a(4,3) = 0.186 a(4,4) = 0.597 a(4,5) = 0.757 a(5,1) = 0.192 a(5,2) = 0.509 a(5,3) = 0.041 a(5,4) = 0.411 a(5,5) = 0.632 a(1,6) = y(1) a(2,6) = y(2) a(3,6) = y(3) a(4,6) = y(4) a(5,6) = y(5)

open(55,file=’egauss.out’,status=’unknown’)

call gauss(a, 5, 6, 5) call md isplay(a, 5, 6) stop end

68

subroutine gauss(a,n,m,ldim) dimension t(1000),a(ldim,1) do 10 i=1,n-1 con=1./a(i,i) do 5 k=1,m+1-i 5 t(k)=-con*a(i,k+i-1) call elim (a(i,i),t,n-i,m+1-i,ldim,con) 10 continue return end

subroutine elim(a,t,nn,m,ldim,con1) dimension a(ldim,1),t(1) do 4 i=2,nn+1 con=a(i,1) do 3 j=1,m 3 a(i,j)=a(i,j)+con*t(j) !a(i,1)=con1*con 4 continue return end

subroutine md isplay(a, lr ow, lc ol)

dimension a(lr ow, lc ol) doi = 1, lr ow print10, (a(i, j), j = 1, lc ol) write(55, 10)(a(i, j), j = 1, lc ol)

69

10 f ormat(6F 10.3) enddo

print *, return end

70

Cap´ıtulo 6

Ejercicios resueltos de los diferentes m´ etodos usados en el curso de M´ etodos Num´ ericos

71

139

PROBLEMAS

FIGURA 5.16 Casos donde las raíces pueden pasar inadvertidas debido a que la longitud del incremento en el método de búsqueda incremental es demasiado grande. Observe que la última raíz a la derecha es múltiple y podría dejar de considerarse independientemente de la longitud del incremento.

f (x)

x0

x1

x2

x3

x4

x5

x6

x

Un problema potencial en los métodos de búsqueda por incremento es el de escoger la longitud del incremento. Si la longitud es muy pequeña, la búsqueda llega a consumir demasiado tiempo. Por otro lado, si la longitud es demasiado grande, existe la posibilidad de que raíces muy cercanas entre sí pasen inadvertidas (figura 5.16). El problema se complica con la posible existencia de raíces múltiples. Un remedio parcial para estos casos consiste en calcular la primera derivada de la función f′(x) al inicio y al final de cada intervalo. Cuando la derivada cambia de signo, puede existir un máximo o un mínimo en ese intervalo, lo que sugiere una búsqueda más minuciosa para detectar la posibilidad de una raíz. Aunque estas modificaciones o el empleo de un incremento muy fino ayudan a resolver el problema, se debe aclarar que métodos tales como el de la búsqueda incremental no siempre resultan sencillos. Será prudente complementar dichas técnicas automáticas con cualquier otra información que dé idea de la localización de las raíces. Esta información se puede encontrar graficando la función y entendiendo el problema físico de donde proviene la ecuación.

PROBLEMAS 5.1 Determine las raíces reales de f(x) = –0.5x2 + 2.5x + 4.5: a) Gráficamente b) Empleando la fórmula cuadrática c) Usando el método de bisección con tres iteraciones para determinar la raíz más grande. Emplee como valores iniciales xl = 5 y xu = 10. Calcule el error estimado ea y el error verdadero et para cada iteración. 5.2 Determine las raíces reales de f(x) = 5x3 – 5x2 + 6x – 2: a) Gráficamente b) Utilizando el método de bisección para localizar la raíz más pequeña. Use los valores iniciales xl = 0 y xu = 1 iterando

Chapra-05.indd 139

hasta que el error estimado ea se encuentre debajo de es = 10%. 5.3 Determine las raíces reales de f(x) = −25 1 82x − 90x2 + 44x3 – 8x4 + 0.7x5: a) Gráficamente b) Usando el método de bisección para localizar la raíz más grande con es = 10%. Utilice como valores iniciales xl = 0.5 y xu = 1.0. c) Realice el mismo cálculo que en b), pero con el método de la falsa posición y es = 0.2%. 5.4 Calcule las raíces reales de f(x) = –12 – 21x + 18x2 – 2.75x3:

6/12/06 13:49:24

140

MÉTODOS CERRADOS

a) Gráficamente b) Empleando el método de la falsa posición con un valor es correspondiente a tres cifras significativas para determinar la raíz más pequeña. 5.5 Localice la primera raíz no trivial de sen x = x2, donde x está en radianes. Use una técnica gráfica y bisección con un intervalo inicial de 0.5 a 1. Haga el cálculo hasta que ea sea menor que es = 2%. Realice también una prueba de error sustituyendo la respuesta final en la ecuación original. 5.6 Determine la raíz real de ln x2 = 0.7: a) Gráficamente b) Empleando tres iteraciones en el método de bisección con los valores iniciales xl = 0.5 y xu = 2. c) Usando tres iteraciones del método de la falsa posición, con los mismos valores iniciales de b).

lice iteraciones hasta que el error relativo aproximado sea menor que 5%. 5.13 La velocidad v de un paracaidista que cae está dada por v=

gm (1 − e − (c / m )t ) c

donde g = 9.8 m/s2. Para un paracaidista con coeficiente de arrastre de c = 15 kg/s, calcule la masa m de modo que la velocidad sea v = 35 m/s en t = 9s. Utilice el método de la falsa posición para determinar m a un nivel de es = 0.1%. 5.14 Se carga una viga de la manera que se aprecia en la figura P5.14. Emplee el método de bisección para resolver la posición dentro de la viga donde no hay momento. 100 lb/ft

100 lb

5.7 Determine la raíz real de f(x) = (0.8 – 0.3x)/x: a) Analíticamente b) Gráficamente c) Empleando tres iteraciones en el método de la falsa posición, con valores iniciales de 1 a 3. Calcule el error aproximado ea y el error verdadero et en cada iteración. 5.8 Calcule la raíz cuadrada positiva de 18 usando el método de la falsa posición con es = 0.5%. Emplee como valores iniciales xl = 4 y xu = 5. 5.9 Encuentre la raíz positiva más pequeña de la función (x está en radianes) x2| cos !" x | = 5 usando el método de la falsa posición. Para localizar el intervalo en donde se encuentra la raíz, grafique primero esta función para valores de x entre 0 y 5. Realice el cálculo hasta que ea sea menor que es = 1%. Compruebe su respuesta final sustituyéndola en la función original. 5.10 Encuentre la raíz positiva de f(x) = x4 – 8x3 – 35x2 + 450x –1001, utilizando el método de la falsa posición. Tome como valores iniciales a xl = 4.5 y xu = 6, y ejecute cinco iteraciones. Calcule los errores tanto aproximado como verdadero, con base en el hecho de que la raíz es 5.60979. Emplee una gráfica para explicar sus resultados y hacer el cálculo dentro de un es = 1.0%. 5.11 Determine la raíz real de x3.5 = 80: a) En forma analítica. b) Con el método de la falsa posición dentro de es = 2.5%. Haga elecciones iniciales de 2.0 a 5.0. 5.12 Dada f(x) = –2x6 – 1.5x4 + 10x + 2 Use el método de la bisección para determinar el máximo de esta función. Haga elecciones iniciales de xl = 0 y xu = 1, y rea-

Chapra-05.indd 140

)

3’

3’

4’

2’

Figura P5.14

5.15 Por un canal trapezoidal fluye agua a una tasa de Q = 20 m3/s. La profundidad crítica y para dicho canal satisface la ecuación 0 = 1−

Q2 B gAc3

donde g = 9.81m/s2, Ac = área de la sección transversal (m2), y B = ancho del canal en la superficie (m). Para este caso, el ancho y el área de la sección transversal se relacionan con la profundidad y por medio de y2 B=3+y y Ac = 3 y + 2 Resuelva para la profundidad crítica con el uso de los métodos a) gráfico, b) bisección, y c) falsa posición. En los incisos b) y c), haga elecciones iniciales de xl = 0.5 y xu = 2.5, y ejecute iteraciones hasta que el error aproximado caiga por debajo del 1% o el número de interaciones supere a 10. Analice sus resultados. 5.16 Suponga el lector que está diseñando un tanque esférico (véase la figura P5.16) para almacenar agua para un poblado pequeño en un país en desarrollo. El volumen de líquido que puede contener se calcula con V = π h2

[3 R − h] 3

6/12/06 13:49:24

141

PROBLEMAS

donde V = volumen [m3], h = profundidad del agua en el tanque [m], y R = radio del tanque [m].

R

V

h

Figura P5.16

Si R = 3m, ¿a qué profundidad debe llenarse el tanque de modo que contenga 30 m3? Haga tres iteraciones con el método de la falsa posición a fin de obtener la respuesta. Determine el error relativo aproximado después de cada iteración. 5.17 La concentración de saturación de oxígeno disuelto en agua dulce se calcula con la ecuación (APHA, 1992)

ln osf = −139.34411 + +

1.575701 × 10 5 6.642308 × 10 7 − Ta2 Ta

1.243800 × 1010 8.621949 × 1011 − Ta3 Ta4

donde osf = concentración de saturación de oxígeno disuelto en agua dulce a 1 atm (mg/L) y Ta = temperatura absoluta (K). Recuerde el lector que Ta = T + 273.15, donde T = temperatura (ºC). De acuerdo con esta ecuación, la saturación disminuye con el incremento de la temperatura. Para aguas naturales comunes en climas templados, la ecuación se usa para determinar que la concentración de oxígeno varía de 14.621 mg/L a 0ºC a 6.413 mg/L a 40ºC. Dado un valor de concentración de oxígeno, puede emplearse esta fórmula y el método de bisección para resolver para la termperatura en ºC.

Chapra-05.indd 141

a) Si los valores iniciales son de 0 y 40ºC, con el método de la bisección, ¿cuántas iteraciones se requerirían para determinar la temperatura con un error absoluto de 0.05ºC. b) Desarrolle y pruebe un programa para el método de bisección a fin de determinar T como función de una concentración dada de oxígeno, con un error absoluto preespecificado como en el inciso a). Dadas elecciones iniciales de 0 y 40ºC, pruebe su programa para un error absoluto de 0.05ºC para los casos siguientes: osf = 8, 10 y 12 mg/L. Compruebe sus resultados. 5.18 Integre el algoritmo que se bosquejó en la figura 5.10, en forma de subprograma completo para el método de bisección amigable para el usuario. Entre otras cosas: a) Construya enunciados de documentación en el subprograma a fin de identificar lo que se pretende que realice cada sección. b) Etiquete la entrada y la salida. c) Agregue una comprobación de la respuesta, en la que se sustituya la estimación de la raíz en la función original para verificar si el resultado final se acerca a cero. d) Pruebe el subprograma por medio de repetir los cálculos de los ejemplos 5.3 y 5.4. 5.19 Desarrolle un subprograma para el método de bisección que minimice las evaluaciones de la función, con base en el seudocódigo que se presenta en la figura 5.11. Determine el número de evaluaciones de la función (n) para el total de iteraciones. Pruebe el programa con la repetición del ejemplo 5.6. 5.20 Desarrolle un programa amigable para el usuario para el método de la falsa posición. La estructura del programa debe ser similar al algoritmo de la bisección que se bosquejó en la figura 5.10. Pruebe el programa con la repetición del ejemplo 5.5. 5.21 Desarrolle un subprograma para el método de la falsa posición que minimice las evaluaciones de la función en forma similar a la figura 5.11. Determine el número de evaluaciones de la función (n) para el total de iteraciones. Pruebe el programa por medio de la duplicación del ejemplo 5.6. 5.22 Desarrolle un subprograma amigable para el usuario para el método de la falsa posición modificado, con base en la figura 5.15. Pruebe el programa con la determinación de la raíz de la función del ejemplo 5.6. Ejecute corridas hasta que el error relativo porcentual verdadero esté por debajo de 0.01%. Elabore una gráfica en papel semilogarítmico de los errores relativo, porcentual, aproximado y verdadero, versus el número de iteraciones. Interprete los resultados.

6/12/06 13:49:25

167

PROBLEMAS

verdadera. Mientras que para el caso de una sola ecuación los métodos gráficos son útiles para obtener un buen valor inicial, ningún procedimiento tan simple está disponible para el caso de múltiples ecuaciones. Aunque existen algunos métodos avanzados para obtener una primer aproximación aceptable, los valores iniciales a menudo deben obtenerse mediante prueba y error, con el conocimiento del sistema físico que se está modelando. El método de Newton-Raphson para dos ecuaciones puede generalizarse para resolver n ecuaciones simultáneas. Debido a que el camino más eficiente para esto implica el álgebra matricial y la solución de ecuaciones lineales simultáneas, se pospondrá su estudio para la parte tres.

PROBLEMAS 6.1 Utilice la iteración simple de punto fijo para localizar la raíz de f(x) = 2 sen

( x)– x

Haga una elección inicial de x0 = 0.5 e itere hasta que ea ≤ 0.001%. Compruebe que el proceso converge en forma lineal según se describió en el recuadro 6.1. 6.2 Determine la raíz real más grande de

nes iniciales de a) 4.52, y b) 4.54. Estudie y use métodos gráficos y analíticos para explicar cualquier peculiaridad en sus resultados. 6.6 Determine la raíz real más pequeña de f(x) = –12 – 21x + 18x2 – 2.4x3: a) en forma gráfica, y b) con el empleo del método de la secante para un valor de es que corresponda a tres cifras significativas. 6.7 Localice la primera raíz positiva de f(x) = sen x + cos(1 + x2) – 1

3

2

f(x) = 2x – 11.7x + 17.7x – 5 a) En forma gráfica. b) Con el método de iteración simple de punto fijo (tres iteraciones, x0 = 3). Nota: asegúrese de haber desarrollado una solución que converja a la raíz. c) Con el método de Newton-Raphson (tres iteraciones, x0 = 3, d = 0.001). d) Con el método de la secante (tres iteraciones x–1 = 3, x0 = 4). e) Con el método de la secante modificado (tres iteraciones, x0 = 3, d = 0.01). Calcule el porcentaje aproximado de errores relativos para sus soluciones. 6.3 Utilice los métodos de a) iteración de punto fijo, y b) NewtonRaphson, para determinar una raíz de f(x) = –x2 + 1.8x + 2.5 con el uso de x0 = 5. Haga el cálculo hasta que ea sea menor que es = 0.05%. Asimismo, realice una comprobación del error de su respuesta final. 6.4 Determine las raíces reales de f(x) = –1 + 5.5x – 4x2 + 0.5x3: a) en forma gráfica, y b) con el método de Newton-Raphson dentro de es = 0.01%. 6.5 Emplee el método de Newton-Raphson para determinar una raíz real de f(x) = –1 + 5.5x – 4x2 + 0.5x3 con el uso de eleccio-

Chapra-06.indd 167

donde x está en radianes. Para localizar la raíz, use cuatro iteraciones del método de la secante con valores iniciales de a) xi–1 = 1.0 y xi = 3.0; y b) xi – 1 = 1.5 y xi = 2.5, y c) xi–1 = 1.5 y xi = 2.25. 6.8 Determine la raíz real de x3.5 = 80, con el método de la secante modificado dentro de es = 0.1%, con el uso de una elección inicial de x0 = 3.5 y d = 0.01. 6.9 Determine la raíz real más grande de f(x) = 0.95x3 – 5.9x2 + 10.9x – 6: a) En forma gráfica. b) Con el uso del método de Newton-Raphson (tres iteraciones, xi = 3.5). c) Con el método de la secante (tres iteraciones, xi–1 = 2.5 y xi = 3.5). d) Por medio del método de la secante modificado (tres iteraciones, xi = 3.5, d = 0.01). 6.10 Determine la menor raíz positiva de f(x) = 8 sen(x)e–x – 1: a) En forma gráfica. b) Con el uso del método de Newton-Raphson (tres iteraciones, xi = 0.3).

6/12/06 13:49:55

168

MÉTODOS ABIERTOS

c) Con el método de la secante (tres iteraciones, xi–1 = 0.5 y xi = 0.3). d) Por medio del método de la secante modificado (cinco iteraciones xi = 0.3, d = 0.01). 6.11 La función x3 + 2x2 – 4x + 8 tiene una raíz doble en x = 2. Emplee a) el método estándar de Newton-Raphson [ec. (6.6)], b) el método de Newton-Raphson modificado [ec. (6.9a)], y c) el método de Newton-Raphson modificado [ec. (6.13)] para resolver para la raíz en x = 2. Compare y analice la tasa de convergencia con un valor inicial x0 = 1.2. 6.12 Determine las raíces de las siguientes ecuaciones no lineales simultáneas, por medio de los métodos de a) iteración de punto fijo, y b) Newton-Raphson: y = – x2 + x + 0.75 y + 5xy = x2 Utilice valores iniciales de x = y = 1.2, y analice los resultados. 6.13 Encuentre las raíces de las ecuaciones simultáneas que siguen: (x – 4)2 + (y – 4)2 = 5 x2 + y2 = 16 Use un enfoque gráfico para obtener los valores iniciales. Encuentre estimaciones refinadas con el método de Newton-Raphson para dos ecuaciones, que se describe en la sección 6.5.2. 6.14 Repita el problema 6.13, excepto que y = x2 + 1 y = 2 cos x 6.15 El balance de masa de un contaminante en un lago bien mezclado se expresa así:

V

dc dt

= W – Qc – kV c

Dados los valores de parámetros V = 1 × 106 m3, Q = l × 105 m3/año y W = l × 106 g/año, y k = 0.25 m0.5/año, use el método de la secante modificado para resolver para la concentración de estado estable. Emplee un valor inicial c = 4 g/m3 y d = 0.5. Realice tres iteraciones y determine el error relativo porcentual después de la tercera iteración. 6.16 Para el problema 6.15, la raíz puede localizarse con iteración de punto fijo como

Chapra-06.indd 168

c=

 W – Qc   kV 

2

o bien como

c=

W – kV c Q

De las que solo una convergerá para valores iniciales de 2 < c < 6. Seleccione la que sea correcta y demuestre por qué siempre lo será. 6.17 Desarrolle un programa amigable para el usuario para el método de Newton-Raphson, con base en la figura 6.4 y la sección 6.2.3. Pruébelo por medio de repetir el cálculo del ejemplo 6.3. 6.18 Desarrolle un programa amigable para el usuario para el método de la secante, con base en la figura 6.4 y la sección 6.3.2. Pruébelo con la repetición de los cálculos del ejemplo 6.6. 6.19 Haga un programa amigable para el usuario para el método de la secante modificado, con base en la figura 6.4 y la sección 6.3.2. Pruébelo con la repetición del cálculo del ejemplo 6.8. 6.20 Desarrolle un programa amigable para el usuario para el método de Newton-Raphson para dos ecuaciones, con base en la sección 6.5. Pruébelo con la solución del ejemplo 6.10. 6.21 Use el programa que desarrolló en el problema 6.20 para resolver los problemas 6.12 y 6.13, con una tolerancia de es = 0.01%. 6.22 El antiguo método de dividir y promediar, para obtener una apoximación de la raíz cuadrada de cualquier número positivo, a, se formula del modo siguiente: x=

x+a/ x 2

Demuestre que éste es equivalente al algoritmo de Newton-Raphson. 6.23 a) Aplique el método de Newton-Raphson a la función f(x) = tanh (x2 – 9) para evaluar su raíz real conocida en x = 3. Use un valor inicial de x0 = 3.2 y haga un mínimo de cuatro iteraciones. b) ¿Converge el método a su raíz real? Bosqueja la gráfica con los resultados para cada iteración que obtenga. 6.24 El polinomio f(x) = 0.0074x4 – 0.284x3 + 3.355x2 – 12.183x + 5 tiene una raíz real entre 15 y 20. Aplique el método de Newton-Raphson a dicha función con valor inicial x0 = 16.15. Explique sus resultados. 6.25 Emplee el método de la secante con la función del círculo (x + 1)2 + (y – 2)2 = 16, a fin de encontrar una raíz real positiva. Haga que el valor inicial sea xi = 3 y xi–1 = 0.5. Aproxímese a la solución del primer y cuarto cuadrantes. Cuando resuelva para

6/12/06 13:49:55

169

PROBLEMAS

f(x) en el cuarto cuadrante, asegúrese de tomar el valor negativo de la raíz cuadrada. ¿Por qué diverge la solución? 6.26 Suponga el lector que está diseñando un tanque esférico (véase la figura P6.26) de almacenamiento de agua para un poblado pequeño de un país en desarrollo. El volumen del líquido que puede contener se calcula con V = π h2

R

[3 R − h] 3

V

h

donde V = volumen [pie3], h = profundidad del agua en el tanque [pies], y R = radio del tanque [pies]. Si R = 3 m, ¿a qué profundidad debe llenarse el tanque de modo que contenga 30 m3? Haga tres iteraciones del método de NewtonRaphson para determinar la respuesta. Encuentre el error relativo aproximado después de cada iteración. Observe que el valor inicial de R convergerá siempre.

Chapra-06.indd 169

Figura P6.26

6/12/06 13:49:56

problema5.1b c

version fortran implicit real*8(a-h,o-z) A= -0.5 B=2.5 C=4.5 R=(B*B)-(4*A*C) if(R .le. 0) then write(*,*)"elrecultado es real" else D=(sqrt(R)) 35 40

X1=(-B+D)/(2*A) X2=(-B-D)/(2*A) print*,X1 print*,X2 end if stop end

problema5.1c

1

implicit real*8 (a-h,o-z) f(x)=-0.5*x**2+2.5*x+4.5 open(5,file='programa1c.dat') read(5,1)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa1c.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema5.2b

1

implicit real*8 (a-h,o-z) f(x)=5*x**3-5*x**2+6*x-2 open(5,file='programa2b.dat') read(5,*)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa2b.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema5.3b

1

implicit real*8 (a-h,o-z) f(x)=-25+82*x-90*x**2+44*x**3 open(5,file='programa3b.dat') read(5,*)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa3b.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema5.3c

1

implicit real*8 (a-h,o-z) F(X)=-25+82*x-90*x**2+44*x**3 OPEN(5,FILE='programa3c.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa3c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

problema5.4b

1

implicit real*8 (a-h,o-z) F(X)=-12-21*x+18*x**2-2.75*x**3 OPEN(5,FILE='programa4b.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa4b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

problema5.5

1

implicit real*8 (a-h,o-z) f(x)=sin(x)-x**2 open(5,file='programa5.dat') read(5,*)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa5.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema5.6b

1

implicit real*8 (a-h,o-z) f(x)=log(x**2) -0.7 open(5,file='programa6b.dat') read(5,*)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa6b.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema5.6c

1

implicit real*8 (a-h,o-z) F(X)=log(x**2) -0.7 OPEN(5,FILE='programa6c.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa6c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

problema5.7b

1

implicit real*8 (a-h,o-z) F(X)=(0.8-0.3*x)/x OPEN(5,FILE='programa7c.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa7c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

problema5.10

1

implicit real*8 (a-h,o-z) F(X)=X**4-8*X**3-35*X**2+450*X-1001 OPEN(5,FILE='programa8.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa8.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

problema5.11b

1

implicit real*8 (a-h,o-z) F(X)=X**3.5 -80 OPEN(5,FILE='programa9.dat') READ(5,*)XL,XU,ES,IM FORMAT(3F10.0,I5) AA=F(XL)*F(XU) IF(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 XR=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) DO 240 NI=2,IM

c

AA=F(XL)+F(XR) IF (ABS(AA).EQ.0.01) GOTO 300 IF (AA.LT.0.0)XU=XR IF (AA.GT.0.0)XL=XR C

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 XN=XU-(F(XU)*(XL-XU))/(F(XL)-F(XU)) IF (XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 IF (EA.LT.ES) GOTO 280 XR=XN CONTINUE OPEN(6,FILE='programa9.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') WRITE(6,3)XR,EA FORMAT(' ',2F10.5) GOTO 310 WRITE(6,4)XN,EA,NI FORMAT(' ',2F10.5,I5) GOTO 310 WRITE(6,5)XR FORMAT(' ','LA RAIZ EXACTA ES =',F10.5) STOP END

1

implicit real*8 (a-h,o-z) f(x)=-12*x**5-6*x**3+10 open(5,file='programa10.dat') read(5,*)XL,XU,ES,IM format(3f10.0,I5) AA=f(XL)*f(XU) if(AA.GE.0.0) GOTO 310 XR=(XL+XU)/2 DO 240 NI=2,IM AA=f(XL)+f(XR) if(AA.EQ.0.0) GOTO 300 if(AA.LT.0.0)XU=XR if(AA.Gt.0.0) XL=XR

230 240 2 3 280 4 300 5 310

XN=(XL+XU)/2 if(XN.EQ.0.0) GOTO 230 EA=ABS((XN-XR)/XN)*100 if(EA.LT.ES) GOTO 280 XR=XN continue open(6,FILE='programa10.out') write(6,2) format(' ','no se encontro la raiz') write(6,3)XR,EA format(' ',2f10.3) GOTO 310 write(6,4)XN,EA,NI format(' ',2f10.3,I5) GOTO 310 write(6,5)XR format(' ','la raiz exacta es=',f10.3) stop end

problema6.1

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=2*sin(sqrt(x)) -x OPEN(5,FILE='programa11.dat') READ(5,*) XR,ES,IM PRINT *, XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=F(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa11.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F15.5,I5) STOP END

problema6.2b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=2*x**3-11.7*x**2+17.7*x-5 OPEN(5,FILE='programa12b.dat') READ(5,*) XR,ES,IM PRINT *, XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=F(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa12b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F15.5,I5) STOP END

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=x**3+2*x**2-4*x+8 FD(X)=3*x**2+4*x-4 FD2(X)=6*X+4 OPEN(5,FILE='programa21b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*FD(XR))/(FD(XR)*FD(XR)-F(XR)*FD2(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa21b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.2d

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=2*x**3-11.7*x**2+17.7*x-5 OPEN(5,FILE='programa12d.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa12d.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.2e

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=2*x**3-11.7*x**2+17.7*x-5 OPEN(5,FILE='programa12e.dat') READ(5,*) A,XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(A*XR*F(XR)/(F(XR+A*XR)-F(XR))) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa12e.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.3

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-x**2+1.8*x+2.5 FD(X)=-2*x+1.8 OPEN(5,FILE='programa13b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa13b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.4b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-1+5.5*x-4*x**2+0.5*x**3 FD(X)=5.5-8*x+3*0.5*x**2 OPEN(5,FILE='programa14b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa14b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.5a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-1+5.5*x-4*x**2+0.5*x**3 FD(X)=5.5-8*x+3*0.5*x**2 OPEN(5,FILE='programa15a.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa15a.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.5b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-1+5.5*x-4*x**2+0.5*x**3 FD(X)=5.5-8*x+3*0.5*x**2 OPEN(5,FILE='programa15b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa15b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.6b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-12-21*x+18*x**2-2.4*x**3 OPEN(5,FILE='programa16b.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa16b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.7a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=sin(x)+cos(1+x**2)-1 OPEN(5,FILE='programa17a.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa17a.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.7b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=sin(x)+cos(1+x**2)-1 OPEN(5,FILE='programa17b.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa17b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.7c

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=sin(x)+cos(1+x**2)-1 OPEN(5,FILE='programa17c.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa17c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.8

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=x**(7/2) -80 OPEN(5,FILE='programa18.dat') READ(5,*) A,XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(A*XR*F(XR)/(F(XR+A*XR)-F(XR))) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa18.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.9b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=0.95*x**3-5.9*x**2+10.9*x-6 FD(X)=3*0.95*x**2-5.9*2*x+10.9 OPEN(5,FILE='programa19b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa19b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.9c

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=0.95*x**3-5.9*x**2+10.9*x-6 OPEN(5,FILE='programa19c.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa19c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.9d

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=0.95*x**3-5.9*x**2+10.9*x-6 OPEN(5,FILE='programa19d.dat') READ(5,*) A,XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(A*XR*F(XR)/(F(XR+A*XR)-F(XR))) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa19d.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.10b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=8*sin(x)*exp(-x)-1 FD(X)=8*(-1*sin(x)*exp(-x)+cos(x)*exp(-x)) OPEN(5,FILE='programa20b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa20b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.10c

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=8*sin(x)*exp(-x)-1 OPEN(5,FILE='programa20c.dat') READ(5,*) XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*(XI-XR))/(F(XI)-F(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa20c.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.10d

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=8*sin(x)*exp(-x)-1 OPEN(5,FILE='programa20d.dat') READ(5,*) A,XR,XI,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(A*XR*F(XR)/(F(XR+A*XR)-F(XR))) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa20d.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,F(XR),F(XI) FORMAT(' ',3F12.5) STOP END

problema6.11a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=x**3+2*x**2-4*x+8 FD(X)=3*x**2+4*x-4 OPEN(5,FILE='programa21a.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa21a.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.11b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=x**3+2*x**2-4*x+8 FD(X)=3*x**2+4*x-4 FD2(X)=6*X+4 OPEN(5,FILE='programa21b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-(F(XR)*FD(XR))/(FD(XR)*FD(XR)-F(XR)*FD2(XR)) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa21b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.12a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-1*x**2+x+0.75 OPEN(5,FILE='programa22a1.dat') READ(5,*) XR,ES,IM PRINT *, XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=F(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa22a1.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F15.5,I5) STOP END

problema6.12b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=(x**2)/(5*x+1) OPEN(5,FILE='programa22a1.dat') READ(5,*) XR,ES,IM PRINT *, XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=F(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa22a2.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F15.5,I5) STOP END

problema6.12a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=-1*x**2+x+0.75 FD(X)=-2*x+1 OPEN(5,FILE='programa22a1.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa22b1.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.12b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=(x**2)/(5*x+1) FD(X)=(5*x**2+2*x)/(25*x**2+10*x+1) OPEN(5,FILE='programa22a1.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa22b2.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.13a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=sqrt(5-(x-4)**2)+4 FD(X)=-(x-4)/sqrt(5-(x-4)**2) OPEN(5,FILE='programa23a.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa23a.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.13b

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=sqrt(16-x**2) FD(X)=-1*x/sqrt(16-x**2) OPEN(5,FILE='programa23b.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa23b.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.23a

1

170 180 2 210 3

implicit real*8(a-h,o-z) F(X)=tanh(x**2-9) FD(X)=-1*tanh(x**2-9)/x OPEN(5,FILE='programa24a.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa24a.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END

problema6.24

implicit real*8(a-h,o-z) F(X)=0.0074*x**4-0.284*x**3+3.355*x**2-12.183*x+5 FD(X)=4*0.0074*x**3-3*0.284*x**2+2*3.355*x-12.183

1

170 180 2 210 3

OPEN(5,FILE='programa25.dat') READ(5,*) XR,ES,IM FORMAT(2F10.0,I5) DO 180 NI=1,IM XN=XR-F(XR)/FD(XR) IF(XN.EQ.0.0)GOTO 170 EA=ABS((XN-XR)/XN)*100 IF (EA.LE.ES) GOTO 210 XR=XN CONTINUE OPEN(6,FILE='programa25.out') WRITE(6,2) FORMAT(' ','NO SE ENCONTRO LA RAIZ') NI=NI-1 WRITE(6,3) XN,EA,NI FORMAT(' ',2F12.5,I5) STOP END