5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Capıtulo 1
Introduccion al Calculo Numerico
1. Introduccion al Caculo Numerico
a ) Introduccion al estudio del error
Error absoluto
Error relativo
Notacion decimal en coma flotante (numeros maquina)
Propagacion del error
Condicionamiento
Estabilidad
b) Condicionamiento de un problema
c) Estabilidad - algoritmo inestable.
2. Normas vectoriales y matriciales
1
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
2
1.1. Introduccion al Calculo Numerico
El Analisis Numerico es una herramienta fundamental en el campo de las Ciencias
Aplicadas.
Objetivo: Disenar metodos numericos de calculo que aproximen, de modo eficiente, la
solucion de problemas practicos previamente formulados matematicamente.
Algoritmo: Secuencia finita de operaciones algebraicas y logicas que producen una
solucion aproximada del problema matematico.
AN =⇒ diseno de algoritmos y estudio de su eficiencia
eficiencia - requerimientos de memoria
- tiempo de calculo (rapidez)- estimacion del error (precision)
1.1.1. Introduccion al estudio del error
El error nos proporciona la precision del metodo.
errores de
entrada
(en medidas)
+errores de
almacenamiento
+errores de
algoritmo a analizar
=
errores
de
salida
Sobre los errores de entrada nada podemos decir.
Antes de comenzar recordemos algunso conceptos que en la terminologıa estandar de
los errores se suelen usar.
Definicion 1 Error absoluto: Sea z la soluci´ on exacta del problema y
z la soluci´ on
aproximada, se define el error absoluto como
z − z
donde · es una norma (se definir´ a m´ as adelante).
(Pensar en · como el m´ odulo).
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
3
Desde el punto de vista de las aplicaciones, resulta mucho m as relevante el
error relativo z − zz
Nota 1 Si z = 0 s´ olo se trabaja con errores absolutos.
1.1.2. Notacion decimal en coma flotante
Es una forma de representacion de numero que contiene la informacion relevante com-
puesta por:
signo + fraccion + signo para el exponente + exponente.
Por ejemplo, -618.45 se puede expresar como:
−618,45 + 0 − 61,845 + 1 − 6184,5 − 1 − 6,1823 + 2
que significa:
−618,45 × 100 − 61,845 × 10 − 61845 × 10−1 − 6,1845 × 102
De todas ellas, se llama notacion decimal en coma flotante normalizada a aquella en
la que la fraccion esta comprendida entre 0 y 10.
±m
±E 1
≤m < 10, E
∈N
∪ {0
}Los numeros representables de forma exacta en el ordenador se llaman numeros maquina
Nota 2 Los ordenadores almacenan la informaci´ on en posiciones de memoria (o bit ≡Binary Digit), que s´ olo toman valores 0/1, encendido/apagado, positivo/negativo . . . . . .
luego utilizan la representaci´ on binaria.
Progagacion del error: Los errores anteriores se propagan a traves de los calculos, debido
a la estructura propia del algoritmo. Para estudiar esta propagaci on, y por tanto el errorfinal, atendemos dos conceptos:
- Condicionamiento
- Estabilidad
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
4
Condicionamiento: mide la influencia que tendrıan los errores en los datos en el caso en que
se puede trabajar con aritmetica exacta (⇒ no depende del algoritmo, sino del problema
en si).
Estabilidad: Esta relacionada con la influencia que tienen en los resultados finales la
acumulacion de errores que se producen al realizar las diferentes operaciones elementales
que constituyen el algoritmo.
1.1.3. Condicionamiento de un problema:
Diremos que un problema esta mal condicionado cuando pequenos cambios en los
datos dan lugar a grandes cambios en las respuestas.
Para estudiar el condicionamiento de un problema se introduce el llamado numero
de condicion de dicho problema, especıfico del problema, que es mejor cuanto mas cerca de
1 (el problema esta bien condicionado) y peor cuanto mas grande sea (peor condicionado).
objetivo: Definir el numero de condicion de un problema.
La gravedad de un problema mal condicionado reside en que su resolucion puede
producir soluciones muy dispares en cuanto los datos cambien un poco (algo muy frecuente
en las aplicaciones).
Ejemplo: Tenemos el siguiente sistema lineal:
Sistema 1:
10x1 + 7x2 + 8x3 + 7x4 = 32,
7x1 + 5x2 + 6x3 + 5x4 = 23,
8x1 + 6x2 + 10x3 + 9x4 = 33,
7x1 + 5x2 + 9x3 + 10x4 = 31
=⇒
x1 = 1
x2 = 1
x3 = 1
x4 = 1
.
Sistema 2:
10x1 + 7x2 + 8x3 + 7,2x4 = 32,
7,08x1 + 5,04x2 + 6x3 + 5x4 = 23,
8x1 + 5,98x2 + 9,89x3 + 9x4 = 33,
6,99x1 + 4,99x2 + 9x3 + 9,98x4 = 31
=⇒
x1 = −81
x2 = 137
x3 = −34
x4 = 22
.
Sistema 3:
10x1 + 7x2 + 8x3 + 7,2x4 = 32,
7x1 + 5x2 + 6x3 + 5x4 = 22,9,
8x1 + 6x2 + 10x3 + 9x4 = 32,98,
7x1 + 5x2 + 9x3 + 10x4 = 31,02
=⇒
x1 = 7,28
x2 = −9,36
x3 = 3,54
x4 = −0,5
.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5
Como se ve, pequenos cambios en los datos (del orden de 2 centesimas) en algunos
elementos, producen grandes cambios en las soluciones: 136 unidades del sistema 1 al
sistema 2.
Lo mismo ocurre al perturbar el segundo miembro del sistema: cambios de aproxima-
damente 1 decima producen cambios en la solucion de aproximadamente 13 unidades.
Lo anterior se debe a que el sistema esta mal condicionado.La gravedad de un problema mal condicionado reside en que su resolucion puede
producir soluciones muy dispares en cuanto los datos cambien un poco, cosa muy frecuente
en las aplicaciones.
1.1.4. Estabilidad
Todo algoritmo que resuelve un problema numericamente produce en cada paso un
error numerico.
Un algoritmo se dice inestable cuando los errores que se cometen en cada etapa del
mismo van aumentado de forma progresiva, de manera que el resultado final pierde gran
parte de su exactitud.
Un algoritmos es estable cuando no es inestable (controlado).
Condicionamiento y estabilidad =⇒Permiten estudiar la precision
de un algoritmo para
un problema concreto
1.2. Normas vectoriales y matriciales
Objetivo: Introducir herramientas de medicion de la variacion de los resultados:
1.2.1. Normas vectoriales
Definicion 2 Una aplicaci´ on · : Rn −→ R+ ∪ {0} es una norma si:
a) v = 0 ⇔ v = 0 ∈ Rn
b) λ v = |λ|v, ∀λ ∈ R, v ∈ Rn
c) u + v ≤ u + v ∀ u, v ∈ Rn (desigualdad triangular)
Ejemplo:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
6
Dado v ∈ R3,
v1 =n
i=1
|vi|,
v2 =
n
i=1|vi|2,
y
v∞ = max1≤i≤n
|vi|son normas vectoriales.
v = (1, −1, 3)
v1 = |1| + | − 1| + |3| = 5, v2 =√
1 + 1 + 9 =√
11, v∞ = max{1, 1, 3} = 3.
Ejercicio: Comprobar que · 1, · 2 y · ∞ son normas vectoriales.
Definicion 3 Dos normas vectoriales son equivalente · y · si existen constantesc1, c2 > 0 tales que:
c1v ≤ v ≤ c2v ∀v ∈ Rn.
En la practica esto significa que cuando · esta acotada, tambien · y viceversa.
1.2.2. Norma de una matriz
Definicion 4 Una norma matricial es una aplicaci´ on | · | : Mn −→ R+ ∪ {0} que
verifica las siguientes propiedades:
a) |A| = 0 ⇔ A ≡= 0
b) |λ A| = |λ||A|, λ ∈ R, A ∈ Mn
c) |A + B |≤ |A| + |B|, A, B ∈ Mn.
d) |A · B |≤ |A||B|.
Las tres primeras propiedades garantizanq ue es una normal vectorial y la ultima que
es compatible con el producto de matrices.
Un resultado que nos da una forma simple de construir una norma matricial es el
siguiente:
Definicion 5 Sea · una norma en Rn, se define la norma | · | : Mn → R+ ∪ {0}como
|A| = supv= 0
Avv = sup
v=1
Av
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
7
Cuando una norma matricial se define de la forma anterior (a traves de una norma
vectorial), se dice que es una norma matricial subordinada a la norma vectorial.
Tenemos las siguiente normas matriciales subordinadas a las vectoriales:
a)
|A
|1 = sup
v=0
Av1v1
b)
|A|2 = supv=0
Av2v2
c)
|A|∞ = supv=0
Av∞v∞
Nota 3 No todas las normas matriciales son normas matriciales subordinadas a normas
vectoriales.
Algunas propiedades de las normas matriciales subordinadas
1) Av ≤ |A|v, A ∈ Mn, v ∈ Rn
2) Existe un vector v ∈ Rn para el que se da la igualdad, es decir,
Av = |A| v
3) |I | = 1 (Por I denotamos la matriz identidad).
Veamos como calcular las normas matriciales subordinadas a las normas vectoriales
anteriores:
Teorema 1 Sea A = (ai,j)ni,j=1 ∈ Mn, se tiene:
a)
|A|1 = max1≤ j≤n
ni=1
|ai,j|
b) |A|2 =
ρ(A∗ A) =
ρ(A A∗) = |A∗|2donde ρ(A∗ A) es el radio espectral de A∗ A, que es el m´ aximo de los valores absolutos
de A∗ A. (Por A∗ denotamos a la matriz adjunta de A, la cual coincide con la matriz
transpuesta cuando todos sus elementos son reales).
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
8
c)
|A|∞ = max1≤i≤n
n j=1
|ai,j|
Nota 4 | · |1 y | · |∞ se calculan a partir de los elementos de la matriz, | · |2 no.
Es inmediato observar que
|AT
|∞ =
|A
|1.
Ejemplo:
Sea
A =
1 0 −7
0 2 2
−1 −1 0
.
Calculamos la norma uno, dos e infinito:
|A|1 = max{2, 3, 9} = 9
|A|∞ = max{8, 4, 2} = 8
AT A =
2 0 −7
1 5 4
−7 4 53
|A|2 = 7,3648
Norma Frobenius:
Es una norma matricial no subordinada a ninguna norma vectorial. Viene dada por:
|A|F =
ni,j=1
|ai,j|2
Como vemos, se calcula a traves de los elementos de la matriz.
Ejemplo:
|A|F =√
1 + 49 + 4 + 4 + 1 + 1 =√
60
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Capıtulo 2
Resolucion numerica de sistemas de
ecuaciones lineales
1. Resolucion numerica de sistemas de ecuaciones lineales
a ) Condicionamiento de sistemas. Precondicionadores.
b) Metodos directos: Gauss (con pivote parcial)
c) Metodos iterativos
Jacobi
Gauss-Seidel
9
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
10
2.1. Condicionamiento de sistemas
Objetivo: Resolver sistemas lineales de ecuaciones mediante metodos numericos de calculo.
La existencia de un sistema mal condicionado es una fuente de posibles errores y
dificultades a la hora de resolver un sistema lineal mediante metodos numericos.
El primer problema que se plantea es como definir y cuantificar el condicionamientode un sistema.
Supongamos que tenemos que resolver el sistema lineal Ax = b, donde A es una matriz
de coeficientes, b es el termino independiente y x es la solucion exacta del sistema,
que llamaremos u (el vector x recibira diferentes nombres a lo largo del tema):
A x = b ⇒ x = u
Si modificamos el termino independiente mediante una perturbacion δb, entonces te-
nemos que resolver el sistema Ax = b + δb, que tendra una nueva solucion (distinta de la
solucion exacta) que llamaremos u + δu:
A x = b + δb ⇒ x = u + δu
El sistema esta bien condicionado si cuando δb es pequena, δu tambien lo es. Ob-
servemos que:
Au + A(δu) = b + δb
Au = b
⇒
A(δu) = δb
δu = A−1(δu)
Usando ahora la propiedad para normas matriciales, obtenemos que:
δu ≤ |A−1|δb (2.1)
De la solucion exacta, b ≤ |A|u, lo que implica que:
1
u ≤ |A|b (2.2)
De (2.1) y (2.2), obtenemos:
δuu ≤ |A||A−1|δb
b , (2.3)
donde δu
u representa el error relativo en los resultados, y δb
b el error relativo enlos datos. De la relacion (2.3), parece deducirse que el numero |A||A−1| es el factor
determinante de la relacion, ya que si es pequeno tenemos el efecto deseado, y si no, ocurre
lo contrario. Si |A||A−1| = 1 yδbb = 0,01, entonces
δuu ≤ 0,01.
Parece entonces natural la siguiente definicion:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
11
Definicion 6 Sea | | · | una norma matricial subordinada y A una matriz invertible.
Llamamos numero de condicion de la matriz A respecto de la norma | · | a la
expresi´ on:
cond(A) = |A||A−1|.
De (2.3) deducimos entonces que:
δuu ≤ cond(A)
δbb .
En el caso en el que las perturbaciones se produzcan en la matriz del sistema A, la
matriz se transforma en A+∆A, llamaremos u+∆u a la solucion aproximada del sistema:
(A + ∆A)(u + ∆u) = b.
Usando que Au = b, obtenemos:
(∆A)(u + ∆u) + A∆u = 0 ⇒ ∆A(u + ∆u) = −A∆u,
luego
∆u = −A−1(∆A)(u + ∆u) ⇒ ∆u ≤ |A−1||∆A|u + ∆u,
y∆u
u + ∆u ≤ |A−1||A||∆A||A| = cond(A)
|∆A||A| (2.4)
2.1.1. Propiedades de cond(A)
Proposicion 1 Para cualquier norma subordinada | · |, se verifica:
1. cond(I ) = 1.
2. cond(A) ≥ 1.
3. cond(A) = cond(A−1),
4. cond(k
·A) = cond(A),
∀k
∈R
\{0
}.
Demostracion:
1. cond(I ) = I I −1| = I I = 1 · 1 = 1.
2. 1 = |I | = |AA−1| ≤ |A||A−1| = cond(A),
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
12
3. cond(A−1) = |A−1||(A−1)−1| = |A−1||A| = cond(A),
4. cond(kA) = |kA||(kA)−1| = |k||A||k|−1|A−1| = cond(A).
2
Ejemplo: Estudiar el condicionamiento del sistema Ax = b con A = 1 1 + ε
1 − ε 1 siendo ε > 0 en la norma | · |∞.
Solucion: A−1 =1
ε2
1 −ε − 1
ε − 1 1
, |A|∞ = 2 + ε y |A−1|∞ =
2 + ε
ε2, luego
cond(A, ∞) =(2 + ε)2
ε2>
4
ε2. Si ε ≤ 0,01, entonces cond(A, ∞) > 40000. Esto indica que
una perturbacion de los datos de 0,01 puede originar una perturbacion de la solucion del
sistema de 40000.
Ejercicio: Repetir el ejemplo anterior con la norma | · |1.
Nota 5 Un sistema lineal est´ a bien condicionado si la matriz A est´ a bien condicionada.
Nota 6
Como cond(A) ≥ 1 cuanto m´ as cerca este su valor de 1, mejor condicionado est´ a el
sistema.
2.1.2. Precondicionadores
En los casos en los que tenemos que resolver el sistema Au = b y cond(A) sea muy
grande, podemos intentar alterar el sistema mediante un precondicionador que rebaje
cond(A). Como cond(kA) = cond(A), no vale multiplicar el sistema por un escalar k.
Normalmente, lo que se hace es multiplicar a izquierda por una matriz P invertible, de
modo que:
A = P A,
tenga cond( A) pequeno, para despues resolver el sistema:Au = b, donde b = P b.
Observamos que al ser P una matriz no singular, la solucion de Au = b y P Au = P b
es la misma.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
13
La idea que se utiliza para definir un precondicionador es muy sencilla. Ya que la
matriz mejor condicionada es la identidad (cond(I ) = 1), si elegimos P = A−1 obtenemos
cond(P A) = cond(A−1A) = cond(I ) = 1.
Pero, por otro lado, no queremos calcular la matriz inversa de A, pues como veremos, es
la forma mas cara computacionalmente de resolver un sistema lineal. Es decir, la solucion
del sistema lineal
A x = b
siendo A una matriz no singular es
x = A−1 b.
Entre los muchos problemas de dar la solucion del sistema de esta forma se encuentra en
que, calcular A−1
es equivalente a resolver tantos sistemas lineales como su dimension.Ademas de que se pueden introducir muchos errores de redondeo (recuerda que para
calcular A−1 es necesario dividir por su determinante, por lo que si es pequeno, aunque
no sea cero, introduce errores de calculo el ordenador).
En resumen, tenemos que para mejorar el condicionamiento de P A, lo ideal es tomar
P = A−1, pero no es posible en la practica. Por lo tanto, para definir P lo que hare-
mos es simplemente dar una aproximacion de A−1. Esta claro que cuanto mejor sea la
aproximacion, mejor condicionada estara la matriz P A.
Precondicionadores: Son muchos los posibles precondicionadores que se pueden definir,y aun hoy en dia son objeto de estudio, debido a su importancia para resolver de una
forma satisfactoria sistemas lineales, los cuales vienen asociados a problemas reales: c alculo
de estructuras, diseno de alcantarillados, predicciones meteorologicas, diseno de aviones,
movimientos de aguas en rios y oceanos, etc.
Por ahora nos vamos a conformar con buscar precondicionadores que sean matrices
diagonales. Consideramos el siguiente precondicionador:
P 1 = diag1
et
1A2,
1
et
iA2, . . . ,
1
etnA2.
Donde por ei denotamos el vector definido por ceros y un 1 en la componente i. En
definitiva, etiA es la fila i-esima de la matriz A. Por lo que et
iA2 es la norma dos de la
fila i de la matriz A. Por lo tanto, el precondicionador anterior se define como una matriz
diagonal cuyos elementos son la inversa de la norma de la correspondiente fila.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
14
Observemos que si A es una matriz diagonal (con todos sus elementos positivos),
entonces P = A−1.
Ejemplo:
A =
1 2 0 0
2 2 3 0
0 3 4 10
0 0 10 5
=⇒ P =
1/√
5 0 0 0
0 1/√
17 0 0
0 0 1/√125 0
0 0 0 1/√
125
Se tiene que cond(A) = 67, mientras que cond(P A) = 31, 36. Observamos que se ha
reducido pr´ acticamente a la mitad el n´ umero de condici´ on para ambas definiciones del
precondicionador.
Existen muchas tecnicas para buscar precondicionadores, por ahora nos vamos a con-
formar con esta. Por supuesto, cuanto mas sofisticada sea la definicion del precondiciona-
dor, mejores resultados obtendremos.
2.2. Metodos directos e iterativos
Hay dos tipos de metodos para la resolucion de sistemas lineales:
Directos: proporcionan la solucion exacta (salvo errores de redondeo) en un numero
finito de pasos: Gauss.
Iterativos: proporcionan una sucesion {xk} que aproxima, o converge, a la solucionexacta:
{xk} −→ x.
El calculo se detiene cuando se alcanza un cierto nivel de precision.
2.3. Metodos directos
Consideramos la resolucion numerica de un sistema lineal Au = b, donde A es una
matriz invertible.
El principio de los metodos directos que vamos a estudiar reside en determinar una
matriz M invertible, tal que la matriz MA sea triangular superior. Tenemos que resolver
entonces el sistema lineal:
MAu = Mb,
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
15
por lo que llamaremos el metodo de remontada.
Este principio es la base del metodo de Gauss para la resolucion de sistemas lineales
con matrices cualesquiera.
2.3.1. Observaciones concernientes a la resolucion de sistemas
lineales
Contrariamente a lo que se piensa, la resolucion de un sistema lineal no es equivalente
al calculo de la matriz inversa del sistema A−1 y despues calcular A−1b. El calculo de la
matriz inversa es equivalente a la resolucion de n sistemas lineales (donde n es el orden
de la matriz A):
Au j = e j 1 ≤ j ≤ n,
donde e j es el n-esimo vector de la base Rn. Pasamos ası de resolver un sistema lineal a
resolver n sistemas lineales y multiplicar A−1 por el termino independiente b.
Los metodos que vamos a estudiar estan basados en el siguiente hecho: si tuviesemos
una matriz triangular superior, la resolucion numerica de un sistema lineal Au = b es
inmediata. Tendrıamos, en forma matricial:
a1,1 . a1,n−1 an
0 . . .0 . . .
0 . . .
0 . an−1,n−1 an−1,n
0 . . an,n
u1
.
.
.
un−1
un
=
b1
.
.
.
bn−1
bn
y en forma de ecuaciones:
a1,1u1 + ... + a1,n−1un−1 + a1,nun = b1
... .
... .
... .
an−1,n−1un−1 + an−1,nun = bn−1
an,nun = bn
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
16
El determinante es det(A) = a11a22...ann = 0, luego el sistema se resuelve:
un = a−1nnbn
un−1 = a−1n−1,n−1(bn−1 − an−1,nun)
...
...
...
u1 = a−111 (b1 − a12u2 − ...a1,n−1un−1 − a1nun)
De ese modo, cada componente ui se escribe como combinacion lineal de las bi, bi+1,...,bn,
luego estamos resolviendo un sitema lineal u = Cb donde C es una matriz triangular
superior. Este metodo se conoce como metodo de remontada. Dicho metodo necesita
un total de
1 + 2 + ... + n − 1 =
n(n
−1)
2 sumas
1 + 2 + ... + n − 1 =n(n − 1)
2multiplicaciones
n divisiones
2.3.2. El metodo de Gauss
El metodo de Gauss es un metodo general de resolucion de un sistema lineal de la
forma Au = b donde A es una matriz invertible. Se compone de tres etapas:
1. procedimiento de eliminacion, que equivale a determinar una matriz invertible M
tal que la matriz MA sea una matriz triangular superior,
2. calculo del vector Mb,
3. resolucion del sistema lineal MAu = Mb, por el metodo de remontada.
a) Etapa de eliminacion:
1) Al menos uno de los elementos de la primera columna de A, ai,j , 1 ≤ i ≤ n es diferente
de cero (o det(A) =0), ai,j, que llamaremos el primer pivote de la eliminacion.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
17
2) Intercambiamos la lınea donde esta el pivote con la primera lınea, lo que equivale a
multiplicar a izquierda por la matriz de trasposicion
T (i0, i1) =
1 | |. . . | |
−− −−0
−−1
−− −−}i0
| |−− −− 1 −− 0 −− −−} i1
| | . . .
| i0
| i1
1
con det(T (i0, i1)) = −1.
P = I si a1 1 = 0
T (1, i) si ai 1, i = 1 es el pivote, det(P ) = −1
y P A = (αi j ) tal que α1 1 = 0.
Multiplicamos por combinaciones lineales adecuadas de la primera lınea de P A con
las otras lıneas de P A, de forma que se anulan todos los elementos de la primera columna
situados debajo de la diagonal. Es decir, multiplicamos E P A, donde
E =
1 0 . . . 0
−
α2 1
α1 1
1
... . . .
−αn 1
α1 10 . . . 1
det(E ) = 1
Nota 7 Cuando se aplican tecnicas de eliminaci´ on, el coste (n´ umero de operaciones ar-
titmeticas requeridas) suele ser proporcional al n´ umero de coeficientes no cero de la matriz,
puesto que se usan tecnicas especiales para evitar los c´ alculos asociados a los coeficientes
nulos de la matriz. Adem´ as, el coste de calcular A−1 suele ser del orden de n3.
Ejemplo del metodo de Gauss
A =
0 1 2 1
1 2 1 3
1 1 −1 1
0 1 8 12
b =
1
0
5
2
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
18
→
1 2 1 3 | 0
0 1 2 1 | 1
1 1 −1 1 | 5
0 1 8 12 | 2
→
1 2 1 3 | 0
0 1 2 1 | 1
0 −1 −2 −2 | 5
0 1 8 12 | 2
P 1 =
0 1 0 01 0 0 0
0 0 1 0
0 0 0 1
, E 1 =
1 0 0 00 1 0 0
−1 0 1 0
0 0 0 1
→
1 2 1 3 | 0
0 1 2 1 | 1
0 0 0 −1 | 6
0 0 6 11|
1
→
1 2 1 3 | 0
0 1 2 1 | 1
0 0 6 11 | 1
0 0 0−
1|
6
E 2 =
1 0 0 0
0 1 0 0
0 1 1 0
0 −1 0 1
, P 3 =
1 0 0 0
0 1 0 0
0 0 0 1
0 0 1 0
luego M = P 3 E 2 E 1 P 1
En la practica, la matriz “de paso” M no se calcula explıcitamente, sino que se obtienen
directamente M A y M b.
Nota 8 La matriz M verifica
det (M ) =
1 si Λ es par
−1 si Λ es impar
Λ es el n´ umero de matrices de permutaci´ on P distintas de la identidad.
A u = b ⇔ M A u = M b.
Nota 9 En el caso de que la matriz A no es invertible, el metodo sigue siendo v´ alido, ya
que en ese caso se considerarıa P = E = I .
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
19
Comparacion de metodos:
Cramer
(n + 1) (n! − 1) sumas
(n + 1) (n − 1) n! productos
n divisiones
Cuadro operacion operaciones totales
n Gauss Cramer
10 805 399167999
100 681500 10162
1000 668165500
102573
El metodo de Gauss es el utilizado mas comunmente para resolver sistemas lineales
cuyas matrices no poseen propiedades particulares (sobre todo llenas).
Todo lo anterior nos permite enunciar el siguiente resultado (que no demostraremos):
Teorema 2 Sea A una matriz cuadrada, invertible o no. Existe al menos una matriz
invertible M tal que MA es triangular superior.
2.3.3. Eleccion practica del elemento de pivote
Como vimos anteriormente, el primer paso en el metodo de Gauss es elegir un pivote
no nulo (pues vamos a dividir por este elemento), lo cual corresponde a intercambiar
ecuaciones en el sistema lineal. En la practica no solo se busca un elemento no nulo, sino
que elegimos aquel que tenga mayor valor absoluto.
Ejemplo: Consideramos el sistema de matriz A =
2−26 1
1 1
, termino indepen-
diente b =
1
2
, y solucion u =
u1
u2
, donde la solucion exacta es
u1 = 2 − u2 = 1,00000001490116
u2 =1 − 2−25
1 − 2−26≈ 0,99999998509884
luego la solucion aproxima a u1 = u2 = 1.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
20
Resolvemos ahora con pivote de dos formas distintas:
a) Tomando 2−26 como pivote:
A =
2−26 1
1 1
→
2−26 1
0 1 − 226
≈
2−26 1
0 −226
b = 12 → 1
2 − 226 ≈ 1−226
Luego estamos resolviendo el sistema:
−226u2 = −226 → u2 = 1
2−26u1 + u2 = 1 → u1 = 0
b) Tomando 1 como pivote:
A =
2−26 1
1 1
→
1 1
2−26 1
→
1 1
0 1
−2−26
≈
1 1
0 1
b =
1
2
→
2
1
→
2
1 − 227
≈
2
1
luego el sistema que resolvemos es:
u1 + u2 = 2
u2 = 1
⇒ u1 = 1
Este ejemplo pone de manifiesto que los errores de redondeo con efecto desastroso
provienen de la division por pivotes “muy pequenos”.
Ejemplo: Aplicamos la estrategia de pivote (pivote parcial) a la matriz ampliada(matriz que posee como ultima columna el termino independiente): −1 13 2 | 1
2 0 −1 | 2
1 9 −2 | 0
Si aplicamos la estrategia de pivote parcial, obtenemos:
2 0 −1 | 2
−1 13 2
|1
1 9 −2 | 0
→
2 0 −1 | 2
0 13 3/2
|2
0 9 −3/2 | −1
→
2 0 −1 | 2
0 13 3/2
|2
0 0 −33/13 | −31/13
luego el sistema asociado tiene como solucion:
x =97
66, y =
1
22, z =
31
33.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
21
2.4. Metodos iterativos
Comenzamos enunciando algunos problemas que pueden presentar los metodos direc-
tos, como motivacion a la introduccion de los metodos iterativos que veremos a continua-
cion.
Problemas de los metodos directos para la resolucion de (SL).
1. Cuando el tamano de la matriz A es grande (n >> 100), la propagacion del error
de redondeo es tambien grande, y los resultados obtenidos pueden diferir bastante
de los exactos.
2. Muchas de las matrices que aparecen en (SL) son de gran tamano ( 100 000) pero
la mayorıa de sus elementos son nulos (Esto ocurre, por ejemplo, en la resolucion
de problemas mediante Elementos Finitos). Estas matrices reciben el nombre de
matrices vacıas o huecas, y se dan cuando numero de elementos no nulos es de orden
n.
a ) Si los elementos no nulos estan distribuidos alrededor de la diagonal principal,
son de aplicacion todavıa los metodos directos que conservan la estructura
diagonal, como LU .
b) Si no ocurre lo anterior, la matriz se dice que es dispersa, y al aplicarle los
metodos directos se produce un fenomeno conocido como rellenado (elementos
que eran nulos en la matriz A, ahora ya no lo son). Entonces, si no se realiza una
adaptacion de los metodos directos al caso de matrices dispersas los resultadosno van a ser, en general, buenos (No vamos a estudiar esa adaptacion).
Los metodos iterativos no tienen esos problemas porque se basan en la resolucion
(reiteradas veces) de sistemas diagonales o triangulares (por puntos o por bloques). Lo
que se intenta es que en cada iteracion el numero de operaciones sea mınimo.
2.4.1. Estudio general del Metodo Iterativo
Supongamos un (SL) A u = b, buscamos una matriz B ∈ Mn y un vector c ∈ RN de
forma que la matriz I − B sea inversible y que la unica solucion del sistema lineal
u = B u + c (I −B) u=c
(2.5)
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
22
es la solucion A u = b.
Tenemos entonces que encontrar la solucion de un sistema lineal se puede ver como un
problema de punto fijo (Problema (2.5)). En este caso, la forma de construir un metodo
iterativo es la siguiente:
Considermos u0 ∈ RN un vector arbitrario, se construye una sucesion de vectores
{uk
}∞k=0 dada por
uk+1 = Buk + c, k ∈ N ∪ {0} (2.6)
y se pretende que la sucesion {uk}k converja a la solucion del sistema lineal.
Definicion 7 El metodo iterativo (2.6) es convergente si existe un vector u ∈ RN tal
que:
lımk→+∞
uk = u
para cualquier vector inicial u0 ∈ RN . En ese caso,
u = B u + c
El error en cada iteracion se puede medir, por tanto como:
ek = uk − u
Se tiene que,
ek = uk − u = (Buk−1 + c) − (Bu + c) = B(uk−1 − u) = Bek−1 = . . . = Bke0
⇒ ek = Bke0.
De ese modo, el error en las iteraciones depende de las potencias sucesivas de la matriz B,
lo que nos dara el criterior para la convergencia del Metodo Iterativo.
Supondremos que A es no singular en el sistema.
Teorema 3 Sea B la matriz del Metodo Iterativo. Son equivalente:
1. El metodo iterativo es convergente.
2. ρ(B) < 1.
3. Existe una norma matricial | · | (que se puede tomar subordinada) tal que
|B| < 1
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
23
Recordemos que
ρ(B) = max1≤i≤n
{|λi(A)| : λi(A) ∈ sp(A)}es el radio espectral de B.
Teorema 4 Sea B ∈ Mn. Si exsite una norma matricial (subordinada o no) tal que
|B| < 1
entonces I + B es inversible y
|(I + B)−1| ≤ |I |1 − |B| .
Nota 10 1) Si I + B es una matriz singular, |B| ≥ 1, para toda norma matricial.
2) Si | · | es una normal matricial subordinada,
|(I + B)−1| ≤ 11 − |B| .
Ejemplo: La sucesion de vectores
vk =
2
k3, 1 − 1
k2, e1/k
T
∈ R3
es convergente al vector v = (0, 1, 1)T , v = lımk→+∞ vk.
El estudio de Metodos Iterativos reposa sobre la solucion de los dos problemas siguien-
tes:
1) Dada la mariz del metodo iterativo B, determinar si el Metodo Iterativo es conver-
gente.
2) Dados dos metodos iterativos convergentes, compararlos: el metodo iterativo mas
rapido es aquel cuya matriz tiene menor radio espectral.
Nota 11 1. Sea B una matriz real cualquiera y | · | una norma matricial (subordi-
nada o no). Entonces,
ρ(B) ≤ |B|
2. Dada una matriz B y > 0, existe al menos una norma matricial subordinada tal
que
|B| ≤ ρ(B) +
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
24
2.4.2. Metodos de Jacobi y Gauss-Seidel
Introducimos dos metodos iterativos clasicos para resolver el Sistema Lineal A u = b.
Todos se basan en descomponer la matriz del sistema A como suma de submatrices:
A = M − N,
donde M es una matriz inversible facil de invertir (en el sentido de que el sistema asociado
sea facil de resolver). Se verifica entonces
A u = b ⇔ (M − N ) u = b ⇔ M u = N u + b ⇔
u = B u + c con
B = M −1N
c = M −1b
Consideramos entonces el Metodo Iterativo:
u
0
∈ R
N
uk+1 = B uk + c, k ∈ N ∪ {0}
Como N = M − A, entonces B = M −1N = M −1(M − A) = I − M −1A. Luego,
I − B = M −1A
que es inversible, y por lo tanto el sistema (I − B) u = c tiene solucion unica.
En la practica, en vez de resolver uk+1 como uk+1 = Buk + c, se resuleve el
sistema
Muk+1 = Nuk + M c b
. (2.7)
Nota 12 La construcci´ on de la matrices M y N parece que “toma” elementos de A.
Parece que cuanto m´ as parecida sea la matriz M a A, m´ as se acercar´ a la soluci´ on a la
exacta (si N = 0 entonces M = A). Pero esto va en contra de que el sistema asociado a
la matriz M (2.7) sea f´ acil de resolver.
Introducimos la siguiente notacion:
Dada A = (ai,j )ni,j=1 ∈ Mn con ai i = 0, consideramos la siguiente descomposicion de
la matriz
A =
−F
D
−E
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
25
A = D − E − F , donde
D = diag(a1 1, a2 2, . . . , an n) E = (ei j )ni,j=1 F = (f i j)n
i,j=1
con
ei j = −ai j si i > j
0 si i ≤ j , f i j = −ai j si i < j
0 si i ≥ j .
La llamamos descomposicion D − E − F por puntos de la matriz A.
2.4.3. Metodo de Jacobi
Consiste en tomar M = D, N = E + F .
Ası pues,
A u = b ⇔ D u = (E + F ) u + b ⇔ u = D−1(E + F ) u + D−1 b
lo que conduce al metodo de Jacobi iterativo por puntos::u0 ∈ RN
uk+1 = D−1(E + F ) uk + D−1b k ∈ N ∪ {0}o equivalentemente
u0 ∈ RN
Duk+1 = (E + F )uK + b k ∈ N ∪ {0}La matriz de Jacobi (por puntos) es B = D−1(E + F ) = I − D−1A
Observemos que si
D =
a1 1
a2 2
. . .
an n
, entonces D−1 =
1/a1 1
1/a2 2
. . .
1/an n
,
luego,
D−1A =
1 a1 2/a1 1 . . . a1n/a1 1
a2 1/a2 2 1 . . . a2n/a2 2...
...
an 1/an n an 2/an n . . . 1
.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
26
Por tanto, queda
xk+11
xk+12...
xk+1
n
=
0 −a1 2/a1 1 . . . −a1 n/a1 1
−a2 1/a2 2 0 . . . −a2 n/a2 2...
...
−an 1/an n
−an 2/an n . . . 0
xk1
xk2...
xk
n
+
b1/a1 1
b2/a2 2...
bn/an n
.
de donde:
xk+1 j =
b j − a j 1xk1 − a j 2xk
2 − . . . − a j j−1xk j−1 − a j j+1xk
j−1 − . . . a j nxkn
a j j j = 1, . . . , n
Observemos que las n componentes del vector xk+1 se calculan simultaneamente a par-
tir de las componente de xk. Por eso el metodo de Jacobi tambien se conoce como el
metodo de iteraciones simultaneas.
La primera cuestion que nos planteamos es la convergencia del metodo. Observamosque
I − D−1A∞ = max1≤i≤n
n j=1,j=i
ai j
ai i
.Ejemplo 1:
A =
2 −2 0
2 3 −1
α 0 2
, α ∈ R.
D =
2 0 0
0 3 0
0 0 2
, E =
0 0 0
−2 0 0
−α 0 2
, F =
0 2 0
0 0 1
0 0 0
,
J = D−1(E + F ) =
1/2 0 0
0 1/3 0
0 0 1/2
0 2 0
−2 0 1
−α 0 0
=
0 1 0
−2/3 0 1/3
−α/2 0 0
.
Definicion 8 Una matrix A = (ai j)ni,j=1 ∈ Mn se dice que es diagonal estrictamente
dominante si |ai i| >
n j=1 j=i
|ai,j |, i = 1, . . . , n .
Teorema 5 Si A ∈ Mn es una matriz diagonal estrictamente dominante, el metodo
iterativo de Jacobi por puntos es convergente.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
27
2.4.4. Metodo de Gauss-Seidel
Una estrategia adecuada para mejorar la convergencia del metodo de Jacobi serıa
utilizar en el paso de calculo de la componente
uk+1i
las componentes calculadas hasta el momento:
{uk+11 , uk+1
2 , . . . , uk+1i−1 } en vez de {uk
1, uk2, . . . , uk
i−1}
Es decir, consiste en reemplazar el sistema correspondiente al metodo de Jacobi:
ai iuk+1i = bi −
i−1 j=1
ai j uk j −
n j=i+1
ai juk j , (2.8)
porai iu
k+1i = bi −
i−1 j=1
ai juk+1 j −n
j=i+1
ai j uk j , (2.9)
Matricialmente, las ecuaciones anteriores se escriben como:
D uk+1 = b + E uk+1 + F uk,
es decir,
(D − E ) uk+1 = F uk + b
Luego,
M = D − E, N = F
el metodo iterativo de Gauss-Seidel por puntos se escribe:
u0 ∈ Rn
uk+1 = (D − E )−1F uk + (D − E )−1 b, k ∈ N ∪ {0}
o equivalentemente u0 ∈ Rn
(D − E ) uk+1 = F uk + b k ∈ N ∪ {0}
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
28
La matrix de Gauss-Seidel (por puntos) es entonces:
B = (D − E )−1 F = I − (D − E )−1 A.
En este metodo, para calcular las componentes del vector uk+1, necesitamos tanto
las componentes de uk+1 ya calculadas, como las restantes del vector uk, por lo que se
denomina metodo de las aproximaciones sucesivas. Dicho metodo sera mas rapido ya que
la matriz M contiene mas elementos de A.
Ejemplo: Sea
A =
2 −2 0
2 3 −1
α 0 2
, α ∈ R.
GS = (D − E )−1 F = 2 0 0
2 3 0
α 0 2
−1 0 2 0
0 0 1
0 0 0
= 0 1 0
0 −2/3 1/3
0 −α/2 0
.
La siguiente tabla nos da los radios espectrales de las matrices de los metodos de
Jacobi y Gauss-Seidel para valores concretos del parametro α.
α ρ(GS ) ρ(J )
−1 0,848656 0,860379
−3 0,97263 1,11506
−5 1,08264 1,305158
Teorema 6 Si A es diagonal estrictamente dominante, entonces el metodo de Gauss-
Seidel es convergente.
Hay muchas otras condiciones de convergencia particulares para los distintos metodos.
Entre ellas, la mas significativa es:
Teorema 7 Si A es simetrica y definida positiva, el metodo de Gauss-Seidel es conver-
gente.
Definicion 9 Una matriz A se dice definida positiva si:
vtav > 0, ∀v ∈ Rn ∪ {0}
o equivalentemente si sp(A) ⊂ R+, todos sus autovalores son positivos (no nulos).
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
29
Problema:
Estudiar la convergencia de los metodos de Jacobi y Gauss-Seidel por puntos para las
matrices:
A =
1 2 −2
1 1 1
2 2 1
A =
2 −1 1
2 2 2
−1
−1 2
Problema: Metodo iterativo para el calculo de la inversa de una matriz.
Se consideran las sucesiones de matrices
An = An−1(I + E n + E 2n), E n = I − A An−1
siendo A ∈ Mn inversible y A0 ∈ Mn una matriz arbitraria.
1. Demostrar que E n = (E 1)3n−1
.
2. Probar que si ρ(E 1) < 1, entonces
lımn→+∞
An = A−1
3. Mostrar que si se toma A0 = A∗/tr(A A∗), entonces
lımn→+∞
An = A−1.
2.4.5. Test de parada de las iteraciones
Cuando un metodo iterativo es convergente, la solucion del sistema lineal A u = b
se obtiene como lımite de la sucesion {uk}∞k=0 de iteraciones. Ante la impoisibilidad de
calcular todas las iteraciones, se plantea el problema de determinar un valor k ∈ N ∪ {0}para el cual podemos considerar que uk sea una buena aproximacion de u. Es decir, si se
desea que el error relativo sea inferior a una cierta cantidad prefijada , se debe cumplir,
uk − u < u
para alguna norma vectorial.Sin embargo, como u es desconocido, no se puede trabajar con esa cantidad.
Si calculamos:
uk+1 − uk < uk+1puede que uk+1 no este proximo a u.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
30
Una condicion de parada de las iteraciones adecuada viene dada a partir del vector
residuo
rk = b − A uk = A(u − uk), k ∈ N ∪ {0}.
Entonces, si uk u y Auk b,
rk
b = A(uk
− u)Au <
Es decir, buscamos k ∈ N ∪ {0} tal que
rk < b.
Debe procurarse que la comprobacion de los tests de parada no incremente en exceso
el numero de operaciones necesarias para realizar una iteracion. Reescribiendo los calculos
de forma adecuada obtenemos:
1. Metodo de Jacobi:
D uk+1 = b + (E + F ) uk = b + (−A + D) uk =
= b − A uk + D uk = rk + D uk
Por tanto,
D (uk+1 − uk) = rk.
De ese modo, el metodo de jacobi se implementa de la siguiente forma:
1) Se calcula rk como rk = b − A uk.
2) Se resuelve el sistema D dk = rk.
3) uk+1 = uk + dk
De esta forma, se calcula el valor uk+1 a partir de uk y se aprovechan los calculos
intermedios, en concreto rk, para el test de parada. Tenemos el esquema:
rki = bi −n
j=1 ai j uk
j
dki =
rki
ai i
uk+1i = uk
i + dki
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
31
2. Metodo de Gauss-Seidel Analogamente, la implementacion del metodo de Gauss-
Seidel, se realiza en las siguiente etapas:
rki = bi −i−1
j=1 ai j uk+1 j −n
j=i+1 ai j uk j
dki =
erki
ai i
uk+1i = uk
i + dki
Nota 13 Las normas vectoriales que suelen usarse con mas frecuencia en los test de
parada son · 2 y · ∞.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
32
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Capıtulo 3
Resolucion numerica de ecuaciones
diferenciales ordinarias
El objetivo de este tema sera la resolucion numerica de ecuaciones diferenciales ordi-narias (e.d.o.). Nuestro problema de partida sera determinar de forma aproximada, me-
diante el uso de Metodos Numericos de Calculo, una solucion de una e.d.o. de primer
orden, conociendo el valor de la curva solucion en un punto.
Otros problemas que se estudiaran seran la resolucion numerica de ecuaciones diferen-
ciales de orden superior y problemas de valores frontera.
Interes: Las ecuaciones diferenciales se usan de forma habitual para construir modelos
matematicos en una amplia variedad de problemas de la ciencia y la ingenierıa. En dichos
problemas se buscan los valores de ciertas funciones desconocidas a traves de lo unico quesomos capaces de medir: como los cambios de una variable afectan a otra. Cuando esta
relacion de cambios se traducen a un modelo matematico, el resultado es una ecuacion
diferencial.
Ejemplo: Consideremos la temperatura de un objeto y(t) que se enfrıa. Podrıamos con-
jeturar que la velocidad del cambio de la temperatura del cuerpo esta relacionada con
la diferencia entre su temperatura y la del medio que lo rodea: los experimentos lo con-
firman y la ley del enfriamento de Newton establece que dicha velocidad de cambio es
directamente proporcional a la diferencia de estas temperaturas. Si denotamos por y(t) latemperatura del cuerpo en el instante t, y A la temperatura del medio que lo rodea,
∂y
∂t= −k(y − A).
donde k es una constante positiva, y el signo negativo indica que la temperatura decrece
33
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
34
Figura 3.1: Solucion de la e.d.o. (3.1)
si la temperatura del cuerpo es mayor que la del medio.
Si conocemos la temperatura del cuerpo, y0, en el instante t = 0, que se denomina
condicion inicial, entonces incluımos esa informacion en el enunciado del problema, de
manera que resolvemos:
∂y
∂t= −k(y − A) con y(0) = y0. (3.1)
La solucion se calcula a traves de la tecnica de separacion de variables, obteniendo:
y = A + (y0 − A) e−k t.
Cada eleccion de y0 nos da una solucion distinta. En la Figura 3.1 se muestran varias
soluciones del problema. Se observa que, cuando t crece, la temperatura del cuerpo se
aproxima a la temperatura ambiente. Si y0 < A el cuerpo se calienta, si y0 > A el cuerpo
se enfrıa.
Definicion 10 Problemas de valor inicial
Una soluci´ on de un problema de valor inicial
y(t) = f (t, y), con y(t0) = y0
en un intervalo [t0, t1] es una funci´ on derivable y = y(t) tal que
y(t0) = y0, y(t) = f (t, y(t)), para todo t ∈ [t0, t1].
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
35
Campo de direcciones o pendientes
En cada punto (t, y) del rectangulo R = {(t, y); a ≤ t ≤ b, c ≤ y ≤ d }
a b
c
d
t
y
la pendiente m de la solucion y = y(t) (derivada) se puede calcular mediante la formula
implıcita
m = f (t, y(t)).
Por tanto, cada valor mi,j = f (ti, y j), calculado para cada punto del rectangulo representa
la pendiente de la recta tangente a la solucion que pasa por el punto (ti, y j).
Un campo de direcciones o campo de pendientes es una grafica en la que se representan
las pendientes {mi,j} en una coleccion de puntos del rectangulo, y puede usarse para ver
como se va ajustando una solucion a la pendiente dada: Calculamos la pendiente en el
punto inicial (t0, y0) → f (t0, y0) para determinar en que direccion debemos movernos.
Damos un paso horizontal desde t0, t0 + h y nos desplazamos verticalmente una distancia
apropiada h f (t0, y0), llegando al punto (t1, y1) de manera que el desplazamiento total que
resulta tenga la inclinacion requerida. Una vez en el punto (t1, y1) se repite el proceso
a lo largo de la solucion. Como solo podemos dar un numero finito de pasos, el metodo
reproduce una aproximacion de la solucion.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
36
3.1. El metodo de Euler
Sea [a, b] el intervalo en el que queremos hallar la solucion de un problema de valor
inicial
y(t) = f (t, y(t))
y(t0) = y0
con f lipschitziana.
Construimos un conjunto finito de puntos {(tk, yk)} que son aproximaciones de la solucion
(y(tk) ≈ yk).
Problema: Construir dichos puntos verificando una ecuacion diferencial. Para f (t, y), t
seran las abcisas, y las ordenadas de los puntos (t, y). Dividimos el intervalo [a, b] en M
subintervalos (de igual tamano) en la siguiente particion:
tk = a + k h, k = 0, 1, . . . , M, h =b − a
M
.
El valor de incremento de h se llama tamano del paso.
Resolucion aproximada: De y(t) = f (t, y(t)) en [t0, tM ] con y(t0) = y0.
Suponiendo que y(t), y(t), y(t) son continuas y usando el Teorema de Taylor para
desarrollar y(t) alrededor del punto t = t0, para cada punto t existe un punto c1 entre t0
y t tal que
y(t) = y(t0) + y(t0)(t − t0) + y(c1)(t − t0)2
2.
Al sustituir y
(t0) = f (t0, y(t0)), h = t1 − t0, obtenemos una expresion para y(t1):
y(t1) = y(t0) + h f (t0, y(t0)) + y(c1)h2
2.
Si el tamano de paso es suficientemente pequeno, h2 se puede considerar despreciable
y
y(t1) ≈ y1 = y0 + h f (t0, y0)
que es la aproximacion de Euler.
Repitiendo el proceso, generamos una sucesion de puntos que se aproximan a la grafica
de la solucion y = y(t). El paso general del metodo de Euler es:
tk+1 = tk + h, yk+1 = yk + h f (tk, yk), k = 0, 1, . . . , M − 1.
Ejemplo:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
37
Usamos el metodo de Euler para hallar una solucion aproximada del problema de valor
inicial
y(t) = R y(t), t ∈ [0, 1], y(0) = y0, R = constante.
Debemos:
1. Elegir el tamano de paso h.
2. Usar la formula para calcular las ordenadas y(t), que se llama ecuacion en diferencia.
yk+1 = yk + h R yk = yk (1 + h R) = yk−1(1 + h R)2 =
= . . . = y0(1 + h R)k+1, k = 0, 1, . . . , M − 1.
En la mayorıa de los casos no se puede hallar una formula explıcita para determinar las
aproximaciones, pero este es un caso especial. Concretamente, es la
formula para calcular el interes compuesto a partir de un deposito inicial.
Ejemplo: Supongamos que se depositan 1000 euros durante 5 anos a un interes compuesto
del 10 %. Cual es el capital acumulado al cabo de esos 5 anos?.
y(t) = 0,1 y en [0, 5]
y(0) = 1000
y(t) = 1000 e0,1 t
yk = y0(1 + h R)k = 1000(1 + 0,1 h)k
Tamano paso Numero iteraciones M yM ≈ y(5)
1 5 1000(1 + 0,1)5
= 1610,511/12 60 1000
1 + 0,1
12
60= 1645,31
1/360 1800 100
1 + 0,1360
1800= 1648,61
La solucion exacta es y(5) = 1648,72 = 1000e0,5.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
38
Descripcion geometrica
3.1.1. Tamano del paso frente al error
Los metodos que presentamos se llaman metodos de difererencias o metodos de varia-
ble discreta. En ellos la solucion se aproxima en un numero finito de puntos llamados
nodos. Son de la forma:
yk+1 = yk + h φ(tk, yk) (3.2)
donde φ se llama funcion incremental, y es de paso simple (o de un solo paso) porque en
el calculo del nuevo punto solo interviene el punto inmediatamente anterior.
Cuando usamos metodos de variable discreta para resolver de manera aproximada unmetodo de valor inicial, existen 2 fuentes de error:
- discretizacion
- redondeo
Definicion 11 Error de discretizacion
Supongamos que {(tk, yk)}M k=0 es un conjunto finito de aproximaciones a la ´ unica so-
luci´ on y = y(t) de un problema de valor inicial.
Se define el error de truncamiento global o error de discretizaci´ on global ek como:
ek = y(tk) − yk, k = 0, 1, . . . , M .
que es la diferencia entre la soluci´ on exacta y la calculada con el metodo en el nodo
correspondiente.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
39
Se llama error de consistencia (o error de truncamiento local) εk+1 a:
εk+1 = y(tk+1) − y(tk) − h φ(tk, y(tk)), k = 0, 1, . . . , M − 1.
y es el error que se comete en un solo paso (el que lleva desde el nodo tk al tk+1).
En el metodo de Euler, en cada paso se desprecia un termino
y(ck)h2
2.
Si ese fuera el unico error que se comete en cada paso, al llegar al extremo superior del
intervalo (dar M pasos) el error acumulado serıa:
M k=1
y(ck)h2
2≈ My(c)
h2
2=
h M
2y(c)h =
=(b − a)
2
y(c) h = O(h).
Podrıa haber otros errores, pero esta estimacion es la que predomina.
Teorema 8 ( Precision del metodo de Euler)
Sea y(t) la soluci´ on de y(t) = f (t, y)
y(t0) = y0
Si y ∈ C 2([t0, b]) y {(tk, yk)}M k=0 es la sucesi´ on de aproximaciones generada por el metodo
de Euler, entonces
|ek| = |y(tk) − yk| = O(h)|εk+1| = |y(tk+1) − y(tk) − h f (tk, y(tk))| = O(h2)
El error al final del intervalo, llamado error global final, viene dado por
E (y(b), h) = |y(b) − yM | = O(h).
Nota 14 El error global final se usa para estudiar el comportamiento del error para ta-
ma˜ nos de paso diferentes, luego nos da una idea del esfuerzo computacional a realizar
para obtener las aproximaciones deseadas
E (y(b), h) ≈ c h
E (y(b),h
2) ≈
c
2h ≈
1
2E (y(b), h)
luego si reducimos a la mitad el tama˜ no de paso en el metodo de Euler, el error global
final tambien se reduce a la mitad.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
40
Figura 3.2: Aproximacion de Euler con paso h = 1 para resolver la e.d.o. (3.1)
Figura 3.3: Aproximacion de Euler con tamanos de paso diferentes para resolver la e.d.o.
(3.1)
La siguiente tabla compara las soluciones obtenidas con el metodo de Euler, con dife-
rentes tamanos, para y(t) = (t − y)/2 en [0, 2] con condicion inicial y(0) = 1:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
42
Tenemos que el error que comete es e c h p, y queremos determinar el valor de la
pontencia p. Entonces:
e1 = c h p1
e2 = c h p2→
e1e2
=h p1h p2
⇒ ln
e1e2
= p ln
h1
h2
⇒ p =
ln(e1/e2)
ln(h1/h2) .
Estimacion de la precision:
Si yh e yh/2 coinciden en n dıgitos, se puede suponer que esos n dıgitos son exactos.
Definicion 12 Se dice que un metodo es de orden mayor o igual a p ( p > 0) si para toda
soluci´ on y(t) del problema en [a, b], si y ∈ C p+1([a, b]) y existe k > 0 (que s´ olo depende
de y, φ, pero es independiente de h) tal que:
M −1n=0
|εn+1| ≤ k h p
Consistencia y estabilidad del metodo
Definicion 13 El metodo es consistente al problema y(t) = f (t, y) en [a, b]
y(a) = y0
si
lımh→0
M −1n=0
|y(tn+1) − y(tn) − h φ(tn, y(tn); h)| = 0,
es decir, seg´ un la definici´ on de εn+1, el error de consistencia,
lımh→0
M −1n=0
|εn+1| = 0
Definicion 14 Un metodo se dice estable si existe c > 0 (constante de estabilidad)
independiente de h, tal que
yk+1 = yk + h φ(tk, yk; h)zk+1 = zk + h φ(tk, zk; h) + ρk
0 ≤ k ≤ N, ∀yk, zk, ρk
max0≤k≤M
|yk − zk| ≤ c{ |y0 − z0| error inicial
+
0≤k≤N −1
|ρk|
suma de errores en cada etapa
}
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
43
Definicion 15 Un metodo es convergente a la soluci´ on del problema si
lımh→0
max0≤k≤M
|yk − y(tk)| = 0
Teorema 9 El metodo es consistente si y s´ olo si φ(t, y; 0) = f (t, y) ∀t ∈ [a, b], ∀y ∈ R.
Teorema10 Si el metodo es estable y consistente, entonces es convergente.
3.2. El metodo de Heun
La siguiente tecnica introduce una nueva idea en la construccion del algoritmo para
resolver el problema de valor inicial:
y(t) = f (t, y(t)) en [a, b],
y(t0) = y0
Para obtener el punto (t1, y1), usamos el Teorema fundamental del Calculo, integrando
y(t) en [t0, t1] de manera que: t1t0
f (t, y(t)) dt =
t1t0
y(t) dt = y(t1) − y(t0),
⇒ y(t1) = y(t0) +
t1
t0
f (t, y(t)) dt.
Usamos ahora un metodo de integracion para aproximar la integral que aparece a la
derecha de la expresion anterior. Usando la regla del trapecio con incremento h = t1 − t0,
el resultado es:
y(t1) ≈ y(t0) +h
2(f (t0, y(t0)) + f (t1, y(t1))) .
Ahora bien, como y(t1) (que es lo queremos aproximar) aparece tambien en la expresion
de la aproximacion, entonces aplicamos el metodo de Euler para aproximar la y(t1) que
aperece en f (t1, y(t1)):
y(t1) ≈ y(t0) + h f (t, y(t0)),
y obtenemos:
y(t1) ≈ y0 +h
2(f (t0, y0) + f (t1, y0 + h f (t0, y0)))
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
44
Repitiendo el proceso en cada intervalo [tk, tk+1], en cada paso la aproximacion dada
por el metodo de Euler se usa como una prediccion del valor que queremos calcular, y
luego la regla del trapecio se usa para hacer una correccion y obtener el valor definitivo.
El paso general del metodo de Heun es entonces el siguiente:
Metodo de Heun o Metodo predictor-corrector
Para tk+1 = tk + h,
pk+1 = yk + h f (tk, yk),
yk+1 = yk +h
2(f (tk, yk) + f (tk+1, pk+1))
La siguiente figura sintetiza el efecto del metodo anterior:
Tamano de paso frente al error
El termino del error de la regla del trapecio usada es:
−y(2)(ck)h3
12.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
45
Si el unico error que se cometiese en cada paso fuese el anterior, entonces despues de M
pasos del metodo de Heun el error acumulado serıa:
−M k=1
y(2)(ck)h3
12≈ M y(2)(c)
b − a
M
h2
12= y(2)(c) (b − a)
h2
12= O(h2)
El siguiente teorema establece la relacion entre el error final y el tamano de paso, y sepuede usar para dar una idea del esfuerzo computacional requerido para el metodo si
queremos obtener una precision fijada de antemano.
Teorema 11 (Precision del metodo de Heun)
Supongamos que y(t) es la soluci´ on de un problema de valor inicial
y(t) = f (t, y(t)),
y(t0) = y0.
Si y(t) ∈ C 3([t0, b]), y {(tk, yk)}N k=0 es la sucesi´ on de aproximaciones dadas por el metodo
de Heun, entonces:
|ek| = |y(tk) − yk| = O(h2),
|εk| = |y(tk+1) − y(tk) − h φ(tk, y(tk))| = O(h3),
donde φ(tk, y(tk)) =1
2(f (tk, y(tk)) + f (tk+1, y(tk) + h f (tk, y(tk)))). En particular, el error
global final en el extremo derecho del intervalo verifica:
E (y(b), h) = |y(b) − yM | = O(h2).
Los siguientes ejemplos ilustran el Teorema de la precision del metodo:Si usasemos tamanos de paso h y h/2, obtendrıamos:
E (y(b), h) ≈ C h2,
E (y(b), h2 ) ≈ C
h2
4=
1
4C h2
≈1
4E (y(b), h).
Entonces, si en el metodo de Heun el tamano de paso se reduce a la mitad, el
error global final se reduce a su cuarta parte.
Ejemplo: Usamos el metodo de Heun para resolver el problema:
y(t) =t − y
2, en [0, 3], con y(0) = 1,
y comparamos las soluciones obtenidas para tamano de paso h = 1, h = 12
, h = 14
, y
h = 18
. La solucion exacta es y(t) = 3 e−t/2 + t − 2.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
46
Para h = 14 , t0 = 0, t1 = 1
4 , y0 = 1, luego:
f (t0, y0) =0 − 1
2= −0,5,
p1 = y0 + h f (t0, y0) = 1 −1
4
1
2
= 0,875,
f (t1, p1) =0,25 − 0,875
2= −0,3125,
y1 = y0 +h
2(f (t0, y0) + f (t1, y1)) = 1 +
1
8(−0,5 − 0,3125) = 0,8984375.
Iterando en cada nodo, obtenemos:
y(3)≈
y12 = 1,672269.
La siguiente figura nos da idea de la aproximacion obtenida con el metodo de Heun
con tamano h = 1 y h = 1/2 para y(t) = (t − y)/2 en [0, 2] con condicion inicial y(0) = 1:
La siguiente tabla compara las soluciones obtenidas con el metodo de Heun, con dife-
rentes tamanos, para y(t) = (t − y)/2 en [0, 2] con condicion inicial y(0) = 1:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
48
3.3. Metodos de Runge-Kutta
Para disenar estos metodos, la pendiente usada es un promedio entre los valores de la
pendiente en el lımite izquierdo del intervalo [tk, tk+1] y en otros puntos intermedios:
yk+1 = yk + h f uncion promedio
donde funcion promedio = a f k + b f k con a y b pesos a elegir de la forma:
f k = f (tk, yk),
f k = f (tk + α h , yk + β h f k)
donde los parametros α, β especifican la posicion del punto intermedio.
Runge y Kutta disenaron el algoritmo eligiendo 4 parametros a, b, α y β de manera
que el resultado fuese lo mas preciso posible. Los parametros son independientes entre sı.
Las restricciones se obtienen desarrollando en serie de Taylor la funcion f en (t, y):a + b = 1
α b = β b = 12 ,
que es un sistema de 3 ecuaciones con 4 incognitas. Si elegimos a = 0, b = 1, α = β = 12 ,
obtenemos:
Euler modificado
yk+1 = yk + h f (tk + h2
, yk + h2
f (tk, yk))
Si a = b = 12 , α = β = 1, se obtiene el metodo de Heun o Euler mejorado. Dicho
metodo y el anterior son metodos de segundo orden.
3.3.1. Algoritmos de Runge-Kutta de cuarto orden
Si incluimos mas puntos en el muestreo dentro del intervalo, el metodo basico de
Runge-Kutta puede mejorarse para que el error de truncamiento global (final) sea propor-
cional a h4 (es decir, se trata de un metodo de cuarto orden). Despues de hacer desarrollos
en serie de Taylor, se obtiene:
yk+1 = yk +h
6(f 1 + 2 f 2 + 2 f 3 + f 4) ,
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
49
donde
f 1 = f (tk, yk)
f 2 = f (tk + h2 , yk + h
2 f 1)
f 3 = f (tk + h2 , yk + h
2 f 2)
f 4 = f (tk + h, yk + h f 3)
Nota 16 En general, los metodos de Runge-Kutta son de la forma:
y0 dado
yk+1 = yk + h
r
i=1
ωi ki
,
donde ki = f (tk + ci h, yk + h
i−1 j=1
aij k j
), con c1 = 0, ci ∈ [0, 1], ωi ∈ R, aij ∈ R.
Tamano de paso frente al error
El termino del error de la regla de Simpson con incremento h/2 es:
−y(4)(c1)h5
2880.
Si el unico error que apareciese en cada paso fuera el anterior, entonces en N pasos el
error acumulado al llevar a cabo el metodo serıa:
−
N k=1
y(4)(ck) h
5
2880 ≈ b − a5760 y(4)(c) h4 ≈ O(h4).
Teorema 12 (Precision del metodo de Runge-Kutta)
Supongamos que y(t) es la soluci´ on del problema de valor inicial
y(t) = f (t, y(t)),
y(0) = y0.
Si y(t) ∈ C 5([t0, b]) y {(tk, yk)}N k=0 es la sucesi´ on de aproximaciones generada por el meto-
do de Runge-Kutta de cuarto orden, entonces:
|ek| = |y(tk) − yk| = O(h4),
|εk| = |y(tk+1) − yk − h T r(tk, yk)| = O(h5),
donde T r(tk, yk) es la funci´ on promedio del metodo. En particular, el error global final
verifica:
E (y(b), h) = |y(b) − yM | = O(h4).
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Capıtulo 3
Parte II (e. d. o. s): Resolucion
numerica de ecuaciones diferenciales
ordinarias
3.4. Sistemas de Ecuaciones Diferenciales
Consideramos el problema de valor inicial:
(S )
dx
dt(t) = f (t, x(t), y(t))
dy
dt
(t) = g(t, x(t), y(t))
con
x(t0) = x0
y(t0) = y0
Definicion 16 Una soluci´ on de (S ) es un par de funciones derivables x(t), y(t), tales
que cuando t, x(t), y(t) se sustituyen en f (t, x(t), y(t)), g(t, x(t), y(t)) el resultado es igual
a las derivadas x(t), y(t) respectivamente, es decir:
x(t) = f (t, x(t), y(t))
y(t) = g(t, x(t), y(t))
con
x(t0) = x0
y(t0) = y0
Ejemplo: Consideramos el sistema
dx
dt(t) = x + 2y
dy
dt(t) = 3x + 2y
con
x(0) = 6
y(0) = 4
51
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
52
La solucion exacta es:
x(t) = 6e4t + 2e−t
y(t) = 6e4t − 2e−t
3.4.1. Resolucion numerica
Podemos encontrar una solucion numerica del sistema (S ) en un intervalo [a, b] consi-
derando los diferenciales:
d x = f (t,x,y) dt, dy = g(t,x,y) dt
Reescribiendo los incrementos dt = tk+1 − tk, dx = xk+1 − xk, dy = yk+1 − yk, es facilimplementar el metodo de Euler:
xk+1 − xk ≈ f (tk, xk, yk) (tk+1 − tk)
yk+1 − yk ≈ g(tk, xk, yk) (tk+1 − tk)
Dividiendo el intervalo [a, b] en N subintervalos de anchura h =b − a
N , y tomando los
nodos tk+1 = tk + h, obtenemos las formulas correspondientes al metodo de Euler:
Metodo de Euler
tk+1 = tk + h
xk+1 = xk + h f (tk, xk, yk)
yk+1 = yk + h g(tk, xk, yk)), para k = 0, 1,...,N − 1.
Sin embargo, para obtener un grado de precision razonable, es necesario utilizar un
metodo de orden mayor. Por ejemplo, las formulas para el metodode Runge-Kutta de
orden 4 son:
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
53
Metodo de Runge-Kutta de orden 4
xk+1 = xk +h
6(f 1 + 2 f 2 + 2 f 3 + f 4),
yk+1 = yk + h6
(g1 + 2 g2 + 2 g3 + g4) ,
donde
f 1 = f (tk, xk, yk), g1 = g(tk, xk, yk),
f 2 = f (tk + h2 , xk + h
2f 1, yk + h2g1), g2 = g(tk + h
2 , xk + h2f 1, yk + h
2g1),
f 3 = f (tk + h2 , xk + h
2f 2, yk + h2g2), g3 = g(tk + h
2 , xk + h2f 2, yk + h
2g2),
f 4 = f (tk + h, xk + h f 3, yk + h g3), g4 = g(tk + h, xk + h f 3, yk + h g3).
Ejemplo: Aplicamos al sistemax = x + 2y,
y = 3x + 2ycon
x(0) = 6,
y(0) = 4
el metodo de Runge-Kutta de orden 4 en el intervalo [0, 0,2] tomando 10 subintervalos de
paso h =0,2
10= 0,02.
Para obtener el primer punto t1 = 0,02, las operaciones intermedias necesarias paraobtener x1 e y1 son:
f 1 = f (t0, x0, y0) = f (0, 6, 4) = 14, g1 = g(t0, x0, y0) = g(0, 6, 4) = 26
x0 + h2
f 1 = 6,14 y0 + h2
= 4,26
f 2 = f (0,01, 6,14, 4,26) = 14,66, g2 = g(0,01, 6,14, 4,26) = 26,94
x0 + h2
f 2 = 6,1466 y0 + h2
g2 = 4,2694
f 3 = f (0,01, 6,1466, 4,2694) = 14,68,54, g3 = g(0,01, 6,1466, 4,2694) = 26,9786
x0 + h f 3 = 6,293708 y0 + h g3 = 4,539572
f 4 = f (0,02, 6,293708, 4,539572) = 15,372852 g4 = g(0,02, 6,293708, 4,539572) = 27,960268
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
55
representa un sistema mecanico en el que un muelle, cuya constante de recuperacion es
k, atado a una masa m, que ha sido separado de su posicion de equilibrio a la que
tiende a volver. Se supone que la amortiguacion debida al rozamiento es proporcional a
la velocidad, que g(t) es una fuerza externa, y que se conocen la posici on inicial x(t0) y
la velocidad inicial x(t0). Despejando la derivada segunda, el problema de valor inicial se
puede escribir como:
(E )
x(t) = f (t, x(t), x(t))
x(t0) = x0,
x(t0) = y0
Si llamamos y(t) = x(t), la e. d. o. de segundo orden se puede reescribir como un problema
de valor inicial para sistemas de primer orden con 2 ecuaciones:
(S )
dx
dt (t) = y(t),dy
dt(t) = f (t, x(t), y(t)),
x(t0) = x0,
y(t0) = y0
Al resolver el sistema (S ) con un metodo numerico de Runge-Kutta de orden 4, se generan
dos sucesiones {xk}, {yk}, siendo {xk} la sucesion de (E ).
Ejemplo 1: Movimiento armonico amortiguado.
x(t) + 4 x(t) + 5 x(t) = 0, x(0) = 3, x(0) = −5.
a) Reescritura como sistema equivalente:
x(t) = −4 x(t) − 5 x(t),
x(t) = y(t)
y(t) = −5 x(t) − 4 y(t)con
x(0) = 3
y(0) = −5
b)En la tabla siguiente mostramos los resultados de RK4 en el intervalo [0 , 5], con N = 50
y h = 0,1, y la comparacion con la solucion exacta x(t) = 3 e−2t cos(t) + e−2t sen(t):
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
57
La desviacion es similar a la de una viga en voladizo. Se puede usar la siguiente
ecuacion diferencial, basada en las leyes de la mecanica, para calcular la deflexion:
d2
dz2y(z) =
f
2 E I (L − z)2 (3.1)
donde E es el modulo de elasticidad, L es la altura del mastil e I es el momento de inercia.
En z = 0, dydz = 0. Calcular la deflexion en el tope del mastil en donde z = L usando
metodos analıticos y numericos. Supongase que el casco no gira.
Solucion: La solucion exacta es y(z) =f
24 E I (L − z)4 + C 1 z + C 2. Su derivada es:
y(z) =−f
6 E I (L − z)3 + C 1 ⇒ y(0) =
−f
6 E I L3 + C 1 = 0 C 1 =
f L3
6 E I
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
58
Para calcular la constante C 2, usamos la hipotesis logica de que y(0) = 0, es decir, el
mastil no se mueve en el sitio de union con el casco del barco, lo que hace que:
C 2 =f L4
24 E I ,
de manera que la solucion es:
y(z) =f
6 E I
1
4(L − z)4 + z L3 −
L4
4
,
y entonces,
y(L) =f L4
8 E I .
El modelo anterior es valido siempre que el intervalo de integracion [0, L] sea pequeno,
y la desviacion del mastil (que acabamos de calcular) tambien. Los valores de f y E se
basan en datos experimentales variables y difıciles de medir exactamente.
Reescribimos la ecuacion (3.1) como un sistema de e. d. o. de primer orden que resol-
vemos usando el metodo de Euler:
dy
dz= u,
du
dz=
f
2 E I (L − z)2
Usamos f = 50 libras/pie, L = 30 pies, E = 1,5 × 108 libras/pie, I = 0,06 pies4, yobtenemos que la desviacion el el extremo superior del mastil es y(30) = 0,5625 pies.
Concretamente, para distintos valores de h:
Tamano de paso de Euler y(30)
1.0 0.5744
0.1 0.5637
0.05 0.5631
Los resultados se pueden usar para propositos de diseno. Esto es valioso en el caso en que
la fuerza no es constante sino que varıa de forma complicada en funcion de la altura sobre
la cubierta del velero.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
59
3.5. Problemas de contorno
Otro tipo de ecuaciones diferenciales son de la forma:
x(t) = f (t, x(t), x(t)) para a ≤ t ≤ b,
con condiciones de contorno (o frontera) x(a) = α,
x(b) = β
Observemos que hemos sustituido las condiciones iniciales para x(t0) = x0, x(t0) = y0
por dos condiciones para x(t). Esto se conoce como problema de contorno o problema de valores frontera.
En este caso, antes de implementar el metodo numerico es necesario garantizar que el pro-
blema de contorno posee solucion. Para ello, usamos el siguiente resultado:
Teorema 13 (Problema de contorno) Supongamos que f (t,x,y) es una funci´ on con-
tinua en la regi´ on:
R = {(t,x,y) : a ≤ t ≤ b, −∞ < x < ∞, −∞ < y < ∞},
con derivadas parciales∂f
∂x,
∂f
∂ycontinuas en R.
Si f x(t,x,y) > 0 para todo (t,x,y) ∈ R y existe una constante M > 0 tal que
|f y(t,x,y)| ≤ M para todo (t,x,y) ∈ R, entonces el problema de contorno
x(t) = f (t, x(t), x(t))
x(a) = α, x(b) = β
tiene soluci´ on ´ unica x = x(t) en a ≤ t ≤ b.
Nota 17 Observemos que se ha usado que y = x(t) para la notaci´ on del teorema anterior.
3.5.1. Caso lineal
Corolario 1 (Problemas de contorno lineales) Supongamos que la funci´ on f del
Teorema anterior es lineal y se puede escribir de la forma:
f (t,x,y) = p(t) y + q(t) x + r(t)
y que sus derivadas parciales∂f
∂x= q(t),
∂f
∂y= p(t) son continuas en R (lo que garantiza
que | p(t)| = |f y| ≤ M = max [a, b]{| p(t)|}). Si q(t) > 0 para todo t ∈ [a, b], entonces el
problema de contorno lineal:x(t) = p(t) x(t) + q(t) x(t) + r(t),
x(a) = α, x(b) = β
tiene soluci´ on ´ unica en a ≤ t ≤ b.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Capıtulo 4
El metodo de Diferencias Finitas
La idea del metodo de Diferencias Finitas consiste en aproximar las derivadas que
aparecen en el problema de ecuaciones diferenciales ordinarias (e.d.o.) de forma que se
reduzca a resolver un sistema lineal. Una vez definido el sistema lineal se estudiara teniendoen cuenta los resultados de los Temas 1 y 2.
Comenzamos viendo el metodo de Diferencias Finitas para un problema de contorno
de segundo orden lineal. En concreto, consideramos la ecuacion lineal:
x(t) = p(t) x(t) + q(t)x(t) + r(t) t ∈ [a, b] (4.1)
con
x(a) = α, x(b) = β.
Recordamos que, segun vimos en el tema anterior, suponemos que
q(t) > 0 y que p(t) esta acotada,
bajo estas condiciones el problema de contorno tiene solucion.
Hagamos una particion de [a, b], donde:
a = t0 < t1 < . . . < tN = b,
h =
b− a
N t j = a + j h j = 0, 1, . . . , N
Usando las formulas de diferencias centradas para aproximar las derivadas tenemos (En
el Apendice de este tema se explica como se deduce la siguiente igualdad):
x(t j) =x(t j+1) − x(t j−1)
2 h+ O(h2) (4.2)
65
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
66
x(t j) =x(t j+1) − 2x(t j) + x(t j−1)
h2+ O(h2) (4.3)
Reemplazando (4.2) y (4.3) en (4.1), aproximando x j ≈ x(t j) ∀ j, obtenemos:
x j+1 − 2 x j + x j−1h2
+ O(h2) = p(t j)
x j+1 − x j−1
2 h+ O(h2)
+ q(t j)x j + r(t j) (4.4)
Eliminando los terminos de orden O(h2) en (4.4) e introduciendo la notacion p j = p(t j),
q j = q(t j), r j = r(t j), obtenemos la ecuacion en diferencias:
x j−1 − 2 x j + x j−1h2
= p jx j+1 − x j−1
2 h+ q j x j + r j (4.5)
que se usa para calcular aproximaciones numericas a la solucion de la ecuacion diferencial
(4.1).
Reagrupando, teniendo en cuenta que las incognitas son x j j = 1, . . . , N , tenemos el
sistemas lineal de ecuaciones:
−h
2p j − 1
x j−1 +
2 + h2q j
x j +
h
2 p j − 1
x j+1 = −h2 r j j = 1, 2, . . . , N − 1
x0 = α
xN = β
(4.6)
El sistema (4.6) es un sistema tridiagonal de N − 1 ecuaciones y N − 1 incognitas,
x1, . . . , xN −1 (pues x0 y xN son datos, las condiciones de contorno del problema). En
notacion matricial podemos escribirlo como
A x = b,
donde x =
x1...
xN −1
,
A =
2 + h2q1 p1h/2 − 1 0 . . . . . . 0
− p2h/2 − 1 2 + h2
q2 p2h/2 − 1 . . . . . . 0. . .
− p jh/2 − 1 2 + h2q j p jh/2 − 1. . .
− pN −1h/2 − 1 2 + h2qN −1
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
67
Ademas, si denotamos e0 = ( p1h/2 + 1)α, eN = (− pN −1h/2 + 1)β , tenemos
b =
−h2r1 + e0
−h2r2...
−h
2
r j...
−h2rN −1 + eN
Aplicamos ahora algunos de los resultados que hemos estudiado para resolucion numeri-
ca de sistemas lineales. Por ejemplo, veamos si podemos utilizar los metodos iterativos de
Jacobi y Gaus-Seidel. Uno de los criterios para ver si ambos metodos son convergentes es
probar que la matriz es estrictamente diagonal dominante. Para ello, tenemos que ver que
|2 + h2
q j | > |1 + p jh/2| + |1 − p jh/2|, j = 1, . . . , N − 1.
Por un lado tenemos que h es el tamano de paso en la discretizacion, por lo que lo podemos
tomar suficientemente pequeno, de forma que
1 − p jh/2 ≥ 0, 1 + p jh/2 ≥ 0 ∀ j
Como ademas, tenemos que q > 0, por hipotesis (lo suponıamos al principio del tema,
para poder asegurar que el problema de contorno tiene solucion), llegamos a que
|1 + p jh/2| + |1 − p jh/2| = 1 + p jh/2 + 1 − p jh/2 = 2
y
|2 + h2q j | > 2 (pues q j > 0).
Por lo tanto, la matriz del sistema es estrictamente diagonal dominante. Entonces, pode-
mos aplicar los metodos iterativos de Jacobi y Gaus-Seidel para resolver el sistema lineal
(no os parece apasionante, como justo la hipotesis que se necesita para asegurar la exis-
tencia de solucion en el problema de contorno (q > 0) es la que asegure que el sistema es
estrictamente diagonal dominante?, acabamos de relacionar dos problemas en apariencia
distintos).
En definitiva, de esta forma, resolviendo el problema lineal, con tamano de paso h,
conseguimos una aproximacion numerica: un conjunto finito de puntos {(t j, x j)}N −1 j=1 .
Si se conoce la solucion exacta x(t j), entonces podemos comparar x j con x(t j).
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
69
Se puede probar que las soluciones numericas tienen un error de O(h2): Veamos que si
reducimos el tamano del paso a la mitad, el error disminuye una cuarta parte (Recordamos
que si un metodo es de orden p, entonces e(h/2) = e(h)/2 p):
Vamos a fijarnos en la tabla anterior, por ejemplo en t = 2, vemos que el error que
cometemos en este punto, para h1 = 0,2 es 0,022533. Dividiendo esta cantidad por cuatro
y ası sucesivamente tenemos
0,022533 0,00563325 0,0014083125 0,0035207
mientras que los valores de la tabla son
0,022533 0,005588 0,001394 0,000348
Es un pequeno ejemplo, comprobando que el esquema es de orden dos.
Teorema 14 Si la soluci´ on del problema de contorno es suficientemente regular, x ∈C 4([a, b]), entonces el error e para el metodo de diferencias finitas satisface:
e∞ ≤ C h2,
es decir, es un metodo de segundo orden.
Realmente habra ocasiones en las que nos interesa conseguir metodos de mayor orden.
Con esta finalidad estudiamos el esquema de mejora de Richardson.
4.1. Esquema de mejora de Richardson
Vamos a mejorar la precision de las aproximaciones numericas anteriormente obtenidas
usando el esquema de mejora de Richardson para extrapolar los valores {x j,1}N j=0, {x j,2}
N j=0,
{x j,3}N j=0 {x j,4}
N j=0, correspondiente a los pasos h1, h2, h3 y h4, donde
h2 =h1
2, h3 =
h2
2, h4 =
h3
2,
con la finalidad de obtener 6 cifras de precision.Para ello, primero eliminamos los terminos de orden O(h2) de {x j,1}
N j=0 y O((h
2)2) de
{x j,2}N j=0, generando los valores
{z j,1} =
4x j,2 − x j,1
3
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
71
4.2. Otras condiciones de frontera.
Las condiciones de frontera pueden ser distintas. Concretamente, en lugar de tomar
x(a) = α, podemos tener x(a) y x(b) en lugar de x(b).
Supongamos por ejemplo que:
x(a) = α
x(b) = β o
x(a) = α
x(b) = β o
x(a) = α
x(b) = β
Ahora x0 y xN son tambien incognitas, luego debemos tener 2 ecuaciones mas en el sistema
resultante.
Supongamos que tenemos la condicion x(a) = α, podemos aproximar esta condicion
mediante:x1 − x−1
2h= α (4.7)
donde x−1 es un punto ficticio de la busqueda de x0. Entonces, escribiendo la ecuacion
aproximada por diferencias finitas para el nodo x0:
x1 − 2x0 + x−1h2
= r0 + q0x0 + p0x1 − x−1
2h
que reordenamos como:
(−1 −1
2hq0)x−1 + (2 + h2 p0)x0 + (−1 +
1
2hq0)x1 = −h2r0,
de donde, teniendo en cuenta que por (4.7)
x−1 = x1 − 2hα
La ecuacion queda
(2 + h2 p0)x0 − 2x1 = −(2 + h q0)h α − h2r0. (4.8)
Razonando asimismo en el extremo b, tenemos que si imponemos x(b) = β la ecuacion
resultante es
−2xN −1 + (2 + h2qN )xN = −h2rN + (2 − h pN )h β. (4.9)
De esta forma, las ecuaciones (4.8) y (4.9) se anadirian al sistema (4.6), obteniendo un
problema de (N + 1) incognitas con (N + 1) ecuaciones.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
72
Apendice: Aproximacion de la derivada de una funcion
a partir del desarrollo de Taylor
En esta seccion vamos a introducir aproximaciones de y (x) e y(x), a partir del desa-
rrollo de Taylor. Estas aproximaciones es la base de los metodos de diferencias finitas.
Consideramos en primer lugar el desarrollo de Taylor de tercer orden, para una funci ony(x). Entonces,
y(t) = y(t0) + y(t0)(t − t0) + y(t0)(t − t0)2
2+ y(c)
(t − t0)3
6
de donde tomando t0 = x y t = x + h, tenemos
y(x + h) = y(x) + y(x)h + y(x)h2
2+ y(c1)
h3
6
Tomando ahora, t0 = x, t = x − h, obtenemos:
y(x− h) = y(x) − hy(x) + y(x)h2
2− y(c2)
h3
6.
Restando ambas expresiones tenemos:
y(x + h) − y(x− h) = 2 h y(x) +h3
6(y(c1) + y(c2))
De esta ultima aproximacion deducimos
y
(x) =
y(x + h) − y(x− h)
2 h +
h2
12(y
(c1) + y
(c2))
Por lo tanto, tenemos una aproximacion de segundo orden a la derivada a la funcion en
un punto:
y(x) ≈y(x + h) − y(x − h)
2 h
Aproximacion de y(x)
Usando el desarrollo de Taylor de cuarto orden, llegamos a que
y(t) = y(t0) + y(t0)(t − t0) + y(t0) (t − t0)2
2+ y(t0) (t − t0)3
6+ yiv(c)(t − t0)4
4!
Si t = x + h, t0 = x
y(x + h) = y(x) + y(x)h + y(x)h2
2+ y(x)
h3
6+ yiv(c1)
h4
4!
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
73
Si t = x − h, t0 = x
y(x + h) = y(x) − y(x)h + y(x)h2
2− y(x)
h3
6+ yiv(c2)
h4
4!
Sumando ambas expresiones y despejando el valor de y (x) tenemos:
y(x) = y(x + h) − 2y(x) + y(x− h)h2
+ h2
4!(yiv(c1) + yiv(c2))
Luego usamos la aproximacion de segundo orden de y (x):
y(x) ≈y(x + h) − 2y(x) + y(x − h)
h2.
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Metodos Numericos de Calculo.
Resumen Metodo de Elementos Finitos
En este texto se aplica el metodo de elementos finitos al problema de contorno conuna sola variable
y(x) + q(x)y(x) = f (x), x ∈ [a, b]y(a) = 0, y(b) = 0.
(1)
El metodo de los Elementos Finitos se basa en sustituir en la formulacion variacionaldel problema la funcion incognita por una funcion interpolada y tomar como funcionestest las funciones base del espacio de interpolacion. A continuacion desarrollamos cadauno de estos pasos. En la Etapa 1, se deduce la formulacion variacional del problema.En la Etapa 2, se define el espacio de interpolacion, las funciones base y la funcioninterpolada. En la Etapa 3 se aplica el metodo de los Elementos Finitos para deducirun sistema lineal. Finalmente en la Etapa 4 reescribimos el sistema en forma matricial.
ETAPA 1 Formulacion variacional del problema.
Sea v(x) una funcion regular definida en [a, b] y que verifica las condiciones de
contorno del problema, es decir,
v(a) = 0, v(b) = 0.
para buscar la formulacion variacional del problema basta multiplicar la ecuacion (1)por v(x), integrar en [a, b] e integrar por partes. Tenemos entonces
−
ba
y(x)v(x)dx + ba
q(x)y(x)v(x)dx = ba
f (x)v(x)dx. (2)
Donde hemos aplicado integracion por partes y las condiciones de contorno que verificala funcion test (v(a) = v(b) = 0). Es decir,
ba
y(x)v(x)dx = y(x)v(x)ba
−
ba
y(x)v(x)dx = −
ba
y(x)v(x)dx
ETAPA 2 Interpolacion.
Fijamos un espacio de interpolacion. Consideramos en este resumen el espacio defunciones lineales a trozos. A continuacion construimos una funcion lineal a trozos g(x)que verifique
g(xk) = yk
para ciertos puntos xk del intervalo [a, b].
1
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
Paso 1 : Construimos una particion del intervalo [a, b]. Consideramos una particionregular de paso h = (b− a)/M , siendo M el numero de subintervalos. Entonces,
xk = a + k ∗ h k = 0, . . . , M
Observamos que x0 = a y xM = b.
Paso 2 : Definimos las funciones base gk(x): funciones lineales a trozo que valen 1 sobreel nodo xk y cero en cualquier otro nodo. Entonces, la funcion viene definida por
gk(x) =
0 si x ≤ xk−1
(x − xk−1)/h si xk−1 ≤ x ≤ xk
(xk+1 − x)/h si xk ≤ x ≤ xk+1
0 si x ≥ xk+1
(3)
Comprueba que efectivamente gk(x) verifica que
gk(x j) =
1 si j = k0 si j = k
Paso 3 : Definimos la funcion g(x) tal que g(xk) = yk k = 1, . . . , M − 1. Basandonos enlas funciones base gk(x) es sencillo construir g(x) ya que gk(x) vale 1 en x = xk
y 0 en cualquier otro nodo de la partici on. Tenemos entonces
g(x) =M −1k=1
ykgk(x). (4)
Observamos que el sumatorio es desde k = 1 hasta k = M − 1, no se incluyenlos valores k = 0 y k = M . Esto es debido a las condiciones de contorno. Elproblema sera aproximar los valores de la funcion incognica y(x) sobre los puntosde la particion xk k = 0, . . . , M . A la aproximacion de y(xk) la denotamos poryk. Pero observamos que x0 = a y como y(a) es un dato entonces y0 = y(a) = 0.Analogamente xM = b, por lo que ym = y(b) = 0.
ETAPA 3 Reduccion a un sistema lineal de ecuaciones.
El metodo de los Elementos Finitos se basa en sustituir en la formulacion variacionaldel problema (2) la funcion y(x) por una funcion que la aproxime en un espacio deinterpolacion g(x). Ademas, como funciones test v(x) se tomaran las funciones basedel espacio de interpolacion gk(x).
Es decir, en primer lugar sustituimos en la formulacion variacional (2) y(x) porg(x):
−
ba
g(x)v(x)dx + ba
q(x)g(x)v(x)dx = ba
f (x)v(x)dx. (5)
teniendo en cuenta la definicion de g(x) (4) tenemos
−M −1
k=1
yk b
agk(x)v(x)dx +
M −1
k=1
yk b
aq(x)gk(x)v(x)dx =
b
af (x)v(x)dx. (6)
5/11/2018 C lculo num rico y diferencias finitas-V zquez -UPS - slidepdf.com
http://slidepdf.com/reader/full/calculo-numerico-y-diferencias-finitas-vazquez-ups
En segundo lugar, como funciones test tomamos v(x) = g j(x) para j = 1, . . . , M − 1.Tenemos entonces
−M −1k=1
yk
ba
gk(x)g j(x)dx+M −1k=1
yk
ba
q(x)gk(x)g j(x)dx = ba
f (x)g j(x)dx, j = 1, . . . , M −1.
(7)Observamos que (7) es un sistema lineal de M −1 incognitas (y1, . . . , yM −1) con M −1ecuaciones ( j = 1, . . . , M − 1).
ETAPA 4 Sistema lineal.
Comenzamos sacando factor comun en la incognita yk. Tenemos
M −1k=1
−
ba
gk(x)g j(x)dx+ ba
q(x)gk(x)g j(x)dx
yk = ba
f (x)g j(x)dx, j = 1, . . . , M −1.
(8)Por lo tanto, los elementos que definen la matriz del sistema son:
−
b
agk(x)g j(x)dx y
b
aq(x)gk(x)g j(x)dx
Observamos que el producto gkg j o gkg j es casi siempre cero. En concreto es distintode cero solo para j = k − 1, j = k y j = k + 1, debido a la definicion de gk (3). Porlo tanto, en el sumatorio que aparece en (8) solo tres terminos son no nulos, los queinvolucran a yk−1, yk y yk+1. Por lo tanto, obtenemos que la matriz del sistema estridiagonal.
Denotemos
p jk = − b
a
g j(x)gk(x)dx, q jk = b
a
q(x)g j(x)gk(x)dx, f j = b
a
f (x)g j(x)dx.
Entonces tenemos el sistema lineal
Ay = b,
donde el vector de incognitas es
y =
y1..
yM −1
La matriz del sistema es
A =
p11 + q11 p12 + q12 p21 + q21 p22 + q22 p23 + q23
. . . . . . . . . pM −1,M −2 + qM −1,M −2 pM −1,M −1 + qM −1,M −1
y el termino independiente es
b =
b1..
bM −1
Top Related