METODOS NUMERICOS Y

APLICACIONES A LA FISICA Y
LA INGENIERIA CON C
Impresión Noviembre 2022 :
dell:REPO/[Link]
BY

LUIS CARDON

SALTA
9 de Noviembre 2022
PUBLIQUE O PEREZCA
2
Contents

1 Búsqueda de raı́ces 1
1.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Ejemplo: un problema de Mecánica . . . . . . . . . . . . . . . . . . . 2
1.2 Raı́ces. Fundamentos matemáticos . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Fı́sica de la caı́da libre con fricción. OMITIR. EN ELABORACION. . . . . . 6
1.3.1 REVISAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Métodos numéricos para la búsqueda de raı́ces . . . . . . . . . . . . . . . . 10
1.5 Método de la bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.1 Ventajas y desventajas . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6 Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6.1 Convergencia del método de la bisección . . . . . . . . . . . . . . . . 13
1.7 Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.7.1 Exactitud y precisión. . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7.2 Ejemplo: raı́z de f (x) = x6 − x − 1 . . . . . . . . . . . . . . . . . . . 17
1.8 Resolución computacional. C . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.8.1 Ejemplo: cálculo del coeficiente de arrastre. . . . . . . . . . . . . . . 26
1.9 Resolución computacional con Python . . . . . . . . . . . . . . . . . . . . . 30
1.9.1 Funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.9.2 Definción de funciones y módulos . . . . . . . . . . . . . . . . . . . . 31
1.9.3 Gráfico de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.9.4 Trabajando desde la terminal interactiva . . . . . . . . . . . . . . . . 38
1.9.5 El método de la Bisección en Python . . . . . . . . . . . . . . . . . . 39
1.10 El método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . 44
1.10.1 Algoritmo de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1.10.2 Ventajas y desventajas . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1.10.3 Resolución computacional . . . . . . . . . . . . . . . . . . . . . . . . 46
1.11 Los métodos de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

i
ii

1.11.1 Implementación computacional: método de la secante sencillo . . . . 49

1.11.2 Método de la Falsa Posición . . . . . . . . . . . . . . . . . . . . . . . 50
1.12 Métodos iterativos de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . . 58
1.12.1 Interpretación geométrica . . . . . . . . . . . . . . . . . . . . . . . . 58
1.12.2 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
1.13 Ceros de sistemas de ecuaciones acopladas . . . . . . . . . . . . . . . . . . . 60
1.14 Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
1.14.1 Problema 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
1.14.2 Problema 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
1.14.3 Problema 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
1.14.4 Problema 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
1.15 Una biblioteca personal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
1.16 Bibliotecas numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
1.16.1 Un ejemplo de aplicación con la GNU Scientific library . . . . . . . . 68

2 Búsqueda de raı́ces: aplicaciones 59

2.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
2.1.1 Ejemplo: cálculo del factor de fricción en conductos . . . . . . . . . . 59
2.1.2 Dimensionamiento de conductos . . . . . . . . . . . . . . . . . . . . . 64
2.2 REVISAR, TERMINAR o SACAR . . . . . . . . . . . . . . . . . . . . . . . 69
2.3 Relaciones psicrométricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
2.3.1 Problema. Punto de rocı́o . . . . . . . . . . . . . . . . . . . . . . . . 70
2.3.2 Programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.4 Problemas. REVISAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.4.1 Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.4.2 Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.4.3 Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.5 Método de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

3 Métodos iterativos: aplicaciones 75

3.1 Modelos cero dimensionales . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.2 La analogı́a térmico-eléctrica. Caso estacionario . . . . . . . . . . . . . . . . 76
3.2.1 Resistencia térmica según el modo de la transferencia . . . . . . . . . 78
3.2.2 Resistencia convectiva . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.2.3 Resistencia radiativa . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
A-cursos/cursoMN/[Link] - August 22, 2023 iii

3.3 Problema 0. Flujo de calor en una pared compuesta . . . . . . . . . . . . . . 80

3.4 Problema 1. Radiación entre dos placas paralelas . . . . . . . . . . . . . . . 80
3.5 Problema 2. Radiación entre placas paralelas . . . . . . . . . . . . . . . . . . 80
3.5.1 Planteo de la solución . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.5.2 Solución numérica 1: bisección . . . . . . . . . . . . . . . . . . . . . . 83
3.5.3 Solución numérica: Newton . . . . . . . . . . . . . . . . . . . . . . . 83
3.5.4 Aproximaciones sucesivas 1 . . . . . . . . . . . . . . . . . . . . . . . 84
3.5.5 Aproximaciones sucesivas 2 . . . . . . . . . . . . . . . . . . . . . . . 84
3.6 Problema 3. Convección y radiación combinada . . . . . . . . . . . . . . . . 86
3.7 Cálculo de las resistencias convectivas . . . . . . . . . . . . . . . . . . . . . . 87
3.8 Convección interna y externa . . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.8.1 Radiación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.9 Como resolver por aproximaciones sucesivas . . . . . . . . . . . . . . . . . . 88
3.10 Ejemplo: cálculo de la resistencia térmica global de un colector solar. . . . . 94
3.11 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.11.1 Una cubierta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.11.2 Dos cubierta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

4 Sistemas de funciones no lineales 111

4.1 Funciones lineales y no lineales . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.2 Sistemas de ecuaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.3 El método de Newton-Raphson para sistemas de ecuaciones no lineales . . . 112

5 Integración de funciones 115

5.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.2 La integral definida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.3 Cuadraturas: la evaluación de integrales definidas . . . . . . . . . . . . . . . 117
5.3.1 Fórmulas cerradas clásicas . . . . . . . . . . . . . . . . . . . . . . . . 118
5.4 La regla del trapecio compuesta . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.4.1 Refinamiento automático . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.4.2 Fórmula de Euler-Maclaurin . . . . . . . . . . . . . . . . . . . . . . . 127
5.4.3 Métodos avanzados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.4.4 Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.5 Regla de Simpson compuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
5.6 Regla de Simpson compuesta de 3/8 . . . . . . . . . . . . . . . . . . . . . . 134
iv

5.6.1 Programa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

6 Cuadraturas. Aplicaciones 137

6.1 Ejemplo. La ley de Plank . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
6.1.1 Programa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
6.1.2 La constante de Boltzman . . . . . . . . . . . . . . . . . . . . . . . . 139
6.1.3 Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

7 Ecuaciones Diferenciales Ordinarias 143

7.1 Introducción: ecuaciones diferenciales . . . . . . . . . . . . . . . . . . . . . . 143
7.2 Ecuaciones diferenciales ordinarias: EDO . . . . . . . . . . . . . . . . . . . . 146
7.2.1 Ecuaciones diferenciales ordinarias de orden superior . . . . . . . . . 147
7.2.2 Problemas de valor inicial . . . . . . . . . . . . . . . . . . . . . . . . 148
7.2.3 Problemas de valor de contorno a dos puntos . . . . . . . . . . . . . . 149
7.2.4 Algunos problemas de origen fı́sico . . . . . . . . . . . . . . . . . . . 150
7.2.5 Caı́da libre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
7.2.6 Descarga de un tanque de agua . . . . . . . . . . . . . . . . . . . . . 151
7.2.7 Descarga de un capacitor . . . . . . . . . . . . . . . . . . . . . . . . . 152
7.2.8 Enfriamiento de un cuerpo. . . . . . . . . . . . . . . . . . . . . . . . 152
7.3 El Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
7.3.1 Discretización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
7.3.2 Aproximación en diferencias finitas . . . . . . . . . . . . . . . . . . . 157
7.3.3 Interpretación Geométrica . . . . . . . . . . . . . . . . . . . . . . . . 157
7.3.4 Aproximación mediante integrales . . . . . . . . . . . . . . . . . . . . 158
7.4 El método de Euler: resolución numérica. . . . . . . . . . . . . . . . . . . . . 158
7.5 Resolución numérica del problema de enfriamiento de un cuerpo. . . . . . . 161
7.5.1 Sobre el dominio de la variable independiente y el paso de integración 165
7.5.2 Influencia del paso de integración . . . . . . . . . . . . . . . . . . . . 165
7.5.3 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
7.6 Uso de estructuras en la programación de funciones . . . . . . . . . . . . . . 167
7.7 Métodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
7.7.1 El método de Euler modificado. Método de Heum . . . . . . . . . . . 171
7.7.2 El método de Euler mejorado: punto medio . . . . . . . . . . . . . . 173
7.7.3 Métodos de Runge-Kutta de segundo orden . . . . . . . . . . . . . . 174
7.7.4 El método del punto medio modificado . . . . . . . . . . . . . . . . . 175
A-cursos/cursoMN/[Link] - August 22, 2023 v

7.7.5 El método de Runge Kutta de tercer orden . . . . . . . . . . . . . . . 176

7.7.6 El método de Runge Kutta de cuarto orden . . . . . . . . . . . . . . 176
7.8 Métodos embebidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
7.8.1 El método de Runge Kutta Fehlberg de quinto orden . . . . . . . . . 177
7.8.2 Runge Kutta. Notación vectorial . . . . . . . . . . . . . . . . . . . . 179
7.9 El método de Euler: derivación formal . . . . . . . . . . . . . . . . . . . . . 180
7.10 Errores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
7.10.1 Concepto de orden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
7.10.2 Error de truncamiento . . . . . . . . . . . . . . . . . . . . . . . . . . 183
7.10.3 Error global . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
7.11 Ajuste adaptativo del paso de integración . . . . . . . . . . . . . . . . . . . . 185
7.11.1 Estimación del error de truncamiento en métodos embebidos. . . . . . 185
7.11.2 Ajuste de paso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
7.11.3 Control adaptativo en el método de RKF . . . . . . . . . . . . . . . . 187
7.12 Duplicación de paso. REVISAR . . . . . . . . . . . . . . . . . . . . . . . . . 188
7.12.1 Tolerancias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
7.13 Propiedes de los métodos de integración de ecuaciones diferenciales. HACER 189
7.14 Librerı́as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189

8 Aplicaciones a modelos térmicos 191

8.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
8.2 La analogı́a térmico eléctrica . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
8.2.1 La resistencia convectiva y radiativa . . . . . . . . . . . . . . . . . . . 194
8.2.2 La pared simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
8.2.3 La pared compuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
8.2.4 Las leyes de Kirchhoff . . . . . . . . . . . . . . . . . . . . . . . . . . 196
8.2.5 La pared convectiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
8.3 *Calentamiento o enfriamiento por inmersión . . . . . . . . . . . . . . . . . . 198
8.4 Modelo de capacitancia concentrada. . . . . . . . . . . . . . . . . . . . . . . 203
8.4.1 Calentamiento del fluido . . . . . . . . . . . . . . . . . . . . . . . . . 204
8.4.2 Solución analı́tica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
8.4.3 Forma adimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
8.4.4 *Aplicabilidad del modelo de capacidad concentrada o decaimiento
exponencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
8.4.5 Enfriamiento por inmersión: analogı́a eléctrica . . . . . . . . . . . . . 210
vi

8.5 Enfriamiento por inmersión. Solución analı́tica . . . . . . . . . . . . . . . . . 213

8.5.1 Solución dimensional y adimensional . . . . . . . . . . . . . . . . . . 213
8.5.2 Resolución numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . 216

9 Sistemas de ecuaciones diferenciales ordinarias 219

9.1 Introducción. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219

10 Aplicaciones a los modelos térmicos 223

10.1 Introducción. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
10.1.1 Enfriamiento convectivo y calentamiento del recipiente . . . . . . . . 223
10.2 Cubos acoplados por conducción. . . . . . . . . . . . . . . . . . . . . . . . . 226
10.2.1 Tres y más cubos. Generalización . . . . . . . . . . . . . . . . . . . . 232
10.2.2 Anatomı́a de la implementación mediante arreglos . . . . . . . . . . . 233
10.3 Runge Kutta de mayor orden para sistemas . . . . . . . . . . . . . . . . . . 242
10.3.1 De nuevo dos cubos en contacto . . . . . . . . . . . . . . . . . . . . . 243
10.3.2 cubos-contacto-n3-corregido.c . . . . . . . . . . . . . . . . . . . . . . 244
10.3.3 Runge Kutta para sistemas. Notación vectorial . . . . . . . . . . . . 246
10.3.4 Conducción y generación de calor interna . . . . . . . . . . . . . . . . 254
10.4 Simulación de una aleta disipadora de calor . . . . . . . . . . . . . . . . . . 256
10.5 Una biblioteca de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 258

11 Otras aplicaciones 259

11.1 El problemas de presa y rapaz, las ecuaciones de Lotka-Volterra . . . . . . . 259

12 Ecuaciones diferenciales ordinarias de orden superior 261

12.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
12.2 Las ecuaciones del movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . 262
12.3 Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262
12.4 Reducción a un sistema de ecuaciones . . . . . . . . . . . . . . . . . . . . . . 263

13 Aplicaciones a la mecánica 265

13.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 265
13.1.1 Movimiento armónico simple . . . . . . . . . . . . . . . . . . . . . . . 265
13.2 Como un sistema de ecuaciones de primer orden . . . . . . . . . . . . . . . . 267
13.2.1 Resolución numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . 267
13.2.2 Resolución por Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 268
13.2.3 Resolución por Runge Kutta . . . . . . . . . . . . . . . . . . . . . . . 271
A-cursos/cursoMN/[Link] - August 22, 2023 vii

13.3 Movimiento amortiguado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275

13.3.1 Sobreamortiguado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275
13.3.2 Crı́ticamente amortiguado . . . . . . . . . . . . . . . . . . . . . . . . 275
13.3.3 Subamortiguado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275
13.4 Movimiento forzado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 278
13.5 Resortes acoplados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 278
13.6 El problema de Fermi, Pasta, Ulam, Tsingou (FPUT) . . . . . . . . . . . . . 282

14 Aplicaciones a la hidráulica 285

14.0.1 Descarga de un tanque de agua . . . . . . . . . . . . . . . . . . . . . 285
14.0.2 Descarga de un tanque a la atmósfera. tanqueAa . . . . . . . . . . . 285
14.0.3 Descarga de un tanque con suministro intermitente. tanqueAb . . . . 287
14.0.4 Descarga de tanques acoplados. tanquecisternaAb . . . . . . . . . . . 289
14.1 ODEs no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 291

15 sobras 293
15.1 Procedimiento de Adams y Rogers . . . . . . . . . . . . . . . . . . . . . . . 293
15.2 Una biblioteca de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 295

16 Algebra. INCIPIENTE 301

16.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 301
16.1.1 Sistemas de ecuaciones lineales . . . . . . . . . . . . . . . . . . . . . 301
16.1.2 El proceso de ortogonalización de Gram-Schmidt . . . . . . . . . . . 303
16.1.3 QR decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305

17 Métodos de discretización para ecuaciones diferenciales parciales 307

17.1 Métodos de discretización . . . . . . . . . . . . . . . . . . . . . . . . . . . . 307
17.2 Aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310

18 El método de las diferencias finitas 313

18.1 Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313
18.2 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313
18.2.1 La discretización espacial . . . . . . . . . . . . . . . . . . . . . . . . . 314
18.3 La red de discretización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315
18.4 Las ecuaciones de discretización . . . . . . . . . . . . . . . . . . . . . . . . . 317
18.4.1 Fórmulas en diferencias para la derivada primera . . . . . . . . . . . 317
18.4.2 Diferencia finita hacia atrás . . . . . . . . . . . . . . . . . . . . . . . 319
viii

18.4.3 Diferencia finita centrada . . . . . . . . . . . . . . . . . . . . . . . . . 319

18.4.4 Discretización de la derivada primera. Resumen . . . . . . . . . . . . 319
18.4.5 Fórmulas en diferencias para la derivada segunda . . . . . . . . . . . 320
18.5 Discretización de la ecuación modelo . . . . . . . . . . . . . . . . . . . . . . 320
18.6 Resolución de las ecuaciones de discretización . . . . . . . . . . . . . . . . . 323
18.6.1 Algoritmo de Thomas . . . . . . . . . . . . . . . . . . . . . . . . . . 323
18.6.2 Implementación del TDMA en C . . . . . . . . . . . . . . . . . . . . 327
18.7 Otra manera de plantera el problema . . . . . . . . . . . . . . . . . . . . . . 328

19 Aplicaciones 329
19.1 Conducción unidimensional en una varilla delgada . . . . . . . . . . . . . . . 329
19.1.1 Caso de generación nula . . . . . . . . . . . . . . . . . . . . . . . . . 330
19.1.2 Caso de generación uniforme . . . . . . . . . . . . . . . . . . . . . . . 330
19.2 Derivación de la ecuación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 331
19.3 Resolución computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
19.3.1 Caso S = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 339
19.3.2 Caso S 6= 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 340
19.3.3 Fuente no homogénea . . . . . . . . . . . . . . . . . . . . . . . . . . . 342
19.4 Propiedades no homogéneas. HACIENDO. TERMINAR . . . . . . . . . . . 342
19.4.1 Discretización A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343
19.4.2 Discretización B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344
19.5 Flujo unidimensional en conductos . . . . . . . . . . . . . . . . . . . . . . . 345
19.5.1 Perfil de velocidad en un canal . . . . . . . . . . . . . . . . . . . . . . 346

20 Diferencias finitas. Integración temporal 349

20.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 349
20.1.1 Discretización temporal . . . . . . . . . . . . . . . . . . . . . . . . . 350
20.1.2 La ecuación de discretización . . . . . . . . . . . . . . . . . . . . . . 351
20.1.3 El método explı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352
20.1.4 Difusión en un sólido semi infinito . . . . . . . . . . . . . . . . . . . . 356
20.1.5 Método totalmente implı́cito. PULIR EL PROGRAMA EJEMPLO
PARA DF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 358
20.1.6 Método de Crank Nicolson . . . . . . . . . . . . . . . . . . . . . . . . 362
20.1.7 Método θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363
20.2 Otras condiciones de borde . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364
20.2.1 Condición de borde de Neumann . . . . . . . . . . . . . . . . . . . . 364
A-cursos/cursoMN/[Link] - August 22, 2023 ix

20.2.2 Condición de borde de Robin . . . . . . . . . . . . . . . . . . . . . . 366

20.2.3 Condiciones de borde periódicas. TERMINAR . . . . . . . . . . . . . 367

21 Diferencias finitas en dos dimensiones 369

21.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 369
21.2 Diferencias finitas en dos dimensiones. . . . . . . . . . . . . . . . . . . . . . 370
21.2.1 Ecuación de discretización, forma matricial. . . . . . . . . . . . . . . 371
21.3 Métodos de solución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371
21.3.1 Algoritmos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . 373
21.3.2 Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374
21.3.3 Método de Gauss Seidel. . . . . . . . . . . . . . . . . . . . . . . . . . 383
21.3.4 Control de convergencia. Normas . . . . . . . . . . . . . . . . . . . . 388
21.3.5 SOR: Sobre y bajo Relajación Sucesiva . . . . . . . . . . . . . . . . . 391
21.3.6 Razón de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . 392
21.3.7 Mejores prácticas de implementación . . . . . . . . . . . . . . . . . . 392
21.3.8 Criterio de convergencia para Gauss Seidel. REVISAR y MEJORAR 395
21.3.9 Aceleración de Chebyshev . . . . . . . . . . . . . . . . . . . . . . . . 395
21.4 Paralelización de Gauss Seidel con Rojo-Negro. MEJORAR . . . . . . . . . 396

22 La ecuación de Helmholtz 399

22.1 Discretización por diferencias finitas . . . . . . . . . . . . . . . . . . . . . . . 400

23 Diferencias finitas en dos dimensiones. Integración temporal 403

23.1 El método explı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 406
23.1.1 Resolución numérica. Método explı́cito . . . . . . . . . . . . . . . . . 407
23.1.2 Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 407
23.1.3 Visualización con gnuplot . . . . . . . . . . . . . . . . . . . . . . . . 410
23.1.4 Otras condiciones de borde. HACER . . . . . . . . . . . . . . . . . . 411
23.2 Método totalmente implı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . 411
23.2.1 Propiedades del método implı́cito . . . . . . . . . . . . . . . . . . . . 415
23.2.2 Método Crank Nicolson. REVISAR y COMPLETAR . . . . . . . . . 416
23.3 El método θ. HACER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416
23.4 Resolución numérica con métodos implı́citos . . . . . . . . . . . . . . . . . . 417
23.5 Resolución iterativa de sistemas de ecuaciones algebraicas . . . . . . . . . . . 417
23.5.1 El método de Jacobi en la resolución de las ecuaciones de discretización419
x

23.5.2 El método de Gauss-Seidel en la resolución de las ecuaciones de dis-

cretización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 422
23.5.3 La función SOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424
23.5.4 La función gaussseidelRN . . . . . . . . . . . . . . . . . . . . . . . . 425

24 El método de los volúmenes de control 429

24.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 429
24.2 La ecuación de balance integral sobre un volumen de control . . . . . . . . . 432
24.2.1 La ecuación de conducción . . . . . . . . . . . . . . . . . . . . . . . . 434
24.3 La red de discretización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435
24.3.1 Redes de discretización uniformes . . . . . . . . . . . . . . . . . . . . 436
24.4 Obtención de la ecuación de discretización . . . . . . . . . . . . . . . . . . . 437
24.5 Conductividad en la interfaz . . . . . . . . . . . . . . . . . . . . . . . . . . . 443
24.6 Condiciones iniciales y de borde . . . . . . . . . . . . . . . . . . . . . . . . . 450
24.7 Práctica A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453
24.7.1 Condiciones de borde de Dirichlet . . . . . . . . . . . . . . . . . . . . 453
24.7.2 Condiciones de Neumann y Robin . . . . . . . . . . . . . . . . . . . . 454
24.7.3 Flujo constante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454
24.7.4 Flujo convectivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455
24.8 Práctica B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 457
24.8.1 Condiciones de borde de Dirichlet . . . . . . . . . . . . . . . . . . . . 457
24.8.2 Condiciones de Neumann y Robin . . . . . . . . . . . . . . . . . . . . 457
24.9 Práctica B. Implementación generalizada . . . . . . . . . . . . . . . . . . . . 461
24.9.1 Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 462
24.10Transitorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 462
24.11Cosmética . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463
24.12Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463
24.13Aplicaciones. Problemas de difusión . . . . . . . . . . . . . . . . . . . . . . . 464
24.13.1 El programa vara1Dt-VC.c . . . . . . . . . . . . . . . . . . . . . . . 465
24.13.2 Como modificar vara1Dt-VC.c . . . . . . . . . . . . . . . . . . . . . . 470
24.13.3 Ejemplos 1: material homogéneo, ausencia de fuente, CB de Dirichlet 472
24.13.4 Ejemplos 2. Término fuente constante . . . . . . . . . . . . . . . . . 473
24.13.5 Ejemplos 3. Material con conductividad por zonas . . . . . . . . . . . 474
24.13.6 Ejemplos 4. Borde derecho adiabático . . . . . . . . . . . . . . . . . . 475
24.13.7 Ejemplo 5. Aleta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476
A-cursos/cursoMN/[Link] - August 22, 2023 xi

24.13.8 Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476

24.13.9 Ejemplo 6. Aleta calentada en sector un extremo . . . . . . . . . . . 478
24.13.10Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 479
24.13.11Poiseuille . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 481
24.14Reglas de discretización de Patankar . . . . . . . . . . . . . . . . . . . . . . 483
24.14.1 Flujos consistentes en las caras de los volúmenes de control . . . . . . 483
24.14.2 Transmitancias (coeficientes) positivos . . . . . . . . . . . . . . . . . 484
24.14.3 Término fuente con dependencia lineales de la temperatura y pendiente
negativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485
24.14.4 Suma de los coeficiente vecinos. . . . . . . . . . . . . . . . . . . . . . 485
24.15La aleta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485
24.15.1 Teorı́a de aletas, distintas condiciones de borde . . . . . . . . . . . . 485
24.15.2 Comparación con cubos . . . . . . . . . . . . . . . . . . . . . . . . . 485

25 Volúmenes de control en dos y tres dimensiones 487

25.1 La ecuación de Laplace en 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 487
25.2 La derivación de las ecuaciones de discretización en dos y tres dimensiones . 488
25.3 Aplicación de las condiciones de borde. REVISAR-NO ES para LABORA-
TORIO 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 491
25.3.1 Condiciones de borde de Dirichlet . . . . . . . . . . . . . . . . . . . . 492
25.3.2 Condiciones de borde de Neumann . . . . . . . . . . . . . . . . . . . 493
25.4 Condiciones de borde de tipo flujo . . . . . . . . . . . . . . . . . . . . . . . . 494
25.4.1 Flujo conocido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494
25.4.2 Flujo convectivo en el borde . . . . . . . . . . . . . . . . . . . . . . . 495
25.5 El término fuente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496
25.6 Redes no uniformes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496
25.7 Resolución de las ecuaciones de discretización . . . . . . . . . . . . . . . . . 496

26 Aplicaciones a la conducción de calor 497

26.1 Geometrı́as complejas en redes de discretización cartesianas . . . . . . . . . . 497
26.2 sobras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503

27 Volúmenes de control. Integración temporal 505

27.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 505
27.1.1 Integración temporal del término de acumulación . . . . . . . . . . . 506
27.1.2 Integración temporal del término de difusión (espacial) . . . . . . . . 507
xii

27.1.3 Esténciles o patrones . . . . . . . . . . . . . . . . . . . . . . . . . . . 509

27.2 El método explı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 509
27.2.1 Criterio de estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . 509
27.2.2 Criterio de estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . 509
27.3 Crank Nicolson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 510
27.3.1 Criterio de realismo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 510
27.4 Totalmente implı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 510

28 NO CORRESPONDEAQUI: ESto tiene que ver con multimod. DEjo hasta

fin de año para no cambiar la numeración 511
28.1 Arreglo de tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . 514
28.2 Frases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515
28.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 516
28.4 Gauss Seidel y Sobrerelajación sucesiva . . . . . . . . . . . . . . . . . . . . . 517

29 Conveccion difusión 519

29.0.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 519
29.1 Derivación de la ecuación de transporte 1D . . . . . . . . . . . . . . . . . . . 522
29.1.1 Derivación a partir de la ecuación general de transporte . . . . . . . . 525
29.2 Discretización mediante volúmenes de control . . . . . . . . . . . . . . . . . 525
29.2.1 Evaluación del flujo convectivo mediante diferencias finitas centradas 527
29.2.2 Solución analı́tica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 529
29.2.3 Resolución con diferencias finitas centradas . . . . . . . . . . . . . . . 532
29.3 Discretización mediante el método upwind . . . . . . . . . . . . . . . . . . . 532
29.3.1 Génesis del método de upwind . . . . . . . . . . . . . . . . . . . . . . 532
29.3.2 Génesis del método de volúmenes de control . . . . . . . . . . . . . . 534
29.4 El esquema exponencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535
29.4.1 Esquema exponencial: ecuación de discretización . . . . . . . . . . . 536
29.4.2 Esquema hı́brido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 537
29.4.3 Comportamiento analı́tico: casos lı́mite . . . . . . . . . . . . . . . . . 538
29.4.4 Esquema hı́brido casos lı́mite . . . . . . . . . . . . . . . . . . . . . . 538
29.5 Esquema Ley de la Potencia . . . . . . . . . . . . . . . . . . . . . . . . . . . 539
29.6 Formulación generalizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 539
29.7 Ecuación de discretización generalizada: 2D . . . . . . . . . . . . . . . . . . 540
29.7.1 Difusión falsa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 541
Chapter 1

Búsqueda de raı́ces

1.1 Introducción
Las raı́ces o ceros de una función f (x) son los valores de x tales que satisfacen la ecuación
y = f (x) = 0 (1.1)
Designemos por α a la raı́z, es decir el valor de x que tal que f (x = α) = 0. Las raı́ces
o ceros de la función no son importante por el valor cero mismo. Si el mismo problema se
presenta como z = f (x) + x = x, que podemos describir con la ecuación g(x) = x, ahora
la búsqueda de ceros se transforma en la búsqueda de los valores de x que satisfagan la
igualdad, es decir en la resolución de la ecuación. Estos valores tienen la misma importancia
conceptual que los ceros o raı́ces del problema original. De manera que resolver los ceros o
las raı́ces de una función es equivalente a resolver la ecuación. Las ecuaciones de interés se
presentan aisladas o formando sistemas y describen problemas de interés en múltiples ramas
de la ciencia. En las secciones siguientes describiremos algunos de estos problema, los más
simples y nos introduciremos en los métodos para su resolución numérica, necesarios cuando
la resolución analı́tica es difı́cil o aún imposible.
Cuando podemos invertir explı́citamente la función f , es decir cuando podemos despejar x
explı́citamente en función de y de manera de tener otra función g = f −1 tal que x = f −1 (y)
o x = g(y), la raı́z se obtienen inmediatamente como α = g(0). Decimos que hemos resuelto
la ecuación en forma exacta o analı́tica. Esto no siempre es posible y depende de si la
función en cuestión es algebraica o trascendental.
Las funciones algebraicas son aquellas que se expresan mediante una suma finita de expre-
siones algebraicas, en las que aparecen solo los operadores suma, resta, producto, cociente
y potencia de números fraccionarios. Con mayor o menor dificultad pueden invertirse en
forma analı́tica. El ejemplo siguiente muestra una función fácil de invertir
x2 − 4 (x − 2)(x + 2)
y= = = x − 2 entonces x = y − 2 (1.2)
x+2 x + 2)
Las función trascendentales, en cambio, son aquellas que trascienden las reglas del
álgebra, requieren algo diferente para su cálculo e interpretación. En ellas aparecen fun-
ciones exponenciales, logarı́tmicas, trigonométricas, hiperbólicas u otras funciones especiales
de la fı́sica matemática, tales como la función error, gamma, etc., y no pueden ser repre-
sentadas por expresiones algebraicas. Estas funciones no pueden invertirse y por lo tanto

1
2

la obtención de sus ceros debe realizarse en forma numérica. Los siguientes son ejemplos
sencillos de esta clase de funciones
y = exp(−x) − x y = cos(x) y = x1/x y = log(x) (1.3)

La resolución explı́cita de problema planteado por f (x) = 0, que aparece con mucha frecuen-
cia en la resolución de problemas de matemática, ciencias e ingenierı́a, puede resultar muy
dificultosa para algunas funciones algebraicas e imposible para funciones trascendentales.
Estos casos conducen al problema de búsqueda de raı́ces que debe hacerse en forma
numérica.
El problema tiene antecedentes tempranos en la historia de la ciencia, aunque en sus primeras
manifestaciones, aun traducidos a términos modernos, tal vez no los reconocerı́amos como
problemas de búsqueda de raı́ces. Uno de estos problemas fue la obtención de la raı́z cuadrada
de un número a, problema que puede expresarse como la resolución de la ecuación f (x) =
x2 − a = 0. Este problema apareció ya en tiempos babilónicos ([Link]
org/wiki/Methods_of_computing_square_roots), y llegó a nuestros dı́as tratado por Isaac
Newton (1642-1726) y por Joseph Raphson (c1646-c1715). Con sus nombres reconocemos
hoy a uno de los métodos más eficaces para la búsqueda de raı́ces y que estudiaremos en la
sección 1.10 en la forma que le dió Thomas Simpson (1710-1760) ya basada en el cálculo.

1.1.1 Ejemplo: un problema de Mecánica

Para mostrar un ejemplo de como surge un problema de búsqueda de raı́ces en el ámbito
de la Mecánica consideremos el problema de la caı́da libre con fricción1 El estudiante de
fı́sica encontrará un excelente tratamiento de este problema en el libro Classical Mechanics
de John R. Taylor, [15], y en el artı́culo de Ooi, [12]. Como todo problema de dinámica está
gobernado por la ley de Newton
F = ma (1.4)
donde a es la aceleración del cuerpo, m es su masa, F es la fuerza que actúa sobre aquel2 .
En el caso de la caı́da libre, las fuerzas actuantes son el peso del cuerpo, mg, donde g es la
aceleración de la gravedad, que tiene en este caso la dirección del movimiento, y la fuerza
de arrastre que actúa en la dirección opuesta y que podemos considerar proporcional a una
potencia de la velocidad v. La suma de ambas fuerzas es
ma = F = mg − kv n (1.5)
La escribimos de nuevo anotando el significado de ambos términos
ma = F = mg − kv n
|{z} (1.6)
|{z}
fuerza gravitacional fuerza de arrastre
Lo hicimos para recalcar que el estudiante debe aprender a reconocer los distintos términos
de una ecuación y los efectos que cada una de ellos produce.
1
Estrictamente, al menos para algunos textos, la caı́da libre se refiere a la caı́da de un cuerpo sólo sometido
a la fuerza de gravedad, por eso empezamos diciendo caı́da libre con fricción para referirnos al caso que nos
ocupa, pero establecido el tema, usaremos caı́da libre para acortar.
2
F , a, g, v son, por supuesto, vectores, Para los signos de sus magnitudes consideramos un sistema
coordenado con origen en el punto de lanzamiento, orientado en la dirección del movimiento. En dicho
sistema los vectores son colineales y prescindiremos de la notación vectorial
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 3

El término de arrastre kv n se debe a la resistencia que el fluido ofrece al desplazamiento del

objeto. El coeficiente k se denomina coeficiente de arrastre. El exponente n define una
dependencia potencial con la velocidad que depende del factor dominante en la fuerza de
arrastre: la inercia o la fricción. En el primer caso la fuerza depende del área transversal del
cuerpo y la relación con la velocidad es cuadrática, es decir k = 2. En el segundo la fuerza
de fricción depende del área de la superficie expuesta al rozamiento con el aire y el arrastre
es lineal, es decir k = 1. Este último caso ocurre para objetos de tamaño normal3 en fluidos
como el agua o el aire, y además es el más sencillo de resolver analı́ticamente, por eso lo
utilizaremos como ejemplo en esta sección.
Antes de seguir adelante hagamos incapı́e en notar que la aceleración y la velocidad del objeto
que cae son funciones del tiempo, entonces podemos escribir a(t) y v(t), y por supuesto,
la fuerza también los será, entonces también podremos escribir F (t). Entonces, como la
aceleración es la derivada temporal de la velocidad, v(t), a = dv(t)/dt, la ecuación que
describe el movimiento es una ecuación diferencial ordinaria ( ODE en inglés o EDO) de
primer orden.
dv(t)
m = mg − kv(t)n (1.7)
dt
que podemos reescribir, para n = 1, como una ecuación diferencial ordinaria lineal
dv k
=g− v (1.8)
dt m
y cuya solución analı́tica es fácil de obtener. Suponiendo que la condición inicial, en el
instante t = 0 es v(t = 0) = v0 se tiene que la solución resulta
gm

v(t) = v0 + 1 − e−kt/m (1.9)
k
Conocidos los parámetros contantesg, m, k, y v0 es fácil obtener la velocidad en función del
tiempo. Por supuesto, la ecuación 1.7 y otras más complejas se puede resolver numéricamente,
pero dejaremos el asunto para el capı́tulo 7. Por ahora nos interesaremos en otra clase de
problema.
Supongamos ahora que se trata del diseño de un objeto que será lanzado desde cierta
altura y que se espera que alcance una cierta velocidad terminal más o menos rápidamente.
El parámetro que controla el diseño es el valor del coeficiente de fricción k. El problema
puede plantearse entonces en los siguientes términos: conociendo la masa del objeto m y la
aceleración de la gravedad g, se quiere saber el valor del coeficiente de fricción k tal que,
cuando el objeto es lanzado con velocidad inicial v0 = 0, alcance la velocidad predeterminada
vf transcurrido cierto tiempo tf desde el lanzamiento.
Podemos reescribir el problema dado por la ecuación 1.9 como un problema de búsqueda de
raı́ces de la siguiente manera
gm

1 − e−ktf /m − vf = 0 (1.10)
k
el análisis de esta función muestra que no es posible despejar el valor incógnita k, es decir,
tenemos una función implı́cita de k de la forma

f (k) = 0 (1.11)
3
Objetos no muy pequeños.
4

donde el coeficiente k es la raı́z buscada de la función y la solución de nuestro problema.

Para resolver este tipo de problemas recurriremos a métodos numéricos especı́ficos para la
búsqueda de raı́ces.
El problema asi enunciado, como veremos en seguida, es fáci de resolver. No obstante, la
simplicidad aparente de este problema fı́sico no debe engañar al lector. El problema de
la caı́da libre en un fluido es un problema de gran importancia técnica y gran actualidad.
Cuando escribı́ por primera vez esta sección, una de las aplicaciones más resonantes fué el
uso de paracaı́das para aminorar la velocidad de entrada de los vehı́culos espaciales Viking
y Pathfinder en la atmósfera de Marte. Paul Whiters, [16], utiliza conceptos de fı́sica básica
para ampliar la visión de estos problemas.
En el año 2021, se ha hecho conocida entre nosotros la cientı́fica argentina Clara O´Farrell,
aerodinamicista, quien trabajó en la Nasa en el diseño y desarrollo de paracaı́das que inter-
vendrı́a en el descenso del Perseverance en Marte. En sus palabras, en una charla TedXRi-
oDeLaPlata, ella decı́a Mi especialidad es la aerodinámica y fui parte del equipo que diseñó
y probó el paracaı́das supersónico utilizado para el aterrizaje de Perseverance. ...Cuando se
encontró con la atmósfera de Marte, Perseverance viajaba a 5.5 kilómetros por segundo. Por
lo tanto, tenı́amos 7 minutos para desacelerar y aterrizar. Para eso usamos el paracaı́das
que desplegamos cuando viajaba a casi dos veces esa velocidad. De él cuelga Percy. Imaginen
estar en mi lugar, en el centro de control el dı́a del aterrizaje. Imaginen la adrenalina y la
sensación cuando vi los videos y, por primera vez, el suelo de Marte. ¿No es de pelı́cula?,
[11].

1.2 Raı́ces. Fundamentos matemáticos

Sea una función f (x), continua de valores reales en el intervalo cerrado [a, b] y tal que

f (a)f (b) < 0 (1.12)

Esto implica que los valors de f (a) y f (b) deben tener distinto signo y la gráfica de la función
f (x) debe cortar el eje x un número impar de veces o por lo menos una vez.
El teorema del valor medio o teorema de Bolzano, establece que debe haber por lo
menos una raı́z, cero o solución,que designaremos α, en el intervalo abierto (a, b), tal que
f (α) = 0, es decir
Si f (a)f (b) < 0 ∃ α ∈ (a, b) : f (α) = 0 (1.13)

La raı́z es única en el intervalo [a, b] si la derivada, f ′ (x) ∈ [a, b] existe y preserva su signo.
Desde el punto de vista práctico, cuando hay más de una raı́z, el proceso de separación de
raı́ces implica determinar el signo de la derivada de la función en puntos de una partición
del intervalo [a, b] suficientemente fina.
Como ejemplo consideremos la función seno entre [π/2, 3π/2] mostrada en la Figura 1.1. Es
evidente allı́ (siendo f (3π/2) = −1 y f (π/2) = 1) que la condición f (a)f (b) < 0 garantiza la
existencia de al menos un cero, o un número impar de ellos, entre los extremos del intervalo.
En el intervalo señalado la función es siempre negativa, por lo que solo habrá una sola raı́z.
Para encontrar las raı́ces cuando no existe una solución explı́cita o cerrada para f (x) = 0
(por ejemplo, cuando la función es trascendental), debe recurrirse a métodos numéricos
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 5

f (x) = sin(x) f (x) = sin(x)

1 1
0.8 f(x) f(x)
0.6 0.5
0.4
0.2
0 0
-0.2
-0.4 -0.5
-0.6
-0.8
-1 -1
π/2 π 3π/2 π/2 π 3π/2

Figure 1.1: La función seno en el intervalo [π/2, 3π/2].

especı́ficos que por lo general caen en la categorı́a de métodos iterativos. Estos métodos
se caracterizan por ir aproximando la solución en una sucesión de pasos repetitivos.
Presentaremos aquı́ algunos de ellos. La búsqueda iterativa de raı́ces, si bien es conceptual-
mente muy simple, puede complicarse en muchos casos (raı́ces múltiples, etc.). Si se enfrenta
con estos problemas, el estudiante deberá recurrir a la bibliografı́a pertinente ([4], [5], [6],
[7], [13]) y deberá hacer uso de su ingenio para resolverlos.
6

1.3 Fı́sica de la caı́da libre con fricción. OMITIR. EN

ELABORACION.
El problema de la caı́da libre fue analizado desde un punto de vista filosófico por Aristóteles
(384-322 AC) en el contexto de las causas del movimiento de los cuerpos. Le atribuı́a al medio
la propiedad de resistirse al paso de un cuerpo extraño, propiedad que llamaba tenuidad y que
relacionaba en forma imprecisa con la densidad. En el aire, un medio tenue, el movimiento
era más fácil, mientras que en el agua, un medio más denso, un objeto adquirirı́a menor ve-
locidad. La velocidad adquirida por un objeto era inversamente proporcional a la tenuidad.
Aristóteles sabı́a que para poner en movimiento un objeto se requerı́a una fuerza. Pero
consideraba que para que el movimiento persitiera, luego de removido el primer motor, nece-
sitaba del medio para mantenerse ası́. El medio, mediante ondas, vibraciones o remolinos,
que el propio cuerpo introducı́a sostenı́a el movimiento. Esto estaba en contraposición a
sus ideas sobre la tenuidad, porque en el caso lı́mite del vacı́o, un objeto se detendrı́a en su
movimiento
Las ideas aristotélicas recién fueron rebatidas en la edad media por el clérigo francés Jean
Buridan (1300-1358), quien introdujo el concepto de inercia, que hoy identificamos con la
resistencia de un cuerpo a cambiar el movimiento. El concepto de inercia no requiere de un
medio para soportar el movimiento. Buridan consideraba que mientras el ı́mpetus (cantidad
de movimiento en términos modernos) fuera mayor que la resistencia, el cuerpo seguirı́a
moviéndose. Hoy entendemos que una resistencia, una fuerza externa, cambiarı́a la cantidad
de movimiento y dentendrı́a finalmente el objeto.
Galileo (1564-1642) en su estudio de las leyes del movimiento, se interesó también por el
efecto de la resistencia del aire, o tal vez más correctamente, en el efecto de la ausencia de
un medio. Usó péndulos4 para medir el efecto de la resistencia del aire y concluyó que esta
dependı́a de la velocidad del cuerpo.
Las leyes de la mécánica recién se establecieron por completo luego de que Newton las
enunciara, lo que dió lugar a que el efecto de la resistencia del aire en el movimiento se pudiera
estudiar desde el punto de vista de una teorı́a completa y coherente. Hoy aceptamos las tres
leyes de Newton como suficientemente adecuadas como modelo mecánico para describir el
movimeinto de los cuerpos. Distingimos la segunda ley como ley fundamental de conservación
y ecuación gobernante del movimiento, y distinguimos otras leyes particulares que describen,
segun sus causas, los distintos tipos de fuerzas que se ejercen sobre los objetos. La segunda
ley de Newton que describe la ecuación de movimiento se escribe

ma = mv̇ = mr̈ = F (1.14)

donde F es la fuerza que causa la aceleración del objeto. En nuestro caso, un objeto en caı́da
libre, o aún en el caso más general de un objeto en vuelo libre sobre la superficie de la tierra,
está sometido a la fuerza de la gravedad y a la fuerza de arrastre o de resistencia que ejerce
el medio en que se mueve.
F = mg + f(v)v̂ (1.15)

El problema del arrastre tomó importancia con el comienzo de la aeronáutica, y cientı́ficos,

4
Recientemente, Marı́n Mersenne ([? ]), puso en duda que estas conclusiones pudieran obtenerse de los
experimentos descriptos por Galileo
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 7

ingenieros y pilotos se ocuparon del tema en las primeras decadas del siglo XIX. Se estudió
en forma experimental mediante el uso de péndulos, caı́da libre, y tuneles de viento.
Eiffel, desde un laboratorio montado en la famosa torre que lleva su nombre (1905-1911)
para estudiar el de la resistencia al movimiento que ejerce el aire.
No obstante, no se dispuso una teorı́a cientı́fica que explicara el fenómeno hasta que Prandlt
en 1904 propuso su teorı́a de la capa lı́mite, que comenzó con el desarrollo de la Mecánica
de Fluidos como ciencia moderna. En este ámbito, las fuerzas que el flujo de un fluido ejerce
al objeto (o la que el objeto ejerce al fluido, opuesta a la primera) se denomina fuerza de
arrastre, y su estudio es importante en un gran número de asuntos prácticos. Piense en
cualquier vehı́culo de transporte, desde una bicicleta, su automóvil, un camión o el tren o
incluso un cohete espacial en su vuelo ascencional o en su regreso a tierra, todos ellos, son
objetos que se mueven en un fluido, y mucha de la potencia que empleamos en esos vehı́culos
se usa para contrarrestar el efecto de la perdidad de energı́a ocasionada por el arrastre.

Figure 1.2: Para flujo laminar la separación se produce en algún punto en la parte frontal
de la esfera y la estela es grande. La sección transversal de la capa lı́mite separada aumenta
la sección transversal de la esfera, aumentando ası́ la diferencia de presión causada por el
objeto. Para flujo turbulento, la capa lı́mite permanece adherida a la esfera y la separación
se produce bastante atrás, y la estela es delgada. Su contribución a la fuerza de arrastre de
presión es menor comparada con la fricción. A medida que el número de Reynolds crece, la
separación se produce más atrás y el coeficiente de arrastre disminuye.

El arrastre tiene origen en dos fenómenos, por una lado en la diferencia de presión que se
produce entre barlovento (parte de un objeto de cara el viento) y sotavento (parte de un
objeto a espaldas del viento) de un cuerpo. Este efecto es predominante en los objetos
romos y es áltamente dependiente de la geometrı́a del cuerpo, sobre todo el área de sección
transversal normal a la dirección del flujo medio. Por otro lado, en la fuerza de fricción que
se produce por efecto de la viscosidad en la capa lı́mite superficial del cuerpo, y por ende,
dependiente de su área.
8

1.3.1 REVISAR
Consideremos las ecuaciones gobernantes del movimiento de un cuerpo romo, tal como una
esfera, por ejemplo, que por simplicidad se mueve en una dimensuón. La segunda ley de
Newton establece
d2 x
= mF (1.16)
dt2
Cualquiera sea la forma funcional de la fuerza, se puede como una serie de potencias. Por
ahora consideremos solo los términos de segundo orden

F = av + bv 2 (1.17)

En el caso de una esfera, para número de Reynolds moderado, tiene lugar un régimen de
capa lı́mite, que si bien se produce por efecto de la viscosidad, su contribución a la fuerza
de arrastre es menor. No obstante, en este régimen ocurre la denominada separación de la
capa lı́mite, que a los efectos de la presión, modifica la forma del cuerpo y en consecuencia
la fuerza de arrastre de forma.
Para este régimen, se puede suponer que la viscosidad, y en consecuencia el número de
Reynolds, no tiene efecto en el fenómeno. Existen dos formas de predecir de que depende la
fuerza de arrastre

• Aplicando la ecuación de Bernoullı́ a una lı́nea de corriente que termine en el punto de

estancamiento
1 2
ρv = ∆p (1.18)
2
y como despreciando los efectos viscosos FD = ∆pA resulta que

FD = CρAv 2 (1.19)

donde C es un coeficiente experimental.

• Mediante un análisis dimensional. La única forma en que las variables FD , ρ, A y v

pueden estar relacionadas entre si en forma dimensionalmente correcta, es mediante
una relación como la precedente
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 9

El exponente n depende del número de Reynolds, el número adimensional mas importante

de la mecánica de fluidos. Este número caracteriza el régimen de flujo, y esta dado por
Dv
Re = (1.20)
ν
donde D es una longitud caracterı́stica del la región de flujo, v es velocidad y ν la viscosidad
cinemática. En el caso de una esfera que cae con velocidad v en un fluido estanco, o en el
caso que fluido fluye con velocidad v alrededor de una esfera estacionaria, dos situaciones
semejantes5 , la longitud caracterı́stica es el diámetro de la esfera D.
Esta configuración de flujo da origen a numerosos regı́menes de flujo que ocurren a dis-
tinto número de Reynolds. Solo consideraremos los regı́menes que ocurren en los tres casos
siguientes

• Para Re → 0 se tiene el denominado flujo reptante o flujo de Stokes.

• Para Re más altos, pero no tanto como para desencadenar el régimen turbulento.

• Para Re mayores a 10, es decir en el régimen turbulento.

Flujo reptante

Stokes, [14], demostró que para Re → 0, situación que se obtiene para objetos muy pequeños
o para fluidos muy viscosos, la fuerza de arrastre es proporcional a la velocidad

FD = 3πµDv (1.21)

es decir k = 3πµD y n = 1.
La situación tiene aplicación en el meteorologı́a, en el análisis de la velocidad de las gotas
de lluvia, en geologı́a, en el análisis de la dispersión de polvo y cenizas volcánicas, en los
problemas de sedimentación, en biologı́a en el estudio de organismos acuáticos y aéreos
pequeños, etc., todas ellas situaciones en las que se trata de partı́culas muy pequeñas.

Para régimen laminar. COMPLETAR

Cuando el número de Reynolds es grande pero no mayor a 1,

5
Transformación de Galileo
10

1.4 Métodos numéricos para la búsqueda de raı́ces

Consideraremos los métodos básicos de búsqueda de raı́ces. Estos tienen importancia propia,
pero además, son un ejemplo de procedimientos iterativos que aparecen con frecuencia en la
solución de problemas cientı́ficos.
En los métodos iterativos, por oposición a los métodos directos, se busca una solución como
lı́mite de una secuencia de valores cada vez más cercanos a la solución.
Los algoritmos de búsqueda de raı́ces se pueden clasificar en dos categorı́as:

• Encorchetado de la raı́z. Proceden achicando progresivamente un entorno que

contiene a la raı́z. Tienen convergencia garantizada.

• Mejoramiento de la raı́z. Funcionan mejorando progresivamente un valor de prueba.

Son más rápidos pero no siempre convergen.

Algunos métodos utilizan ambas ideas simultáneamente. En las secciones siguientes pre-
sentaremos los siguientes:

• Método de encorchetado. En ellos, se busca la raı́z en una serie de intervalos, que

siempre la contienen, pero que se hacen cada vez menores en cada iteración. Se basan
en el teorema del valor medio y el método prototı́pico es el método de la bisección
o búsqueda binaria.

• Newton Raphson. Este método se elige un valor de x0 denominado semilla, y se

reemplaza la función por la por la recta tangente a ella en f (x0 ), es decir se reemplaza
la función por una función lineal que aproxima la función original en su cercanı́a. Luego
se obtiene la raı́z de esta función lineal y se utiliza este valor como nueva semilla para
comenzar el ciclo de nuevo.

• Método de la secante. Usa una aproximación parecida a la del método de Newton,

pero usa la secante en vez de la tangente.

• Método de la falsa posición. Este método asegura la raı́z dentro de un intervalo

cerrado, aunque su lógica deriva del método de la secante.

• Otros métodos. El método de Ridders se basa en el método de la falsa posición y

utiliza una exponencial, el método de Muller es similar al método de la secante pero
usando tres puntos y pasa por ellos una parábola, el método de Brent combina el
método de las secante y el de la bisección.

• Métodos de iterativos de punto fijo o aproximaciones sucesivas.

1.5 Método de la bisección

Consideremos el problema de encontrar una raı́z de y = f (x) en el intervalo cerrado [a, b].
Si se cumple que f (a)f (b) < 0, quiere decir que uno de los lı́mites a o b es negativo y el otro
positivo, entonces como la función es continua, pasa por cero por lo menos una vez.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 11

En la búsqueda numérica de la raı́z, el intervalo [a, b] debe elegirse de manera que contenga
sólo una raı́z. Si hubiera dos raı́ces, o cualquier otro número par de ellas, el algoritmo no
detectará ninguna raı́z en el intervalo. Si aún se aplica, el algoritmo encontrará una raı́z,
pero no sabremos cuantas de ellas hay en el intervalo y cual de ellas encuentra. Si hubiera
un número impar de raı́ces, el método de la bisección convergerá siempre a alguna de ellas,
pero no sabremos a cual. Por eso es importante haber visualizado el comportamiento de la
función antes de comenzar la busqueda de raı́ces.
El método de la bisección o búsqueda binaria es el método de encorchetado, bracketing,
más conocido. Consiste en dividir el entorno [a, b] por la mitad. La abscisa que corresponde
al promedio de a y b se denomina c y será nuestra primera aproximación a la raı́z,
a+b
c= (1.22)
2
El valor de c determina dos nuevos intervalos [a, c] y [c, b] y la raı́z α estará incluida en
alguno de ellos, como se muestra en la figura 1.3. La distancia o longitud entre los lı́mites
de cualquiera de estos intervalos, será siempre mayor que la distancia |c − α| entre la raı́z y
c.
|b − c|
|α − c| |α − c| < |b − c| < |b − a|

a c α b
Figure 1.3: Bisección del intervalo de búsqueda [a, b]

Podemos usar la longitud de los segmentos definidos por los intervalos [a, c] y [c, b], que son
respectivamente |c − a| y |b − c|, para acotar la diferencia entre la raı́z α y el valor c, |c − α|,
de la siguiente manera
|c − α| ≤ |b − c| = |c − a| (1.23)
siendo |b − c| o |c − a| una cota superior. Si esta cota es suficientemente pequeña, tomamos
el valor de c como una aproximación al valor de la raı́z y hacemos c ≈ α. La diferencia
entre el valor aproximado y el valor de la raı́z es el error de la solución numérica. En el
caso del método de la bisección no es posible conocer el error, pero como hemos visto si es
posible acotarlo, es decir es posible determina si es menor que una cantidad determinada
suficientemente pequeña.
Cuán pequeño es suficientemente pequeño depende de la precisión requerida para el resul-
tado. La precisión nos indica cuán cercano está el valor aproximado del verdadero valor de
la raı́z. El valor mı́nimo de la precisión aceptable se denomina tolerancia y se designa con
ε. Es un número pequeño, arbitrario, predeterminado, que da una medida de la pequeñez
que debe tener la cota del error para que nos quedemos satisfechos con la solución.
Entonces, exigiremos que la distancia |b − c| o |c − a| sea menor que cierta tolerancia ε, un
número pequeño, arbitrario.
En concreto, c será una aproximación de α en el intervalo [a, b] si
|c − α| ≤ |b − c| = |c − a| < ε (1.24)
Si en caso contrario, no se satisface el criterio |b − c| = |c − a| < ε, podemos tomar como
nuevo intervalo de búsqueda al segmento [a, c] o [c, b]. De entre los dos elegimos aquel que
12

contenga o encorchete la raı́z, es decir, para el cual se cumple la condición que f (a)f (c) < 0
o f (c)f (b) < 0. El procedimiento se repite obteniéndose una sucesion de valores {ci } =
{c0 , c1 , c2 , ..., cN } hasta que aguno de ellos satisface la tolerancia establecida. El algoritmo
se puede describir de la siguiente manera.

Método de la bisección: algoritmo

Algoritmo

1. Datos. Definir la función f (x), los valores de a, b y la tolerancia ε.

2. Determinar si existe al menos una raı́z en el intervalo: f (a)f (b) < 0. Si no la

hay, producir un mensaje y terminar.

3. Calcular el valor de c con

a+b
c= (1.25)
2
4. Si |b − c| ≤ ε, tomar c como raı́z, imprimir, terminar.

5. Definir el próximo intervalo de búsqueda. Si f (c)f (b) < 0 hacer a = c, en caso

contrario hacer b = c

6. Volver al paso 2.

El algoritmo producirá una serie de valores {ci } que como hemos visto aproximan la raı́z
cada vez con mayor precisión. Decimos que la sucesión o serie de valores convege a la raı́z.

1.5.1 Ventajas y desventajas

El método de la bisección tiene las siguientes ventajas

1. Es muy simple, fácil de programa, sólo requiere calcular el promedio de dos números.
2. Al ser un método de encorchetado siempre converge.
3. Si hay más de una raı́z, siempre converge a una.
4. Se puede calcular rigurosamente una cota para el error.
5. La razón de convergencia es predecible.
6. A mayor número de iteraciones mayor presición.

Presenta las siguientes desventajas

1. Necesita un intervalo de búsqueda que contenga una sola raı́z.

2. La convergencia es lenta (razón de convergencia lineal).
3. No tiene en cuenta los lı́mites de precisión de la computadora.
4. Si la solución tiene polos o singularidades, la solución puede converger a ellos.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 13

1.6 Convergencia
Los métodos iterativos producen una sucesión o un serie de resultados que es de esperar
que se acerquen progresivamente a la solución buscada. Esto no siempre es ası́. Si la serie
tiende a un valor lı́mite decimos que la serie converge. En caso contrario decimos que
la serie diverge. El cómputo de los procesos iterativos divergentes termina en valores no
computables (overflow o NaN, Not a Number en C) y, de no tomar medidas preventivas, en
la terminación del programa.
Para algunos algoritmos, la convergencia está siempre garantizada, para otros, ésta depende
de ciertas condiciones que varı́an con el algoritmo. Una propiedad importante de los algorit-
mos iterativos de búsqueda de raı́ces es la rapidez de la convergencia, que puede determinarse
analı́ticamente.
Antes de pasar a los detalles de la convergencia en el método de la bisección, demos primero
una mirada práctica, experimental, del asunto. En la figura 1.4 se muestra gráficamente, para
el ejemplo estudiado, los valores aproximados de la raı́z en función del número de iteración.
La tabla de datos para el gráfico se produce en la lı́nea 41 del programa biseccion-1.c.
Se observa que las primeras aproximaciones a la raı́z difieren de ella por más o por menos.

f (x) = x6 − x − 1
1.26
dat : 1:2
1.24

1.22

1.2

1.18

1.16

1.14

1.12
0 2 4 6 8 10 12 14 16 18 20

Figure 1.4: Convergencia hacia la raı́z.

La diferencia, grande al principio, decrece paulatinamente hasta que en el gráfico parece

alcanzar un valor constante. Este comportamiento es evidencia de que el algoritmo de la
bisección, aplicado a esta función, está convergiendo a la raı́z.

1.6.1 Convergencia del método de la bisección

Los métodos iterativos en general producen una secuencia que esperamos converja a la
solución, tenemos que
α = limn→∞ cn (1.26)
14

En la práctica, solo podemos alcanzar un valor aproximado. ¿Hasta cuándo es aconsejable

seguir con la iteración? o ¿Cuándo es aconsejable detenerla?
Llamamos tolerancia a un valor ǫ, pequeño, que será mayor que el error que admitimos en
la solución
|f (cn ) − f (α)| < ǫ (1.27)
como f (α) = 0, solo tenemos que chequear que |f (cn )| ≈ 0 con cierta tolerancia.
En el caso particular del método de la bisección es posible determinar de antemano, con
independencia de la función, cuantas iteraciones son necesarias para que nuestra solución se
obtenga dentro de una tolerancia al error determinada.
Hemos visto que el método mantiene siempre encorchetada la raı́z. También vimos que
podemos usar la longitud del intervalo, del corchete, como cota del error. Si queremos que
esa cota sea menor que cierta tolerancia, ǫ, podemos escribir
|cn − α| < |bn − cn | = |an − cn | < |(b − a)n | < ǫ (1.28)
Ahora, siendo que que la longitud del intervalo de búsqueda se reduce a la mitad en cada
iteración tenemos que
|(b − a)n |
|(b − a)n+1 | = (1.29)
2
Con cada iteración obtendremos los valores de la sucesión para la longitud del intervalo, que
escrita en términos del intervalo inicial será
|(b − a)| |(b − a)| |(b − a)|
|(b − a)1 | = , |(b − a)2 | = , |(b − a)3 | = , ... (1.30)
2 4 8
generalizado para n iteraciones se obtiene
|(b − a)|
|(b − a)n | = (1.31)
2n
Para obtener una determinada precisión será suficiente hacer
|(b − a)|
|cn − α| ≤ <ǫ (1.32)
2n
de donde se obtiene el número de iteraciones necesarias. El requisito para una determinada
precisión puede escribirse como
(b − a)
< 2n (1.33)
ǫ
y tomando logaritmos en ambos miembros y teniendo en cuenta que loga xn = nloga x, resulta
que el número de iteraciones necesario para obtener una determinada exactitud es
 
ln(|b − a|/ǫ) |b − a|incial
n> o n = log2 (1.34)
ln2 ǫ

1.7 Error
Al obtener una solución aproximada de un problema es importante estimar el error que se
comete, que se puede expresar como error absoluto o error relativo. El error absoluto
es la diferencia entre el valor aproximado y el valor de la raı́z.
ea = |αk − α| (1.35)
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 15

mientras que el error relativo se define como sigue

|αk − α|
er = (1.36)
|α|
En la figura 1.5 se muestra en un gráfico logarı́tmico como desciende el error hasta alcanzar
la precisión establecida. Para el presente caso se observa que aproximadamente cada 5
iteraciones se baja el error absoluto en dos órdenes de magnitud.
El análisis de la figura 1.4 y 1.5 muestra que, si se grafican los resultados en una escala
adecuada (que represente bien el intervalo de búsqueda de la solución, ası́ como los valores
de esta en ese intervalo), después de 15 iteraciones no puede apreciarse cambios en la solución.
Esto se obtiene ya para una precisión de 105 . 6

1
♦ bidat : 1:4 ♦
0.1 ♦
♦
♦
0.01 ♦
♦
♦
♦
0.001 ♦
error ♦
♦
0.0001 ♦
♦
♦
1e-05 ♦
♦
♦
♦
1e-06 ♦
♦
♦
1e-07 ♦
0 5 10 15 20 25
iteración
Figure 1.5: Evolución de error

1.7.1 Exactitud y precisión.

En ciencia e ingenierı́a, los términos exactitud (accuracy) y precisión (presition) tienen sig-
nificados bien definidos que es necesario aclarar.
En los métodos numéricos, la precisión se refiere a cuan bién la implementación computa-
cional representa al número que quiere representar (Hoffman, pp. 4). En C, por ejemplo,
podemos aumentar la precisión de un número utilizando distintas clase de números de punto
flotante. La precisión con que se representan los números reales no es homogénea.
Si hablamos de una medida experimental, por ejemplo, exactitud se refiere a la diferencia
entre el valor real, verdadero, de la medida y el que podemos obtener con nuestro instrumento,
que puede incluir errores sistemáticos y observacionales. El valor verdadero de la medida
puede definirse en términos estadı́sticos en los que no incursionaremos aquı́.
En el caso de los métodos numéricos que nos ocupan, la exactitud se refiere a la diferencia
entre el verdadero valor de una solución, que en algunos pocos casos podemos obtener en
6
Ver Atkinson, ENA pp 64, o NR
16

forma analı́tica, y el valor aproximado que podemos calcular con nuestro algoritmo numérico.
En los métodos numéricos, la exactitud está medida por el error.

f (x) = x6 − x − 1 f (x) = x6 − x − 1
1.5 2
a a
1.4 c 1.8 c
1.3 b 1.6 b
1.2 1.4
1.1 1.2
1 1
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20

Figure 1.6: Evolución del entorno de búsqueda y de la aproximación de la raı́z: a) bisección,

b) Régula Falsi. Observe que en el segundo caso, el algoritmo trabaja como la secante, sin
cambio de punto de referencia.

Razón de convergencia

El análisis numérico puede determinar la velocidad con los distintos métodos proceden hacia
la convergencia. De ello surge una clasificación que los distingue en métodos de convergencia
lineal, cuadrática y superlineal.
Una secuencia {cn }∞
n ≥ 0 se dice que converge con orden p ≥ 1 a un punto α si se cumple
que
|α − cx+1 | ≥ K|α − cx+1 |n con n ≥ 0 y K > 0. (1.37)
otra forma de escribirlo es
|α − cx+1 |
limn7→∞ = =K (1.38)
|α − cx+1 |n
Si p = 1 se habla convergencia lineal y K, la razón de convergencia debe ser menor que
uno.
Si p > 1 se habla convergencia superlineal. Para p = 2 y p = 3 se habla convergencia
cuadrática y cúbica.
Para convergencia lineal se puede inducir que

|α − cx+1 | ≥ K n |α − c0 | n≥0 (1.39)

resultado que hemos usado más arriba.

Criterios prácticos de terminación

Existen varios criterios de terminación para un algoritmo iterativo. Estos pueden tener
convergencia lenta, no obstante por lo general se preestablece un máximo de iteraciones,
alcanzado el cual el algoritmo se detiene. Si la terminación se hizo mediante este criterio
debiera analizarse si los resultados obtenidos son suficientemente convergidos.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 17

Un criterio más acorde con la idea de convergencia consiste en detener la iteración cuando la
diferencia entre dos resultados sucesivos resultan suficientemente pequeña, es decir |cn+1 −
cn )| < ǫ. Este criterio puede ser engañoso si la convergencia es muy lenta o cuando la raı́z
es cercana a cero. La forma relativa de este criterio es
|cn+1 − cn |
< ǫr (1.40)
|cn+1 |

Esta forma se prefiere cuando no se conoce los valores que tomarán las raı́ces o la función. Los
valores de tolerancia recomendable depende de los valores de la raı́z. Se pueden especificar en
forma absoluta o relativa. Si la raı́z esta cercana a cero, la forma relativa no es conveniente.
Si la raı́z es de orden 1, una tolerancia absoluta de ǫ = 10−6 será aceptable. Si la raı́z ocurre
para valores más altos, la tolerancia deberá ser menor, en cuyo caso conviene el uso de un
criterio relativo.
Teniendo en cuenta que las raı́ces hacen cero la función, aún también podrı́a usarse como
criterio de convegencia |f (xn+1 )| < ǫ. La también tolerancia puede elegirse en relación a la
precisión de la máquina, véase por ejemplo la sección Bisection Methods en el libro Numerical
Recipies, [13].

1.7.2 Ejemplo: raı́z de f (x) = x6 − x − 1

Consideremos de la función
f (x) = x6 − x − 1 (1.41)
en el intervalo [1, 2], y cuya raı́z, en ese intervalo es α = 1.13472413840152.
Es fácil explorar la función con gnuplot de la siguiente manera

f(x)=x*x*x*x*x*x-x-1
set xrange[-10:10]
set yrange[-10000:10000]
plot f(x)

para un primer intento se puede omitir la especificación del rango. La figura 1.7, realizada
con gnuplot, muestra el comportamiento general de la función, y en la cercanı́a de la raı́z.

f (x) = x6 − x − 1 f (x) = x6 − x − 1
10000 1
f(x) f(x)
5000 0.5
0 0
-5000 -0.5
-10000 -1
-10 -5 0 5 10 1 1.04 1.08 1.12 1.16 1.2

Figure 1.7: La función f (x) = x6 − x − 1 en los intervalos de x, [−10, 10] y [1, 1.2]
18

1.8 Resolución computacional. C

Resolución computacional

Por supuesto, existe varias maneras de implementar un algoritmo en un lenguaje de progra-

mación. En esta y en las siguintes secciones, mostraremos algunas variantes.
El algoritmo presentado en la sección 1.5 se ha programado en el código C biseccion-1.c
reproducido en el listado 1.1. La función dada por la ecuación 1.41 se ha programado como
una función de C, cuya definición se da en las lı́neas 49 a 52. El intervalo de búsqueda está
definido en las lı́neas 17 y 18. En la lı́nea 32 se hace una primera evaluación de c. Luego se
comienza un ciclo iterativo while en la lı́nea 33. El ciclo continua mientras el error |b − c|
sea mayor que la tolerancia establecida. El núcleo del algoritmo, que elige y modifica el
intervalo de búsqueda se programa entre las lı́neas 34 a 44.

Listing 1.1: biseccion-1.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4

5 float ff ( float R ) ;
6

7 int main ()
8 {
9 // biseccion -1. c
10 // o b t i e n e una raiz de la f u n c i o n
11 // f ( x ) = x ^6 - x -1
12 // en el i n t e r v a l o [ a : b ] = [ 1 : 2 ]
13
14 int i =0;
15 float x ,a ,b ,c , e p s i l o n =1e -6;
16
17 a =1;
18 b =1.2;
19

20 // C a l c u l a la ff f ( x ) vs . x
21 // para e x p l o r a r g r a f i c a m e n t e
22

23 // for ( x = a ;x <= b ; x = x +0.01)

24 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
25 // return 0;
26
27 // existe la raiz en el i n t e r v a l o
28 if ( ff ( a ) * ff ( b ) >0)
29 { printf (" No hay raiz en el i n t e r v a l o\ n ") ;
30 return 0;}
31

32 c =( a + b ) /2.;
33 while ( fabs (b - c ) > e p s i l o n)
34 {
35 if ( ff ( b ) * ff ( c ) < 0. )
36 a=c;
37 else
38 b=c;
39
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 19

40 c =( a + b ) /2.;
41
42 i = i +1;
43 printf ("%2 d %1.6 f %1.2 g \ n " ,i ,c , fabs (b - c ) ) ;
44 }
45

46 return 0;
47 }
48

49 float ff ( float x )
50 {
51 return pow (x ,6.) -x -1.;
52 }

La salida del programa, producida por la lı́nea 43 se muestra en el listado siguiente

./bi
1 1.150000 0.05
2 1.125000 0.025
3 1.137500 0.013
4 1.131250 0.0063
5 1.134375 0.0031
6 1.135938 0.0016
7 1.135156 0.00078
8 1.134766 0.00039
9 1.134570 0.0002
10 1.134668 9.8e-05
11 1.134717 4.9e-05
12 1.134741 2.4e-05
13 1.134729 1.2e-05
14 1.134723 6.1e-06
15 1.134726 3.1e-06
16 1.134725 1.4e-06
17 1.134724 7.2e-07
Figure 1.8: Salida del programa biseccion-1.c

Se observa que en la iteración 17 se ha superado la precisión establecida previamente en

1 × 10−6 , y el ciclo se interrumpe. Evidentemente el procedimiento converge a la solución.
En las sección ?? analizaremos en más detalle la convergencia y el error.

Resolución computacional. Alternativa I

El mismo algoritmo de la bisección se puede incluir en una función de C como se muestra

en el programa biseccion-2.c reproducido en el listado 1.2. La función biseccion puede
ser utilizada por otros procedimientos, como por ejemplo, para resolver problemas de raı́ces
múltiples, o pude ser reemplazada por otra función que implemente otro algoritmo.
Listing 1.2: biseccion-2.c
1 # i n c l u d e < stdio .h >
20

2 # i n c l u d e < stdlib .h >

3 # i n c l u d e < math .h >
4 float ff ( float R ) ;
5 float b i s e c c i o n( float a , float b ) ;
6 float ff1 ( float x ) ;
7

8 int main ()
9 {
10 // biseccion -2. c
11
12 // o b t i e n e una raiz de la f u n c i o n
13 // f ( x ) = x ^6 - x -1
14 // en el i n t e r v a l o [a - b ]=[1 -2]
15

16 int i =0;
17 float a ,b , c ;
18 a =1.;
19 b =2.;
20 // C a l c u l a la ff f ( x ) vs . x
21 // para e x p l o r a r g r a f i c a m e n t e
22
23 // for ( x =0.; x <= b ; x = x +0.1)
24 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
25 // return ;
26

27 c = b i s e c c i o n (a , b ) ;
28 printf (" Raiz =% f F ( raiz ) =% f \ n " ,c , ff ( c ) ) ;
29 return 0;
30 }
31

32 float ff ( float x )
33 {
34 return pow (x ,6.) -x -1.;
35 }
36 float b i s e c c i o n( float a , float b )
37 {
38 float c ;
39 int i =0 , imax =100;
40 float e p s i l o n =1.e -7;
41 // existe la raiz en el i n t e r v a l o
42 if ( ff ( a ) * ff ( b ) >0)
43 f p r i n t f( stderr , "# No hay raiz en el i n t e r v a l o \ n ") ;
44

45 c =( a + b ) /2.;
46 while ( fabs (b - c ) > e p s i l o n && i < imax )
47 {
48 if ( ff ( b ) * ff ( c ) < 0. )
49 a=c;
50 else
51 b=c;
52

53 c =( a + b ) /2.;
54 i = i +1;
55 // f p r i n t f( stderr , "%2 d %1.5 f %1.2 g \ n " ,i ,c , fabs (b - c ) ) ;
56 }
57 return c ;
58 }
59
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 21

60 float ff1 ( float x )

61 {
62 return pow (x ,3.) -x -1.;
63 }

Un defecto de este programa es que el nombre de la función de trabajo, en este caso ff,
esta embebido en la función biseccion, lo que es poco conveniente. Para trabajar con
otra función, con otro nombre, habrı́a que modificar la función biseccion, cosa que por lo
general no querremos hacer. La función biseccion podrı́a estar en una librerı́a que no debe
modificarse. Se podrı́a cambiar la función misma ff lo que resultarı́a propenso a errores.
Para solucionar esta clase de problemas, las funciones de C admiten como argumento el
nombre de funciones. La alternativa III a continuación muestra, sin mucha explicación,
como se usa.

Resolución computacional. Alternativa II

Una caracterı́stica de C algo más avanzada permite pasarle a una función el nombre de otra.
El prototipo de una función que recibe como argumento otra función es

void func ( float (*f)(float) );

aquı́ float (*f)(float) es un puntero a una función que retorna un flotante y que toma
como argumento otro flotante. En nuestro caso, el prototipo de la función biseccion será

float biseccion(float a, float b , float (*f)(float) );

{...}

La función se llamarán como se muestra ahora

c=biseccion(a, b , sin);
c=biseccion(a, b , cos);
c=biseccion(a, b , mifuncion);

La definición de la función será, para el caso general

float func ( float (*f)(float) )

{ float x;
...
return (*f)(x);
}

Para el caso particular

float biseccion(float a, float b, float (*f)(float) )

{ float x;
...
22

if( (*f)(a)*(*f)(b);
....
}

El programa biseccion-3.c ser reproduce en listado 1.3. La función biseccion se muestra

por separado7 .

Listing 1.3: biseccion-3.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4

5 float f1 ( float x ) ;
6
7 float b i s e c c i o n( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r) ;
8
9 int main ()
10 {
11 // biseccion -3. c
12 // o b t i e n e una raiz de la f u n c i o n
13 // f ( x ) = x ^6 - x -1
14 // en el i n t e r v a l o [a - b ]=[1 -2]
15

16 int i =0 , i m p r i m a t u r;
17 float a ,b , c ;
18 a =1.;
19 b =2.;
20 // C a l c u l a la ff f ( x ) vs . x
21 // para e x p l o r a r g r a f i c a m e n t e
22
23 // for ( x =0.; x <= b ; x = x +0.1)
24 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
25 // return ;
26

27 c = b i s e c c i o n (a ,b , f1 ,1) ;
28 printf (" Raiz =% f f ( raiz ) =% g \ n " ,c , f1 ( c ) ) ;
29 return 0;
30 }
31

32 float f1 ( float x )
33 {
34 return pow (x ,6.) -x -1.;
35 }

Ahora la función biseccion recibe como argumento el nombre de la función de trabajo. Por
supuesto esta debe estar declarada y definida dentro del alcance de la llamada.
La función biseccion se dispone en el archivo lib-biseccion.c para posibilitar su reuti-
lización y se muestra en el listado 1.4.

Listing 1.4: lib-biseccion.c

7
Esta función y otras para obtener raı́ces debieran guardarse en un único archivo que hace las veces de
librerı́a, o mejor dicho de biblioteca. Aquı́ se las ha mantendio en un archivo individual con el propósito de
facilitar la numeración automática desde 1 que hace el procesador de texto al imprimirlas.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 23

1 # i n c l u d e < stdio .h >

2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4

5 float b i s e c c i o n( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r )

6 {
7 float c ;
8 int i =0 , imax =100;
9 float e p s i l o n =1.e -4;
10 // existe la raiz en el i n t e r v a l o
11 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i a b c
fabs (b - c ) \ n ") ;
12
13 if ((* f ) ( a ) * (* f ) ( b ) >0)
14 { f p r i n t f( stderr , "# No hay raiz en el i n t e r v a l o\ n ") ;
15 return 0;}
16 c =( a + b ) /2.;
17 while ( fabs (b - c ) > e p s i l o n && i < imax )
18 {
19 if ( (* f ) ( b ) *(* f ) ( c ) < 0. )
20 a=c;
21 else
22 b=c;
23
24 c =( a + b ) /2.;
25 i = i +1;
26 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i ,a ,b ,c ,
fabs (b - c ) ) ;
27 }
28 return c ;
29 }

Ahora, la compilación debe hacerse mediante la siguiente orden:

gcc -o bi biseccion-3.c lib-biseccion.c -lm

La salida del programa biseccion-3.c, para ǫ = 10−4 , es la siguiente

Resolución computacional. Alternativa IV

Burden y Faires, [7], sugiere que, para determinar qué subintervalo de [a, b] contiene la raı́z es
mejor utilizar la función signo. En C99 dispone la función copysign(float x, float y)
de la librerı́a tgmath, que compone un valor con el signo de y y el valor de x. La función
devolverá 1 o −1. Da el mismo resutado pero evita overflows en la multiplicación f (a)f (b).
Esta alternativa se ha implementado en la función biseccion6, en el archivo
lib-biseccion-6.c, listada en 1.5
Listing 1.5: lib-biseccion-6.c
1 # include < stdio .h >
2 # include < stdlib .h >
3 # include < math .h >
4 # include < tgmath .h >
5
24

$ ./bi
# i a b c fabs(b-c)
1 1.00000 1.50000 1.25000 0.25
2 1.00000 1.25000 1.12500 0.12
3 1.12500 1.25000 1.18750 0.062
4 1.12500 1.18750 1.15625 0.031
5 1.12500 1.15625 1.14062 0.016
6 1.12500 1.14062 1.13281 0.0078
7 1.13281 1.14062 1.13672 0.0039
8 1.13281 1.13672 1.13477 0.002
9 1.13281 1.13477 1.13379 0.00098
10 1.13379 1.13477 1.13428 0.00049
11 1.13428 1.13477 1.13452 0.00024
12 1.13452 1.13477 1.13464 0.00012
13 1.13464 1.13477 1.13470 6.1e-05
# Raiz=1.134705 f(raiz)=-0.000201099
$
Figure 1.9: Salida del programa biseccion-3.c

6 float b i s e c c i o n 6
7 ( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r)
8 {
9 float c ;
10 int i =0 , imax =100;
11 float e p s i l o n =1.e -4;
12 // existe la raiz en el i n t e r v a l o
13 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i a b c
fabs (b - c ) \ n ") ;
14
15 if ( c o p y s i g n (1.0 ,(*f ) ( a ) ) * c o p y s i g n (1.0 ,(*f ) ( b ) ) >0)
16 { f p r i n t f( stderr , "# No hay raiz en el i n t e r v a l o\ n ") ;
17 return 0;}
18 c =( a + b ) /2.;
19 while ( fabs (b - c ) > e p s i l o n && i < imax )
20 {
21 if ( c o p y s i g n (1.0 ,(*f ) ( b ) ) * c o p y s i g n (1.0 ,(*f ) ( c ) ) <0)
22 a=c;
23 else
24 b=c;
25
26 c =( a + b ) /2.;
27 i = i +1;
28 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i ,a ,b ,c ,
fabs (b - c ) ) ;
29 }
30 return c ;
31 }
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 25

La función bisección de Numerical Recipies

A efectos de mostrar que cada algoritmo se puede desarrollar de muchas maneras diferentes,
presentamos aquı́ en una una versión libre, ajustada a nuestra notación, la función zbrac
del libro Numerical Recipies de Press et al. Se recomienda al lector usar la implementación
original de los autores.

Listing 1.6: biseccion-5.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4 # define MAX 20
5 float ff ( float R ) ;
6
7 main ()
8 {
9 // biseccion -5. c
10 // o b t i e n e una raiz de la f u n c i o n
11 // f ( x ) = x ^6 - x -1
12 // en el i n t e r v a l o [a - b ]=[1 -2]
13

14 int i =0;
15 float a ,b ,c , e p s i l o n =1 e -5;
16 float f , fc ,r , dx ;
17
18 a =1.;
19 b =2.;
20

21 // C a l c u l a la ff f ( x ) vs . x
22 // para e x p l o r a r g r a f i c a m e n t e
23

24 // for ( x =0.; x <= b ; x = x +0.1)

25 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
26 // return ;
27
28 // existe la raiz en el i n t e r v a l o
29

30 fc = ff ( b ) ;
31 f = ff ( a ) ;
32

33 if ( fc *f >0)
34 printf (" No hay raiz en el i n t e r v a l o \ n ") ;
35
36 // Se ordena la d i r e c c i o n de b u s q u e d a
37 if (f <0)
38 { r=a;
39 dx =b - a ;}
40 else
41 { r=b;
42 dx =b - a ;}
43
44 // B i s e c c i o n
45 do
46 {
47 dx = dx /2.;
48 c = r + dx ;
26

49 fc = ff ( c ) ;
50 if ( fc <0) r = c ; // m a n t i e n e a c o t a d a la raiz
51 i ++;
52 f p r i n t f( stderr , "%2 d %1.5 f %1.2 g \ n " ,i ,r , fabs ( dx ) ) ;
53 }
54 while ( fabs ( dx ) > e p s i l o n || i < MAX ) ;
55

56 printf ("%2 d %1.5 f %1.2 g \ n " ,i ,r , fabs ( dx ) ) ;

57

58 }
59

60 float ff ( float x )
61 {
62 return pow (x ,6.) -x -1.;
63 }

Esta implementación difiere de las anteriores en que primero se determina (lı́neas 37 − 42)
cual de los extremos del intervalo corresponde a la función con valor negativo y toma este
valor, que se denomina r en el programa, como punto de referencia. La longitud del intervalo
es dx, y el valor medio del intervalo se calcula sumando dx/2. al valor de la referencia dx.
El el algoritmo de bisección se desarrolla en las lı́neas 45 − 54. El listado siguiente muestra
los primeros 10 valores de la salida

$ head dat
1 1.00000 0.5
2 1.00000 0.25
3 1.12500 0.12
4 1.12500 0.062
5 1.12500 0.031
6 1.12500 0.016
7 1.13281 0.0078
8 1.13281 0.0039
9 1.13281 0.002
10 1.13379 0.00098
$

1.8.1 Ejemplo: cálculo del coeficiente de arrastre.

Retomemos el ejemplo presentado en la sección 1.1.1 sobre el cálculo del coeficiente de
arrastre para un objeto en caı́da libre con velocidad terminal especificada. El primer paso
en la resolución de un problema de este tipo es escribir la función de la forma

g(k) = 0 (1.42)

que para este caso es

gm

1 − e−ktf /m − vf = 0 (1.43)
k
la función g(k) podrá programarse como una función de C (lı́neas 53 a 57 del listado).
Para conocer la forma de la función g(k) y para oner una aproximación al resultado podemos
obtener una gráfica de g(k) vs. k. Para ello se calcula g(ki ) para valores discretos de k,
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 27

que llamaremos ki , pertenecientes al intervalo de k que se quiere investigar, se imprime los

resultados en forma de tabla (lı́neas 16 a 18, comentadas) y se grafican con gnuplot. El
gráfico de la función se muestra en la figura 1.10 y del mismo se puede obtener un valor
aproximado para la raı́z, α ∼ 5. El programa arrastre-b2.c reproducido en el listado 1.7

80
”[Link]”
70
60
50
40
f (c)
30
20
10
0
-10
0 1 2 3 4 5 6 7 8 9 10
c
Figure 1.10: La función f (k) vs. k.

es una modificación del programa biseccion-3.c y usa la funcón biseccion. Solo se ha

cambiado la definición de la función ff y se ha provisto en ella parámetros necesarios que se
muestran en la tabla 1.1.

m g t vt ǫ k
2
kg m /seg seg m/seg 1 kg/seg
5 10 9 10 1 × 10−6 4.99938

Table 1.1: Coeficiente de arrastre. Parámetros y resultados del problema.

Listing 1.7: arrastre-b2.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4
5 float b i s e c c i o n( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r) ;
6 float ff ( float k ) ;
7
8 int main ()
9 {
10
11 // arrastre - b2 . c
12 // busca el c o e f i c i e n t e de a r r a s t r e
13 // para a l c a n z a r una cierta v e l o c i d a d predeterminada
14 // en un lapso de tiempo dado
15 // aplica el metodo de la b i s e c c i o n
16

17 int i =0;
28

18 float a ,b , k ;
19 // C a l c u l a la ff f ( x ) vs . x
20 // para e x p l o r a r g r a f i c a m e n t e
21 // for ( k =0.; k <=10.; k = k +0.1)
22 // printf ("% f % f \ n " ,k , ff ( k ) ) ;
23 // return ;
24

25 // i n t e r v a l o de b u s q u e d a
26 a =0.;
27 b =6.;
28

29 k = b i s e c c i o n (a ,b , ff ,1) ;
30
31 printf (" Raiz =% f f ( raiz ) =% g \ n " ,k , ff ( k ) ) ;
32 return 0;
33 }
34 float ff ( float k )
35 {
36 float v ,t ,m , g ;
37

38 g =10.; // a c e l e r a c i o n de la garvedad , m / s
39 t =9.; // tiempo final , s
40 m =5.; // mass , kg
41 v =10.; // v e l o c i d a d final , m / s
42

43 return ( g * m / k ) *(1. - exp ( - k * t / m ) ) -v ;

44 }

El listado siguiente muestra los resultados de las últimas iteraciones antes de alcanzar la
convergencia de donde se obtiene que, dentro de la precisión deseada k ∼ 4.99938:

$ ./bi
$ i a b c fabs(b-c)
1 3.00000 6.00000 4.50000 1.5
2 4.50000 6.00000 5.25000 0.75
3 4.50000 5.25000 4.87500 0.38
4 4.87500 5.25000 5.06250 0.19
5 4.87500 5.06250 4.96875 0.094
6 4.96875 5.06250 5.01562 0.047
7 4.96875 5.01562 4.99219 0.023
8 4.99219 5.01562 5.00391 0.012
9 4.99219 5.00391 4.99805 0.0059
10 4.99805 5.00391 5.00098 0.0029
11 4.99805 5.00098 4.99951 0.0015
12 4.99805 4.99951 4.99878 0.00073
13 4.99878 4.99951 4.99915 0.00037
14 4.99915 4.99951 4.99933 0.00018
15 4.99933 4.99951 4.99942 9.2e-05
Raiz=4.999420 f(raiz)=-7.58616e-05

La figura 1.11 muestra gráficamente la aproximación convergente desde el comienzo de la

iteración.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 29

5.3
”[Link]”
5.2

5.1

c 4.9
4.8

4.7

4.6

4.5
0 5 10 15 20 25
iteracion

Figure 1.11: Convergencia a la raı́z

30

1.9 Resolución computacional con Python

El primer paso en la busqueda de raı́ces debe ser estudiar la función que tenemos entre manos.
La forma más instructiva de hacerlo es mediante un gráfico. Para ello es conveniente:

• Definir una función de Python que represente nuestra función matemática.

• Dibujar la función.
• Obtener la raı́z en forma gráfica u obtener el rango de abscisas en que se encuentra la
raı́z.

1.9.1 Funciones
El script [Link] mostrado en el listado 1.8 es un ejemplo muy pequeño y simple
de la estructura de un programa en Python.

✞ Listing 1.8: [Link]

1 # f u n c i o n e s 0 0 . py
2 import numpy as np
3

4 # =================================
5 # D e f i n i c i o n de l a f u n c i o n
6 def fun1 (x):
7 return x∗∗3−5∗x−9
8

9 # =================================
10 # Programa p r i n c i p a l
11 # Entrada de d a t o s
12 a=1
13 b=2
14

15 # Impresion
16 rango = [Link] (a, b, 5)
17 print (rango )
18 print (fun1 (rango ))
✝ ✆

El script tiene tres partes, aunque sin ninguna discontinuidad entre ellas, salvo la barra
horizontal que hemos incluido como comentario para resaltar lo dicho. En la primera usual-
mente se importan los módulos con los que se quiere trabajar. Recordemos que los módulos
son archivos terminados en .py que contienen colecciones de código relacionado, funciones,
clases, variables, etc. Volveremos a ello en seguida. En la segunda usualmente se definen
varias funciones que luego usará el programa principal. En este script se han definido una
función de Python fun1 que en este caso implementa la función matemática
y = f (x) = x3 − 5x − 9 (1.44)
Para cada valor del argumento x calculará el valor de la función y. En seguida, después de
graficar esta función, obtendremos alguna de sus raı́ces.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 31

En la tercera parte del script, que hemos denominado el programa principal, se programa
la tarea que se requiere llevar a cabo. Se hace la entrada de datos y se ejecuta la tarea y si
es necesario se hace el posprocesamiento. En nuestro caso se ha establecido los datos dando
valores a las variables a y b, se hace uso de la función linspace del módulo numpy para crear
un rango de valores entre a y b, se imprime estos valores y se calcula e imprime la imágen
de estos valores a través de la función fun1
La ejecución del programa resulta en la siguiente salida 1.12

python [Link]
[1. 1.25 1.5 1.75 2. ]
[-13. -13.296875 -13.125 -12.390625 -11. ]
Figure 1.12: Salida del programa [Link]

En general, el programa principal de un problema de apenas mayor complejidad se com-

pondrá de una secuencia de funciones que con toda seguridad llaman a otras y estas a su vez
muchas más. Python provee varios métodos para organizar y agrupar nuestras funciones de
acuerdo a su funcionalidad de manera de facilitar su posterior reusabilidad. Una de estas
formas es el uso de módulos. Muchos de estos ya están programados, como el módulo NumPy,
(Numerical Python), que como se presenta a si mismo es es una biblioteca para trabajar
con datos numéricos y es el núcleo de los ecosistemas cientı́ficos de Python. Pero lo más
interesante es que podremos programar nuestros propios módulos para satisfacer nuestras
las particulares necesidades de nuestros programas.

1.9.2 Definción de funciones y módulos

Para adentrarnos en el uso de la programación de módulos, supongamos, como es natural en
este capı́tulo, que queremos calcular las raı́ces de varias funciones distintas. Anticipándonos,
probablemente querramos probar varios métodos de solución diferentes o varias formas de
programar el mismo algoritmo. Parece lógico entonces que, para empezar, agrupemos la
definición de estas funciones en un único módulo. Como dijimos, un módulo es un archivo
cuyas funciones, parámetros, etc. Se puede incluir o importar en cualquier otro módulo,
script o programa. Las funciones del módulo se hacen accesible al script que lo incorpora.
El archivo [Link] mostrado en el listado 1.9 define varias funciones (solo dos por
ahora) de ejemplo. Pero como verá, hemos dejado un par de lı́neas de código que usan estas
funciones. Por supuesto este programa no sabe de nuestra intención de usralo como módulo
y se puede ejecutar como cualquier otro.

✞ Listing 1.9: [Link]

1 # f u n c i o n e s 0 . py
2 import numpy as np
3

4 def fun1 (x):

5 return x∗∗3−5∗x−9
6

7 def fun2 (k,t,v,m,g):

8 ””” Caida l i b r e con f r i c c i o n
32

9 k c o e f i c i e n t e de f r i c c i o n
10 t ti e m p o
11 v velocidad ini cial
12 m masa
13 g a c e l e r a c i o n de l a g r a v e d a d ”””
14

15 return (g∗m/k)∗(1−[Link](−k∗t/m))−v
16

17

18 print (fun1 (5))

19 print (fun2 (5 ,3 ,10 ,5 ,9.8))
✝ ✆

Es bastante común que cuando estemos programando la solución de un problema no solo

definamos varias funciones a medida que avanzamos sino que las vayamos usando. Asi se ha
hecho en el script [Link] De alguna manera esto nos irá diciendo si progresamos
en la dirección correcta. Esta forma de trabajar es bastante común en un ambiente en
donde quienes programan no son profecionales informáticos. Veremos que es posible adoptar
algunas buenas prácticas en seguida. No obstante, nos encontraremos entonces que junto con
nuestras funciones tenemos parte o todo el programa principal. Supongamos por ahora que
en nuestro archivo de funciones han quedado las dos lı́neas de código que imprimen algunos
resultados y hacen uso de las funciones definidas.
Este archivo, como cualquiera terminado en .py puede ser incorporado como módulo en
cualquier otro programa haciéndole accesible sus funciones. Solo debe introducirse la sigu-
iente sentencia en la primera parte de nuestro script,
✞
1 import funciones as fun
✝ ✆

El programa [Link] que se lista más adelante es igual al programa en sus efectos al
[Link] salvo que ya no definimos en él a la función fun1 ni ninguna otra.

✞ Listing 1.10: [Link]

1 # e j e m p l o 1 . py
2 import numpy as np
3 import funciones0 as f
4

5 # =================================
6 # Programa p r i n c i p a l
7 # Entrada de d a t o s
8 a=1
9 b=2
10

11 # Impresion
12 rango = [Link] (a, b, 5)
13 print (rango )
14 print (f.fun1 (rango ))
✝ ✆

La salida del programa es la siguiente ¡Ahora atención! La salida contiene la que pretendemos
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 33

python [Link]
91
-0.6879132700050654
[1. 1.25 1.5 1.75 2. ]
[-13. -13.296875 -13.125 -12.390625 -11. ]
Figure 1.13: Salida de [Link]

en nuestro programa principal, [Link], y la de nuestro módulo [Link]. ¡Esto

no es lo que queremos!
En efecto, la incorporación de este módulo en otro programa no solo habilitará el uso de las
funciones allı́ definidas sino que además también ejecutará su código principal. El programa
que llama al módulo tendrá su propio código principal y se ejecutarán ambos. Esto no es
lo que queremos. Para que el módulo pueda correr su programa principal si se lo llama a
ejecución autónoma desde el intérprete python y a la ves no lo haga si se lo incluye como
módulo en otro script, debe modificarse como en el script funciones1 mostrado en el listado
??.

✞ Listing 1.11: [Link]

1 # f u n c i o n e s 1 . py
2 import numpy as np
3

4 def fun1 (x):

5 return x∗∗3−5∗x−9
6

7 def fun2 (k,t,v,m,g):

8 ””” Caida l i b r e con f r i c c i o m
9 k v a r i a b l e i n d e p e n d i e n t e c o e f i c i e n t e de a r r a s t r e
10 parametros :
11 t ti e m p o
12 v velocidad ini cial
13 m masa
14 g a c e l e r a c i o n de l a g r a v e d a d ”””
15 return (g∗m/k)∗(1−[Link](−k∗t/m))−v
16

17 # Programa p r i n c i p a l
18 # Se l l a m a a e j e c u c i o n l a s f u n c i o n e s d e f i n i d a s
19 # E s t e programa s e e j e c u t a r a como programa autonomo .
20 # Sus f u n c i o n e s ta m b i e n pueden s e r l l a m a d a s d e s d e o t r o
programa
21

22 if name == " m a i n ":

23 print (fun1(−10))
24 print (fun2 (5, 3 ,10 ,5 ,9.8) )
✝ ✆

Cambiando en el scrip [Link] la sentencia =import funciones0 as f por la senten-

cias import funciones1 as f solo se ejecutará el programa principal del propio script.
34

1.9.3 Validación de funciones

Ahora, antes de proseguir, nos detengamos todavı́a en nuestro nuevo módulo [Link].
Es recomendable (Langtangen 2016, [? ]) verificar que las funciones programadas funcio-
nen correctamente. La forma más obvia de hacerlo es comparar sus resultados con los que
podemos calcular a mano o con una calculadora. Python provee varios medios para facili-
tar la prueba de un código. Uno de ellos es la denominada prueba afirmativa, edecuada
para hacer una comprobación de unidad ([Link]
#unit-tests-vs-integration-tests). Básicamente en esta prueba comparamos el resul-
tado de nuestro script con tres resultados correctos que hemos calculado por algún otro
medio. Por cada función definida podemos hacer una función de prueba. Para el caso de la
función fun1 en el programa [Link], la función test_fun1() hace la comparación
entre la salida de la función para un grupo de datos especificado con los resultados esperados
calculados a mano.
En la función test_fun1() se crea un arreglo, con tres valores de entrada, datos y se crea
otro arreglo con los resultados de aplicar la función a esos datos. En un tercer arreglo se
almacena los valores calculados manualmente, con calculadora, con symolab [Link]
[Link]/ u otro recurso online confiable, a partir de los datos. Luego se calcula el valor
absoluto de la diferencia entre unos y otros y de estas diferencias se toma el valor máximo
y se asigna a la variable exito. Esta es una variable boleana que puede tomar valor True o
False Si la variable exito resulte falsa, la sentencia assert exito interrumpirá la ejecución
del programa. En caso contrario no hará nada. El ejemplo mostrado, sólo demostrativo de
la estructura de la función de prueba es muy sencillo. Por lo general nuestras funciones serán
más complicadas y se podrán progamar de distintas maneras. No obstante, los resultados
deberı́a ser los mismos. Una función de prueba como la mostrada ayudarı́a a comprobar si
las modificaciones no alteran los resultados.
Con estos arreglos nuestro módulo funciones queda

✞ Listing 1.12: [Link]

1 # f u n c i o n e s . py
2 import numpy as np
3

4 # d e f i n i c i o n de f u n c i o n e s
5

6 def fun1 (x):

7 return x ∗∗3 − 5∗ x − 9
8

9 def fun2 (x):

10 return x ∗∗3 − 5∗ x
11

12 # B l o c k de p r u e b a
13

14 def test fun1 ():

15 ””” Compara l a s a l i d a de l a f u n c i o n con r e s u l t a d o s
manuales
16 c a l c u l a d o s p a ra a=−10, b =0 , c= 10 ”””
17 tolerancia = 1.E−10
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 35

18 datos =[Link] ([−10,0,10])

19 exacta =[Link] ([−959, −9, 941])
20 numericos =fun1 (datos )
21 diferencia = abs(numericos−exacta ).max ()
22 excito = diferencia < tolerancia
23 assert excito
24

25

26 def test fun2 ():

27 ””” Compara l a s a l i d a de l a f u n c i o n con r e s u l t a d o s
manuales
28 c a l c u l a d o s p a ra a=0 , b=−np . s q r t ( 5 ) , c=np . s q r t ( 5 ) ”””
29 tolerancia = 1.E−15
30 datos =[Link] ([−10,0,10])
31 exacta =[Link] ([−950, 0, 950])
32 numericos =fun2 (datos )
33 diferencia = abs(numericos−exacta ).max ()
34 excito = diferencia < tolerancia
35 assert excito
36

37

38 # Programa p r i n c i p a l
39 # Si se llama a ejecucion e l s c r i p t
40 # s e e j e c u t a r a como programa autonomo .
41 # S i l a s f u n c i o n e s no pasan l a p r u e b a l a
42 # ejecucion se interrumpe
43 # S i s e i m p o r t a no e j e c u t a e l programa p r i n c i p a l
44 #
45 if name == " m a i n ":
46 test fun1 ()
47 test fun2 ()
48 print ("Prueba pasada con exito ")
49 print ("fun1 ", fun1(−10))
50 print ("fun2 ", fun2 (0))
51 print ("fun2 ", fun2 ([Link] (5)))
52 print ("fun2 ", fun2(−[Link] (5)))
✝ ✆

En los ejemplos sucesivos cada nueva función se programará en el módulo [Link] y

se hará su correspondiente función de prueba.

1.9.4 Gráfico de funciones

Python tiene numerosas opciones de graficación. Aquı́ usaremos las más básicas del módulo
matplotlib. El archivo [Link] listado en 1.13 muestra un ploteo sencillo. No
obstante la sencillez, requiere de los módulos numpy y matplotlib y la definición de las
funciones que queremos graficar. Entonces, creamos dos vectores, uno que denominaremos
x, con los valores discretos de las absisas y otro, que denominaremos y, con la imagen cor-
36

respondiente a través de la función. El vector de abscisas igualmente espaciadas se crea con

linspace mientras que el vector de ordenadas y se crea aplicando x a la función previamente
definida.
El archivo [Link] listado en 1.13 muestra una función adecuada para facilitar el
ploteo de funciones.

✞ Listing 1.13: [Link]

1 # p l o t e o B a s i c o . py
2 import numpy as np
3 import matplotlib .pyplot as plt
4 import funciones as f
5

6 # Dato : i n t e r v a l o de c a l c u l o
7 a=1
8 b=10
9 # Crea un ra n g o de v a l o r e s de x y de y
10 x = np. linspace(a, b, 1000) # a b s c i s a s d i s c r e t a s
11 y = f.fun1 (x) # coordenadas d i s c r e t a s
12 # P l o t e o s i m p l e de l o s d a t o s
13 [Link] (x,y)
14 [Link] ()
✝ ✆

El archivo [Link] listado en 1.14 muestra una scrit mejorado para facilitar
y adornar el ploteo de funciones:

✞ Listing 1.14: [Link]

1 # p l o t e o I n t e r m e d i o . py
2 import numpy as np
3 import matplotlib .pyplot as plt
4 import funciones as f
5

7 # I n t e r v a l o de c a l c u l o
8 a=−5
9 b=5
10 # Crea un ra n g o de v a l o r e s de x y de y
11 x = np. linspace(a, b, 50) # a b s c i s a s d i s c r e t a s
12 y = f.fun1 (x) # coordenadas d i s c r e t a s
13

14 # P l o t e o d e c o r a d o de l o s d a t o s
15 plt. xlabel ( ' x ' )
16 plt. ylabel ( ' y ' )
17 [Link] ( ' fun1 ' )
18 plt. axhline(y=0, c="blue ")
19 [Link] (x,y, '−ob ' , label ="f(x)=x∗∗3−5∗x−9")
20 plt. legend ()
21 plt. savefig( ' fig−fun1 .eps ' )
22 [Link] ()
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 37

✝ ✆

El programa, además, crear un archivo con la figura 1.14 que se muestra a continuación.

fun1
100
f(x)=x**3-5*x-9

75

50

25

0
y

−25

−50

−75

−100

−4 −2 0 2 4
x

Figure 1.14: Gráfica de la función fun1 creada con el programa [Link] y

archivada como [Link].

Por último, podemos crear una función para hacer este trabajo, es decir el ploteo cada vez
que necesitemos hacerlo. La función se incorporará al modulo [Link]. [Link]
se puede ejecutar de forma autónoma o se puede llamar como un módulo desde otro script.
La función se muestra en el listado 1.15. La ejecución autónoma produce la figura que se ve
en pantalla y el archivo [Link] que se ha incorporado aquı́ como figura 1.15.

✞ Listing 1.15: [Link]

1 # p l o t e o . py
2 import numpy as np
3 import matplotlib .pyplot as plt
4 import funciones as f
5

6 def miploteo(ff , a,b,n, titulo ,leyenda , archivo):

7 ””” Funcion de p l o t e o
8 f f funcion a plotear
9 x , y v a l o r e s d i s c r e t o s de l a f u n c i o n
10 a , b intervalo
11 t i t u l o t i t u l o o nombre de l a f u n c i o n
12 leyenda : la funcion a grafica
13 a r c h i v o : nombre d e l a r c h i v o a g r a f i c a r
14 ”””
15 # P l o t e o d e c o r a d o de l o s d a t o s
16 # Crea un ra n g o de v a l o r e s de x y de y
17 x = [Link] (a, b, n) # a b s c i s a s d i s c r e t a s
18 y = ff(x) # coordenadas d i s c r e t a s
19
38

20 plt. xlabel ( ' x ' )

21 plt. ylabel ( ' y ' )
22 [Link] (titulo )
23 plt. axhline(y=0, c="blue ")
24 [Link] (x,y, '−ob ' , label = leyenda)
25 plt. legend ()
26 plt. savefig( archivo)
27 [Link] ()
28

29

30 if name == " main ":

31

32 miploteo(f.fun1 ,−10,10,50,"fun1 ","t=x^3 −5x−9","fig−fun1−

test ")
✝ ✆

fun1
1000
t=x^3 -5x-9

750

500

250

0
y

−250

−500

−750

−1000

−10.0 −7.5 −5.0 −2.5 0.0 2.5 5.0 7.5 10.0

x

Figure 1.15: Gráfica de la función fun1 creada con el programa [Link] y archivada
como [Link].

Notará que en módulo [Link] no hemos includo un bloque de prueba. Las pruebas
para una función como la presente son más complejas. Primero hay que determinar que se
quiere probar. El estudiante interesado puede acceder a la web para mayor información.

1.9.5 Trabajando desde la terminal interactiva

Hasta ahora hemos corrido los programas llamando a python como si fueran autónomos.
Podemos usar la terminal interactiva ipython para probar nuestra funciones con distintos
datos. Las lineas siguientes muestran como.

ipython
Python 3.9.12 (main, Apr 5 2022, 06:56:58)
Type ’copyright’, ’credits’ or ’license’ for more information
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 39

IPython 8.4.0 -- An enhanced Interactive Python. Type ’?’ for help.

In [1]: import funciones as f

In [2]: import miploteo as mip

In [3]: [Link](f.fun1,-4,4,50,"fun1", "x^3 -5x -9", "fig-fun1-1")

La salida de este sentencia es la misma que la de anates, salvo que hemos cambiado los datos.

1.9.6 El método de la Bisección en Python

Siendo Python un lenguaje de alto nivel, existen muchas maneras distintas de programar
cualquier método. Muchas de estas han sido incorporadas a los módulos estandar de Python,
por ejemplo en ScyPy. Ya veremos como usar sus funciones, pero el propósito aqui es
utilizar estos algoritmos para familiarizarnos con los métodos numéricos y su programación.
Ya tendremos tiempo de comprar todo hecho.
El programa [Link] listado en 1.16 es una implementación muy simple del método
de la bisección, sigue a rajatabla el algoritmo descripto en la sección 1.5. El script tiene dos
partes, aunque no discontinuidad entre ellas, salvo la barra horizontal que hemos incluido
como comentario para resaltar lo dicho. En la primera se han definido las funciones que
compondrán nuestro script, en este caso, la función de Python que calculará la función
matemática en cuestión, fun1 y la función de Python que implementa el algoritmo de la
bisección, biseccion-1. En la segunda parte, que hemos denominado el programa principal,
se establecen los datos y se llama a ejecución la función la función biseccion-1. Cualquier
programa de mayor complegidad tendrá estas dos partes. El programa listado en 1.16 es
una implementación muy simple.

✞ Listing 1.16: [Link]

1 # b i s e c c i o n −1. py
2 import matplotlib .pyplot as plt
3 import numpy as np
4 import funciones as f
5

6 # Metodo de l a b i s e c c i o n
7 def biseccion (a,b,tol):
8 ””” Metodo de l a b i s e c c i o n ”””
9 paso = 1
10 condicion = True
11 while condicion :
12 c = (a + b)/2
13 print ( ' Iteracion %d, c = %0.6 f and f.fun1 (c) = %0.6 f '
% (paso , c, f.fun1 (c)))
14

15 if f.fun1 (a) ∗ f.fun1 (c) < 0:

16 b = c
40

17 else :
18 a = c
19

20 paso = paso + 1
21 condicion = abs(f.fun1 (c)) > tol
22

23 print ( ' \ nLa raiz es: %0.8 f ' % c)

24

25 # Programa p r i n c i p a l
26 print ("#biseccion −[Link]:")
27 ”””
28 # Entrada de d a t o s
29 a = i n p u t ( ' I n t e r v a l o , e x tre m o a : ' )
30 b = i n p u t ( ' I n t e r v a l o , e x tre m o b : ' )
31 t o l = input ( ' Tolerancia : ')
32

33 # Convecion a f l o t a n t e s
34 a = float (a)
35 b = float (b)
36 tol = float ( tol )
37 ”””
38 a=1
39 b=5
40 tol =1.E−6
41

42

43 #Nota : R e d a c c i o n a l t e r n a t i v a
44 # a = f l o a t ( i n p u t ( ' I n t e r v a l o , e x tre m o a : ' ) )
45 # b = f l o a t ( i n p u t ( ' I n t e r v a l o , e x tre m o b : ' ) )
46 # t o l = f l o a t ( input ( ' Tolerancia : ') )
47

48

49 # Chequeo d e l i n t e r v a l o i n i c i a l
50 if f.fun1 (a) ∗ f.fun1 (b) > 0.0:
51 print ( ' El intervalo dado no contiene la raiz o contiene un
numero par de ellas . ' )
52 print ( ' Eliga un nuevo intervalo . ' )
53 else :
54 biseccion (a,b,tol)
✝ ✆

La salida del programa es la siguiente:

A-cursos/cursoMNu/raı́[Link] - August 22, 2023 41

python [Link]
#[Link]:
Iteracion 1, c = 3.000000 and f.fun1(c) = 3.000000
Iteracion 2, c = 2.000000 and f.fun1(c) = -11.000000
Iteracion 3, c = 2.500000 and f.fun1(c) = -5.875000
Iteracion 4, c = 2.750000 and f.fun1(c) = -1.953125
Iteracion 5, c = 2.875000 and f.fun1(c) = 0.388672
Iteracion 6, c = 2.812500 and f.fun1(c) = -0.815186
Iteracion 7, c = 2.843750 and f.fun1(c) = -0.221588
Iteracion 8, c = 2.859375 and f.fun1(c) = 0.081448
Iteracion 9, c = 2.851562 and f.fun1(c) = -0.070592
Iteracion 10, c = 2.855469 and f.fun1(c) = 0.005297
Iteracion 11, c = 2.853516 and f.fun1(c) = -0.032680
Iteracion 12, c = 2.854492 and f.fun1(c) = -0.013700
Iteracion 13, c = 2.854980 and f.fun1(c) = -0.004204
Iteracion 14, c = 2.855225 and f.fun1(c) = 0.000546
Iteracion 15, c = 2.855103 and f.fun1(c) = -0.001829
Iteracion 16, c = 2.855164 and f.fun1(c) = -0.000641
Iteracion 17, c = 2.855194 and f.fun1(c) = -0.000048
Iteracion 18, c = 2.855209 and f.fun1(c) = 0.000249
Iteracion 19, c = 2.855202 and f.fun1(c) = 0.000101
Iteracion 20, c = 2.855198 and f.fun1(c) = 0.000027
Iteracion 21, c = 2.855196 and f.fun1(c) = -0.000011
Iteracion 22, c = 2.855197 and f.fun1(c) = 0.000008
Iteracion 23, c = 2.855196 and f.fun1(c) = -0.000001
Iteracion 24, c = 2.855197 and f.fun1(c) = 0.000003
Iteracion 25, c = 2.855197 and f.fun1(c) = 0.000001
Iteracion 26, c = 2.855197 and f.fun1(c) = -0.000000
La raiz es: 2.85519654
42

Este programa tiene el inconveniente que solo trabaja con la función fun1. Un pequeño
cambio en la entrada de datos de la función biseccion habilita la entrada del nombre de la
función como su argumento. La nueva versión del programa, ahora listado en 1.17 muestra
este cambio. También se ha usado el módulo [Link], se ha realizado una función
de prueba test_biseccion.py y se ha redactado como módulo, de manera que ahora este
script se puede ejecutar en forma autónoma o se puede llamar desde otro script.

✞ Listing 1.17: [Link]

1 # p r u e b a b i s e c −1. py
2 import matplotlib .pyplot as plt
3 import numpy as np
4 import funciones as f
5

6 # Metodo de l a b i s e c c i o n
7 def biseccion (ff ,a,b,tol):
8 ””” Metodo de l a b i s e c c i o n ”””
9 paso = 1
10 condicion = True
11 while condicion :
12 c = (a + b)/2
13 #p r i n t ( ' I t e r a c i o n %d , c = %0.6 f and f f ( c ) = %0.6 f ' % (
paso , c , f f ( c ) ) )
14

15 if ff(a) ∗ ff(c) < 0:

16 b = c
17 else :
18 a = c
19

20 paso = paso + 1
21 #c o n d i c i o n = a b s ( a−c )> t o l #a b s ( f f ( c ) ) > t o l
22 condicion = abs(ff(c)) > tol
23

24 #p r i n t ( ' \ nLa r a i z e s : %0.8 f ' % c )

25 return c
26

27 def test biseccion ():

28 ””” Compara l a r a i z c a l c u l a d a p o r b i s e c c i o n con l a r a i z
exacta
29 p a ra t r e s f u n c i o n e s ”””
30 tolerancia = 1.E−5
31 numericos =[Link] ([ biseccion (f.fun1 ,1,5,1. E−15),biseccion
(f.fun2 ,−5,−0.9,1.E−10)])
32 exacta =[Link] ([2.85519,− [Link] (5) ])
33 diferencia = abs(numericos−exacta ).max ()
34 excito = diferencia < tolerancia
35 assert excito
36

37
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 43

38 if name == " main ":

39

40 test biseccion ()
41 print ( biseccion (f.fun1 ,1,5,1. E−10))
42 print ( biseccion (f.fun2 ,−5,−0.9,1.E−10))
✝ ✆
44

1.10 El método de Newton-Raphson

El método de Newton utiliza una estrategia diferente. Aproxima la función f (x) en un punto
x0 por la recta tangente a la función en ese punto. El punto x0 debe ser un valor cercano a la
raı́z que denominamos semilla o valor de prueba. La pendiente de la recta tangente está
dada por la derivada de la función, f ′ (x0 ). Luego, se puede obtener un valor aproximado de la
raı́z tomando la intersección de esta recta con el eje de las abscisas. Entonce si la intersección
ocurre para x1 , este valor es la primera aproximación a la raı́z α. Luego el procedimiento
se repite tomando el valor x1 como nueva semilla, calculando la nueva pendiente f ′ (x1 ) y
usando la recta definida por ella para obtener una segunda aproximación x2 . El proceso
se repite para x3 , x4 y ası́ hasta xn y se sigue hasta que la diferencia entre los dos últimos
valores consecutivos de la serie de aproximaciones, xn y xn+1 sea suficientemente pequeña.
Para derivar el algoritmo podemos comenzar haciendo el desarrollo en series de Taylor alrede-
dor de x + ∆x, que hacemos solo hasta los términos lineales
(x + ∆x) − x ′
f (x) = f (x + ∆x) + f (x) + · · · (1.45)
1!
de donde surge la aproximación en diferencias hacia adelate para la derivada primera
f (x + ∆x) − f (x)
f ′ (x) ≈ (1.46)
(x + ∆x) − x
Ahora, como dijimos, reemplazamos la función f (x) por la recta tangente g(x), y podemos
escribir
f (x + ∆x) − f (x) g(x + ∆x) − f (x)
f ′ (x) ≈ ≈ (1.47)
(x + ∆x) − x (x + ∆x) − x
La igualdad es aproximadamenta válida si ∆x es pequeño. Nos interesa ahora elegir ∆x
tal que x + ∆x = β, donde β es la raı́z de la recta tangente (esta será nuestra primera
aproximación a la raı́z de f (x)). Esta se obtiene en la intersección de g(x) con las abscisas,
en donde g(x + ∆x) = 0 y la ecuación 1.47 ser reduce a
0 − f (x)
f ′ (x) ≈ (1.48)
(x + ∆x) − x
reordenando se tiene una expresión que nos permite calcular una aproximación a la raı́z
f (x)
x + ∆x = x − (1.49)
f ′ (x)
Teniendo en cuenta que repetiremos este cálculo varias veces, ahora es conveniente cambiar
a una notación más adecuada: x pasará a ser xn y nuestra primera aproximación x + ∆x se
designará con xn+1 , entonces tendremos

f (xn )
xn+1 = xn − (1.50)
f ′ (xn )
El algoritmo consiste en tomar un punto x0 como estimación inicial para la raı́z. Prefer-
entemente este valor debe ser próximo a la raı́z. Se denomina semilla o valor de prueba.
Se aplica la fórmula 1.50 para obtener un valor mejorado x1 . Ahora sucesivamente se utiliza
el último valor de x para reiniciar el proceso iterativo.
La figura 1.16 esquematiza el procedimiento.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 45

f (x) gi (x) = f (xi ) + f ′ (xi )x

g0 (x)
g1 (x)
∗
x3 x2 x1 x0

Figure 1.16: Esquema del proceso de solución con el algoritmo de Newton.

1.10.1 Algoritmo de Newton

Algoritmo

1. Datos: Establecer la función f (x), f ′ (x), ǫ, y un valor de prueba o semilla, x0 .

2. Si f ′ (x0 ) = 0, división por cero, salir.

3. Calcular
f (x0 )
x1 = x0 − (1.51)
f ′ (x0 )

4. Si |x1 − x0 | < ǫ, tomar α = x1 , salir.

5. Hacer x0 = x1

6. Volver al paso 2.

Obsérvese que no hace falta guardar todos los valores intermedios xn . Solo hace falta tener
registro de dos de ellos a los efectos de poder establecer la convergencia. No obstante, con
frecuencia, en las funciones y procedimientos de librerı́as matemáticas, como por ejemplo las
de Numerical Recipies ([13]), se guardan los valores de xi , f (xi ) en un arreglo para los fines
de impresión posterior, ya que la tarea de enviar un resultado a pantalla puede consumir
demasiado tiempo.

1.10.2 Ventajas y desventajas

El método de Newton tiene las siguientes ventajas

1. La convergencia es muy rápida (cuadrática) si se dispone de una semilla cercana a la

solución

2. Se puede calcular una cota para el error.

3. Se puede generalizar a problemas multidimensionales.

46

4. Se puede aplicar a funciones complejas.

Presenta las siguientes desventajas

1. No garantiza la convergencia, es decir el algoritmo puede diverger.

2. No tiene en cuenta los lı́mites de precisión de la computadora.

3. Es necesario conocer y calcular la derivada f ′ (x) en forma explı́cita. Esto puede ser
costoso computacionalmente o incluso imposible.

4. La convergencia deviene en lineal en la cercanı́a de raı́ces múltiples.

5. No existe cota práctica o rigurosa para el error.

Otros algoritmos relacionados Hausholder, Halley.

1.10.3 Resolución computacional

El programa newton-1.c reproducido en el listado 1.18 implementa el método de Newton-
Raphson para el caso de la función f (x) = x6 − x − 1.
Listing 1.18: newton-1.c
1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4

5 float ff ( float x ) ;
6 float dff ( float x ) ;
7

8 int main ()
9 {
10 // newton -1. c
11

12 // o b t i e n e una raiz de la f u n c i o n
13 // f ( x ) = x ^6 -x -1
14 // en el i n t e r v a l o [a - b ]=[1 -2]
15

16 int i =1;
17 float x0 =0. , x1 , x2 , fx1 , e p s i l o n =1 e -6 , sol = 1 . 1 3 4 7 2 4 1 3 8 4 0 1 5 2 L ;
18
19 x1 =2.;
20 while ( fabs ( x1 - x0 ) > e p s i l o n)
21 {
22 i = i +1;
23 x0 = x1 ;
24 x1 = x0 - ff ( x0 ) / dff ( x0 ) ;
25 printf (" % d % f % f \ n " , i , x1 , sol ) ;
26 }
27 return 0;
28 }
29 float ff ( float x )
30 {
31 return pow (x ,6.) -x -1.;
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 47

32 }
33 float dff ( float x )
34 {
35 return 6.* pow (x ,5.) -1;
36 }

La siguiente es la salida del programa anterior.

cardon@cardo:~/A-cursos/cursoMN/programas/raı́ces> ./new
2 1.680628 1.134724
3 1.430739 1.134724
4 1.254971 1.134724
5 1.161538 1.134724
6 1.136353 1.134724
7 1.134731 1.134724
8 1.134724 1.134724
cardon@cardo:~/A-cursos/cursoMN/programas/raı́ces>
Figure 1.17: Salida del programa newton-1.c

Para alcanzar la misma precisión el método de la bisección requiere 17 iteraciones.

1.11 Los métodos de la secante

Cuando no se puede obtener la derivada de la función para aplicar el método de Newton,
o se desconoce la función en forma explı́cita y sus valores se obtienen numéricamente, se
puede reemplazar la derivada f ′ (x) por la secante8 . Considere el intervalo [x0 , x1 ] en donde
se supone la existencia de la raı́z. La pendiente de la función en ese intervalo se puede
aproximar por la pendiente de la recta secante entre los puntos (x0 , f (x0 )) y (x1 , f (x1 ))
f (x1 ) − f (x0 )
f ′ (x) ≈ (1.52)
x1 − x0
además, la recta secante es colineal con la recta definida por (x1 , f (x1 )) y (x2 , f (x2 )). Este
último punto es donde la recta tangente interseca al eje de las abscisas (nótese que, por ello,
f (x2 ) = 0). Resulta la igualdad
f (x1 ) − f (x0 ) f (x1 ) − 0
= (1.53)
x1 − x0 x1 − x2
que nos permite obtener el algoritmo deseado
x1 − x0
x2 = x1 − f (x1 ) (1.54)
f (x1 ) − f (x0 )
Obtenida la aproximación x2 , podemos dejar de lado el valor de x0 para repetir el mismo
cálculo para el intervalo [x1 , x2 ], obteniéndose una nueva aproximación x3 . En general puede
repetirse para cualquier intervalo [xn , xn+1 ], por lo que, para este último caso se tiene
f (xn ) − f (xn−1 ) f (xn ) − 0
f ′ (x) ≈ = (1.55)
xn − xn−1 xn − xn+1
8
Esto implica reemplazar la derivada por una diferencia.
48

El algoritmo resultante es el siguiente

xn − xn−1
xn+1 = xn − f (xn ) (1.56)
f (xn ) − f (xn−1 )

El método de la secante puede interpretarse de la siguiente manera. Se toma dos valores

de prueba, x0 y x1 . Luego se traza la recta secante definida por los puntos correspondien-
tes, (x0 , f (x0 )) y (x1 , f (x1 )). La intersección de ella con las abscisas define el nuevo valor
aproximado para la raı́z, α ∼ x2 .
Para el próximo paso de la iteración se debe actualizar el intervalo de búsqueda. Una forma
muy elemental de hacerlo es la siguiente. Manteniéndose fijo x1 , se reemplaza x0 por x2 y se
repite el procedimiento hasta que la serie x2 , x3 , x4 ....xn converge hasta la precisión deseada,
es decir hasta cuando la diferencia entre dos de sus valores consecutivos es menor que la
tolerancia preestablecida. Este procedimiento se muestra esquemáticamente en la figura
1.18. El algoritmo es el siguiente.
•

f (x)

x0 x2 x3 x1
•

Figure 1.18: Esquema del proceso de solución con el algoritmo de la secante.

Algoritmo simple

Algoritmo

1. Datos. Establecer la función f (x), ε, y el intervalo inicial de búsqueda [x0 , x1 ]

2. Calcular
x1 − x0
x2 = x1 − f (x1 ) (1.57)
f (x1 ) − f (x0 )

3. Si |x2 − x0 | < ε, tomar α = x1 , salir.

4. Hacer x0 = x2

5. Volver al paso 2.

En este procedimiento se mantuvo fijo el punto x1 , f (x1 ), como se muestra en la figura 1.18
marcado con un punto negro. El otro extremo del intervalo, originalmente en x0 , se actualizó
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 49

por x2 , x3 , etc. Dadas las caracterı́sticas de la función, la raı́z estuvo siempre en el interior
del intervalo de búsqueda. Veremos en seguida que este procedimiento no es siempre el más
conveniente.

Convergencia

limn→∞ |α − cx+1 | ∼ k|α − c0 |1.618 (1.58)

1.618 la razón áurea9 .

1.11.1 Implementación computacional: método de la secante sen-

cillo
Programa 1.19
Listing 1.19: secante-b1.c
1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4 float f u n c i o n( float x ) ;
5 float s e c a n t e 1( float x1 , float x0 , float (* f u n c i o n) ( float ) , int
i m p r i m a t u r) ;
6
7 int main ()
8 {
9 // secante - b1 . c
10 int i =1;
11 float x0 , x1 , c ;
12 x1 =2.;
13 x0 =1.;
14 c = s e c a n t e 1( x1 , x0 , funcion ,1) ;
15 printf (" Raiz =% f f ( raiz ) =% g \ n " ,c , f u n c i o n( c ) ) ;
16 return 0;
17 }
18 float f u n c i o n( float x )
19 {
20 return pow (x ,6.) -x -1.;
21 }

La subrutina secante1 en el archivo lib-secante-1.c se muestra en el listado 1.20

Listing 1.20: lib-secante-1.c
1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4
5 float s e c a n t e 1( float x1 , float x0 , float (* f ) ( float ) , int i m p r i m a t u r)
6 {
7 int i =1;
9
Atkinson AINA, pp52
50

8 float x2 , fx1 , fx0 , e p s i l o n =1 e -4;

9
10

11 fx1 =(* f ) ( x1 ) ;
12 fx0 =(* f ) ( x0 ) ;
13

14 x2 = x1 - fx1 * ( x1 - x0 ) /( fx1 - fx0 ) ;

15 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i x1 x0 x2
fabs ( x2 - x0 ) \ n ") ;
16
17 while ( fabs ( x2 - x0 ) > e p s i l o n)
18 {
19 i = i +1;
20 x0 = x2 ;
21 fx0 =(* f ) ( x0 ) ;
22 x2 = x1 - fx1 * ( x1 - x0 ) /( fx1 - fx0 ) ;
23 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i , x0 , x1 , x2
, fabs ( x2 - x0 ) ) ;
24 }
25 return x2 ;
26 }

La siguiente es la salida del programa anterior:

./bi
# i x1 x0 x2 fabs(x2-x0)
2 1.01613 2.00000 1.03067 0.015
3 1.03067 2.00000 1.04372 0.013
4 1.04372 2.00000 1.05535 0.012
...
34 1.13379 2.00000 1.13392 0.00014
35 1.13392 2.00000 1.13404 0.00012
36 1.13404 2.00000 1.13414 0.0001
Raiz=1.134140 f(raiz)=-0.00600566

Observe que el método converge algo más lento que el método de Newton. Se hubiera
podido mantener fijo el otro extremo, el punto x0 , para luego de obtenida la aproximación
x2 , actualizar el intervalo de búsqueda reemplazando x1 por x2 . Este procedimiento se
esquematiza en la figura 1.19. Puede observarse que en este caso que, dadas las caracterı́sticas
de la función, el intervalo de búsqueda no siempre encierra la raı́z buscada. Aparece como una
propiedad deseable que la raı́z se encontrara siempre contenida en el intervalo de búsqueda.
Esto asegurarı́a una convergencia monótona a la raı́z. Para ello basta elegir adecuadamente,
en cada paso, el nuevo extremo del dominio de búsqueda. Podemos usar para ello el mismo
algoritmo que en el método de la bisección. La combinación del método de la tangente con
una selección del intervalo de búsqueda que mantenga encerrada a la raı́z se denomina el
método de la Regula Falsi.

1.11.2 Método de la Falsa Posición

El método de la Posición Falsa, o Regula Falsi, es básicamente el método de la secante en
el que se tiene cuidado de mantener la raı́z encorchetada entre los lı́mites del intervalo de
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 51

f (x)

x0 x2 x4 x3 x1
•

•
Figure 1.19: Esquema del proceso de solución con el algoritmo de falsa posición.

búsqueda. Consideraremos dos modos de hacerlo, en el primero, prácticamente construire-

mos el mismo algoritmo que en el caso de la bisección aunque cambiaremos la fórmula de
aproximación la raı́z. En el otro describiremos el algoritmo como lo hace Press et al. en
Numerical Recipies.
Para describirlo utilizaremos la notación empleada para describir el método de la bisección.
[a, b] será el intervalo de búsqueda y c nuestra aproximación, que ahora calcularemos según

(a − b)
c = b + fb (1.59)
(fa − fb )
que no es otra cosa que la misma ecuación 1.56

Algortimo: Regula Falsi y bisección

El algoritmo es el siguiente
Algoritmo

1. Datos: Definir la función f (x) y los valores de a, b y ε.

2. Calcular el valor de c con

(a − b)
c = b − fb (1.60)
(fa − fb )

3. Si |b − c| ≤ ε, tomar c como raı́z, imprimir, terminar.

4. Si f (c)f (b) < 0 hacer a = c, en caso contrario hacer b = c

5. Volver al paso 2.

El programa regulafalsi-b1.c reproducido en el listado 1.23 implementa la forma simple

del método e la Falsa Posición.
Listing 1.21: regulafalsi-b1.c
1 # i n c l u d e < stdio .h >
52

2 # i n c l u d e < stdlib .h >

3 # i n c l u d e < math .h >
4 # define MAX 50
5

6 float f u n c i o n( float x ) ;
7 float r e g u l a f a l s i 1 ( float x , float B , float (* f ) ( float ) , int i m p r i m a t u r )
;
8
9 int main ()
10 {
11 // regulafalsi -1. c
12 // o b t i e n e una raiz de la f u n c i o n
13 // f ( x ) = x ^6 - x -1
14 // en el i n t e r v a l o [a - b ]=[1 -2]
15

16
17 float a ,b , c ;
18 a =1.;
19 b =2.;
20

21 // C a l c u l a la ff f ( x ) vs . x
22 // para e x p l o r a r g r a f i c a m e n t e
23

24 // for ( x =0.; x <= b ; x = x +0.1)

25 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
26 // return ;
27

28 c = r e g u l a f a l s i 1 (a ,b , funcion ,1) ;
29 printf (" Raiz =% f f ( raiz ) =% g \ n " ,c , f u n c i o n( c ) ) ;
30 return 0;
31 }
32
33 float f u n c i o n( float x )
34 {
35 return pow (x ,6.) -x -1.;
36 }

El programa usa la funcion regulafalsi1 de la biblioteca lib-raices.c mostrada en el

listado 1.22

Listing 1.22: regulafalsi1

1 float r e g u l a f a l s i 1 ( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r )
2 {
3 int i =0;
4 float c , fa , fb , fc , e p s i l o n =1 e -4;
5 fa =(* f ) ( a ) ;
6 fb =(* f ) ( b ) ;
7 if ( fa * fb >0) printf (" No hay raiz en el i n t e r v a l o \ n ") ;
8 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i a b c
fabs (b - c ) \ n ") ;
9
10 do
11 {
12 c = b - fb *( b - a ) /( fb - fa ) ;
13 fc =(* f ) ( c ) ;
14 if ( fb * fc < 0. )
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 53

15 {a=c;
16 fa =(* f ) ( a ) ;
17 } else
18 {b=c;
19 fb =(* f ) ( b ) ;
20 }
21 i = i +1;
22 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i ,a ,b ,c ,
fabs (b - c ) ) ;
23 }
24 while ( fabs ((* f ) ( c ) ) > e p s i l o n && i < MAX ) ;
25 return c ;
26 }

Implementación de Numerical Recipies

A efectos de comparación presentamos aquı́ una versión ajustada a nuestra notación y algo
simplificada de la implementación de Press et al., [13]. Se recomienda al lector usar la
implementación de los autores.
Se parte del intervalo inicial de búsqueda es [x0 , x1 ]. tal que f (x0 )f (x1 ) < 0. Se hace
a = min(f (x0 ), f (x1 ) y b = max(f (x0 ), f (x1 )), con lo que se consigue que en el nuevo
intervalo de búsqueda [a, b], f (a) sea siempre el valor negativo y f (b) el positivo. En cada
iteración se cambiará a o b por c, de manera que se mantenga ...
El método de la Posición Falsa o Regula Falsi, se caracteriza por mantener la raı́z acorralada
entre los lı́mites del intervalo de búsqueda.
Para implementarlo es conveniente usar la siguiente notación. Las coordenadas que definen
el intervalo, que fueron denominadas x0 y x1 será designadas ahora por xa o xb según el valor
que tome la función f () sea alto o bajo (mayor que cero o menor que cero)). Si xi , con i
igual a 1 o 0, es tal que f (xi ) > 0 , designaremos xa = xi y f (xa ) = fa . El otro valor será
designado con el subı́ndice b.
Con esta notación, el algoritmo representado por la ecuación 1.54 se puede escribir como

(xa − xb )
x = xb + fb (1.61)
(fa − fb )

El algoritmo de solución es el siguiente

 54

Regula Falsi. Algoritmo

Algoritmo

1. Datos: Establecer la función f (x), el intervalo [x0 , x1 ], y la tolerancia ε.

2. Determinar si en el intervalo [x0 , x1 ] existe una raı́z. Si no existe, abortar.

3. Determinar cual de los valores x0 y x1 corresponden a xa y xb

4. Calcular
(xa − xb )
x = xb + fb (1.62)
(fa − fb )

5. Si |x − xb | < ε, tomar α = x, salir.

6. Actualizar el intervalo de búsqueda. Si f (x)f (xa ) < 0 hacer xb = x, caso

contrario hacer xa = x

7. Volver al paso 4.

Para determinar cual de los valores x0 y x1 corresponden a xa y xb implementamos el siguiente

algoritmo
Algoritmo

1. Suponemos fa = f (x0 ) y fb = f (x1 ) y calculamos fa y fb .

2. Si fb < 0 la elección fue correcta.

3. Si fb > 0 intercambiar xa por xb y recalcular las funciones fa y fb .

En la siguiente implementación se ha evitado recalcular las funciones llamadas en la lı́nea 1

y 2, por ello el intercambio (swap) de las lı́neas 12 a 14

1 //
2 fb=funcion(x1);
3 fa=funcion(x2);
4 if(fb<0)
5 {
6 xb=x1;
7 xa=x2;
8 }
9 else
10 {
11 xb=x2;
12 xa=x1;
13 swap=fb;
14 fb=fa;
15 fa=swap;
16 }
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 55

Ventajas y desventajas

Ventajas

• No requiere el cálculo de la derivada.

• Es rápido. Razón de convergencia Q = 1.6.

• Convergencia más que lineal cerca de la raı́z.

• Convergencia lineal cerca de raı́ces múltiples.

Desventajas

• La iteración puede diverger.

• No hay cota práctica para el error.

Convergencia

f (x) = x6 − x − 1 f (x) = x6 − x − 1
1.5 2
a a
1.4 c 1.8 c
1.3 b 1.6 b
1.2 1.4
1.1 1.2
1 1
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20

Figure 1.20: Evolución del entorno de búsqueda y de la aproximación de la raı́z: a) bisección,

b) Régula Falsi. Observe que en el segundo caso, el algoritmo trabaja como la secante, sin
cambio de punto de referencia.

Implementación computacional. Falsa Posición.

El programa regulafalsi-b1.c listado en 1.23, implementa el método de la régula falsi a

la manera de Press et al. ([13]). Para el caso mostrado, ejecuta las mismas operaciones que
el programa simplificado precedente.
Listing 1.23: regulafalsi-b1.c
1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4 # define MAX 50
5

6 float f u n c i o n( float x ) ;
56

7 float r e g u l a f a l s i 1 ( float x , float B , float (* f ) ( float ) , int i m p r i m a t u r )

;
8

9 int main ()
10 {
11 // regulafalsi -1. c
12 // o b t i e n e una raiz de la f u n c i o n
13 // f ( x ) = x ^6 - x -1
14 // en el i n t e r v a l o [a - b ]=[1 -2]
15
16

17 float a ,b , c ;
18 a =1.;
19 b =2.;
20

21 // C a l c u l a la ff f ( x ) vs . x
22 // para e x p l o r a r g r a f i c a m e n t e
23
24 // for ( x =0.; x <= b ; x = x +0.1)
25 // printf ("% f % f \ n " ,x , ff ( x ) ) ;
26 // return ;
27

28 c = r e g u l a f a l s i 1 (a ,b , funcion ,1) ;
29 printf (" Raiz =% f f ( raiz ) =% g \ n " ,c , f u n c i o n( c ) ) ;
30 return 0;
31 }
32

33 float f u n c i o n( float x )
34 {
35 return pow (x ,6.) -x -1.;
36 }

La función que implementa el algoritmo, regulafasi1 se incluye en el archivo lib-regulafalsi-1.c

listada en 1.24

Listing 1.24: lib-regulafalsi-1.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4 # define MAX 50
5
6 float r e g u l a f a l s i 1 ( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r )
7 {
8 int i =0;
9 float c , fa , fb , fc , e p s i l o n =1 e -4;
10 fa =(* f ) ( a ) ;
11 fb =(* f ) ( b ) ;
12 if ( fa * fb >0) printf (" No hay raiz en el i n t e r v a l o \ n ") ;
13 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i a b c
fabs (b - c ) \ n ") ;
14
15 do
16 {
17 c = b - fb *( b - a ) /( fb - fa ) ;
18 fc =(* f ) ( c ) ;
19 if ( fb * fc < 0. )
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 57

20 {a=c;
21 fa =(* f ) ( a ) ;
22 } else
23 {b=c;
24 fb =(* f ) ( b ) ;
25 }
26 i = i +1;
27 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i ,a ,b ,c ,
fabs (b - c ) ) ;
28 }
29 while ( fabs ((* f ) ( c ) ) > e p s i l o n && i < MAX ) ;
30 return c ;
31 }

La salida del programa es la siguiente:

./bi
# i a b c fabs(b-c)
1 1.01613 2.00000 1.01613 0.98
2 1.03067 2.00000 1.03067 0.97
3 1.04372 2.00000 1.04372 0.96
4 1.05535 2.00000 1.05535 0.94
...
48 1.13464 2.00000 1.13464 0.87
49 1.13465 2.00000 1.13465 0.87
50 1.13466 2.00000 1.13466 0.87
Raiz=1.134660 f(raiz)=-0.000662124

Es interesante notar cuan rápido converge cuando se intercambian los lı́mites de la búsqueda.
58

1.12 Métodos iterativos de punto fijo

Todos los métodos para encontrar raı́ces proceden mediante repetición o iteración del al-
goritmo básico. Algunos de estos métodos pertenecen a una clase especial denominados
métodos iterativos de punto fijo o también métodos de aproximaciones sucesivas.
Se denomina punto fijo a las solución x = p que satisface la ecuación x = f (x). Consider-
emos la ecuación
x = g(x) (1.63)
cada solución, o valor de x que satisfaga la ecuación se denomina un punto fijo de la función
g(x). El problema de encontrar sus soluciones se denomina un problema de punto fijo.
En los casos precedentes, nuestro problema fue encontrar los ceros de la ecuación f (x) = 0.
En muchas ocasiones, es posible y conveniente reescribir la ecuación f (x) de manera que
el problema de encontrar sus ceros pueda formularse como un problema de punto fijo. Por
ejemplo, para encontrar los ceros o raı́ces de la función f (x) = x2 − a = 0, con a < 0, se
pueden plantear los siguientes problemas de punto fijo

1. x = g1 (x) = x2 + x − a
2. x = g2 (x) = x/a
3. x = g3 (x) = 1/2(x + x/a)

Los problemas de punto fijo se pueden resolver por métodos iterativos de aproximaciones
sucesivas. En general un método de aproximaciones sucesivas es aquel que permite encon-
trar las solución por aplicación repetida del siguiente algoritmo
xn+1 = g(xn ) (1.64)
Se comienza calculando g(x0 ), donde x0 es un valor arbitrario pero preferentemente cercano
a la solución. Este valor se denomina usualmente semilla o valor de prueba. Se obtiene
x1 = g(x0 ). Luego se reemplaza x0 por x1 y se repite el cálculo. Se genera una secuencia
{xn }∞
n que bajo ciertas condiciones es convergente y conduce a encontrar la solución, es decir
el punto fijo.
Se puede demostrar10 (véase [6], [8]) que para una función g(x), continua en el intervalo [a, b],
diferenciable en (a, b), dada una semilla x0 , la aplicación de un proceso iterativo de aproxi-
maciones sucesivas como el indicado por la ecuación 1.64, conducirá a una serie convergente
si se cumple que en la cercanı́a del punto fijo es
Maxa≤x≤b |g ′(x)| < 1 (1.65)

1.12.1 Interpretación geométrica

El principio en que se basan estos métodos se puede mostrar convenientemente siguiendo el
cálculo en forma gráfica. Consideremos el caso el sistema de ecuaciones
y = g(x)
(1.66)
y=x
10
Caso particular del teorema de punto fijo de Banach
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 59

Evidentemente es equivalente a escribir x = g(x).

Podemos reescribir nuevamente como las siguientes dos funciones
y = g(x) (1.67)
cuya gráfica será en general, una curva cualquiera, y
y=x (1.68)
cuya gráfica es una recta. La solución se encontrará en la intersección de ambas gráficas,
como se muestra en la figura 1.21.
a) b) y = g(x)
y=x y=x
y = g(x)

x0 x1 x2 x3 x3 x2 x1 x0

Figure 1.21: Algoritmo iterativo de aproximaciones sucesivas a) 0 < g ′(α) < 1, convergente.
b) g ′(α) > 1, divergente
a) b)
y=x y=x

y = g(x) y = g(x)

x0 x2x3 x1 x3 x0 x1

Figure 1.22: Algoritmo iterativo de aproximaciones sucesivas a) −1 < g ′ (α) < 0, convergente.
b) g ′(α) < −1, divergente

Para obtener el valor de la intersección iterativamente podemos seguir el procedimiento

siguiente. Se calcula la ecuación 1.67 para un valor de prueba x0 , obteniéndose
y0 = g(x0 ) (1.69)
60

Ahora se introduce el valor de y0 en la inversa de la ecuación 1.68 para obtener x1 ,

x1 = y0 (1.70)

Luego se introduce y1 en la 1.67 y se recomienza el ciclo. La figura 1.21 muestra el procedi-

miento en forma gráfica. El mismo queda expresado por

xn+1 = g(xn ) (1.71)

donde n es un ı́ndice que indica el número de iteraciones. El algoritmo puede darse por
convergido cuando
|xn+1 − xn | < ε (1.72)
o cuando
|xn+1 − xn |
<ε (1.73)
|xn |

1.12.2 Algoritmo
El algoritmo puede escribirse de la siguiente manera
Algoritmo
1. Se propone un valor de prueba x0

2. Se obtiene una aproximación mejorada con

x1 = g(x0 ) (1.74)

3. Si |x1 − x0 | < ε, terminar

4. Se hace x0 = x1

5. Se vuelve al paso 2

1.13 Ceros de sistemas de ecuaciones acopladas

El algoritmo pude aplicarse para la resolución de problemas en donde se tiene varias ecua-
ciones acopladas.
x1 = F (x2 , x3 , ..., xN , otros parámetros) (1.75)
x2 = G(x1 , , x3 , ..., xN , otros parámetros) (1.76)
x3 = H(x1 , , x3 , ...xN , otros parámetros) (1.77)
Este tipo de problemas es frecuente en la resolución de circuitos eléctricos, circuitos térmicos,
redes hidráulicas entre muchos otros problemas de ingenierı́a.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 61

1.14 Problemas resueltos

1.14.1 Problema 1
Las ecuaciones paramétricas del vuelo de un proyectil son:

A) para el caso ideal (sin resistencia del aire):

x = vx0 t (1.78)
1
y = vy0 t − gt2 (1.79)
2
B) Para el caso de resistencia lineal:
m k
x= vx0 (1 − e− m t ) (1.80)
k
mgt m mg k
y=− + (vy0 + )(1 − e− m t ) (1.81)
k k k
C) Para el caso de resistencia cuadrática:
1
x= (ln(1 − kvx0 t) (1.82)
kvx0
g 1 gt2 gt
y = (vy0 − ) ln(1 − kvx0 t) − − (1.83)
2kvx0 kvx0 4 2kvx0

Datos:
m = 1kg
g = 9.81ms−2
v0 = 9.81ms−1 , velocidad inicial
θ = 20o
k = 0.1m−1

Consigna:

• Programe funciones de C para calcular las coordenadas x e y de la trayectoria de

ambos modelos de vuelo de proyectil, A, B o C. Genere una tabla tal que permita
graficar con gnuplot las dos trayectorias simultáneamente, desde el origen, en (x = 0,
y = 0) hasta que el proyectil cae al suelo en (x = R =?,y = 0).

• Redireccione la tabla, desde la shell, a un archivo [Link] y haga un dibujo con

gnuplot de ambas trayectorias superpuestas.

• Imprima el alcance del proyectil, en cada caso.

• Utilice los dos canales disponibles en C, stderr y stdout, para imprimir la tabla y los
alcances por vı́as diferentes. Use el redireccionamiento de la shell para enviar la tabla
al archivo [Link] y el alcance al archivo [Link].
62

Resolución

El programa 1.25 resuelve el problema propuesto.

Listing 1.25: biseccion-parcial-1-2014-1.c

1 # i n c l u d e < stdio .h >
2 # i n c l u d e < stdlib .h >
3 # i n c l u d e < math .h >
4
5 float ff ( float R ) ;
6

7 int main ()
8 {
9 // biseccion - parcial -1 -2014 -1.c
10 // a l c a n c e de un p r o y e c t i l
11 // r e s i s t e n c i a c u a d r a t i c a
12 // Parker , 1997 , W a r b u r t o n et al . , 2010
13
14 int i =0;
15 float a ,b ,c ,R , e p s i l o n =1e -6;
16

17 a =1.;
18 b =20.;
19

20 // C a l c u l a la ff f ( c ) vs . c para g r a f i c a r
21 // for ( c =0.; c <= b ; c = c +0.1)
22 // printf ("% f % f \ n " ,c , ff ( c ) ) ;
23 // return ;
24

25 // existe la raiz en el i n t e r v a l o
26 if ( ff ( a ) * ff ( b ) >0)
27 printf (" No hay raiz en el i n t e r v a l o \ n ") ;
28

29 c =( a + b ) /2.;
30 while ( fabs (b - c ) > e p s i l o n)
31 {
32 if ( ff ( b ) * ff ( c ) < 0. )
33 a=c;
34 else
35 b=c;
36

37 c =( a + b ) /2.;
38 i = i +1;
39 printf ("%2 d , %1.5 f %1.2 g \ n " ,i ,c , fabs (b - c ) ) ;
40 }
41 return 0;
42
43 }
44

45 float ff ( float R )
46 {
47 float v , vx0 , vy0 ,g , k ;
48 v =9.81; // modulo de la v e l o c i d a d i n i c i a l
49 k =0.1; // c o e f i c i e n t e de a r r a s t r e
50 g =9.81; // a c e l e r a c i o n de la g r a v e d a d
51 theta =20.:
52 vx0 = v * cos ( theta * 3 . 1 4 1 5 / 1 8 0 . ) ;
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 63

53 vy0 = v * sin ( theta * 3 . 1 4 1 5 / 1 8 0 . ) ;

54
55 return ( vy0 + g /(2.* k * vx0 ) ) *( R / vx0 ) -( g /(4* k * k * vx0 * vx0 ) ) *( exp (2* k * R ) -1.) ;
56 }

1.14.2 Problema 2
El alcance R de un proyectil cuando se considera resistencia del aire cuadrática, ha sido
aproximado por Parker (1977) y recientemente por Warburton et al. (2010), para el caso
especial en que vx0 >> vy0 , es decir trayectorias de ángulo bajo, por la siguiente ecuación
trascendental  
g R g 2kR


vy0 + − e − 1 =0 (1.84)
2kvx0 vx0 4(kvx0 )2

Datos:
Use mismos coeficientes que en el problema anterior.

Consigna: Obtenga el alcance R:

Utilice el método de la bisección. Programe de tal manera que la función sólo se calcule una
vez por ciclo. Utilice un lazo do while.

Resolución

El programa 1.26 resuelve el problema propuesto.

Listing 1.26: biseccion-parcial-1-2014-2b.c

1 # i n c l u d e " lib - raices . h "
2

3 float ff ( float R ) ;
4 int main ()
5 {
6 // biseccion - parcial -1 -2014 -2b . c
7

8 int i =0;
9 float a ,b ,c , R ;
10

11 a =1.;
12 b =20.;
13

14 // C a l c u l a la ff f ( c ) vs . c para g r a f i c a r
15 // for ( c =0.; c <= b ; c = c +0.1)
16 // printf ("% f % f \ n " ,c , ff ( c ) ) ;
17
18 // return ;
19 c = b i s e c c i o n 7(a ,b , ff ,1) ;
20 printf (" Raiz =% f f ( raiz ) =% g \ n " ,c , ff ( c ) ) ;
21 return 0;
22 }
23

24 float ff ( float R )
64

25 {
26 float v , vx0 , vy0 ,g ,k , theta ;
27 v =9.81;
28 k =0.1;
29 g =9.81;
30 theta =20;
31 vx0 = v * cos ( theta * 3 . 1 4 1 5 / 1 8 0 . ) ;
32 vy0 = v * sin ( theta * 3 . 1 4 1 5 / 1 8 0 . ) ;
33 return ( vy0 + g /(2.* k * vx0 ) ) *( R / vx0 ) -( g /(4.* k * k * vx0 * vx0 ) ) *( exp (2.* k * R )
-1.) ;
34 }

El programa hace uso de la función biseccion7 incluida en la librerı́a lib-raices

Listing 1.27: biseccion7

1 float b i s e c c i o n 7( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r)
2 {
3 float c ;
4 int i =0 , imax =100;
5 float e p s i l o n =1.e -4;
6 float fa , fb , fc ;
7

8 // existe la raiz en el i n t e r v a l o
9 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "# i a b c
fabs (b - c ) \ n ") ;
10
11 if ((* f ) ( a ) * (* f ) ( b ) >0)
12 { f p r i n t f( stderr , "# No hay raiz en el i n t e r v a l o\ n ") ;
13 return 0;}
14

15 fb =(* f ) ( b ) ;
16 fa =(* f ) ( a ) ;
17 c =( a + b ) /2.;
18 do
19 {
20 fc =(* f ) ( c ) ;
21
22 if ( fb * fc < 0. )
23 {a=c;
24 fa = fc ;
25 }
26 else
27 {b=c;
28 fb = fc ;}
29
30 c =( a + b ) /2.;
31 i = i +1;
32 if ( i m p r i m a t u r ==1) f p r i n t f( stderr , "%4 d %9.5 f %9.5 f %9.5 f %9.2 g \ n " ,i ,a ,b ,c ,
fabs (b - c ) ) ;
33 }
34 while ( fabs (b - c ) > e p s i l o n) ;
35 return c ;
36 }
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 65

1.14.3 Problema 3
Obtenga la raı́z de sen(x) = 0 mediante el método de Newton Raphson y como problema de
punto fijo (aproximaciones sucesivas). Ayuda: El algoritmo de Newton Raphson es

f (xn )
xn+1 = xn − (1.85)
f ′ (xn )

1.14.4 Problema 4
Un objeto cae verticalmente en aire sujeto a la fuerza de la gravedad y a la fuerza viscosa.
La altura del objeto, h, en función del tiempo, t, está dada por

mg m2 g
h(t) = h0 − t + 2 (1 − e−kt/m ) (1.86)
k k
Datos g = 9.8m/s2 , m = 0.1kg, h0 = 100m, k = 1.5kgs/m.
Calcule cuánto tarda el cuerpo en caer al suelo con una exactitud de 0.01s.
Elija usted el método de resolución.
66

1.15 Una biblioteca personal

Las distintas funciones que se han presentado hasta ahora, lib-biseccion, lib-secante,
lib-newton, etc., se han agregado al archivo lib-raices.c que hará las veces de biblioteca11.
Este archivo viene acompañado por su archivo de encabezados lib-raices.h y ambos se
pueden utilizar para poner a disposición cualquiera de las funciones desarrolladas. Para ello,
debe incluirse en el programa principal el archivo de encabezados,

#include "lib-raices.h"

que reemplaza a los prototipos de la funciones individuales que incluı́amos antes, ya que el
archivo lib-raices.h hará la inclusión de todos los prototipos de todas las funciones. Eso
es suficiente para hacer accesibles al programa todas las funciones del archivo lib-raices.c
El programa debe compilarse ahora como

gcc -o no nombre.c lib-raices.c -lm

El archivo nombre.c, como siempre, debe estar en el directorio de trabajo asi como también
los dos nuevos archivos, lib-raices.c y lib-raices.h que componen nuestra biblioteca12.
El listado del archivo de encabezados lib-raices.h, listado 1.28, sirve como manual de uso
de las funciones de la biblioteca personal. Cada prototipo indica el retorno de la función y
los datos o argumentos que requiere y su tipo:
Listing 1.28: lib-raices.h
1 # include < stdio .h >
2 # include < stdlib .h >
3 # include < math .h >
4 # include < tgmath .h >
5
6 # define MAX 500
7

8 float b i s e c c i o n( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r ) ;

9 float b i s e c c i o n 6( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r) ;
10 float b i s e c c i o n 7( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r) ;
11
12 float r e g u l a f a l s i 1 ( float a , float b , float (* f ) ( float ) , int i m p r i m a t u r) ;
13 float s e c a n t e 1( float x1 , float x0 , float (* f ) ( float ) , int i m p r i m a t u r) ;

C dispone de herramientas para hacer de estos archivos una verdadera biblioteca (library
en inglés), de tipo estática o dinámica. Las bibliotecas estáticas o dinámicas precompilan
las funciones del archivo y las ponen a disposición para ser usadas en distintas etapas del
ensamblado de un programa o en la ejecución del mismo.
Las bibliotecas estáticas, cuyos archivos se distinguen con la extensión .a, consisten en un
archivo objeto que se ensambla con el objeto resultante del programa fuente en la etapa de
11
El término en inglés es library, que significa biblioteca. Por asociación sonora con el término español
librerı́a, con frecuente, y es probable que lo encuentre ası́ en este mismo libro, usamos este vocablo para
referirnos a una biblioteca.
12
Estos archivos, como se describe más abajo, nos son técnicamente bibliotecas de C, no obstante cumplen
la misma función y pueden transformarse en ellas.
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 67

compilación para producir un único ejecutable o binario. Las funciones usadas de la librerı́a
pasan a formar parte del binario ası́ creado. Los ejecutables son algo más grandes que los
de los programas ensamblados con librerı́as dinámicas.
Los objetos de las bibliotecas dinámicas, cuyos archivos se distinguen con la extensión .so,
se unen al programa del usuario recién en la etapa de ejecución. Se denominan también
bibliotecas compartidas por que existe una única versión objeto de la misma, reutilizable
por diferentes programas a la vez. Las bibliotecas dinámicas se vinculan en dos etapas. En
una primera etapa, durante la compilación, el ensamblador hace diversas verificaciones, en
particular que las funciones que el programa usa estén en el programa o en la biblioteca. Los
objetos de funciones de librerı́a no se ensamblan al programa. En una segunda etapa, durante
la ejecución un cargador sube estos objetos a memoria y los hace disponibles al programa
ejecutable. Tienen la ventaja de que existe un único objeto que puede ser compartido por
diversos programas.

1.16 Bibliotecas numéricas

El propósito de las páginas precedentes ha sido introducirnos en la programación en C y su
aplicación a la implementación de algoritmos básicos. Los programas dados como ejemplos
muestran una secuencia casi literal de la formulación matemática del algoritmo, cercana a
la primera interpretación que puede realizar un estudiante. Por supuesto no pretenden ser
su mejor implementación.
Para resolver problemas de producción es conveniente, en todo caso que sea posible, utilizar
funciones de biblioteca de reconocida calidad, que las hay en el dominio público y comerciales.
Comentaremos aquı́ dos de estas bibliotecas:

• Numerical Recipies: Esta es la biblioteca que acompaña a la serie de libros Numer-

ical Recipies, de Press et al., [13]. Estos libros, con sus correspondientes biblioteca
se han escrito para varios lenguajes de programación, C, C++, Fortran 77, Fortran
90, Java, Basic, Modula 2 y Lisp. Las funciones de Numerical Recipies C++ pueden
llamarse desde Python. La última versión del libro, la tercera, corresponde a C++,
fue publicada por Cambridge University Press en 2007.

• GNU Scientific Library: La GNU Scientific Library, [9], es una biblioteca en el

dominio público, abierta, mantenida por la Fundación GNU. Provee una amplia gama
de rutinas matemáticas, como generadores de números aleatorios, funciones especiales
y ajuste de mı́nimos cuadrados. Hay más de 1000 funciones en total con un amplio
conjunto de pruebas. Tiene implementados varios métodos iterativos para la búsqueda
de raı́ces.

Existen varias otras librerı́as para C, C++, un listado puede encontrarse en la Wikipedia,
[2]. No debe perderse de vista que existen importantes librerı́as para álgebra numérica en
otros lenguajes, en particular Fortran. Una fuente importante de información puede hallarse
en la Wikipedia, [2]. Las librerı́as comerciales más conocidas son HLS, antes Hardwell,
[3], NAG, [10], IMSL, [1] que disponen versiones para los lenguajes de interés cientı́fico
principales: fortran, C, Java, Python, etc..
68

1.16.1 Un ejemplo de aplicación con la GNU Scientific library

La GNU Scientific library, o GSL es una colección de software para cálculos numéricos en
matemática aplicada y ciencia. Esta escrito en C y existen wrapers que permiten llamarlas
desde otros lenguajes. Es parte del Proyecto GNU y se distribuye bajo la Licencia Pública
General de GNU. La página oficial de la librerı́a es [Link]
Una descripción más detallada se dará en el apéndice ??.
Para la obtención de raı́ces, la librerı́a contiene varios métodos

• Newton: gsl_root_fdfsolver_newton

• Secante: gsl_root_fdfsolver_secant

• Steffenson: gsl_root_fdfsolver_steffenson

• Brent: gsl_root_fsolver_brent

• Bisección: gsl_root_fsolver_bisection

El uso de esta librerı́a requiere conocimientos algo más avanzados de C que los utilizados
hasta aquı́. En particular el uso de estructuras y punteros
Para compilar se necesita el archivo de encabezados demo_fn.h que se puede bajar del sitio
web de GSL:
Listing 1.29: demofn.h
1 // d e m o _ f n. h
2 struct q u a d r a t i c _ p a r a m s
3 {
4 double a , b , c ;
5 };
6

7 double q u a d r a t i c ( double x , void * params ) ;

8 double q u a d r a t i c _ d e r i v ( double x , void * params ) ;
9 void q u a d r a t i c _ f d f ( double x , void * params ,
10 double *y , double * dy ) ;

El programa demofn.c implementa la definición de la función a integrar.

Listing 1.30: demofn.c
1 // d e m o _ f n. c
2 double q u a d r a t i c ( double x , void * params )
3 {
4 struct q u a d r a t i c _ p a r a m s * p
5 = ( struct q u a d r a t i c _ p a r a m s *) params ;
6
7 double a = p - >a ;
8 double b = p - >b ;
9 double c = p - >c ;
10

11 return ( a * x + b ) * x + c ;
12 }
13
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 69

14 double q u a d r a t i c _ d e r i v ( double x , void * params )

15 {
16 struct q u a d r a t i c _ p a r a m s * p
17 = ( struct q u a d r a t i c _ p a r a m s *) params ;
18
19 double a = p - >a ;
20 double b = p - >b ;
21
22 return 2.0 * a * x + b ;
23 }
24

25 void q u a d r a t i c _ f d f ( double x , void * params ,

26 double *y , double * dy )
27 {
28 struct q u a d r a t i c _ p a r a m s * p
29 = ( struct q u a d r a t i c _ p a r a m s *) params ;
30

31 double a = p - >a ;
32 double b = p - >b ;
33 double c = p - >c ;
34
35 *y = (a * x + b) * x + c;
36 * dy = 2.0 * a * x + b ;
37 }

El programa brent.c encuentra los ceros de resuelve x2 −5 = 0 con gsl_root_fsolver_brent

Listing 1.31: brent.c
1 # include < stdio .h >
2 # include < gsl / g s l _ e r r n o.h >
3 # include < gsl / g s l _ m a t h.h >
4 # include < gsl / g s l _ r o o t s.h >
5

6 # i n c l u d e " d e m o _ f n. h "
7 # i n c l u d e " d e m o _ f n. c "
8
9 int main ( void )
10 {
11 int status ;
12 int iter = 0 , m a x _ i t e r = 100;
13 const g s l _ r o o t _ f s o l v e r _ t y p e * T ;
14 gsl_root_fsolver *s;
15 double r = 0 , r _ e x p e c t e d = sqrt (5.0) ;
16 double x_lo = 0.0 , x_hi = 5.0;
17 gsl_function F;
18 struct q u a d r a t i c _ p a r a m s params = {1.0 , 0.0 , -5.0};
19
20 F . f u n c t i o n = & q u a d r a t i c;
21 F . params = & params ;
22
23 T = gsl_root_fsolver_brent ;
24 s = gsl_root_fsolver_alloc (T);
25 g s l _ r o o t _ f s o l v e r _ s e t (s , &F , x_lo , x_hi ) ;
26

27 printf (" using % s method \ n " ,

28 gsl_root_fsolver_name (s));
29
70

30 printf ("%5 s [%9s , %9 s ] %9 s %10 s %9 s \ n " ,

31 " iter " , " lower " , " upper " , " root " ,
32 " err " , " err ( est ) ") ;
33

34 do
35 {
36 iter ++;
37 status = g s l _ r o o t _ f s o l v e r _ i t e r a t e ( s ) ;
38 r = gsl_root_fsolver_root (s);
39 x_lo = g s l _ r o o t _ f s o l v e r _ x _ l o w e r ( s ) ;
40 x_hi = g s l _ r o o t _ f s o l v e r _ x _ u p p e r ( s ) ;
41 status = g s l _ r o o t _ t e s t _ i n t e r v a l ( x_lo , x_hi ,
42 0 , 0.001) ;
43

44 if ( status == G S L _ S U C C E S S )
45 printf (" C o n v e r g e d :\ n ") ;
46

47 printf ("%5 d [%.7 f , %.7 f ] %.7 f %+.7 f %.7 f \ n " ,

48 iter , x_lo , x_hi ,
49 r , r - r_expected ,
50 x_hi - x_lo ) ;
51 }
52 while ( status == G S L _ C O N T I N U E && iter < m a x _ i t e r) ;
53
54 gsl_root_fsolver_free (s);
55
56 return status ;
57 }

El programa newton.c resuelve los ceros de resuelve x2 −5 = 0 con gsl_root_fsolver_newton

Listing 1.32: newton.c
1 # include < stdio .h >
2 # include < gsl / g s l _ e r r n o.h >
3 # include < gsl / g s l _ m a t h.h >
4 # include < gsl / g s l _ r o o t s.h >
5

6 # i n c l u d e " d e m o _ f n. h "
7 # i n c l u d e " d e m o _ f n. c "
8

9 int main ( void )

10 {
11 int status ;
12 int iter = 0 , m a x _ i t e r = 100;
13 const g s l _ r o o t _ f d f s o l v e r _ t y p e * T ;
14 gsl_root_fdfsolver *s;
15 double x0 , x = 5.0 , r _ e x p e c t e d = sqrt (5.0) ;
16 g s l _ f u n c t i o n _ f d f FDF ;
17 struct q u a d r a t i c _ p a r a m s params = {1.0 , 0.0 , -5.0};
18
19 FDF . f = & q u a d r a t i c;
20 FDF . df = & q u a d r a t i c _ d e r i v ;
21 FDF . fdf = & q u a d r a t i c _ f d f ;
22 FDF . params = & params ;
23
24 T = gsl_root_fdfsolver_newton ;
25 s = gsl_root_fdfsolver_alloc (T);
A-cursos/cursoMNu/raı́[Link] - August 22, 2023 71

26 g s l _ r o o t _ f d f s o l v e r _ s e t (s , & FDF , x ) ;
27
28 printf (" using % s method \ n " ,
29 gsl_root_fdfsolver_name (s));
30
31 printf ("% -5 s %10 s %10 s %10 s \ n " ,
32 " iter " , " root " , " err " , " err ( est ) ") ;
33 do
34 {
35 iter ++;
36 status = g s l _ r o o t _ f d f s o l v e r _ i t e r a t e ( s ) ;
37 x0 = x ;
38 x = gsl_root_fdfsolver_root (s);
39 status = g s l _ r o o t _ t e s t _ d e l t a (x , x0 , 0 , 1e -3) ;
40

41 if ( status == G S L _ S U C C E S S )
42 printf (" C o n v e r g e d :\ n ") ;
43
44 printf ("%5 d %10.7 f %+10.7 f %10.7 f \ n " ,
45 iter , x , x - r_expected , x - x0 ) ;
46 }
47 while ( status == G S L _ C O N T I N U E && iter < m a x _ i t e r) ;
48

49 gsl_root_fdfsolver_free (s);
50 return status ;
51 }

Compilación con GSL

La librerı́a GSL utiliza para cálculos algebraicos las librerı́a cblas. Esta es una versión en
C de la librerı́a BLAS (Basic Linear Algebra Subprograms). Esta librerı́a, originalmente
en fortran pero ahora en varios otros lenguajes, especifica o estandariza operaciones básicas
de álgebra lineal como ser la multiplicación de vectores, matrices y vectores etc. Es un
estándar pero tiene versiones de referencia desarrolladas mantenidas por [Link]. Por
supuesto existen versiones comerciales de estas librerı́as optimizadas para computadoras
de alta perfomance. Numerosos sistemas de computación numérica mantienen las mismas
interfaces que hacen compatible el uso de sus funciones: entre ellas Armadillo, LAPACK,
LINPACK, GNU Octave, Mathematica, MATLAB, NumPy, R, y Julia.
Si las librerı́as gsl y cblas están bien instaladas en un sistema Linux, se puede puede
compilar con el siguiente comando

gcc -o bb brent.c -L/usr/local/lib -lgsl -lgslcblas -lm

72

Bibliography
[1] IMSL Numerical Libraries. URL [Link]

[2] Wikipedia. List of numerical libraries. URL [Link]

of_numerical_libraries.

[3] The HSL Mathematical Software Library, 2021. URL [Link]

[4] F. S. Acton. Numerical Methods That Works. Harper and Row, 1970.

[5] K. E. Atkinson. Elementary Numerical Analysis. John Wiley and Sons, 1993.

[6] Kendall E. Atkinson. An Introduction to Numerical Analysis. John Wiley and Sons,
second edition, 1989.

[7] R. L. Burden and J. Douglas Faires. Análisis Numérico. Thomson Learning, 1989.

[8] Mc. Cracken and Dorn. Métodos Numéricos y Programación en Fortran. Limusa, 1980.

[9] M. Galassi et al. Gnu Scientific Library Reference Manual, 2018. URL [Link]
[Link]/software/gsl/.

[10] NAG. Nag Library, 2021. URL [Link]

[11] Clara O’Farrel. ¿hubo vida en marte? URL [Link]

universo/hubo-vida-en-marte/. Accesado: 8-12-2021.

[12] Ooi. Titulo. Titulo, 2004.

[13] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipies

in Fortran 77. The Art of Scientifc Computing. Cambridge University Press, 1992.

[14] G. G. Stokes. On the effect of internal friction of fluids on the motion of pendulums.
Transactions of the Cambridge Philosophical Society, 9, Part II:8–106, 1851.

[15] John R. Taylor. Classical Mechanics. University Science Books, 2005.

[16] Paul Whiters. Some cool motion sensor stuff. URL [Link]
v=KyktvC7w7Js.