Raíces de funciones no
lineales
Introducción a matrices y
sistemas lineales
UNSAM ECyT
Métodos Numéricos
[Link]
Raíces
Sea f: [a, b] → |R (o bien f: |R → |R).
Una raíz de f es un r ∈ [a, b] tal que f(r) = 0.
Ejemplos (con funciones lineales y no
lineales).
Idea de invertir la función.
Dificultad en general.
Objetivos
Queremos hallar raíces de funciones reales
f : [a, b] → |R o f : |R → |R (con otras
condiciones a determinar luego).
Nos importa aproximar las raíces.
No nos importará cuántas raíces habrá en
total, sólo aproximar una de ellas.
Asumimos que podemos calcular f(x) para
cualquier x. Esto es importante para los
algoritmos.
Método de bisección
Ideas
Juego de pensar un número.
¿Qué preguntas hago?
Búsqueda en un arreglo ordenado.
¿Cómo minimizamos las comparaciones?
Búsqueda binaria
Otros métodos
Método de bisección
Requisito
Que se pueda evaluar f(x) en cualquier punto real
x deseado dentro del dominio de definición de f.
Ejemplos
x3 + x2 + 4, x sin x, x3 + cos x, √x, log x, 1 / (x4 + 1), etc.
Idea
Ir evaluando la f en distintos puntos, corriendo de
un lado a otro… “hasta que se llegue a algo”.
Teorema de Bolzano
Si f : [a, b] → |R es continua en el intervalo
cerrado real [a, b], tal que
f (a) ≤ 0 ≤ f (b) o bien f (b) ≤ 0 ≤ f (a),
entonces existe r ∈ [a, b] tal que f (r) = 0.
En otras palabras, una función continua en
un intervalo cerrado y acotado tal que f(a)
y f(b) tienen distinto signo, tiene al menos
una raíz dentro de ese intervalo.
Teorema del valor intermedio
(Es equivalente al anterior)
Si f : [a, b] → |R es continua en el intervalo
cerrado real [a, b], e y0 ∈ |R es tal que
f (a) ≤ y0 ≤ f (b)
entonces existe x0 ∈ [a, b] tal que f (x0) = y0.
En otras palabras, una función continua en
un intervalo cerrado y acotado toma todos
los valores entre f(a) y f(b) al menos una vez.
Ejemplos
x + cos x a = 0,
x5 + x4 + 2 x3 – 100 a = -100, b = 100
cos x
x+1
…
f
x0 x2 x3 x1
Algoritmo de bisección
// se tiene una función f : R -> R continua
// input: dos reales a, b tales que a < b y f(a)*f(b) < 0, y un real epsilon > 0
float biseccion() {
while (b-a >= epsilon) {//o bien otras condiciones...
c=(a+b)/2;
if (f(c)==0) return c;
else if (sgn(f(c))*sgn(f(a)) < 0)
b=c;
else
a=c;
}
return c;
}
int sgn(float z) {
if (z==0) return 0;
else if (z>0) return 1;
else return -1;
}
Criterios de parada
Sea x0, x1, x2, … la sucesión de puntos obtenida.
n ≥ N, fijando de entrada N
(nro. de pasos prefijado)
(bn - an) < ε, fijando ε > 0 de entrada
(error absoluto)
|pn – pn-1| / pn < ε , fijando ε > 0 de entrada
(error relativo, sólo utilizable si pn ≠ 0)
|f(pn)| < ε, fijando de entrada ε > 0; pero este criterio solemos
descartarlo… |f(pn)| “chico” no significa estar “cerca de una raíz”.
Pregunta
Supongamos que alguien nos plantea este otro
posible criterio de parada:
|pn - r | < ε, fijando ε > 0 de entrada, donde r es la raíz
buscada
Es decir, comparar pn con r y parar cuando estemos
cerca de r.
¿Qué decimos a esto?
Respuesta: ¡no sirve, ya que no conocemos r!
Parada
Elegimos el criterio de parada estándar:
(bn - an) < ε
¿Cuántos pasos tendremos que dar con el método
de bisección para llegar a esa condición?
Veamos que podemos calcular esto.
Pasos
a0 = a, b0 = b, los dos puntos iniciales
Cada intervalo es la mitad del intervalo del paso
anterior.
|b1 - a1| = ½ |b0 - a0| = ½ |b - a|
|b2 - a2| = ½ |b1 - a1| = ½ (½ |b0 - a0|) = ¼ |b - a|
|b3 - a3| = ½ |b2 - a2| = … = ½ ( ¼ |b - a|) = 1/8 |b - a|
…
|bn - an| = 1/2n |b - a| al cabo de n pasos
Error
|bn - an| = 1/2n |b - a|
Queremos que el error de la aproximación buscada sea < ε, entonces
podemos plantear:
1/2n |b - a| < ε
Como deseamos averiguar el número de pasos, despejamos n:
|b - a| / ε < 2n
Tomando logaritmos, se mantiene la desigualdad.
log2 (|b - a| / ε) < n
Conclusión:
Tomando un n > log2 (|b - a| / ε), se tiene una aproximación con la
precisión buscada, es decir error < ε.
Es decir, si ejecutamos esa cantidad de pasos nos aseguramos un
error < ε. Obsérvese que en general, puede pasar que no
conozcamos el error exacto que se comete con pn. Sólo podemos
asegurar que aquel sea menor que la cota ε.
Observaciones
n > log2 (|b - a| / ε)
En realidad tomaremos como n un entero
mayor que este valor.
Es una cota fácil de calcular.
Es una cota independiente de la función f.
Observar que, cuando |b - a| es grande, n
crece, y cuando ε es chico, n también crece.
Convergencia
Llamemos r a una raíz de f en el intervalo [a, b], suponiendo que allí
hay una sola.
Sea pn el punto dado por el método en el paso n.
Resulta pn = (an + bn) / 2 ya que es el punto medio del intervalo actual.
|pn – r| ≤ |bn – an| ≤ 1/2n (b – a) → 0 cuando n → ∞
Luego, |pn – r| → 0.
Esto dice que pn → r, es decir pn tiende a la raíz r.
Esto prueba la convergencia del método.
(Si, en cambio, hubiera varias raíces en [a,b], se podría explicar, pero
lo omitimos…)
Convergencia
Veamos que tanto an como bn también convergen a
una raíz.
Primero, una observación: |pn - an| ≤ 1/2n (b - a).
Luego,
|r - an| ≤ |r – pn + pn – an|
≤ |r – pn| + |pn – an| (desigualdad triangular)
≤ 1/2n (b – a) + 1/2n (b – a) (por la observación)
= 1/2n-1 (b – a) → 0 cuando n → ∞.
Luego |r – an| → 0, es decir an → r.
Análogamente se prueba que bn → r.
Método de Newton
Idea
En lugar de tomar el punto medio, tomar el
próximo punto siguiendo la tangente a la curva
desde esta evaluada en el punto anterior.
Esto obliga a que la función tenga tangente
definida en cada punto, al menos dentro de cierto
intervalo.
Newton
f(x0)
f(x0)
x1 r x0
x0 - x1
f (x0)
tangente = = f ’(x0)
x0 - x1
Newton
f (x0)
tangente = = f ’(x0)
x0 - x1
f(x0) = f ’(x0) (x0 - x1)
f (x0)
= x0 - x1
f ’(x0)
f (x0)
x1 = x0 -
f ’(x0)
Newton
f (x0)
x1 = x0 -
f ’(x0)
f (xn)
xn+1 = xn -
f ’(xn)
Ejemplo del método de Newton
• Hallar una aproximación de √ 2
• Podemos lograrlo exhibiendo una función
que pueda calcularse con facilidad y que a
la vez tenga como raíz ese valor (por
ejemplo)
• Una función de este tipo es: f(x) = x2 - 2
Ejemplo: aproximar √2
f (x) = x2 – 2
f ’(x) = 2x
x0 = 1 (que sea un valor cercano a la raíz buscada…)
f (xn) xn2 - 2
xn+1 = xn - xn+1 = xn -
f ’(xn) 2 xn
Ejemplo: aproximar √2
xn2 - 2
xn+1 = xn -
2 xn
xn2 2
xn+1 = xn - +
2 xn 2 xn
xn 1
xn+1 = + = ½ (xn + 2/xn)
2 xn
Observaciones
• Notar que se parte de un número x0, que sea
una aproximación a ojo de la raíz que busca,
vale decir, que esté “cerca” de la raíz. Y en los
sucesivos pasos del método, va encontrando
otras aproximaciones, de la misma raíz, cada
vez mejores. A esas otras aproximaciones las
va llamando x1, x2, x3, ... etc. Conocida x0 se
calcula x1. Conocida x1 se calcula x2, y así
siguiendo, cada una usa el valor de la anterior.
Por eso se tendrá una fórmula que es xn+1 = ...
algo que depende de xn es decir de la
aproximación obtenida en el paso anterior.
Observaciones
• Este método no garantiza llegar a una raíz, ni converger
a ella. A veces funcionará y a veces no, dependiendo de
la f, y del punto de partida x0.
• Al igual que pasa con el método secante, el método de
Newton a priori no garantiza nada. Pero lo bueno es
que, cuando converge, ¡converge a una raíz de f más
rápidamente que lo que converge el método de
bisección! Es decir, el error en cada paso va tendiendo a
0 más rápido que el error en cada paso del otro método.
• En la práctica se puede comenzar con un x0 cualquiera
que uno crea que esté cerca de la raíz buscada, según
el caso particular de la función de que se trate o, en
última instancia, midiendo a ojo.
Teorema
Si f es C2 (segunda derivada continua) en un
entorno cerca de una raíz r de f, y f ’(x) ≠ 0 en
un entorno de r, entonces existe ε > 0 tal que, si
tomamos cualquier x0 en el intervalo (r-ε, r+ε), el
método de Newton produce una sucesión xn que
converge a r (cuadráticamente).
Dem. (Ver Zarowski, usa el Teorema del Punto
Fijo de Banach)
Observaciones
• Es importante que el método necesita que la
función f sea derivable. Porque en cada paso,
como verán, se divide por f ' (xn) (es decir la
derivada de f aplicada a la aproximación
anterior). Por eso, el método no sirve
indefectiblemente para cualquier función
continua, como sí servía bisección. La f tiene
que ser derivable (al menos en una región cerca
de la raíz buscada...) Tiene entonces esta
desventaja.
Observaciones
• Para este método hay un teorema principal que dice que cuando la f
tiene una raíz, y además f es derivable con continuidad dos veces
(o sea f tiene derivada segunda y además esta es continua) cerca
de la raíz buscada, y cerca de esa raíz buscada la función f ' no se
hace 0, el método de Newton convergerá cuando se parta de un
valor x0 suficientemente cerca de esa raíz. ¡Interesante!
• Dado que no podemos conocer el valor de la raíz de antemano,
este teorema parecería de poca utilidad...
• Sin embargo se puede desarrollar una estrategia adecuada. Si no
se dispone de un buen valor inicial x0, se puede utilizar el método
de bisección para conseguir un valor suficientemente cerca de la
solución y, posteriormente, continuar resolviendo el problema pero
ahora por el método de Newton. Es decir, seguir iterando desde allí
con este método.
Implementación
• Se podría decir varias cosas.
• Hay que tener forma de calcular f(x) y f '(x)
para los xn que se vayan obteniendo.
• Se puede utilizar cualquier criterio de
parada (de acuerdo al error deseado)
similares a los que vimos para el método
bisección.
Proposición
Supongamos que f : |R → |R es C1 (es decir
derivable con derivada continua en todo
R), y f ’(x) ≠ 0 ∀x∈|R. Si el esquema
iterativo de Newton converge, entonces
converge a una raíz de f.
Dem. f (xn)
xn+1 = xn -
f ’(xn) (sigue)
Demostración
Llamemos x* = lim xn. Tomamos límite a
ambos lados de la igualdad. Como f es
derivable entonces es continua, y
f (xn)
xn+1 = xn -
f ’(xn)
por ser f y f’ continuas, tienden a
f (x*) f (x*)
x* = x* - ⇒ 0= -
f ’(x*) f ’(x*)
⇒ 0 = f(x*)
Proposición (consistencia)
Supongamos que f : |R → |R es C1. En el
esquema iterativo de Newton, si existe n
tal que xn = xn+1, y f ’(xn) ≠ 0, entonces
xn = xn+1 = xn+2 = … es una raíz de f.
Dem. f (xn)
xn = xn+1 = xn -
f ’(xn)
f (xn) f (xn)
xn = xn - ⇒ 0=-
f ’(xn) f ’(xn)
⇒ 0 = f(x*)
Método Secante (breve reseña)
f (xn) (xn - xn-1)
xn+1 = xn –
f (xn) – f (xn-1)
Secante
f (xn) (xn - xn-1)
xn+1 = xn -
f (xn) – f (xn-1)
Equivale a
f (xn)
xn+1 = xn -
f (xn) – f (xn-1)
xn – xn-1
Escrito así se observa que es semejante al método de Newton, ya que el
denominador es una aproximación a la derivada f ’(xn)
Cara o… secante
• Secante no necesita que f sea derivable, sólo
que f sea continua; mientras que Newton sí
necesita que sea derivable.
• Secante hace más operaciones en cada paso.
• Secante necesita dos puntos de partida, x0 y
x1, mientras que Newton necesita uno solo, x0.
• Ninguno de ambos métodos, secante y Newton,
garantizan la convergencia para cualquier
situación.
Orden de convergencia
Supongamos que pn → p, δ > 0.
Si existe λ > 0 tal que para todo n
lim |pn+1 – p|
δ
= λ
n→∞ |pn – p|
se dice que (pn) converge con orden δ.
δ = 1: lineal, δ > 1: superlineal, δ = 2: cuadrático
Ejemplos: 1/n converge a 0 linealmente,
y lo mismo pasa con las sucesiones
2/n, 3/n, 1/n2, 1/2n…
Orden de convergencia
• Bisección: orden lineal (demostración)
• Secante: orden 1+√5 (δ resulta ser el
número de oro). 2
• Newton: orden 2 (cuadrático), siempre que
se cumplan las condiciones del teorema…
Ejercicios “de raíz”
• Defina un esquema que utilice el método
de Newton para aproximar √3.
• El método secante necesita dos
aproximaciones iniciales x0 y x1.
Supongamos que tenemos sólo una, x0.
¿Cómo podemos proceder para tener otra
y poder aplicar el método?
Ejercicios más… lineales
• ¿Qué ocurrirá al aplicar el método de
Newton para hallar la raíz de una función
lineal? Es decir, f(x) = ax + b, con a y b
constantes, a ≠ 0.
• Misma pregunta para el método secante.
Recordar que los métodos antedichos
están para ser usados cuando el despeje
de x sea difícil o imposible.
Sistemas lineales y matrices
• Ecuaciones lineales
• Espacios vectoriales: |Rn y |Rm×n
• Vectores
• Suma y producto por escalar
• Suma, neutro, inverso aditivo
Suma de matrices
en detalle
A = (ai,j) ∈ |Rm×n (ai,j)1≤i≤m,1≤j≤n
B = (bi,j) ∈ |Rm×n (bi,j)1≤i≤m,1≤j≤n
A + B = (ai,j + bi,j) ∈ |Rm×n
(ai,j + bi,j) 1≤i≤m,1≤j≤n
Propiedades de la suma
A+B=B+A
(A + B) + C = A + (B + C)
A+0=A=0+A
λ(A + B) = λA + λB (λ∈R)
A + (-A) = 0 = (-A) + A
Inversa aditiva: -A = (-ai,j) ∈ |Rm×n
Producto de matrices
A = (ai,j) ∈ |Rm×n
B = (bi,j) ∈ |Rn×p
AB = (ci,j)
n
Donde ci,j = Σ (ai,k bk,j) ∈ |Rm×p
k=1
Ejemplos…
Propiedades del producto
(A B) C = A (B C)
A (B+C) = AB + AC, (A+B) C = AC + BC
λ(A B) = (λA) B = A (λB) (λ∈R)
A0=0=0A
AI=A=IA
Para matrices cuadradas A y B:
AB ≠ BA en general
Notaremos An = AA…A (n factores iguales a A)
Inversa
Sea A ∈ |Rn×n. Si existe A’ ∈ |Rn×n tal que
AA’ = A’A = I, a la matriz A’ la llamaremos
inversa de A.
Si A’ existe, es única, y en ese caso la
notamos con A-1.
-1
Ejemplo: 1 2 1 -2/3
0 3 0 1/3
Transpuesta
Si A = (ai,j) ∈ |Rm×n,
se define At = (aj,i) ∈ |Rn×m
Observación: (At)t = A
Una matriz A es simétrica si A = At.
Una matriz A es antisimétrica si A = -At.
Ejemplos…
Ejemplos
• Ejemplos de matrices que tienen inversa:
• Ejemplos de matrices que no tienen inversa:
11 123 00 01
11 369 00 00
112
00 234 012 123
23 234 034 124
103 056 125
Propiedades de la inversa
• Si A es inversible, (A-1)-1 = A
• Si A y B son inversibles, (AB)-1 = B-1A-1
• Si A es inversible y λ ≠ 0, (λA)-1 = λ-1A-1
• Si A es inversible, (At)-1 = (A-1)t
Determinante (detalles en clase)
• Conjunto de matrices inversibles de reales de
n×n: Gl(n, |R).
• El determinante es una función det : |Rn×n → |R
• Los casos 2×2 y 3×3…
Teorema: A ∈ Gl(n,|R) ⇔ det A ≠ 0.
Sistemas lineales
Tener m ecuaciones lineales con n
incógnitas equivale a tener una ecuación
de matrices y vectores de la forma:
Ax=b
donde
A ∈ |Rm×n , b ∈ |Rn, x ∈ |Rn
Ecuaciones lineales
• Métodos conocidos:
– Método de determinantes o regla de
Cramer
– Otros…
• Veremos:
– Método de descomposición LU
– Método iterativo de Jacobi
–…
Ecuaciones lineales
• Otra forma es plantearlo matricialmente:
Ax=b
A-1 (A x) = A-1 b
(A-1 A) x = A-1 b
I x = A-1 b
x = A-1 b
Esto da un método, pero… invertir una
matriz suele ser mucho esfuerzo.
Método de determinantes
Dado el S.E.L. A x = b
con A ∈ |Rn×n inversible, b ∈ |Rn,
el vector solución x ∈ |Rn está dado por
det Ai i = 1, 2, …, n
xi =
det A
donde Ai es en la matriz A pero cambiando la
columna i por el vector b
Ejemplo
Matrices triangulares
Sea A ∈ Rn×n (una matriz cuadrada)
• A = (ai,j) es triangular inferior si
∀1 ≤ i,j ≤ n (i < j ⇒ ai,j = 0)
• A = (ai,j) es triangular superior si ∀1 ≤ i ≤ n
∀1 ≤ i,j ≤ n (i > j ⇒ ai,j = 0)
• A = (ai,j) es estrictamente triangular inferior si
∀1 ≤ i,j ≤ n (i ≤ j ⇒ ai,j = 0)
• A = (ai,j) es estrictamente triangular superior si
∀1 ≤ i,j ≤ n (i ≥ j ⇒ ai,j = 0)
Matrices diagonales y escalares
• A = (ai,j) ∈ Rn×n es diagonal si
∀1 ≤ i,j ≤ n (i ≠ j ⇒ ai,j = 0)
• A = (ai,j) ∈ Rn×n es escalar si
∃λ∈|R ∀1 ≤ i,j ≤ n (ai,j = λ) ∧
∀1 ≤ i, j ≤ n (i ≠ j ⇒ ai,j = 0)
Observaciones
A escalar ⇒ A diagonal (⇐)
A diagonal ⇔ (A triangular inf ∧ A triangular sup)
Propiedades de det
• det At = det A
• Si T = (ai,j) es triangular, det T = Π ai,i
• det (λA) = λndet A ej. en 2×2:
• En general no vale “det (A+B) = det A + det B”
• det (AB) = det A . det B
• det I = 1
• A inversible ⇔ det A ≠ 0
• Si A es inversible, det (A-1) = (det A)-1
Sistemas triangulares
2 0 0 0
A = 0 1 0 b = 3
1 -3 3 -6
Matrices triangulares
¿Cuántos elementos hay desde la diagonal
hacia abajo en una matriz de n×n?
n
Sistemas triangulares
• A x = b, donde A es triangular (inferior o
superior)
• Se despejan los xi más fácilmente que si A
no fuese triangular.
• Puede tener interés contar el número de
operaciones elementales que se necesitan
para despejar los xi.
Idea del método LU
Si es posible descomponer a A como el
producto de dos matrices L y U, entonces el
S. E. L. A x = b se podrá transformar en dos
S. E. L.: L y = b, U x = y. Entonces, una vez
resuelto el primero, la solución (y) se coloca
como vector independiente en el segundo, y
luego se resuelve el segundo para hallar x.
Este x es la solución de A x = b.
Teorema (factorización LU)
Ejemplo con el método LU
x1 + 2x2 = 7
2x2 + 3x3 = 0
x1 + x3 = -2
El sistema lineal se podría dar en forma variada,
en ese caso primero hay que normalizarlo.
Ejemplo con el método LU (sigue)
1 2 0 7
A = 0 2 3 b = 0
1 0 1 -2
Ejemplo con el método LU (sigue)
1 2 0 1 0 0 u1 u2 u3
A = 0 2 3 = l1 1 0 0 u4 u5
1 0 1 l2 l3 1 0 0 u6
Siguiendo las reglas del producto de matrices, al vincular por
ejemplo la 1ª fila de L con la 3ª columna de U, se deduce que u3 = 0.
1) Despejar los li y los ui, en general algo sencillo.
2) Resolver Ly = b con sustitución hacia adelante.
3) Resolver Ux = y con sustitución hacia atrás.
Ejemplo 2
1 4 1
A = 2 8 1 b = …
0 -1 1
¡No se puede aplicar este método (a pesar de
que A es inversible)! Ocurre que no se
satisface la hipótesis del teorema anterior.
Ejemplo 3
0 0 3
A = 2 3 4 b = …
0 0 5
En este caso, A no es inversible. Tampoco
es LU factorizable. (Tampoco se satisface la
hipótesis del teorema anterior.)
Ejercicios
• Tomar sistemas lineales que sepan que
tienen solución única y vean si pueden
aplicar el método LU.
• Contar la cantidad de operaciones que se
hacen, para 2×2, 3×3,…, n×n.
• Probar el corolario siguiente…
Corolario
Si A es LU factorizable, entonces como
consecuencia de lo anterior se tiene que
det(A) = Π ui,i
Dem.
A = LU, con L triangular inferior con 1 en la
diagonal, y U triangular superior. Entonces
Norma matricial
Se llama así a cualquier función
|| . ||* : |Rm×n → |R que cumpla, para toda matriz A,
Ejemplo en R:
Otro ejemplo en R:
Por ejemplo:
Normas más usadas
norma 2, de Fröbenius,
pitagórica o euclidiana
i,
Otra norma:
Propiedades de algunas normas
La norma desigualdad
2 cumple de Hölder
No todas
las normas
cumplen
esto:
Como consecuencia A 2 A 2
de la desigualdad de
Hölder, en norma 2: A n A n
Teorema
D = (di,j) es diagonal ⇒ D-1 = (di,j-1) = (1/di,j)
(D es inversible ⇔ ∀1 ≤ i ≤ n di,i ≠ 0)
d1,1 0 0 … 0
0 d2,2 0 … 0
D = 0 0 d3,3 … 0
…
0 0 0 … 0 dn,n
d1,1-1 0 0 … 0
0 d2,2-10 … 0
⇒ D-1 = 0 0 d3,3-1 … 0
…
0 0 0 … 0 dn,n-1
Método de Jacobi
Es un método iterativo, a diferencia de los
anteriores que son exactos.
Se basa en hacer un esquema iterativo y mejorar
la aproximación que se tiene en cada paso…
Esa aproximación será un vector x(n), para n = 0,
1, 2, …
Siempre se puede expresar una matriz A como
A = L + D + U, donde L es estrictamente triang.
inf., U es estrictamente triang. sup., D es
diagonal.
Método de Jacobi
A=L+D+U
Ax = b
(L + D + U) x = b
(L + U) x + D x = b
D x = -(L + U) x + b
x = D-1 [-(L + U) x + b]
x = -D-1 (L + U) x + D-1b
Observar que invertir D es sencillo porque es una matriz
diagonal. De ahí la conveniencia del método…
Método de Jacobi
La última igualdad anterior equivale a despejar en
ecuación una variable distinta, y en función de ese
despeje, “actualizar” su valor una y otra vez (en los
pasos k = 1, 2, 3, …).
para i = 1, 2, …, n
Método de Jacobi
A partir de la identidad x = -D-1 (L+U) x + D-1b armamos
el esquema iterativo siguiente, que dada una
aproximación x(k) la “mejora” obteniendo x(k+1):
x(0) = un vector inicial cualquiera
x(k+1) = -D-1 (L+U) x(k) + D-1 b
Entonces, una pregunta muy importante es si esta
sucesión converge.
Es una sucesión de vectores. Que converja quiere decir
(para nosotros) que cada componente converja a un
número real en forma independientemente.
Convergencia de Jacobi
Sean H = -D-1 (L+U), v = -D-1 b
El esquema de iteración es
x(0)
x(k+1) = H x(k) + v
Supongamos que ||H|| < 1 (usemos ||.||2).
Convergencia de Jacobi
Llamemos x al vector de raíces exactas del S.E.L.
A x = b, tomemos H = -D-1(L+U), v = D-1 b
x=Hx+v⇒v=x–Hx
x(k+1) = H x(k) + v
e(k) = x(k) – x es el vector de error en el paso k
e(k+1) = x(k+1) – x = H x(k) + v – x = H x(k) + (I – H)x – x
= H x(k) + x – Hx – x = H x(k) – Hx = H (x(k) – x) = H e(k)
Convergencia de Jacobi
e(k+1) = H e(k) = H2 e(k-1) = H3 e(k-2) = … = Hk+1 e(0)
Luego, e(k+1) = Hk+1 e(0)
Entonces, si || H || < 1,
||e(k+1)|| = || Hk+1 e(0) || ≤ ||e(0)|| ||Hk+1||
≤ ||e(0)|| ||H||k+1 (por la desigualdad de Hölder)
= c ||H||k+1, donde c = ||e(0)|| es una constante.
Tenemos que c ||H||k+1 → 0 cuando k → ∞,
luego e(k+1) → 0.
Esto dice que x(k) – x → 0, luego x(k) → x,
es decir x(k) converge al vector de raíces.
Criterio de parada
Para k grande,
x(k+1) - x(k)
= (H x(k) + v) – (H x(k-1) + v) definición
= H (x(k) - x(k-1)) factorizando H
= H (H(x(k-1) - x(k-2))) definición
= H2 (x(k-1) - x(k-2)) asociatividad
= H2 (H(x(k-2) - x(k-3))) definición
= H3 (x(k-2) - x(k-3)) asociatividad
…
= Hk (x(1) - x(0))
Criterio de parada
Tenemos que para todo n ≥ 0
x(k+1) - x(k) = Hk (x(1) - x(0))
y notar que x(1) - x(0) es un vector constante.
Tomando norma 2,
|| x(k+1) - x(k) || = || Hk (x(1) - x(0)) ||
≤ ||Hk|| ||x(1) - x(0)|| por la desigualdad de Hölder
≤ ||H||k ||x(1) - x(0)|| por la desigualdad de Hölder
= ||x(1) - x(0)|| ||H||k → 0 cuando k → ∞ porque ||H|| < 1
Luego, || x(k+1) - x(k) || → 0 cuando k → ∞.
Observaciones
Esto es válido para cualquier esquema matricial
de la forma x(k+1) = H x(k) + v.
La convergencia no depende del vector inicial x(0).
La condición || H || < 1 se aplica a la matriz H (no a
la matriz A). A diferencia del criterio que
veremos después.
(Existen otros criterios mejores a partir de los
autovalores de H, pero más costosos de
verificar.)
Teorema
Si el esquema anterior converge,
entonces converge a una solución x del
S.E.L. A x = b.
Dem.
x(k+1) = -D-1 (L+U) x(k) + D-1b
tienden a
x* = -D-1 (L+U) x* + D-1b
Esta última igualdad equivale a A x* = b.
Teorema
Utilizaremos esto de la forma que sigue a continuación:
Diagonal estrictamente dominante
Además, si se cumple la condición ≤ para todas las filas
y la condición < al menos para una de ellas, entonces
también se asegura la convergencia.
Método iterativo general
Additive splitting de la matriz A:
x(k+1) = B x(k) + f
Se puede escribir en la forma
x(k+1) = P−1 N x(k) + P−1 b (para k ≥ 0)
Aquí, B = P−1 N , f = P−1 b .
donde A = P − N, que es la diferencia de dos
matrices P y N, donde P es inversible (llamada
preconditioning matrix).
En el método de Jacobi, se toma
P = D , N = -(L + U).
Método iterativo general
Para que esto sirva, P debe ser inversible y
calcular su inversa debe ser sencillo.
Alternativamente, se puede expresar así:
x(k+1) = x(k) + P-1 [b − A x(k)]
r(k) residuo en el paso k
Método de Jacobi sobrerrelajado
tomando un número ω tal que 0 < ω ≤ 1. Esto, en la forma matricial,
equivale a tomar la matriz B del método iterativo general como
B = ω H + (1- ω) I
Si ω = 1, en la fórmula superior el sumando derecho se hace 0, y se tiene
el método de Jacobi ya visto.
Complejidad de los algoritmos
iterativos
Eliminación de Gauss precisa O(n3) operaciones
(es decir una cantidad acotada por cte . n3).
Cada ciclo en el algoritmo iterativo lleva O(n2)
operaciones.
Para m iteraciones, son O(m n2) operaciones. Esto
lleva a comparar m con n.
Si A es esparsa con O(n) elementos no nulos, los
métodos iterativos precisan O(m n) operaciones.
Sistemas no lineales
DF(x) es la matriz jacobiana de f en x.
Es el método de Newton para sistemas de
ecuaciones no lineales.
Esta presentación es temporaria e incompleta y se complementa con la explicación oral y las
anotaciones hechas durante la clase. Se recomienda consultar al docente, la bibliografía, y el
blog o el campus virtual, tener en cuenta el contexto de la materia, los mensajes, las clases
anteriores y las próximas.