Métodos_Numéricos_-_Ingeniería_aeronáutica_-_ETSEIAT_-_UPC

Vista previa del material en texto

Andrés Zarabozo 
Métodos Numéricos. Apuntes 
 
 
 
 
 
 
 
 
 
 
 
 
Ingeniería Aeronáutica 
ETSEIAT 
2012 
Andrés Zarabozo Martínez Métodos numéricos 
- 2 - 
 
 
Acerca de estos apuntes 
Estos apuntes se han realizado para cubrir el temario de la asignatura “métodos numéricos”, 
que se imparte en el quinto curso de Ingeniería Aeronáutica en la Escola Tècnica Superior 
d’Enginyeries Industrial i Aeronàutica de Terrassa, de la Universitat Politècnica de Catalunya 
(ETSEIAT – UPC). 
 
 
 
Licencia 
Esta obra está bajo una licencia Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0) de 
Creative Commons. Para ver una copia de esta licencia, visite: 
http://creativecommons.org/licenses/by-sa/3.0/deed.es_ES 
En líneas generales: 
Es libre de: 
 Compartir – Copiar, distribuir y comunicar públicamente la obra. 
 Transformar la obra y crear obras derivadas. 
 Hacer un uso comercial de esta obra. 
Bajo las condiciones siguientes: 
 Reconocimiento — Debe reconocer al autor de la obra original (pero no de una 
manera que sugiera que tiene su apoyo o apoya el uso que hace de su obra). 
 Compartir bajo la Misma Licencia — Si altera o transforma esta obra, o genera una 
obra derivada, sólo puede distribuir la obra generada bajo una licencia idéntica a ésta. 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 3 - 
Índice 
1. Programación en Fortran 
1.1. Tipos de variables .......................................................................................................... 7 
1.2. Expresiones aritméticas ................................................................................................ 7 
1.3. Bucles y condiciones ..................................................................................................... 8 
1.4. Vectores y matrices ....................................................................................................... 9 
1.5. Subrutinas ................................................................................................................... 11 
1.6. Lectura y escritura de datos ........................................................................................ 12 
1.7. Programas ................................................................................................................... 13 
 Mínimo error de un tipo de variable .............................................................. 13 Programa 1.
 Cálculo del número π ...................................................................................... 15 Programa 2.
2. Interpolación 
2.1. SPLINE .......................................................................................................................... 17 
2.2. Interpolación cúbica en un solo intervalo ................................................................... 18 
2.3. Programa ..................................................................................................................... 20 
 Entrada de datos e interpolación ................................................................... 20 Programa 3.
3. Integración 
3.1. Fórmula de Newton-Cotes .......................................................................................... 23 
3.2. Regla del punto medio ................................................................................................ 25 
3.3. Regla del trapecio ........................................................................................................ 25 
3.4. Regla de Simpson ........................................................................................................ 27 
3.5. Extrapolación de Richardson ....................................................................................... 28 
3.6. Programas ................................................................................................................... 30 
 Obtención de la serie de Fourier .................................................................... 30 Programa 4.
 Iteraciones y extrapolación ............................................................................ 33 Programa 5.
4. Resolución numérica de ecuaciones no lineales 
4.1. Método de la bisección ............................................................................................... 36 
4.2. Método de Newton - Raphson .................................................................................... 36 
4.3. Método de la secante .................................................................................................. 37 
4.4. Método de regula falsi ................................................................................................ 38 
4.5. Teorema del punto fijo ................................................................................................ 39 
4.6. Programas ................................................................................................................... 41 
Andrés Zarabozo Martínez Métodos numéricos 
- 4 - 
 Deflexión de un cable ..................................................................................... 41 Programa 6.
 Resolución mediante Newton-Raphson ......................................................... 44 Programa 7.
5. Resolución numérica de sistemas lineales 
5.1. Método de eliminación de Gauss ................................................................................ 45 
5.2. Matrices mal condicionadas ........................................................................................ 46 
5.3. Descomposición LU ..................................................................................................... 47 
5.4. Doolittle ....................................................................................................................... 47 
5.5. Crout ............................................................................................................................ 48 
5.6. Matrices dispersas ....................................................................................................... 49 
5.7. Método del gradiente conjugado................................................................................ 51 
5.8. Ejemplo de aplicación del método del gradiente conjugado ...................................... 56 
5.9. Programas ................................................................................................................... 58 
Programa 1. Resolución de un sistema mediante Gauss .................................................... 58 
Programa 2. Sistema con matriz mal condicionada ............................................................ 61 
Programa 3. Doolitle ........................................................................................................... 64 
Programa 4. Crout ............................................................................................................... 66 
Programa 5. Matriz banda ................................................................................................... 69 
6. Resolución numérica de sistemas no lineales 
6.1. Método de Newton para sistemas no lineales............................................................ 72 
7. Ecuaciones diferenciales, problemas de valor inicial 
7.1. Tipos de ecuaciones diferenciales ............................................................................... 74 
7.2. Método de Euler.......................................................................................................... 74 
7.3. Ejemplo de resolución mediante el método de Euler. ................................................ 75 
7.4. Método implícito de Euler ........................................................................................... 75 
7.5. Método de Heun ......................................................................................................... 76 
7.6. Método de Runge-Kutta de segundo orden ...............................................................77 
7.7. Método de Runge-Kutta de órdenes superiores ......................................................... 78 
7.8. Método adaptativo de Runge-Kutta ........................................................................... 80 
7.9. Esquema de Crank-Nicolson ........................................................................................ 81 
7.10. Método de Runge-Kutta para EDOs de grados superiores ......................................... 81 
8. Método de diferencias finitas 
8.1. Diferencias finitas en una dimensión. ......................................................................... 82 
8.2. Diferencias finitas en dos dimensiones ....................................................................... 83 
Andrés Zarabozo Martínez Métodos numéricos 
- 5 - 
8.3. Ejemplo de problema mecánico .................................................................................. 84 
8.4. Programas ................................................................................................................... 90 
Programa 6. Ejemplo de problema mecánico ..................................................................... 90 
9. Método de elementos finitos 
9.1. Forma fuerte del problema ......................................................................................... 96 
9.2. Método de Rayleigh-Ritz y método de Galerkin ......................................................... 96 
9.3. Ejemplo de formulación fuerte para problema estructural ........................................ 97 
9.4. Funciones de forma ................................................................................................... 100 
9.5. Cuadratura de Gauss ................................................................................................. 101 
9.6. Tipos de elementos en 2D ......................................................................................... 101 
9.7. Ejemplo de problema térmico................................................................................... 105 
9.8. Programa ................................................................................................................... 111 
Programa 7. Radiador espacial .......................................................................................... 111 
10. Apéndice 
10.1. Programación en Fortran .......................................................................................... 123 
 Mínimo error de un tipo de variable ............................................................ 123 Programa 1.
 Cálculo del número π .................................................................................... 124 Programa 2.
10.2. Interpolación ............................................................................................................. 125 
 Entrada de datos e interpolación ................................................................. 125 Programa 3.
10.3. Integración ................................................................................................................ 127 
 Obtención de la serie de Fourier .................................................................. 127 Programa 4.
 Iteraciones y extrapolación .......................................................................... 129 Programa 5.
10.4. Resolución numérica de ecuaciones no lineales ....................................................... 131 
 Deflexión de un cable ................................................................................... 131 Programa 6.
 Resolución mediante Newton-Raphson ....................................................... 132 Programa 7.
10.5. Resolución numérica de ecuaciones lineales ............................................................ 133 
 Resolución de un sistema mediante Gauss .................................................. 133 Programa 8.
 Sistema con matriz mal condicionada .......................................................... 134 Programa 9.
 Doolittle ...................................................................................................... 136 Programa 10.
 Crout ........................................................................................................... 138 Programa 11.
 Matriz Banda .............................................................................................. 140 Programa 12.
10.6. Método de diferencias finitas ................................................................................... 142 
 Ejemplo de problema mecánico ................................................................. 142 Programa 13.
Andrés Zarabozo Martínez Métodos numéricos 
- 6 - 
10.7. Método de elementos finitos .................................................................................... 145 
 Radiador especial ........................................................................................ 145 Programa 15.
10.8. Datos de entrada del Programa 13, radiador espacial .............................................. 154 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 7 - 
1. Programación en Fortran 
1.1. Tipos de variables 
En Fortran se suelen usar cinco tipos básicos de variables: 
 Boolean: 
Es una variable de un solo bit y por lo tanto solo puede tener dos posibles valores 
(.true. o .false.). 
 Integer: 
Solo permite tener un valor enteros, en general ocupa bits. El primer bit contiene el 
signo y por lo tanto el módulo del número ocupa los restantes bits. El valor máximo 
es igual a . 
 Real: 
Es similar al Double usado en otros lenguajes de programación como C. Permite 
números en decimales y en general ocupa bits. De forma similar al Integer el primer 
bit contiene el signo del número. 
 Complex: 
Son números complejos compuestos de dos números Real separados por una coma y 
entre paréntesis. El primer número representa el valor real y el segundo el valor 
complejo. 
 Character: 
Se utiliza para escribir letras. Se suele generar un vector donde cada posición es un 
caracter. 
Para definir una variable se debe declarar primero el tipo y luego el nombre de la variable. Las 
variables tienen que tener un nombre entre y letras o números siendo la primera siempre 
una letra. En Fortran no hay diferencia entre mayúsculas y minúsculas. 
 
real x 
integer i 
 
 
También se puede definir el tamaño en memoria de las distintas variables como por ejemplo: 
 
real*8 x 
integer*8 i 
 
 
1.2. Expresiones aritméticas 
Las cinco operaciones aritméticas básicas son: la suma (+), la resta (-), la multiplicación (*), la 
división (/) y el exponente (**). Fortran también reconoce otro tipo de funciones como por 
ejemplo funciones trigonométricas, logaritmos… Las operaciones aritméticas devuelven un 
resultado (que puede ser un Real o Integer) y debe ser almacenado en una variable. 
Andrés Zarabozo Martínez Métodos numéricos 
- 8 - 
 
x = 1. + 3. 25 * 5. / 3. 
x = cos(1.25) + y**2. 
 
 
El orden de las operaciones es: 
 Primero: 
Todas las operaciones exponenciales y se desarrollan de derecha a izquierda. 
 Segundo: 
Todas las multiplicaciones y divisiones, el orden es de izquierda a derecha. 
 Tercero: 
Todas las operaciones de suma y resta, el orden es de izquierda a derecha. 
En este libro se va a centrar el estudio del uso de Fortran para aplicación en métodos 
numéricos. Los programas se deben escribir de forma óptima para el cálculo. La latencia es el 
tiempo que se consume en una operación. Las operaciones básicas (menos los exponentes) 
tienen latencias muy bajas, en especial la suma, la resta y la multiplicación. La división tiene 
una latencia un poco más alta aunque en los nuevos procesadores como por ejemplo los Intel 
i7 la diferencia en latencia es despreciable, en ordenadores antiguos la latencia puede ser 
hasta cinco veces más lenta. Las operaciones exponenciales, trigonométricas, logarítmicas… 
son demasiado lentas y deben ser evitadas en lamedida de lo posible. Por ejemplo si se quiere 
calcular la siguiente ecuación: 
 ( ) ∑ 
 
 
 
 (1.1) 
 
Se debería evitar hacer la operación exponente y se podría escribir de la siguiente forma (la 
terminología de los bucles en Fortran se explica más adelante). 
 
p = a(n) 
do i .eq. n, 1, -1 
 p = p * x + a (n - 1) 
end do 
 
 
1.3. Bucles y condiciones 
El bucle estándar en Fortran es do. Se le da un valor inicial a una variable Integer y se le 
define el valor final. El código que se ejecuta en el bucle se debe finalizar con end do. 
do i .eq. 1, 10, 2 
 x = x + i 
end do 
Andrés Zarabozo Martínez Métodos numéricos 
- 9 - 
Por defecto cada vez que se ejecuta el bucle se le suma uno a la variable (en el ejemplo i) pero 
se puede cambiar el valor que se añade poniendo después del valor final una coma y el 
sumando. 
La función if ejecuta solo el código si se cumple la condición impuesta. Para poner condiciones 
en Fortran se utiliza la siguiente sintaxis (abreviaciones del inglés): 
- Mayor que: .gt. 
- Menor que: .lt. 
- Mayor o igual a: .ge. 
- Menor o igual a: .le. 
- Igual a: .eq. 
- Diferente de: .ne. 
Se puede complementar la función if utilizando las funciones else (si la condición del if no se 
cumple se ejecuta el código que sigue al else) o elseif (parecido al else pero se le añade una 
condición más para que se ejecute el código). Se pueden poner tantos elseif como se quieran 
en un mismo conjunto de if. Se debe finalizar el código de esta función con endif. 
 
if (a .gt. 2) then 
 b = 1 
else 
 b = 0 
end if 
 
 
La función do while ejecuta repite el código hasta que la condición deje de cumplirse. El 
código dentro del do while se ejecuta hasta que se finaliza (con end do) y vuelve a la línea del 
principio para comprobar si se cumple aún la condición. Si se desea salir del bucle en algún 
punto intermedio se puede utilizar la función exit. 
 
do while (a .lt. 5) 
 a = a + 1 
end do 
 
 
1.4. Vectores y matrices 
Se usan los vectores y matrices para utilizar un solo nombre que refiera a un conjunto de 
direcciones de memoria. Todos los datos almacenados en un vector o una matriz deben ser del 
mismo tipo. 
Para usar un vector o matriz (a partir de ahora solo se referirán a matrices ya que el 
funcionamiento es casi el mismo) se debe primero decláralo. Se puede declarar de forma 
estática o de forma dinámica. Si se declara de forma dinámica no se define la dimensión de la 
matriz, mientras que si se define de forma estática si que se define la dimensión. 
Andrés Zarabozo Martínez Métodos numéricos 
- 10 - 
Para declarar una matriz de forma estática se pone el tipo de las variables que almacena 
seguido del nombre de la matriz y en paréntesis sus dimensiones, se utiliza una coma para 
definir las distintas columnas de la matriz (el valor de la dimensión de cada columna debe ser 
un integro). 
 
real a(4,9) 
integer b(2,n) 
 
 
Para declarar una matriz de forma dinámica se debe utilizar el parámetro allocatable y utilizar 
dos puntos para las columnas que no se quiere definir la dimensión. Una vez se decide dar un 
valor a la dimensión se utiliza la función allocate. 
 
real, allocatable : : a(4, :) 
integer n 
 
print*,"Dimensions of the matrix?" 
read(*,*) n 
 
allocate(a(:,n)) 
 
 
Una matriz creada de forma dinámica puede ser eliminada (para liberar espacio en la 
memoria) con la función dellocate. Para saber si una matriz dinámica está creada o no en un 
momento dado se puede usar la función intrínseca allocatable (la función devuelve un valor 
booleano). 
Las tablas se puede utilizar globalmente, referenciando un conjunto de elementos seguidos o 
de forma individual. Si solo se indica el nombre de la matriz se referencia todos los elementos. 
Para referenciar un valor específico se debe introducir la posición del valor de la matriz. Si se 
quiere referenciar un rango se utiliza el doble punto para delimitar el rango. 
 
real M(2,4), V(3) 
 
V(2) = 1. 
M(2,2:4) = V(2) 
 
 
En el ejemplo anterior se introduce el valor de V(2) (es decir ) en los tres valores de la matriz 
M (M(2,2), M(2,3), M(2,4)). 
Las operaciones que se pueden hacer con matrices son: 
 x * M: Multiplicación de cada elemento de una matriz M y un escalar x. 
 x + M: Suma de cada elemento de una matriz M y un escalar x. 
 M * M: Multiplicación de cada elemento por su homologo elemento entre dos 
matrices de misma dimensión. 
Andrés Zarabozo Martínez Métodos numéricos 
- 11 - 
 traspose(M): Transpuesta de una matriz. 
 matmul(M,M): Multiplicación matricial entre dos matrices. 
 dot_product(V,V): Producto escalar entre dos vectores. 
Las matrices se almacenan en direcciones consecutivas de memoria. En el caso de un vector 
los valores se almacena secuencialmente, empezando por el primero del vector V(1) hasta el 
último valor V(n). 
V(1) V(2) V(3) … V(n-1) V(n) 
Figura 1.1. Orden de los valores de un vector en la memoria. 
En el caso de matrices, los valores se almacenan por columnas. Por ejemplo si se tiene una 
matriz M(2,2,2) el orden en la memoria es: 
M(1,1,1) M(2,1,1) M(1,2,1) M(2,2,1) M(1,1,2) M(2,1,2) M(1,2,2) M(2,2,2) 
Figura 1.2. Orden de los valores de una matriz en la memoria. 
1.5. Subrutinas 
Las subrutinas o subroutine son como funciones en otros lenguajes de programación. Se crea 
una subrutina de la misma forma que el programa principal, con la única diferencia que se 
deben poner las variables de referencia que se le quieren pasar a la subrutina. 
Es recomendable definir al principio de la rutina si la variable es de entrada o de salida. De esta 
forma el compilador puede informar de errores en el programa, si por ejemplo se le intenta 
dar un valor a una variable de entrada. 
En las subrutinas las variables se envían como referencia, es decir, se envía la dirección de 
memoria de la variable. Si la variable es una matriz se envía la dirección del primer valor de la 
matriz. 
 
subroutine areaCircle(A, r) 
 implicit none 
 real, intent (in) :: r 
 real, intent (out) :: A 
 
 real, parameter :: pi = acos(-1) 
 
 A = pi * r * r 
end subroutine 
 
 
La subrutina de ejemplo calcula el área de un círculo. El valor de entrada es el radio y devuelve 
el área. Para llamar a la subrutina se utiliza la función call. 
 
real A1, r1 
r1 = 3. 
call areaCircle(A1, r1) 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 12 - 
1.6. Lectura y escritura de datos 
Cuando se quiere trabajar con un archivo para leer o escribir datos se debe usar la función 
open. Se le asigna un número de unidad al archivo abierto que es el identificador. Dentro de la 
función se debe definir el archivo y el acceso. El acceso se puede definir por ejemplo como 
secuencial, es decir, que cada vez que se escriba o se lea pasa a la siguiente línea. Para el 
archivo se debe introducir la ruta completa o bien si el archivo está en la misma carpeta se 
puede simplemente poner el nombre. 
Para leer un archivo se debe utilizar la función read, seguido de una variable o matriz en 
donde se quiera almacenar los datos de lectura. Para escribir se usa write seguido de las 
variables o los datos que se quieren introducir. 
Una vez se termina de usar el archivo es conveniente cerrarlo usando close. 
 
open(unit = 1, file='fileData', access='SEQUENTIAL') 
 read(1,*) n 
 write(1,*) 2 * n 
close(1) 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 13 - 
1.7. Programas 
 Mínimo error de un tipo de variable Programa 1.
El primer programa sirve de ejemplo para poder ver un poco como funcionan los programas 
escritos en Fortran. 
El objetivo es calcular el mínimo valor para un tipo de variable. Se utilizan variables real ya que 
no tiene sentido hacerlo con integer (el error sería ). Se busca el menor número tal que 
 . En Fortran existe una función intrínseca que calcula lo mismo y que se llama 
epsilon(). 
Como todo programa se debe empezar declarando el programa y lasvariables que se van a 
utilizar. Se necesita al menos una variable donde se pueda almacenar el resultado deseado 
(dándole un valor inicial) y también se necesita una variable para el bucle que se va a realizar. 
program program1 
 implicit none 
 integer i 
 
 Small = 1. 
 
 
end program 
 
 
El bucle que se quiere realizar consiste en ir dividiendo por dos el valor de la variable small e ir 
comprobando que se cumpla la ecuación . Se elige en este ejemplo un bucle tipo do 
pese a que sería más lógico utilizar un bucle while, de esta forma se puede ver como funciona 
también la condición if. 
Dentro de cada iteración se debe comprobar que se cumpla la ecuación descrita, y si se cumple 
se debe dividir el valor del resultado por dos. En caso de que no se cumpla la condición hay 
que salir de bucle utilizando exit. 
 
do i = 1, 1000 
 if (1. + small .gt. 1.) then 
 small = 0.5 * small 
 else 
 exit 
 end if 
end do 
 
 
El resultado aún no es exactamente el valor almacenado en small, el valor bueno es el 
resultado multiplicado por dos. Para poder ver el valor en pantalla se usa la función print*,. Se 
puede también utilizar la función epsilon(small) para poder comparar el resultado. 
Andrés Zarabozo Martínez Métodos numéricos 
- 14 - 
 
small = 2 * small 
 
print*, small 
print*, epsilon(small) 
 
 
Se puede cambiar el tipo de variable cuando se define esa, pudiendo probar con variables con 
más precisión. 
 
real*8 small8 
real*16 small16 
 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 15 - 
 Cálculo del número π Programa 2.
Utilizando solo una cuerda de y un cuadrado de se puede calcular el número . Se 
dibuja un cuarto de círculo de radio un metro dentro del cuadrado y se van lanzando bolas 
sobre el cuadrado. Estadísticamente el número de bolas dentro de cada zona es proporcional 
al área. El área de un cuarto de círculo es: 
 
 
 
 
 
 
 (1.2) 
 
Siendo el número de bolas totales lanzadas y el número de bolas que caen en el cuarto de 
círculo, la relación entre y es: 
 
 
 
 
 
 
 
 ⁄
 
 
 
 
 (1.3) 
 
Por lo tanto el número es igual a: 
 
 
 
 (1.4) 
 
Para obtener un número aleatorio se utiliza la función random_number(). De esta forma se 
simboliza que se lanza una bola y aleatoriamente cae en una posición concreta del cuadrado. 
Como en todos los programas lo primero que se hace es declarar el programa y las variables. 
Para obtener mayor precisión se utilizan variables con double precision. Se necesitan cuatro 
variables. Dos variables e para definir la posición donde cae la bola, se podría también 
utilizar un vector. Y luego dos variables que definan el número de bolas totales y las que han 
caído dentro de cuarto de círculo. 
program program2 
 implicit none 
 real*8 x, y 
 integer*8 total, inside 
 
 total = 0 
 inside = 0 
 
 
Se debe empezar el proceso iterativo mediante un bucle. Se muestra a continuación otra 
forma de hacer bucles. Se utiliza la función continue y la función go to. La segunda debe 
tener un número integro para que el compilador sepa a donde tiene que ir. 
El primer paso del bucle es generar la posición aleatoria y sumarle un valor al número total de 
bolas lanzadas. 
Una vez se lanza la bola se debe comprobar si cae dentro del cuarto de círculo. Como el centro 
del eje de coordenadas se define en el centro del circulo se debe cumplir que . 
Andrés Zarabozo Martínez Métodos numéricos 
- 16 - 
Lo siguiente ya es calcular el valor de . Como las iteraciones se hacen muy rápido se puede 
marcar que muestre el valor de cada cierto número de iteraciones. La función mod(x,y) 
devuelve el resto de la división de íntegros. 
Hay que acordarse de que siempre que se hace un cálculo con variables de diferentes tipos se 
deben convertir al tipo de la variable que se está calculando. Esto puede evitar errores de 
cálculos que no son detectados por el compilador. En este problema, como se han calculado 
las variables con double precision se deben convertir con la función dble(). 
 
10 continue 
 
 !Ball thrown 
 call random_number(x) 
 call random_number(y) 
 total = total + 1 
 
 !Looks where the ball fell 
 if(x*x + y*y .lt. 1.) inside = inside + 1 
 
 !Shows the result every 1000 iterations 
 if (mod(total, 1000) .eq. 0) then 
 print*, total, 4. * dble(inside) / dble(total) 
 end if 
 
go to 10 
 
end program 
 
 
Este programa tiene un bucle que nunca finaliza, hay que tener cuidado con este tipo de 
bucles ya que, pese a que en este programa esto no afecta al resultado, en otro programa más 
complicado puede meterse en un bucle infinito y no calcular nada. Con este tipo de bucles se 
puede salir utilizando la función exit normalmente dentro de un if. 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 17 - 
2. Interpolación 
2.1. SPLINE 
En análisis numérico, un spline es una curva diferenciable definida en porciones mediante 
polinomios. Esto es útil para poder generar curvas en funciones donde se tienen una serie de 
puntos definidos o también para simplificar la función mediante tramos definidos como 
polinomios. 
El caso más sencillo de spline es el polinomio de primer orden, es decir, rectas que unan los 
puntos. Éste genera curvas con puntas donde la pendiente no está definida en los puntos. 
Cuanto mayor es el orden del polinomio de la aproximación, aparecen más oscilaciones, 
generando curvas de peor calidad. 
 
Figura 2.1. Aproximación de primer orden (recta) 
y aproximación de un orden muy superior. 
Se suelen usar polinomios de ordenes bajos, normalmente suelen ser de segundo o tercer 
orden. En general se busca una curva que pase por los puntos definidos y que lo haga con la 
misma pendiente que la función en los puntos. Esto se consigue con polinomios de tercer 
orden. 
Si lo único que se tiene es la distribución de puntos de la función se debe hacer una estimación 
de la pendiente utilizando puntos de alrededor. Si la separación entre los puntos es uniforme 
la pendiente de un punto se puede igualar a la pendiente de la recta entre los dos puntos 
adyacentes. 
 
Figura 2.2. Aproximación de la pendiente mediante dos puntos adyacentes. 
𝑥𝑖− 𝑥𝑖 𝑥𝑖+ 
Andrés Zarabozo Martínez Métodos numéricos 
- 18 - 
Siendo la pendiente: 
 
 
 
|
 
 
 + − 
 + − 
 ( ) (2.1) 
 
Donde el error es ( ) y decrece de forma cuadrática con la disminución del intervalo entre 
los puntos. 
En el caso de estar en los puntos extremos, al faltar uno de los puntos, se aproxima la 
pendiente con el propio punto. Por ejemplo si se está en el último punto no se tendría un 
punto a la derecha y por lo tanto la pendiente sería: 
 
 
 
|
 
 
 − 
 − 
 ( ) (2.2) 
 
Se puede ver que el error ( ) ha aumentado considerablemente con esta aproximación. 
2.2. Interpolación cúbica en un solo intervalo 
Para encontrar la función interpolada de tercer orden en un intervalo, se hace primero un 
cambio de variable. Teniendo una variable [ ] se transforma a una variable definida en 
un intervalo unitario [ ]. El cambio de variable es simplemente: 
 
 
 
 (2.3) 
 
Se usan polinomios de Hermite para encontrar la función interpolada. Se tienen cuatro 
funciones básicas de Hermite, siendo la función interpolada una combinación lineal de estas 
cuatro funciones. 
Las funciones solo pueden empezar en o y deben acabar en o . Las en y 
solo pueden ser o . 
 
Figura 2.3. Las cuatro funciones básicas de Hermite. 
Se busca un polinomio que pase por dos puntos conocidos ( y ) y que además la pendiente 
también coincida ( 
 y 
 ). 
𝐻 
𝐻 
𝐻 
𝐻4 
 
 𝜉 
Andrés Zarabozo Martínez Métodos numéricos 
- 19 - 
El polinomio es por lo tanto: 
 ( ) 
 ( ) ( ) 
 4( ) (2.4) 
 
El valor de las funciones básicas de Hermite se pueden calcularfácilmente ya que son 
polinomios de tercer orden. Se calcula como ejemplo la cuarta función 4. Esta función tiene 
doble raíz en (ya que la pendiente en ese punto es ) y la tercera raíz en . A parte 
se sabe que la pendiente en es . Por lo tanto la función tiene la siguiente forma, donde 
 es una constante que se debe calcular. 
 4 
 ( ) (2.5) 
 
Como la pendiente en es , se puede obtener el valor de la constante. 
 
 4
 
|
 
 ( )| (2.6) 
 
Por lo tanto el cuarto polinomio es: 
 4 
 ( ) (2.7) 
 
Las demás funciones se obtienen de forma similar y son: 
 ( )( )
 ( )
 
 ( ) (2.8) 
 
Para el valor de la pendiente también se tiene que hacer el cambio de variable. 
 
 
 
 
 
 
 
 
 
 (2.9) 
 
Sabiendo del cambio de variable que: 
 
 
 
 
 
 
 (2.10) 
 
Por lo tanto: 
 
 
 
 ( )
 
 
 (2.11) 
 
Finalmente la función interpolada es: 
 
 ( )( )
 
 ( ) ( ) 
 
 ( ) 
 ( ) ( )
 (2.12) 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 20 - 
2.3. Programa 
 Entrada de datos e interpolación Programa 3.
En este programa se ve un ejemplo de programa que puede leer y escribir en un archivo. 
Además se utiliza en el programa una subrutina para hacer la interpolación. 
El programa debe de ser capaz de leer un archivo con una serie de puntos de una función. 
Pregunta cuantos puntos se quieren tener de salida y genera unos puntos utilizando la 
interpolación cúbica. Los nuevos puntos calculados se escriben en un fichero de salida para 
poder comparar la función de entrada y la de salida. 
Se empieza escribiendo la subrutina. Se debe empezar definiendo las variables que se van a 
usar. Esta subrutina calcula un valor de ( ) (yx) para una posición de (x). Para calcular la 
función interpolada se necesitan los valores de la función en los extremos del tramo ( y ) y 
las pendientes en esos puntos ( 
 y 
 ), se utiliza un vector de cuatro posiciones para guardar 
estos valores (v1(4)). Además se necesita saber el valor del tramo y , se guardan también 
en un vector (v2(2)). 
 
subroutine hermite (v1, v2, x, yx) 
 implicit none 
 real, intent(in) :: v1(4), v2(2), x 
 real, intent(out) :: yx 
 
 
Se tiene que hacer el cambio de variable de a definido por la ecuación (2.3). 
 
 real :: chi, span 
 
 span =v2(2) - v2(1) 
 chi = (x - v2(1)) / span 
 
 
Se calcula ahora el valor ( ) con la formula de la interpolación descrita en la ecuación (2.12). 
Se puede utilizar el carácter & al final de una línea para indicar que el código sigue como si no 
hubiese un salto de línea, esto facilita la lectura del código y no afecta al programa. 
 
 yx = v1(1) * (1. + 2. * chi) * (1. - chi) * (1. - chi) + & 
 v1(2) * span * chi * (1. - chi) * (1. - chi) + & 
 v1(3) * chi * chi * (3. - 2. * chi) + & 
 v1(4) * span * chi * chi * (chi - 1.) 
 
end subroutine 
 
 
Se debe hacer ahora el programa principal que es el que lea los datos de entrada calcule los 
nuevos valores y los escriba en un fichero. 
Andrés Zarabozo Martínez Métodos numéricos 
- 21 - 
Lo primero que se debe hacer es definir las variables. Se usan vectores y matrices para guardar 
los datos de entrada. Se define una matriz para guardar el valor de la función y el valor de su 
derivada, y un vector para almacenar la posición en el eje . Como el programa no sabe aún 
cuantos valores de la función se tienen (los debe leer del archivo) se deben definir estas 
matrices como allocatable. 
Como variables adicionales se definen el número de valores de entrada (n) y el número de 
valores de salida (nExit). Además de estas variables se definen otras que ayudarán a la hora de 
programar. 
 
program program3 
 implicit none 
 real, allocatable :: ydy(:,:), x(:) 
 real xInt, yInt, step 
 integer n, i, nExit, last 
 
 
Se debe leer el archive con los datos para saber la dimensión de las matrices. Se pone como 
primer dato el número de valores de la función, y luego, en cada línea los valores de la función. 
Se utiliza la función allocate para dar las dimensiones a las matrices. 
 
open(unit = 1, file='initialData', access='SEQUENTIAL') 
 read(1,*) n 
 allocate(ydy(2,n), x(n)) 
 
 do i = 1, n 
 read(1,*) x(i), ydy(1,i) 
 end do 
 
close(1) 
 
 
Como solo se tiene la información de los puntos y de la función en los puntos se debe 
aproximar la derivada. Se utiliza la ecuación (2.2) para los puntos externos y la ecuación (3.1) 
para el resto. 
Una vez calculadas las derivadas se pide el número de puntos de salida. 
 
ydy(2,1) = (ydy(1,2) - ydy(1,1)) / (x(2) - x(1)) 
do i = 2, n - 1 
 ydy(2,i) = (ydy(1, i+1) - ydy (1, i-1)) / (x(i+1) - x(i-1)) 
end do 
ydy(2,n) = (ydy(1,n) - ydy(1,n-1)) / (x(n) - x(n-1)) 
 
 
print *, 'Introduce number of points:' 
read(*,*) nExit 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 22 - 
Se iguala la variable xInt a x(1). La variable xInt es la posición de la función interpolada. Se 
Calcula también el paso entre los puntos interpolados. 
 
 
 
 (2.13) 
 
Como se va a empezar a escribir resultados se abre el archivo de resultados. El parámetro que 
se le ha puesto en status genera un nuevo documento cada vez que se ejecuta el programa. 
En caso de que el documento ya exista, borra el anterior y crea uno nuevo. 
 
xInt = x(1) 
step = (x(n) - x(1)) / real(nExit - 1) 
open(unit=2, file='dataExit', access='SEQUENTIAL', status='REPLACE') 
 
 
Se empieza ahora el bucle de cálculo. Cada vez que se hace una iteración se calcular los puntos 
que forman el intervalo en el que xInt está. 
 
last = 1 
 
do i = 1, nExit 
 !Find interval 
 do while (xInt .gt. x(last+1)) 
 last = last + 1 
 end do 
 
 
Se llama a la subrutina para obtener el valor de yInt para el punto xInt. Una vez se obtiene el 
valor, se escribe el resultado en el archivo. Finalmente antes de terminar el bucle se toma un 
nuevo valor de xInt y se asegura que todos los valores iniciales de se han tomado. 
 
 call hermite(ydy(:,last:last+1), x(last:last+1), xInt, yInt) 
 
 write(2,*) xInt, yInt 
 xInt = xInt + step 
 xInt = min(xInt, x(n)) 
 end do 
 close(2) 
 
end program 
 
 
Debido a que Fortran únicamente envía del vector la dirección de memoria del primer valor 
asignado, se puede escribir la llamada de la subrutina de la siguiente manera. 
 
call hermite(ydy(1,last), x(last), xInt, yInt) 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 23 - 
3. Integración 
3.1. Fórmula de Newton-Cotes 
Si se tiene una función ( ) con [ ], la integral se define según: 
 Sean los elementos de área por debajo de la función donde 
 Si es suficientemente grande tal que | | y la partición de los elementos es 
regular, la integral entre y de ( ) es igual a un sumatorio de los elementos. 
 ∫ ( ) 
 
 
 
 
∑( − ) ( )
 
 
 [ − ] (3.1) 
 
Si la función está definida en un rango de valores de , podría haber problemas para obtener 
los valores de ( ). 
Una primera forma de calcular la integral es a través de la formula de Newton Cotes. Se puede 
muestrear la función en una serie de puntos y extrapolar los puntos a funciones que se puedan 
integrar como polinomios. Finalmente se hace la integral de los polinomios. 
Como convenio, cuando se escogen los puntos de los intervalos para integrar, se escogen los 
dos extremos y el resto se suelen coger de forma que el intervalo entre puntos es regular 
(siendo muy cómodo que el intervalo sea igual). 
Se considera una función como la que se ve en la Figura 3.1. 
 
Figura 3.1. Función aleatoria dividida en intervalos. 
Si los intervalos son iguales éste se puede calcular fácilmente con la posición de los extremos. 
 
 
 
 (3.2) 
 
Teniendo el valor de ( ) para todos los puntos ( ) se calculan las funcionespolinómicas que aproximan a ( ) para cada intervalo. Esta aproximación polinómica puede 
𝑦 
𝑥 
𝑓(𝑥) 
𝑎 𝑏 
𝑥 𝑥 𝑥 . 𝑥𝑛− 𝑥𝑛 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 24 - 
ser del orden que se desee, pudiendo aproximar con rectas, parábola u otros.La integral es 
entonces la integral de los polinomios (el conjunto se define como ). 
 ∫ ( ) 
 
 
 ∫ ( ) 
 
 
 (3.3) 
 
El polinomio se calcula mediante la interpolación polinómica de Lagrange. Para un intervalo 
que empieza en , el polinomio interpolado en la forma de Lagrange es la siguiente 
combinación lineal. 
 
 ( ) (3.4) 
 
Donde es la delta de Kronecker. 
 {
 
 
 (3.5) 
 
Por lo tanto el polinomio interpolado queda: 
 
 ( ) 
( )( ) ( − )( + ) ( )
( )( ) ( − )( + ) ( )
 (3.6) 
 
Si es muy grande, el polinomio tiene muchas oscilaciones y puede variar mucho respecto a la 
función inicial. Los valores típicos de son y . Al tener un pequeño, el polinomio de 
Lagrange queda muy simplificado. 
El error es la diferencia entre la integral que se quiere obtener y el valor obtenido. 
 ∫ ( ) 
 
 
 (3.7) 
 
En general el error al interpolar una función en un polinomio de orden es del orden de + . 
Por lo tanto el error al integral es del orden de: 
 ∫ + 
 
 
 + (3.8) 
 
Este error se denomina error simple ya que es el error de un solo intervalo. El error producido 
por todos los intervalos se denomina el error compuesto. 
 ∑∫ 
 ( ) 
 + 
 +( − ) 
 
 
 (3.9) 
 
 ∑ + 
 
 
 + (3.10) 
 
Como además se sabe que ( ) ⁄ , el error compuesto es por lo tanto del orden de: 
 ( ) + (3.11) 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 25 - 
3.2. Regla del punto medio 
La regla del punto medio establece que la función en los intervalos es constante, tomando solo 
un punto del intervalo. Los intervalos se normalizan para que quede una función comprendida 
entre y . 
 
Figura 3.2. Normalización de la función. 
Siendo: 
 ∫ ( ) 
 
 
 
 
 
∫ ( ) 
 
− 
 (3.12) 
 
Se aproxima ( ) a una recta cuyo valor es ( ). La integral de esta función es por lo tanto: 
 ∫ ( ) 
 
− 
 ∫ ( ) 
 
− 
 ( ) (3.13) 
 
Pese a que este caso sería el caso en que el error simple no es . Según el teorema de 
Taylor, una función se puede aproximar como un polinomio de la siguiente forma: 
 ( ) ( ) ( ) ( ) (3.14) 
 
La aproximación que se usa en este caso es de orden ( ( ) ( )). Pero el segundo 
término (que en teoría tendría que definir el error) es cero al integrarse, ya que es una función 
antisimétrica entre y . Por lo tanto el error viene definido por el tercer término y es: 
 
 
 (3.15) 
 
3.3. Regla del trapecio 
Se aumenta el orden y ahora se toman dos puntos del intervalo en vez de uno. El polinomio 
que aproxima el intervalo es por lo tanto una recta que pasa por los dos límites del intervalo. 
Se vuelve a normalizar la función para que esté contenida entre y . 
 
Figura 3.3. Polinomio normalizado produciendo un trapecio en el intervalo. 
 
𝑓(𝑥) 
𝑔(𝜉) 
𝜉 𝑥 
 𝑏 𝑎 
𝑔(𝜉) 
𝜉 
 
𝑃(𝜉) 
Andrés Zarabozo Martínez Métodos numéricos 
- 26 - 
La integral del polinomio es igual al área del trapecio. 
 
 ( ) ( )
 
 ( ) ( ) (3.16) 
 
La integral completa es la suma de las integrales de los distintos intervalos. 
 
 ( ) ( )
 
 
 ( ) ( )
 
 
 ( − ) ( )
 
 (3.17) 
 
Agrupando términos iguales se puede simplificar la expresión. 
 
 ( ) ( )
 
 ∑ ( )
 − 
 
 (3.18) 
 
Se define el primer término de la ecuación anterior como . 
 
 ( ) ( )
 
 
 ( ) ( )
 
 (3.19) 
 
Al estar aproximando de forma lineal el error simple es un término de segundo orden que al 
integrar se convierte en tercer orden. El error compuesto es de segundo orden. Los errores son 
del mismo orden de magnitud que los de la regla del punto medio, pero no son iguales. 
 
 
 (3.20) 
 
En general no basta con calcular una vez la integral para obtener un resultado válido. Se deben 
probar distintos valores de intervalos para comprobar si el resultado converge. Para no tener 
que calcular desde el principio cada vez. Es recomendable aprovechar los resultados obtenidos 
en las anteriores integraciones y una forma de hacerlo es dividir por dos el valor del 
intervalo , es decir poner un punto más en el centro de cada intervalo. De esta forma solo se 
tienen que calcular la mitad de puntos. 
Si se hace con este método se empezaría calculando la integral con un intervalo . El 
resultado saldría de hacer el siguiente cálculo: 
 ( ∑ ( )
 − 
 
) 
 
 
 (3.21) 
 
Nunca se debe dar por válido el resultado de la primera integración y por lo tanto hay que 
hacer una segunda. Se tiene ⁄ y por lo tanto , y además se define el 
resultado del sumatorio de la primera integración. 
 ( ∑ ( )
 − 
 
) ( ∑ ( )
 − 
 
) (3.22) 
 
Se ve fácilmente que solo es necesario hacer la mitad de los cálculos al solo tener que calcular 
los intervalos con subíndice impar. De forma similar se pueden ir haciendo iteración y el 
trabajo de cada iteración es el mismo. 
Andrés Zarabozo Martínez Métodos numéricos 
- 27 - 
3.4. Regla de Simpson 
En la regla de Simpson se extrapolan los intervalos en parábolas. Se toman tres puntos de cada 
intervalo ( ), los puntos son los extremos y el punto central del intervalo. De forma similar 
se centra a los estudios anteriores, se centra el intervalo en [ ]. 
 
Figura 3.4. Extrapolación parabólica utilizando tres puntos. 
Para obtener el valor de la integral por del intervalo se puede usar el polinomio interpolante 
de la ecuación (3.6), como se ha dicho antes el orden es . 
 
 ( )
( )( )
( )( )
 ( )
( )( )
( )( )
 ( )
( )( )
( )( )
 ( )
 
 
 ( )
 
 
 ( )
 
 
 
(3.23) 
 
Se integra el polinomio entre y . 
 ∫ 
 
− 
 ( ) 
 
 
 ( ) 
 
 
 ( ) 
 
 
 (3.24) 
 
Para calcular la integral completa se deben sumar las integrales de todos los intervalos. En este 
caso se toman grupos de tres puntos, y por el número de puntos (Figura 3.5) debe ser un 
número par. 
 
Figura 3.5. Límites de los intervalos de integración. 
𝑃(𝜉) 
𝜉 
 
𝑔(𝜉) 
𝑦 
𝑥 
𝑓(𝑥) 
𝑎 𝑏 
𝑥 𝑥 𝑥 . 𝑥𝑛− 𝑥𝑛 
Andrés Zarabozo Martínez Métodos numéricos 
- 28 - 
La integral completa es por lo tanto (siendo como antes ( ) ⁄ ): 
 [
 
 
 ( ) 
 
 
 ( ) 
 
 
 ( )] [
 
 
 ( − ) 
 
 
 ( − ) 
 
 
 ( )] (3.25) 
 
Si se agrupan los términos similares se obtiene: 
 
 [
 
 
 ( ) 
 
 
 ( ) 
 
 
 ( ) 
 
 
 ( − ) 
 
 
 ( − ) 
 
 
 ( )]
 (
 ( ) ( )
 
 
 
 
∑ ( ) 
 
 
∑ ( ))
 
(3.26) 
 
El cálculo del error es similar al de la regla del punto medio. Como es una regla cuadrática se 
podría pensar que el error seria 
4 que viene de integrar . Pero es una función 
antisimétrica y al integrarla entre y da cero. El error simple es por lo tanto 
 y 
el error compuesto es 
4. 
3.5. Extrapolación de Richardson 
Se puede acelerar el cálculo numérico (aumentando por lo tanto el orden de convergencia) 
con métodos como la extrapolación de Richardson. Permite a partir de un orden de 
convergencia , tener una aproximación . Es preferible que sea . 
Supongase que ( ) es una estimación de orden para: 
 
 
 ( ) (3.27) 
 
Al no poder hacer el límite infinito, se hace una aproximación de utilizando ( ). 
 ( )(3.28) 
 
Donde son constantes desconocidas reales y son constantes conocidas tal que 
 . 
Si solo se toma el primer término, la aproximación queda: 
 ( ) 
 ( + ) (3.29) 
 
Dividiendo el paso por dos se escribe como función de ⁄ . 
 (
 
 
) (
 
 
)
 
 ( + ) (3.30) 
 
Utilizando las ecuaciones (3.29) y (3.30) se elimina la constante desconocida . 
 ( ) ( ) (
 
 
) ( + ) (3.31) 
 
 
 
 ( ) (
 
 )
 
 ( + ) 
(3.32) 
Andrés Zarabozo Martínez Métodos numéricos 
- 29 - 
Con este proceso, se obtiene una mejor aproximación de eliminando el mayor término de 
error que es ( ). Este proceso se puede repetir para eliminar términos de error para 
obtener una mejor aproximación. 
El problema es que no se suele saber el valor de . Se suele usar igual al orden de cuadratura 
del método usado o hacer una aproximación del valor haciendo algunas iteraciones. 
Se calcula la aproximación de para el intervalo ⁄ . 
 (
 
 
) (
 
 
)
 
 ( + ) (3.33) 
 
Se elimina la constante desconocida utilizando la ecuación (3.30). 
 ( ) (
 
 
) (
 
 
) ( + ) (3.34) 
 
Si está a punto de converger se podría decir que el hecho de dividir el intervalo por la mitad no 
afecta al resultado. 
 ( ) (
 
 
) (
 
 
) (
 
 
) (3.35) 
 
Pudiendo aislar : 
 ( (
 
 
) (
 
 
)) (
 
 
) ( ) 
 
 (
 (
 
 ) 
( )
 (
 
 ) (
 
 )
) (3.36) 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 30 - 
3.6. Programas 
 Obtención de la serie de Fourier Programa 4.
En este programa se quiere calcular los coeficientes de una serie de Fourier. Los coeficientes 
de la serie de cosenos se obtienen utilizando la fórmula integral. 
 
 
 
∫ ( ) 
 
− 
 
 
 
∫ ( ) ( ) 
 
− 
 (3.37) 
 
De forma similar se pueden obtener los coeficientes de la serie de senos. 
 ∫ ( ) ( ) 
 
− 
 (3.38) 
 
Una vez se tienen los coeficientes de la serie de cosenos y los coeficientes de la serie de senos 
se pueden calcular la amplitud y el desfase. 
 √ 
 
 (
 
 
) (3.39) 
 
Las integrales se hacen utilizando el método del trapecio. Se crea una subrutina que calcule la 
integral de una función. A esta subrutina se le debe dar un vector y(n) donde se almacenen los 
valores de la función. 
 
subroutine trapice(y, b, z, n) 
 implicit none 
 integer, intent (in) :: n 
 real, intent (in) :: y(n), b 
 real, intent (out) :: z 
 
 z = 0.5 * (y(1) + y(n))+ sum(y(2:n-1)) 
 z = z * b / (real(n) - 1.) 
 
end subroutine 
 
 
La subrutina fourier debe calcular las fases y las amplitudes. Debe ir llamando a la subrutina 
trapice para calcular los coeficientes de la serie de cosenos y senos. 
Los valores de las amplitudes y de las fases se almacenan en una matriz llamada 
modes(2,0:m). Se le asigna un rango definido a la matriz (0:m) para mantener el convenio de 
subíndices. 
La variable n es el número de puntos seleccionados de la función para hacer las integrales. La 
variable m es el número de términos de la serie de Fourier que se quieren calcular. 
Se añaden además una serie de parámetros como pi o su inversa invPi. Como se sabe que se 
va a dividir por bastantes veces es recomendable tener una variable que guarde su inversa y 
así solo tener que hacer multiplicaciones. 
Andrés Zarabozo Martínez Métodos numéricos 
- 31 - 
 
subroutine fourier (y, modes, n, m) 
 implicit none 
 integer, intent (in) :: n, m 
 real, intent (in) :: y(n) 
 real, intent (out) :: modes(2, 0:m) 
 
 integer i 
 real z(n), x(n), dx, invPi, a, b 
 real, parameter :: pi = acos(-1.) 
 
 invPi = 1. / pi 
 
 
Se inicializa el vector de los puntos . 
 
 dx = 2. * pi / (real(n) - 1.) 
 x(1) = -pi 
 
 do i = 2, n 
 x(i) = x(i-1) + dx 
 end do 
 
 
Se calculan ahora los coeficientes de las series y luego finalmente las amplitudes y fases. El 
primer coeficiente se calcula fuera del bucle, este valor es ya la amplitud y no hay fase 
porque no multiplica una función trigonométrica. 
 
 !First coeficient 
 call trapice(y, 2.*pi, a, n) 
 modes(1,0) = a * 0.5 * invPi 
 
 do i = 1, m 
 !Cos coeficients 
 z = y * cos(real(i) * x) 
 call trapice(z, 2.*pi, a, n) 
 a = a * invPi 
 
 !Sin coeficients 
 z = y * sin(real(i) * x) 
 call trapice(z, 2.*pi, b, n) 
 b = b * invPi 
 
 !Amplitudes and phases 
 modes(1,i) = sqrt(a*a + b*b) 
 modes(2,i) = -atan(b / a) 
 end do 
 
end subroutine 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 32 - 
Finalmente se escribe la funcion principal del programa. Se utiliza una función de ejemplo que 
ya está escrita con sus amplitudes y fases para poder comprobar si el programa funciona bien. 
La función escogida es: 
 . ( . ) ( . ) (3.40) 
 
El programa debe primero seleccionar una serie de puntos y calcular el valor de la función 
almacenándolos en el vector y(n). El programa llama a la subrutina fourier que es la que se 
encarga de devolver las amplitudes y fases. 
 
program program4 
 implicit none 
 integer, parameter :: n = 300 
 integer i 
 real, parameter :: pi = acos(-1.) 
 
 real y(n), step, x, modes(2, 0:8) 
 
 x = -pi 
 step = 2. * pi / (real(n) - 1.) 
 
 !Discretization 
 do i = 1,n 
 y(i) = 0.5 + 3. * cos(2. * x + 0.75) + 1. * cos(6. * x + 0.25) 
 x = x + step 
 end do 
 
 !Amplitudes and phases 
 call fourier (y, modes, n, 8) 
 !Print results 
 do i = 0, 8 
 print* , modes(1,i), modes(2,i) 
 end do 
 
end program 
 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 33 - 
 Iteraciones y extrapolación Programa 5.
En este programa se calcula la integral de una función mediante el método del trapecio pero 
disminuyendo en cada iteración el tamaño del intervalo entre dos puntos. 
Como ya se ha explicado en el tema 3.3, si se divide el intervalo en dos, en los puntos donde ya 
se ha evaluado la función no hace falta volver a calcularlos y por lo tanto solo hay que hacer la 
mitad de cálculos para cada iteración. 
Para este programa se utilizan variables con double precision. Para ello las funciones 
trigonométricas que se deben usar son dsin y dcos. En ambas se debe introducir una variable 
de double precision. Para transformar el tipo de variable se usa la función dble. 
Se crea una subrutina que calcula la función. La función de ejemplo es: 
 ( ) (3.41) 
 
 
 
subroutine evaluateFunction(x, y) 
 implicit none 
 real* 8, intent(in) :: x 
 real* 8, intent(out) :: y 
 
 y = dcos(5. * x) + 3. * x * x 
end subroutine 
 
 
El programa principal se encarga de ir evaluando la función y sumando esos valores a la 
variable suma. Con ésta se puede calcular la integral en función del valor del intervalo de la 
iteración. 
Se inicializan los valores de las distintas variables y se evalúa la suma de los puntos extremos, 
que si además se divide por dos se obtiene el coeficiente de la ecuación (3.19). El número de 
intervalos viene definido por la variable n. 
 
program program5 
 implicit none 
 integer :: i, j, n 
 real*8, parameter :: a = 0., b = 1. 
 real*8 suma, prevInteg, integ, h, sumEnds, aux, x 
 
 n = 10 
 h = (b - a) / dble(n) 
 suma = 0. 
 
 !Sum of limits divided by two 
 call evaluateFunction(a, aux) 
 sumEnds = aux 
 call evaluateFunction(b, aux) 
 sumEnds = (sumEnds + aux) * 0.5 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 34 - 
Lo siguiente es hacer el bucle iterativo. En cada iteración simplemente hay que sumar los 
valores de la función a la variable suma. Se debe diferenciar la primera iteración de las demás 
ya que en la primera hay que evaluar la función en todos los puntos y en las demás 
simplemente hay que evaluarla cuando j (subíndice del punto ) sea impar. 
 
 do i = 1, 10 
 !Sum of the values of the function 
 if(i .eq. 1) then 
 x = a + h 
 do j = 1, (n - 1) 
 call evaluateFunction (x, aux) 
 suma= suma + aux 
 x = x + h 
 end do 
 
 else 
 x = a + h 
 do j = 1, (n - 1), 2 
 call evaluateFunction (x, aux) 
 suma = suma + aux 
 x = x + 2 * h 
 end do 
 end if 
 
 
La integral se obtiene multiplicando la suma de todos los puntos (dividiendo por dos los 
límites) y multiplicándolo por el intervalo, ecuación (3.18). Una vez obtenido el valor de la 
integral se calcula el tamaño del intervalo y el número de intervalos para la próxima iteración. 
 
 !Integral 
 integ = h * (sumEnds + suma) 
 
 !Number of intervals and size of interval 
 n = 2 * n 
 h = 0.5 * h 
 
 
Se sabe que el resultado analítico de la integral es: 
 ∫ ( ( ) ) 
 
 
 . (3.42) 
 
El resultado que se imprime es la comparación del resultado calculado y el resultado analítico. 
 
 !Results 
 print*, integ - (.2 * dsin(dble(5.)) + 1.) 
 
end program 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 35 - 
El programa ya está acabado. Si se ejecuta se puede ver como el error va disminuyendo a 
medida que se van haciendo iteraciones. Si se aumenta el número de iteraciones (en este 
programa se tiene que cambiar el último valor del bucle do) se puede ver que llega un punto 
en el que el error vuelve a aumentar. Esto es debido a los errores de redondeo. En este 
programa se necesita aumentar mucho el número de iteraciones porque se han utilizado 
variables de mucha precisión justamente para evitar este problema. 
Se puede extrapolar el resultado de la integral utilizando el método e extrapolación de 
Richardson. El valor de la integral extrapolada se calcula utilizando la ecuación (3.32). Como se 
usa el método del trapecio, el orden de convergencia es . 
 
 
 ( ) (
 
 
)
 
 
 ( ) (
 
 
)
 
 
(3.43) 
 
El resultado de cada iteración se puede mostrar junto al resultado utilizando la extrapolación 
de Richardson. 
 
 !Richardson extrapolation 
 if(i .eq. 1) then 
 print*,"No extrapolation", " With extrapolation" 
 analitical = .2 * dsin(dble(5.)) - 1. 
 else 
 aux = (4. * integ - prevInteg) / 3. 
 print*, integ - analitical, aux - analitical 
 end if 
 prevInteg = integ 
 end do 
 
end program 
 
 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 36 - 
4. Resolución numérica de ecuaciones no lineales 
4.1. Método de la bisección 
Los métodos de resolución numérica de ecuaciones no lineales permiten la obtención de 
resultados de para una función tal que ( ) . A la solución se le llama raíz o cero de la 
función. 
El primer método se basa en el teorema de Bolzano. Sea una función real continua en un 
intervalo cerrado [ ] con ( ) y ( ) de signos contrarios. Entonces existe al menos un 
punto del intervalo abierto ( ), con ( ) . 
El teorema como tal no especifica el número de puntos, pero afirma que al menos existe uno. 
 El método consiste en ir tomando el punto medio ( ( ) ⁄ ) de un intervalo [ ] tal 
que ( ) y ( ) de signos contrarios. Se forman dos tramos y se escoge el tramo tal que la 
función en los límites sea de signo contrario. Si por ejemplo ( ) y ( ) se tienen 
dos opciones: 
 
 
 
 {
 ( ) [ ]
 ( ) [ ]
 (4.1) 
 
El proceso se repite en cada iteración. La solución es: 
 
 
 
 (4.2) 
 
Si no hubiese errores de redondeo al hacer infinitas iteraciones se obtendría el resultado. 
 
 
 ̃ (4.3) 
 
La velocidad con la cual la sucesión de iteraciones converge se llama orden de convergencia. Si 
se tiene una secuencia ̃ que converge a un valor , se dice que la sucesión converge con 
orden a . 
 
 
| ̃ + |
| ̃ |
 
 (4.4) 
 
El número es llamado orden de convergencia, y en el caso del método de Bolzano . Es 
por lo tanto un método muy lento. 
En general no se puede calcular el orden de convergencia exacto, pero utilizando unas cuantas 
iteraciones se puede estimar, aunque siempre considerando que a función converge. 
4.2. Método de Newton - Raphson 
El método de Newton - Raphson es un método con un orden de convergencia mucho mayor, 
en general su orden de convergencia es . 
Andrés Zarabozo Martínez Métodos numéricos 
- 37 - 
Para éste método se debe escoger un punto que se aproxime a la solución. La cercanía del 
punto a la solución es importante para que el resultado converja. Se puede aproximar la 
función en una recta en ese punto utilizando la serie de Taylor. 
 ( ) ( ) 
 ( ) ( ) ( )
 (4.5) 
 
Se busca el punto en el cual ( ) . 
 ( ) 
 ( ) ( ) (4.6) 
 
 
 ( )
 ( )
 (4.7) 
 
 
Figura 4.1. Primera iteración del método de Newton-Raphson. 
Se obtiene una regla de recurrencia y por lo tanto se pueden ir haciendo iteraciones hasta 
obtener el resultado deseado. Además es una regla cuadrática ( ). Pero tiene algunos 
problemas: 
 Contiene la derivada de la función, y en ocasiones no es posible obtenerla si por 
ejemplo no se tiene una expresión analítica de la función. 
 Es un método abierto, la convergencia no está garantizada por un teorema de 
convergencia global. El resultado depende del primer valor que se estima. Se pueden 
por ejemplo tener problemas con máximos y mínimos. 
4.3. Método de la secante 
Con este método se evita el primer problema que tiene el método de Newton-Raphson. Se 
aproxima la tangente mediante la secante que pasa por dos puntos. En este método se 
necesitan dos puntos iniciales. La aproximación de la derivada es: 
 ( ) 
 ( ) ( )
 
 (4.8) 
 
Y por lo tanto: 
 ( ) 
 
 ( ) ( )
 (4.9) 
 
𝑥 𝑥 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
Andrés Zarabozo Martínez Métodos numéricos 
- 38 - 
 
 
Figura 4.2. Dos primeras iteraciones del método de la secante. 
En este caso se gana en simplicidad y flexibilidad de ejecución pero se pierde en orden de 
convergencia. Este método no es cuadrático pero su orden de convergencia suele ser entre . 
y . , y por lo tanto no es tan rápido como el método de Newton-Raphson. 
Este método tiene el mismo problema que el método de Newton-Raphson por lo que no se 
puede asegurar la convergencia y ésta depende de los valores iniciales. 
4.4. Método de regula falsi 
El método de regula falsi o falsa posición es una combinación del método de la secante y el 
método de la bisección. El orden de convergencia ronda la unidad (convergiendo más 
lentamente que el método de la secante) pero éste método asegura la convergencia a 
diferencia de por ejemplo el método de la secante. 
Se debe empezar con dos valores iniciales tal que los signos de la función es esos puntos son 
opuestos, garantizando que hay al menos una raíz en el interior del intervalo (teorema de 
Bolzano). Se traza una recta entre los dos puntos y se calcula el cero de esta recta. 
 
Figura 4.3. Dos primeras iteraciones del método regula falsi. 
𝑥 
𝑥 𝑥 𝑥 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
𝑥 
𝑥 𝑥 
𝑥 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
Andrés Zarabozo Martínez Métodos numéricos 
- 39 - 
Se calcula el punto a partir de los puntos y . 
 
 ( ) ( ) 
 ( ) ( )
 (4.10) 
 
Una vez se obtiene el nuevo punto se calcula ( ) y se elige un nuevo intervalo que cumpla 
que la raíz esté dentro de él. En el caso del ejemplo de la Figura 4.3 se toman los números y 
 , ya que ( ) ( ) . Se repite el proceso de forma iterativa. 
Éste método puede resultar lento si por ejemplo se tiene una función como el de la Figura 4.4 
y además no se puede asegurar que el resultado haya convergido correctamente (no se puede 
asegurar que el error sea mayor de lo estimado). Esto es debido a que en esta función las 
iteraciones avanzan solo en una dirección. 
 
Figura 4.4. Dos primeras iteraciones. 
4.5. Teorema del punto fijo 
El último método de estudio consiste en obtener una ecuación del tipo ( ). Se resuelve 
de forma iterativa utilizando un valor inicial + ( ). 
Gráficamenteeste método consiste en encontrar el punto donde se cruza una función 
y otra función ( ). Para ello se empieza con un valor , se calcula el valor de ( ) y se 
busca del valor de que sea igual. 
 
Figura 4.5. Dos primeras iteraciones, convergiendo. 
𝑥 
𝑥 
𝑥 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
𝑥 
𝑥 𝑥 𝑥 
𝑓(𝑥 ) 
𝑓(𝑥 ) 
Andrés Zarabozo Martínez Métodos numéricos 
- 40 - 
Este método puede no converger. Para que converja se debe cumplir que en la solución 
| ( )| . Si no se cumple la condición, la solución diverge alejándose en cada iteración de 
la solución, como se puede ver en la Figura 4.6. 
 
Figura 4.6. Tres primeras iteraciones, divergiendo. 
 
𝑥 𝑥 𝑥 
𝑥 
Andrés Zarabozo Martínez Métodos numéricos 
- 41 - 
4.6. Programas 
 Deflexión de un cable Programa 6.
En muchas películas de superagentes como por ejemplo James Bond se suelen ver como 
atraviesan una calle entre edificios solamente con un cable. Considerando que realmente 
pueden clavar el cable en el edificio de en frente y que se mantiene sujeto, se quiere estudiar 
si realmente el cable se mantiene rígido como muestran normalmente. 
Se tiene una cable sujeto en los extremos entre dos paredes y aguantando un peso en el 
centro, como muestra la Figura 4.7. El cable está caracterizado por sus propiedades y . En 
general las cuerdas que se usan para estos casos se llaman cuerdas de piano y son bastante 
rígidas. Está claro que en ningún caso estas cuerdas se pueden enrollar en un reloj o en el 
cinturón. 
 
Figura 4.7. Diagrama del cable con un peso en el medio. 
La tensión se relaciona con el peso del superagente utilizando el equilibrio de fuerzas 
verticales, siendo el ángulo que forma el cable deflectado con la horizontal. 
 (4.11) 
 
Suponiendo que la deformación es lineal la tensión del cable genera una pequeña 
deformación . 
 (4.12) 
 
La longitud del cable estirado se puede relacionar con la longitud inicial del cable . 
 
 
 
 (4.13) 
 
Por lo tanto la deformación es: 
 
 
 
 
 
 
 (4.14) 
 
Se juntan las ecuaciones (4.11), (4.12) y (4.14) se obtiene una relación entre el peso del 
superagente y el ángulo . 
 (
 
 
 ) (4.15) 
 
𝐴 𝐵 
 𝑙 
𝑃 
𝑇 𝑇 
Andrés Zarabozo Martínez Métodos numéricos 
- 42 - 
Simplificando la ecuación un poco se obtiene la ecuación que se quiere resolver 
numéricamente. En muchos problemas de física es bueno adimensionalizar las ecuaciones ya 
que estas muestran en un parámetro la relación entre propiedades importantes del problema. 
En este caso se relaciona el peso del superagente con la rigidez del cable. 
 
 
 
 (4.16) 
 
Se debe hacer una aproximación del resultado para tener un valor inicial. Se puede considerar 
que el ángulo es pequeño | | . Esto no implica que y ya que la 
ecuación ⁄ no tiene sentido. El problema está en el coseno de la ecuación (4.15). Se 
aproxima con un orden superior: 
 
 
 
 (4.17) 
 
Además se puede incluso hacer otra aproximación para facilitar la resolución de la ecuación. 
 
 
 
 
 
 
 
 
 
 
 
 
 (4.18) 
 
Por lo tanto la solución aproximada es: 
 
 
 
 
 
 
 √
 
 
 
 (4.19) 
 
Se utiliza el método de la secante para resolver la ecuación. Se puede aproximar el resultado 
utilizando un número cercano de la aproximación obtenida, por ejemplo . . 
Se toman los siguientes valores de las propiedades del cable y del superagente: 
 − (4.20) 
 
Lo primero que hay que hacer en este programa es calcular los valores iniciales. Además se 
crea una variable ⁄ . − . 
 
program Program6 
 implicit none 
 integer i, j, n 
 integer, parameter :: itMax = 7 
 real, parameter :: k = 0.25E-3 
 real x1, x2, x3, fx1, fx2, slope 
 
 x1 = (2. * k)**(1./3.) 
 x2 = 1.1 * x1 
 
 
Se empieza ahora el bucle iterativo. En este se debe evaluar la función en el último valor 
calculado de (para el código es x2). Como el primer valor (el primer x1) solo se tiene que 
calcular una vez éste se calcula justo antes de iniciar el bucle. 
Andrés Zarabozo Martínez Métodos numéricos 
- 43 - 
Como protección se pone una salida del bucle en caso de que x1 sea igual a x2. Este caso 
significaría que ya se tiene la solución y se evita además hacer cálculos con divisiones por cero. 
 
 fx1 = k + sin(x1) - tan(x1) 
 print*, 'x= ', x1, 'f(x) = ', fx1 
 
 do i = 1, itMax 
 if (x1.eq.x2) exit 
 
 fx2 = k + sin(x2) - tan (x2) 
 
 print*, 'x = ', x2, 'f(x) = ',fx2 
 
 
 
Para hacer la nueva iteración se calcula la pendiente utilizando los dos resultados anteriores y 
se usa la ecuación (4.9) para calcular el nuevo valor x3. Una vez calculado se vuelven a 
establecer las variables x1 y x2. 
 
 slope = (fx2 - fx1) / (x2 - x1) 
 x3 = x2 - fx2 / slope 
 
 x1 = x2 
 x2 = x3 
 
 fx1 = fx2 
 end do 
 
end program 
 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 44 - 
 Resolución mediante Newton-Raphson Programa 7.
Se resuelve la misma ecuación mediante el método de Newton-Raphson. En este caso no se 
necesitan tener dos valores iniciales pero sí que se necesita la derivada de la función. 
 ( ) 
 
 
 (4.21) 
 
Como se utilizan varias veces las funciones trigonométricas se puede optimizar el programa 
creando dos variables cosx y sinx que calculen solo una vez por iteración la función 
trigonométrica. 
 
program Program7 
 implicit none 
 integer i, j, n 
 integer, parameter :: itMax = 7 
 real, parameter :: k = 0.25E-3 
 real x1, x2, fx1, slope, cosx, sinx 
 
 x1 = (2. * k)**(1./3.) 
 
 do i = 1, itMax 
 
 cosx = cos(x1) 
 sinx = sin(x1) 
 
 fx1 = k + sinx - sinx / cosx 
 slope = cosx - 1. / (cosx * cosx) 
 
 print*, 'x = ', x1, 'f(x) = ', fx1 
 
 if (fx1 .eq. 0) exit 
 
 x2 = x1 - fx1 / slope 
 
 x1 = x2 
 end do 
 
end program 
 
 
 
 
Andrés Zarabozo Martínez Métodos numéricos 
- 45 - 
5. Resolución numérica de sistemas lineales 
5.1. Método de eliminación de Gauss 
El método de eliminación de Gauss es un algoritmo de álgebra lineal para determinar 
soluciones de un sistema de ecuaciones lineales. 
Se obtienen las soluciones del sistema lineal mediante la reducción del sistema dado a otro 
equivalente en el que cada ecuación tiene una incógnita menos que la anterior. 
Se transforma la matriz del sistema de ecuaciones en una matriz triangular superior. 
Obteniendo la solución de la última incógnita y pudiendo ir hacia atrás resolviendo las demás 
incógnitas del sistema. 
 
(
 
 
 
 
 
 
 
 
 )
 
 
 
 
{
 
 
 
 
 
 
 
 
 }
 
 
 
 
 
{
 
 
 
 
 
 
 
 
 }
 
 
 
 
 (5.1) 
 
A un sistema de ecuaciones lineal se le pueden hacer una serie de operaciones llamadas 
elementales. 
a. Multiplicar una ecuación por un escalar no nulo. 
b. Intercambiar de posición dos ecuaciones. 
c. Sumar a una ecuación un múltiplo de otra. 
En el caso del sistema matricial, estás se traducen en: 
a. Multiplicar una fila de la matriz y la fila del vector solución ( ) por un escalar. 
b. Intercambiar dos columnas de la matriz y dos filas del vector incógnitas. 
c. Sumar los términos de una fila a otra de la matriz y del vector solución. 
Para poder obtener la matriz triangular se debe, para cada fila: 
 Se debe ir a la columna que empiece a partir de la diagonal . 
 Si éste valor es nulo se debe cambiar la columna por otra que esté a la derecha y que 
tenga un valor no nulo. Si coincide que se está en la última fila con todos los valores 
nulos, el sistema es indeterminado. 
 Se deben obtener ceros por debajo de sumando múltiplos adecuados de la fila 
superior a las filas de debajo. 
Una vez obtenida la matriz triangular se pueden ir obteniendo los caloresde las resultados 
partiendo del último. 
El buen resultado de este método depende de la calidad de la matriz, que está relacionado con 
el radio espectral de la matriz, definido por el supremo de entre los valores absolutos de los 
elementos de su espectro. 
Andrés Zarabozo Martínez Métodos numéricos 
- 46 - 
5.2. Matrices mal condicionadas 
Existen varias formas de definir una matriz mal condicionada. Una matriz de un sistema de 
ecuaciones lineales ( ) está mal condicionada si, por ejemplo, pequeños errores o 
variaciones en los elementos de o tienen un gran efecto en la solución exacta de . 
Otra forma de definirlo de forma sencilla es que una matriz está mal condiciona si el 
determinante es anormalmente pequeño. Se puede comparar con la norma de la matriz ‖ ‖, 
estando mal condicionada si | | ‖ ‖. La norma de la matriz (RMS) se puede calcular como: 
 ‖ ‖ √∑∑ 
 
 
 (5.2) 
 
Como ejemplo se define el siguiente sistema de ecuaciones siguiente, donde : 
 {
 
 ( ) 
 (5.3) 
No se puede despreciar ya que si se hiciese el sistema no tendría solución (el determinante 
es nulo). Se utiliza el método de eliminación de Gauss. 
 {
 
 
 (5.4) 
 
La solución del sistema es por lo tanto: 
 
 
 
 
 
 
( 
 
 
) (5.5) 
 
Se puede ver que tanto la solución de como la de depende mucho del valor de . Un 
problema de redondeo del programa podría arruinar el resultado. Se puede comprobar que el 
determinante es mucho más pequeño que la norma. 
 | | ‖ ‖ √ ( ) √ (5.6) 
 
Cuando se tiene una matriz mal condicionada se debe evitar que los valores problemáticos 
(números muy pequeños) estén en la diagonal. Se debe por lo tanto reordenar la matriz en 
caso de tener ese problema. Si se quiere ser riguroso en la reordenación de la matriz se debe 
conseguir una matriz lo más parecido a una matriz diagonal dominante. 
Se debe calcular un parámetro de calidad en todas las filas (desde la fila del pívot hasta la 
fila n). Éste es el cociente entre el valor de la columna y los elementos de la fila. 
 |
 
 (| |)
| (5.7) 
 
Donde viene dado por el elemento de la columna ( [ ] donde es la posición del pívot) 
y viene dado por [ ]. 
Andrés Zarabozo Martínez Métodos numéricos 
- 47 - 
Para garantizar que la matriz esté bien condicionada se cambia la fila del pívot por aquella fila 
que tenga mayor. 
Otro método menos riguroso consiste en buscar el máximo valor de los elementos de la 
columna por debajo del pívot y cambiar la fila por la del pívot. Aunque de esta forma no se 
garantiza que no haya números muy grandes en el resto de la fila. 
5.3. Descomposición LU 
La descomposición LU (donde LU viene de las siglas inglesas Lower-Upper) es una forma de 
factorización de una matriz como el producto de una matriz triangular inferior y una superior. 
Esto puede ser útil cuando se debe resolver varios sistemas de ecuaciones donde la matriz 
no varía y lo que varía es el vector . Un ejemplo sería tener una estructura y se deben calcular 
varios estados de carga distintos. 
 
(
 
 
 
 
 
 
 )
 
 
 
(
 
 
 
 
 
 
 )
 
 
 (5.8) 
 
En este libro se estudian dos métodos de descomposición LU. Ambos dependen de qué matriz 
( o ) tenga unos en su diagonal. El método Doolittle introduce unos en la matriz ( ) 
y el método Crout tiene los unos en la matriz ( ). 
Una vez se tienen las matrices y se resuelve el sistema y luego . Como y 
 son diagonales, no hace falta hacer el método de Gauss para obtener la solución. 
 ( ) (5.9) 
 
5.4. Doolittle 
En el siguiente ejemplo se puede ver como haciendo el método Doolittle es análogo a hacer el 
método de Gauss. Se empieza con las dos matrices y de ejemplo. 
 
(
 
 
 
 
 
 4 4 4 )
 
 
 
(
 
 
 4
 4
 4
 44)
 
 
 (5.10) 
 
 
(
 
 
 4
 4 4
 4 4 4
 4 4 4 4 4 4 4 4 4 4 4 4 44)
 
 
 (5.11) 
 
En los elementos de la matriz aparecen las operaciones que se hacen al utilizar el método de 
eliminación de Gauss. 
Andrés Zarabozo Martínez Métodos numéricos 
- 48 - 
Para obtener las matrices con el método Doolittle se hacen las operaciones similares a la 
triangulación de Gauss. Se triangula la matriz dejando en los ceros el multiplicador. 
Como ejemplo se parte de la siguiente matriz : 
 (
 
 
 
) (5.12) 
 
Se empieza con el pívot de la primera columna. En los elementos y se deja el valor del 
multiplicador (valores encuadrados). 
 (
 
 
 
) (5.13) 
 
Se hace lo mismo para el siguiente pívot. 
 (
 
 
 
) (5.14) 
 
Por lo tanto las matrices y son: 
 (
 
 
 
) (
 
 
 
) (5.15) 
 
5.5. Crout 
De forma similar al método Doolittle se multiplican dos matrices y genéricas para ver la 
forma de la matriz . 
 (
 
 
 
 4 4 4 44
) (
 4
 4
 4
 
) (5.16) 
 
 
(
 
 
 4
 4 4
 4 4 4
 4 4 4 4 4 4 4 4 4 4 4 4 4 44)
 
 
 (5.17) 
 
De la primera columna se pueden sacar las variables y utilizando la primera fila se pueden 
obtener las variables . Se reduce el tamaño de la matriz para obtener los demás valores. 
 (
 4 4
 4 4 4
 4 4 4 4 4 4 4 4 4 4 4 4 44
) (5.18) 
Andrés Zarabozo Martínez Métodos numéricos 
- 49 - 
De forma similar en la primera columna se pueden obtener las variables , ya que los 
productos a los que suman estas variables contienen valores conocidos. De la primera fila se 
resuelven las variables utilizando los valores conocidos. Se vuelve a reducir la matriz. 
 (
 4 4 4
 4 4 4 4 4 4 4 4 4 4 44
) (5.19) 
 
Se pueden resolver y 4 y luego 4. Este proceso se va repitiendo con cada submatriz 
utilizando la primera columna para obtener las variables y la primera fila para obtener las 
variables . 
De forma general si se tiene la siguiente matriz modificada para ir obteniendo los valores de 
 y de (hasta el pívot ): 
 
(
 
 
 
 
 
 
 
 
 )
 
 
 
 (5.20) 
 
Los valores (donde son los elementos de la columna a partir de ) se obtienen de la 
siguiente forma: 
 ∑ 
 − 
 
 (5.21) 
 
Los valores de (donde son los elementos de la fila por detrás de ) se obtienen de la 
siguiente forma: 
 
 
 
 ( ∑ 
 − 
 
 ) (5.22) 
 
5.6. Matrices dispersas 
Una matriz dispersa es una matriz donde la mayor parte de los elementos son ceros. Existen 
métodos para agilizar la resolución de sistemas de ecuaciones disminuyendo la carga 
computacional y el tamaño consumido en memoria para almacenar todos los valores. 
 Almacenamiento Skyline 
Se generan vectores en cada fila de la matriz dispersa tomando desde el primero hasta el 
último elemento no nulo de cada fila. Es posible que hayan ceros entre los elementos 
seleccionados. Por ejemplo, se quiere almacenar la siguiente matriz: 
 
(
 
 
 
 
 
 
 )
 
 
 (5.23) 
Andrés Zarabozo Martínez Métodos numéricos 
- 50 - 
Los cuatro vectores que se forman son: 
 { } {
 
 
 
} {
 
 
 
} {
 
 
 
 
} { } (5.24) 
 
Se escriben todos estos valores en un vector contiguo. De esta forma se almacenan solo 
elementos