Cálculo Numérico-[Miguel Pasadas Fernández]

Tema 1 Preliminares Índice 1. Introducci Introducción oń 2. Teoremas eoremas básicos asicos 3. Errores Errores 3.1 Fuentes usuales de error 3.2 Representación o n de n´ umeros umeros 3.3 Tipos de errores 4. Propagaci´ Propagaci´ on on del error 5. Condicionam Condicionamien iento to y Estabilidad Estabilidad 6. Normas Normas de computaci´ computaci´ on on para el curso

1

1

Intr In trod oduc ucci ci´ ´ on on

Cálcu al culo lo nu num´ méric er ico o (An´ (A n´ alis al isis is Nu Num´ m´ eric er ico, o, M´ eto et odos do s Nu Num´ m´ eric er icos os): ): Nace con con el hom hombre, bre, pero pero se desa desarr rrol olla la junt juntoo con con las las m´ aquinas aquinas de • Nace cálculo. alculo.

• Estudia métodos etod os numéricos ericos para resolver problemas problema s matemáticos aticos de las ciencia cien cias, s, técnicas ecni cas e ingenie ing enierr´ıa.

Problemas que aborda:

• Finito-dimensionales – Resolución on de ecuaciones y sistemas de ecuaciones – Interpolación on y Aproximación on – Cálculo alculo de valores y vectores propios

• Infinito-dimensionales – Derivación on num´ num éric er ica. a. – Integración on num´ numéric er ica. a. – Resolución on numérica erica de E.D.O. y E.D.P.

Objetivo: Obtener la solución, on, exacta o con aceptable aproximación, on, con el menor esfuerzo/tiempo posible.

Ejemplo1 Ejem plo1:: M´ etodo eto do num´ erico eric o (m´ ( m´ etodo eto do iterativo) itera tivo) para calcular calc ular



1 5 xn+1 = xn + 2 xn 2



.

√ 5:

Aplicaci´ on on: (con 10 cifras) x0 = x1 = x2 = x3 = x4 =

2, 2.25, 2.236111111, 2.236067978, 2.236067977.

√ El valor exacto de 5 es 2.2360679774997896964.. Prerrequisitos de la materia : ´ ebraa line Alg lineal: al: • Algebr

espaci espacios os vecto vectoria riales les,, base, base, dimens dimensi´ i´ on, on, depende dependenci nciaa lineal, lineal, matrices, matrices, sistemas sistemas de ecuaciones ecuaciones lineales, lineales, valores alores propios. propios.

• Cálculo alculo infinitesimal e integral en una y dos variables.

3

2

Teoremas b´ asicos

Teorema de Bolzano Sea f c (a, b) tal que f (c) = 0.

∈ C [a, b] tal que f (a)f (b) < 0 entonces existe

∈

Teorema del valor intermedio Sea f C [a, b] y L un real comprendido entre f (a) y f (b). Entonces existe c [a, b] tal que f (c) = L .

∈

∈

Teorema de Rolle Sea f c [a, b] tal que f  (c) = 0.

1

∈ C [a, b] tal que f (a) = f (b). Entonces

∈

1

Teorema del valor medio Sea f f (b) f (a) que f  (c) = .

∈ C [a, b].

− b−a

existe

∈ [a, b] tal

Entonces existe c

Primer Teorema del valor medio del calculo integral Sea f entonces existe c (a, b) tal que

∈ C [a, b]

∈

f (c) =

b

1 b

−a



f (x)dx.

a

Segundo Teorema del valor medio del calculo integral Sean f, g C [a, b] tal que g (x) 0, para todo x [a, b]. Entonces existe c (a, b) tal que

≥

∈

∈

b



∈

b

f (x)g (x)dx = f (c)

a



g (x)dx.

a

Teorema de Taylor (formula de Taylor) Sean f C n+1 [a, b] y x 0 [a, b] fijo. Entonces para todo x [a, b] existe c = c (x) comprendido entre x0 y x, tal que f (x) = P n (x) + Rn (x), donde

∈

∈

n

 ( )=

P n x

k=0

f k (x0 ) (x k!

∈

f (n+1) Rn (x) = (x (n + 1)!

k

−x ) , 0

−x ) 0

∞

Suma de series geometricas Si r <

 1 entonces n=0

divergente en caso contrario. 4

cr n =

n+1

c

1

.

− r , siendo

3

Errores

3.1

Fuentes usuales de error

• Idealización:

rozamientos, vientos, atracciones, gravedad, relatividad, efectos ”despreciables”.

• Experimental-incertidumbre:

lectura aparatos, interferencias,

estimaciones estad´ısticas.

• Humano: equivocaciones aritméticas o de propagación. • Discretización: aproximación de un proceso matemático infinito por uno finito.

• Redondeo: 3.2

las máquinas tienen una precisi´ on limitada.

Representaci´ o n de n´ umeros

N´ umeros en notaci´ on binaria 156310 = 1 103 + 5 102 + 6 101 + 3 100 = 1 21 0 + 1 29 + 1 24 + 1 23 + 1 = 110000110112 .

× ×

× ×

× ×

×

×

×2

1

+1

×2

0

Algoritmo de conversi´ on decimal-binario Invértase el orden de la sucesión de restos de dividir por 2 reiteradamente.

N´ umeros no enteros 1563.2510 = 11000011011.012 0.710 = 0.101102 (periódico en binario). Otras bases usuales: 8, 16. 5

Notaci´ on cient´ıfica (Representación en punto flotante) z = σ (0.d1 d2 . . . dn )β e .

donde

• σ = ±1 (signo), • β ∈ N − {0, 1} (base), • e ∈ Ω ⊂ Z, n

• 0.d d . . . d 1 2

n

 = i=1

di β −i , di

∈ N, 0 ≤ d

i

< β (mantisa).

Normalizaci´ on: d1 > 0. β, n y Ω son caracter´ısticas de la m´ aquina.

on depende de n y de β . La precisi´ Ejemplo 2 Sea una máquina con β = 2, n = 4, e

∈ {−3, −2, −1, 0, 1, 2, 3}.

Los n´ umeros representables ser´ıan:

6

3.3

Tipos de errores

Definici´ on 1 Sea x∗ la representaci´ on de x

∈ R en una m´ aquina dada.

Se define el error absoluto de tal representaci´ on como ea = x

∗

| −x |

y el error relativo como ∗

er

|x − x | . = |x|

El error relativo es más intuitivo y da mejor idea de la precisión.

Definici´ on 2 Métodos de representaci´ on/operaci´ on:

• El error de truncatura aparece cuando se prescinde de las cifras de la mantisa a partir de una dada:

(0.d1 . . . dn dn+1 . . .)∗ = 0.d1 . . . dn .

• El

error de redondeo aparece cuando se corta la sucesi´ on de deci-

males de la mantisa mediante el redondeo de la ´ ultima cifra:

 0  ) =  0

(0.d1 . . . dn dn+1 . . .

.d1 . . . dn

si dn+1 <

.d1 . . . dn + β −n

si dn+1

∗

Propiedad 1 El error relativo de repretaci´ on est´ a acotado −n+1

• por truncatura: por β • por redondeo: por β 2

,

−n+1

.

7

β ,

2

≥ β 2 .

Ejemplo 3 Con β = 10 y n = 5, π

∗

π

∗

  • 10 = 0 31415 mediante truncatura, con error 0 0001,   • = 0 31416 mediante redondeo, con error 0 0005. 10

.

< .

.

< .

Definici´ on 3 Si d es el mayor entero para el cual −d+1

∗

|x − x | < β |x| 2

se dice que x∗ aproxima a x con d d´ıgitos significativos.

Ejemplos 1. β = 10, x = 3.141592, x∗ = 3 .14, −2

∗

|x − x | ≈ 0.000507 < 10 2 |x|

luego x∗ aproxima a x con 3 d´ıgitos significativos. 2. β = 10, x = 106 , x∗ = 999996, −5

∗

|x − x | ≈ 0.000004 < 10 2 |x|

luego x∗ aproxima a x con 6 d´ıgitos significativos. 3. β = 10, x = 0.000012, x∗ = 0 .000009, −0

∗

|x − x | ≈ 0.25 < 10 2 |x|

luego x∗ aproxima a x con 1 d´ıgito significativo. −d

| − x | < β 2

Definici´ on 4 Si d es el mayor entero para el cual x que x∗ aproxima a x con d decimales.

∗

se dice

¿Con cuántos decimales se realizan las aproximaciones de los ejemplos anteriores? 8

4

Propagaci´ on del error

Consideremos una má quina en la que β = 10, n = 6, mediante error de truncatura. Por tanto, el error relativo caracter´ıstico de representación es de 10−6+1 = 10−5 . Sean a = 1001 y b = 1000. Entonces a2

con error

2

− b ≈ 1002000 − 1000000 = 2000,

2001 2000 20001

−

−4

≈ 5 × 10

.

Pero a2

2

−b

= (a

− b)(a + b) = a + b = 2001,

con error 0. Luego procesos matemáticos equivalentes pueden no ser computacionalmente equivalentes. Si a = 1, b = 108 y c =

8

−10 , entonces

• a + (b + c) = 1, • (a + b) + c = 0. Por otro lado a + b = b , luego ¿es a = 0?.

Ejercicio n

 Se desea obtener

ai con n muy grande (miles de millones). ¿Cu´ a l es la

i=1

mejor forma de ordenar los calculos? 1. Si los ai son todos del mismo signo y magnitudes parecidas, 2. Si los ai son todos del mismo signo y magnitudes muy diversas, 3. Si los ai son de cualquier signo y magnitudes parecidas, 4. Si los ai son de cualquier signo y magnitudes muy diversas? 9

Advertencia Las operaciones de sumar (absoluta), multiplicar y dividir introducen error relativo de magnitud igual al de representacion, pero la resta (absoluta) puede introducir un error relativo gigantesco.

5

Condicionamiento y Estabilidad

Definici´ on 5 Un proceso est´ a bien condicionado si peque˜ nas variaciones en sus datos de entrada provocan peque˜ nas variaciones en la soluci´ on, y mal condicionado si las mismas condiciones provocan grandes variaciones en la soluci´ on. Un proceso de c´ alculo es estable si los errores de representaci´ on y redondeo introducidos tanto a la entrada como durante las operaciones intermedias no provocan perturbaci´ on importante en los resultados; e inestable en caso contrario. Sólo si se tiene un problema bien condicionado y se resuelve con un proceso estable se puede tener garant´ıa de precisión en el resultado. Por ejemplo, es fácil demostrar por inducci´ on que la sucesi´ on de valores 1 puede generarse indistintamente a partir de los siguientes algorit2 n≥0 mos:

  n

(I) s0 = 1, sn = 12 sn−1 , n (II) s0 = 1, s1 = 12 , sn =

≥ 1.

23 s 2 n−1

−

11 s , 2 n−2

n

≥ 2.

Sin embargo, con el segundo (operando con 6 cifras de precisión) el decimosexto término es s15 = 113, frente al valor 21 0.00031.

−

15



An´ alogamente, la sucesi´ on 31 n≥0 puede generarse a partir del algoritmo 1 10 s0 = 1, s1 = , sn = sn−1 sn−2 , 3 3 n 2, que tambi´ en es inestable.

  n

−

≥

10

6

Normas de computaci´ on para el curso

• Salvo indicación contraria, se debe trabajar con todas las cifras de la calculadora, incluso si se pide poca precisión.

• En particular, si se pide un resultado con 5 cifras decimales de precisión, NO se deben redondear los cálculos intermedios a 5 decimales.

• Para procesos iterativos de aproximaciones sucesivas, se detendrá el proceso cuando se repitan:

– Tantas cifras como la precisión requerida, si el proceso tiene asegurada una convergencia rápida (velocidad superior a la lineal, que ya se verá). – Tantas cifras como la precisión requerida MAS DOS, si el proceso converge lentamente (velocidad lineal).

• Ojo a las funciones trigonomtricas: ladora en modo radianes.

11

hay que poner siempre la calcu-

Tema 2 Resoluci´ on de Ecuaciones No Lineales Índice 1. Introducción 2. Método de Bisección 2.1 Algoritmo del Método de Bisecci´ on 2.2 Análisis de Método de Bisección 3. Método de Regula-Falsi 3.1 Algoritmo del Método de Regula-Falsi 3.2 Análisis de Método de Regula-Falsi 4. Método de la Secante 5. Método de Newton-Raphson 6. Métodos Iterativos 6.1 Algoritmo de los Métodos Iterativos 6.2 Interpretación Gráfica 6.3 Convergencia de los Métodos Iterativos 6.4 Convergencia Global del Método de Newton-Raphson 6.5 Aplicación del Teorema de convergencia global 7. Aceleración de la convergencia 7.1 Aceleración de Aitken 7.2 Aceleración de Steffensen 1

1

Introducci´ on

on amortiguada de una estructura Problema: Oscilaci´

Supongamos que la oscilaci´ on de una estructura, dotada de un sistema de amortiguaci´ on, ante un movimiento oscilatorio, viene dada por la función t

y (t) = 10 e cos2t. 2

¿En qué instante t la posición de la estructura es y(t) = 4? Se trata de resolver la ecuación t

10 e cos2t = 4. 2

de incógnita t. Este problema es imposible de resolver por medios anal´ıticos sencillos. Sea f : Ω ⊂

R

on en una variable → R. Consideraremos la ecuaci´ f (x) = 0.

umero s ∈ Ω se dice una soluci´ Definici´ on 1 El n´ on de la ecuaci´ on si se verifica que f (s) = 0, es decir, si s es una ra´ız de la funci´ on f .

2

2

M´ etodo de Bisecci´ on

2.1

Algoritmo del m´ etodo de Bisecci´ on

El método de Bisección para la resoluci´0on de la ecuación f (x) = 0 se basa en el Teorema de Bolzano que nos asegura la existencia de, al menos, una ra´ız de una función f (x) en un cierto intervalo [a, b], bajo ciertas condiciones.

on continua en [a, b] Teorema de Bolzano Sea f : [a, b] → R una funci´ tal que f (a)f (b) < 0 . Entonces existe c ∈ ( a, b) tal que f (c) = 0 .

Supongamos que f (x) es continua y cambia de signo en los extremos de [ a, b]. Basándonos en el anterior teorema, podemos aproximar una solució n de la ecuaci´ on f (x) = 0 dividiendo el intervalo inicial en dos subintervalos iguales y eligiendo aquel en el que f (x) cambia de signo. Después se repite el proceso hasta que se verifique algún criterio de parada.

Algoritmo del M´ etodo de Bisección 1. a0 = a , b0 = b 2. Para n = 0, 1, . . ., hacer: 1 (an + bn ) 2 ◦ Si f (an )f (mn ) < 0, tomar an+1 = an, bn+1 = mn ; en caso contrario, tomar an+1 = m n, bn+1 = b n .

◦ mn =

Ejemplo Resolver mediante al algoritmo de bisección la ecuación ex − x = 0 en [0, 1]. 3

2.2

An´ alisis del M´ etodo de Bisecci´ on

C´ alculo previo del n´ umero de interaciones Recordemos que

Se define el error absoluto de una aproximación s respecto del valor exacto s como e = | s − s|.

Para garantizar que el error del Método de Bisección sea menor o igual que un cierto valor de tolerancia ε se aplica el siguiente resultado:

Teorema1 (Error absoluto m´ aximo del M´ etodo de Bisecci´ on) Sea f : [a, b] → R una funci´ on continua en [a, b] tal que f (a)f (b) < 0 y f (s) = 0, para alg´ un s ∈ ( a, b). Sea {mn }n=0,1,... la sucesi´ on de aproximaciones de s obtenidas mediante el M´ etodo de Bisecci´ on y en = | s − mn |, para n = 0, 1, . . .. Entonces en ≤

b−a . n+1

2

Esquema de Demostraci´ on

4

en = |mn − s| ≤ m n − an = b n − mn =

1 (bn − an) = 2

1 (bn−1 − an−1 ) = . . . = 22 1 (b0 − a0 ). 2n+1 Luego en ≤

b.a

2n+1

. 

Por tanto, para garantizar que en < ε, se debe verificar que b−a ε − 1. log 2

log n ≥

on aproximada del problema de la oscilaci´ on amortiguada Ejemplo: Resoluci´ de una estructura

Se trata de resolver la ecuación t

f (t) = 10 e cos2t − 4 = 0 . 2

Supongamos que deseamos que en ≤ ε = 10−3 . Como f (0) = 6 > 0 y f (1) = −6.524 < 0 entonces podemos tomar [ a, b] = [0, 1].

5

El n´ umero de iteraciones que debemos realizar para asegurar la tolerancia de error considerada es: 1 10−3 − 1 ≈ 8 .966, log 2

log n ≥

es decir, n = 9 .

3 3.1

M´ etodo de Regula-Falsi Algoritmo del M´ etodo de Regula-Falsi

Se trata de realizar un refinamiento del Método de de Bisección, eligiendo la aproxiamción m a distancias de a y b proporcionales a f (a) y f (b). La ecuación de la recta que pasa por los puntos ( a, f (a)) y ( b, f (b)) es y − f (a) x−1 = , f (b) − f (a) b−a

de donde se tiene que el corte con el eje OX es, haciendo x = 0 y despejando y , el valor m =

af (b) − bf (a) . f (b) − f (a)

6

Se verifica que:

• mn converge más rápidamente a s que en el Método de Bisección. • Un extremo es fijo. • La amplitud de los intervalos no tiene a cero. • No admite acotación del error. Luego el algoritmo es el siguiente

Algoritmo del M´ etodo de Regula-Falsi 1. a0 = a , b0 = b 2. Para n = 0, 1, . . ., hacer:

◦ mn =

an f (bn ) − bn f (an ) f (bn ) − f (an )

◦ Si f (an )f (mn ) < 0, tomar an+1 = an, bn+1 = mn ; en caso contrario, tomar an+1 = m n, bn+1 = b n .

4

M´ etodo de la secante

Se trata de un m´ etodo iterativo en el que, en cada paso, se calcula una aproximaci´ on de la solución en lugar de un intervalo que la contiene. 7

Se parte de x 0 = a y x 1 = b y se calcula, iterativamente para cada n ≥ 1, la intersección de la secante que une los puntos ( xn−1 , f (xn−1 ) y (xn , f (xn )) con el eje de abscisa, obteniéndose la abscisa xn+1 =

xn−1 f (xn ) − xn f (xn−1 ) . f (xn ) − f (xn−1 )

Por tanto el algoritmo es el siguiente:

Algoritmo del M´ etodo de la Secante 1. x0 = a , x1 = b 2. Para n = 1, 2 . . ., hacer xn+1 =

xn−1 f (xn ) − xn f (xn−1 ) f (xn ) − f (xn−1 )

Ejemplo Calculemos mediante 5 pasos del método de la secante una aproximacin de la solución del problema de la oscilación amortiguada de una estructura. Se trataba de calcular la solución de la ecuación t

f (t) = 10 e cos2t − 4 = 0 . 2

8

Tomando x 0 = 0 y x1 = 1, se obtiene la siguiente sucesión de aproximaciones: x2 = x3 = x4 = x5 =

0.479078 0.517905 0.513640 0.513652,

con lo cual podemos afirmar que una aproximación con cuatro decimales exactas de la solució n es 0.5126.

5

M´ etodo de Newton-Raphson

Se trata de llevar el l´ımite el método de la secante y, por tanto, en cada iteración n, considerar la recta tangente a f (x) en (xn, f (xn )) y tomar como siguiente aproximación xn+1 la intersección de dicha tangente con el eje de abscisas. Por tanto, teniendo en cuenta que la ecuaci´ on de la recta tangete a la gráfica de f (x) en el punto ( xn , f (xn )) es y − f (xn ) = f  (xn )(x − xn ),

se tiene que

xn+1 = x n −

f (xn ) . f  (xn )

Luego el algoritmo del Método de Newton-Raphson es

Algoritmo del M´ etodo de Newton-Raphson 1. Dado x0 2. Para n = 1, 2 . . ., hacer xn+1 = x n −

9

f (xn ) f  (xn )

Observaciones sobre el M´ etodo de Newton-Raphson • Es el más rápido =0 • Requiere que f  (s) 

• Elegir x0 puede ser delicado • La acotación del error es complikcada • El criterio de parada más usual es la repetición de cifras. Ejemplos A continuación presentamos la solución aproximada de algunas ecuaciones con el método:

¿Qué ha ocurrido en el tercer caso? Obviamente se ha presentado un caso de divergencia: 10

También se podr´ıa haber presentado un caso de oscilaci´ on como el siguiendo gráfico indica:

11

6

M´ etodos Iterativos

Se trata de transformar la ecuación f (x) = 0 (cálculo de una ra´ız de la función f (x)) en una ecuación del tipo x = g (x) (cálculo de un punto fijo de la función g (x)) de forma que sean equivalentes, es decir, tengan la misma solución.

Ejemplo (Método de Newton-Raphson) Si f (x) es una función de clase C 1 y s una ra´ız (f (s) = 0) tal que f  (s)  = 0, entonces, resolver f (x) = 0 es un problema equivalente a calcular un punto f (x) fijo de la función g (x) = x −  ya que f (s) = 0 si y sólo si g (s) = s , como f (x) se puede comprobar fácilmente.

6.1

Algoritmo de los m´ etodos iterativos

Los métodos iterativos se basan en el cálculo de un punto fijo para una cierta función g(x). El siguiente resultado determina condiciones suficientes para la existencia de un punto fijo para g (x).

Teorema del punto fijo Sea g : [a, b] →

R una

funci´ on derivable verifi-

cando:

• g ([a, b]) ⊆ (a, b), • maxx∈[a,b] |g  (x)| < 1 . Entonces existe un ´ unico s ∈ [a, b] tal que g (s) = s y, adem´ as, para todo on {xn } generada por la iteraci´ on xn+1 = g (xn ) conx0 ∈ [a, b], la sucesi´ verge a s.

Demostraci´ on Sea h(x) = g(x) − x, para todo x ∈ [a, b]. Entonces h(x) es continua y verifica que h(a) > 0 y h(b) < 0, por lo que se verifican las condiciones del Teorema de Bolzano. En consecuencia, existe un s ∈ (a, b) tal que h(s) = 0, es decir g(s) = s . Para demostrar la unicidad, supongamos que existe dos valores s, t ∈ (a, b) tales que g (s) = s y g (t) = t, entonces, por el Teorema del Valor Medio, existe c ∈ ( s, t) tal que g(t) − g(s) = g  (c)(t − s), es decir, g (c) = 1, lo que contradice la segunda condición. 12

Sea ahora {xn } la sucesió n generada a partir de x0 ∈ [a, b] mediante la iteración xn+1 = g (xn , para m ≥ 0, y sea L = maxx∈[a,b] |g  (x)| < 1. Se verifica entonces que en = |xn − s| = |g (xn−1 ) − g(s)| = | g (x)|en−1 ≤ Le n−1 ≤ L 2 en−2 ≤ ·· · ≤ L ne0 . 

Algoritmo de un M´ etodo Iterativo 1. Dado x0 ∈ [ a, b] 2. Para n = 1, 2 . . ., hacer xn+1 = g (xn)

6.2

Interpretaci´ on gr´ afica de los m´ etodos iterativos

Seún sea g (x) y se elija x0 , los métodos iterativos pueden ser convergentes o divergentes y, en ambos casos pueden variar en forma espiral o en escalera, como se indican en los siguientes gráficos.

13

6.3

Convergencia de los m´ etodos iterativos

on { xn } hacia Definici´ on 2 Se define el orden de convergencia de una sucesi´ un valor s como aquel n´ umero real p ≥ 1 tal que em+1 = λ  = 0 , ∞. p n→+∞ en

lim

Se deduce que si {xn } → s con orden de convergenia p ≤ 1 entonces dn+1 ≈ λe pn ,

n → + ∞.

Algunos o´rdenes de convergencia de los métodos estudiados son los siguientes:

14

on de clase Teorema de Convergencia Local Sea g : [a, b] → R una funci´ 1  C en un entorno del punto fijo s. Si |g (s)| < 1, entonces existe δ > 0 tal que en el intervalo [s − δ, s + δ ] se dan las condiciones del teorema del punto fijo.

Demostraci´ on Basta tomar L tal que |g (s)| < L < 1 y δ > 0 tal que  |g (x)| ≤ L , para todo x ∈ [ s − δ, s + δ ].

on Teorema de Convergencia Cuadr´ atica Sea g : [a, b] → R una funci´ 2  de clase C en un entorno del punto fijo s. Si |g (s) = 0, entonces la convergencia {xn } → s es al menos cuadr´ atica.

Demostraci´ on Teniendo en cuenta que en+1 = | xn+1 − s| = | g(xn ) − g (s)|

1 = | g (s)(xn − s) + g (cn )(xn − s)2 | 2 1 0 + g  (cn)e2n; 2 concluimos que 1 1 en+1 = lim |g (cn)| = |g  (s)|  = ∞ . n→∞ en n→+∞ 2 2 lim



15

6.4

Convergencia global del m´ etodo de Newton-Raphson

Teorema de Convergencia LOCAL para el M´ etodo de NewtonRaphson Sea f : [a, b] → R una funci´ on de clase C 2 en un entorno de s tal que = 0. Entonces existe δ > 0 tal que, para todo x0 ∈ f (s) = 0 y f  (s)  [s − δ, s + δ ], la sucesi´ on del m´ etodo de Newton-Raphson converge a s. 3 Adem´ as, si f es de clase C [a, b], la convergencia es al menos cuadr´ atica.

Teorema de Convergencia GLOBAL para el M´ etodo de NewtonRaphson Sea f ∈ C 2 [a, b] verificando 1. f (a)f (b) < 0 , 2. f  (x)  = 0, para x ∈ [ a, b], 3. f  (x)f  (y ) ≥ 0 , para todo x, y ∈ [ a, b],



|f (a)| |f (b)| 4. max , |f  (a)| |f  (b)|



≤ b − a.

Entonces existe un ´ unico s ∈ [ a, b] tal que f (s) = 0 y, para todo x0 ∈ [ a, b], la sucesi´ on del m´ etodo de Newton-Raphson converge a s. Adem´ as, si f ∈ 3 atica C [a, b] la convergencia es al menos cuadr´

Demostraci´ on Por [1.], aplicando el Teoreme de Bolzeno, existe al menos una ra´ız s para f (x); y por [2.], aplicando el Teorema de Rolle, la ra´ız s es u ńica. Asimismo, por [2.] y [3.], f  y f  conserva su signo en [a,b]. Supongamos que f  (x) > 0 y f  (x) ≥ 0, para todo x ∈ [ a, b] (en el resto de los casos se razona análogamente). 16

f (x) f (x)f  (x)  Tomando g (x) = x −  se tiene que g (x) = , y por ser f (x) (f  (x))2 f (x) creciente, se verifica que g(x) es decreciente en [ a, s] y creciente en [ s, a].

Entonces

• si x ∈ [ a, s], s = g (s) ≤ g (x) ≤ g (a) = b por [4.], luego g (x) ∈ [ s, b]; • si x ∈ [s, b], s = g (s) ≤ g(x) ≤ g (b) = b po definició n de g , luego g (x) ∈ [ s, b]. Por tanto, para todo x0 ∈ [ a, b], la sucesión del método de Newton-Raphson {xn } ⊂ [ s, b]. Además xn+‘1 = g (xn) = xn −

f (xn ) on es estric< xn , luego la sucexsi´ f  (xn )

tamente creciente y acotada. En consecuencia tiene l´ımite, es decir, existe s ∈ [ s, b] tal que lim = s  . n→+∞

Finalmente, tomando l´ımite en xn+1 = g (xn), se llega a que s = g (s ) y, por la unicidad del punto fijo s, concluimos que s = s  , es decir, lim xn = s .  n→+∞

6.5

Aplicaci´ on del teorema de convergencia global

Problema: Calcular la ra´ız positiva de c Se trata de hallar la ra´ız positiva de la ecuaci´ on f (x) = x 2 − c = 0. Veamos si se verifican las cuatro condiciones del teorema de convergencia tomando 0 < a < b tales que a2 < c < b2 1. f (a)f (b) = ( a2 − c)(b2 − c) < 0. 2. f  (x) = 2x > 0, para todo x ∈ [ a, b]. 3. f  (x) = 2  = 0, para todo x ∈ [ a, b]. (b + a)(b − a) b2 − c b2 − a2 |f (b)| 4.  = = ≤ ≤ b − a, luego se verifica 2b 2b 2b |f (b)| esta hipótesis. En conclusión, el método de Newton-Raphson aplicado a f (x) = x2 − c en un intervalo 0 < a < b, con a2 < c < b2 , converge a la ra´ız positiva de c. 17

7 7.1

Aceleraci´ on de la convergencia M´ etodo de Aitken

Teorema de aceleraci´ on de Aitken Sea {xn } → s al menos linealmente. Se define

(xn+1 − xn )2 xn = x n − . xn+2 − 2xn+1 + xn



Entonces {xn } → s m´ as r´ apidamente en el sentido





en = 0, n→+∞ en

lim

siendo en = xn − s y en = x n − s.

 

Demostraci´ on Supongamos que Entonces

en+1 = kn con lim kn = λ y |λ| < 1. n→+∞ en

(xn+1 − xn)2 en = xn − s = x n − − s = xn+2 − 2xn+1 + xn ((xn+1 − s) − (xn − s))2 (en+1 − en )2 = e n − = xn − s − (xn+2 − s) − 2(xn+1 − s) + ( xn − s) en+2 − 2en+1 + en e2n (kn − 1)2 kn (kn+1 − kn ) = e n en − , en (kn+1 kn − 2kn + 1) kn+1 kn − 2kn + 1

  luego

7.2



en → 0. en



M´ etodo de Steffensen

Dado x0 , consideremos el método iterativo xn+1 = g (xn ). Sean x1 y x2 las dos primeras aproximaciones.

18

Definimos

(x1 − x0 )2 x0 = x 0 − , x2 − 2x1 + x0 

Y, en general, para cada n >= 0, se define (n+1)

= x n x0 (n+1) (n+1) = g (x0 ), x1 (n+1) (n+1) = g (x1 ), x2 (n+1) xn+1 = x 0 

−

(n+1)

(x1 (n+1)

x2

(n+1) 2

− x0

(n+1)

− 2x1

)

(n+1)

+ x0

.

Entonces { xn} → s más rápidamente que los anteriores métodos. El método etodo de Steffensen. as´ı definido se denomina m´

Ejemplo La siguiente tabla muestra las sucesiones de aproximaciones obtenidas mediante el método de iteración funcional xn+1 = g (xn ), el método de aceleración de Aitken y el de superaceleración de Steffensen. En concreto se considera g(x) = e−x y x0 = 0.5 (que tiene convergencia lineal):

19

LU

n n a11 x1 + a12 x2 + · · · + a1n xn = b 1 , a21 x1 + a22 x2 + · · · + a2n xn = b 2 , ... an1 x1 + an2 x2 + · · · + ann xn = b n ,

A

   =  

a12 · · · a22 · · ·

a1n a2n

an1 an2 · · ·

ann

a11 a21

b

   =  

b1 b2 bn

  

  

  

Ax = b .

{x(m) } n

n

x = A −1 b

x1 , . . . , xn xi =

det Ai , det A

i = 1, . . . , n ,

A

Ai

i

n!n − 1

(n + 1) n! − 1 n

106

n

106

n

5 10 20

≈ 3600 ≈ 4 × 108 ≈ 1 ,02 × 1021

3,6 6 32,4

39

An×n = (aij ) AT = A −1 AT = A

i  = j

aij = 0

|i − j | > 1 i>j

aij = 0 aij = 0

i > j + 1

aij = 0

xt Ax(≥) > 0 ∀ x  =0

| j  =i

aij |(<)|aii |, ∀ i

Ux = b

det U  =

0

Un×n

n

bi xi =

 −

j =i+1

uii

uij x j ,

i = n, n − 1, . . . , 1.

3x1 + 2 x2 − 2x3 + 4 x4 = − 5, 3x2 − 5x3 − 3x4 = 0, 4x3 + x4 = − 3, 2x4 = 6.

6 = 3, 2 3 −3 − x4 = − , x3 = 4 2 5x3 + 3 x4 1 = , x2 = 3 2 −5 − 2x2 + 2 x3 − 4x4 x1 = = − 7. 3 x4 =

Ax = b Ux = c

a21 , a31 , . . . , an1 (2)

(2)

(3)

(3)

a32 , . . . , an2 a43 , . . . , an3

(n−1)

an,n−1 (k)

aik (k)

akk

k i

k

−

(k)

aik

(k)

akk

n

  

2 3 −2 1 0 1 3 2 −1 −2 2 −2 0 1 2 1 1 0 1 −1

2  0  0 0

  

3 −2 1 0 3 3 3 − −2 2 2 2 0 2 −5 5 1 0 1 2 2

   

2  0  0  0 2  0  0  0

3 −2 1 0 3 3 3 − −2 2 2 14 0 12 −5 − 3 13 0 −5 3 3 3 −2 1 3 3 3 − 2 2 0 12 −5 0

0

11 2

   

 0  −2  14   − 3   43  − 18

 14   − 33   − 334   23   3386  33

4n3 + 3n2 − 7n = O c(Gaussn ) = 6

106

n3

2  3

.

106

n

5 10 20 100 100 0

90 705 5510 6 71 55 0 6 67

90 0,7 5,5 0,67 11

Ax = b Dx = c

k k = 1, . . . , n n

D (k )

(k)

(k)

(k)

a1k , . . . , ak−1,k , ak+1,k , . . . , ank n3 + n 2 − 2 n

n

det 0 det A  =0

= 

 0 0001 1 1  1 1 2  1 00010  ,

, ... 0,99990...

 0 0001 1 0 −10000 0 ,

1 −10000



1

|aik | = máx |ark |

i > k

k≤r ≤n

k

 0 0001



1 1 2 0,0001 1 1

 1

1 1 1 1 2

,



∼

exacta

−→

 1 00010 

, 1 1 0,0001 0,0001 0,0002

(2)

...

0,99990

1 2 0 1 1

 0 0001 (1)

,

 1 p.p.

−→

1

 0 p.p.

−→

1

(n)

.

máx |aij | = máx |a2 j | = · · · = máx |anj | .

.

(i, j ) ( )

k ax |a(rsk) | , |aij | = m´ ≤ ≤ k

r

n

k ≤s≤n

i

 10  7 8

k

7 8 7 5 6 5 6 10 9 7 5 9 10

 10  7 8

7 8 7 5 6 5 6 10 9 7 5 9 10

j

32 23 33 31

32,1 22,9 33,1 30,9

  

  

sol.exacta

−→

k

1  1  1

,

1

sol.exacta

−→

 92   −12 6   45  , , , −1,1

.

LU LU A = LU

L

U Ax =

b Ly = b , Ux = y

2n2

 1  −2  3

0 2 0 −5 2

 100  −2 2 0  302 −5 2 1  2 3 −3  0 3 2 0 0 1 0 0

0 0 2 1

0 0 0 1

 2   0  0

3 −3 0 3 2 2 0 1 1 0 0 0 −2

  −7   =⇒ =  1    2 0   1 −7  −1  1  = = ⇒  2 2

0 −7 0 16 0 −17 1 39

0 2 1 0 −2

  

0

.

39

.

y

x

0

x

 −7   16  =   −17 

.

LU

LU k = 1, . . . , n k−1

kk ukk = a kk

 −  −

kr urk ;

r =1 k−1

aik

ir urk

r =1

ik =

ukk

,

i = k + 1 , . . . , n;

,

j = k + 1 , . . . , n .

k −1

akj

 −

kr urj

r =1

ukj =

kk

4n3 + 2n 2n3 = O( ) 6 3

n2 + n2

1

L

1

U

LU A = LU

aij = 0

|i − j | > 1

k = 1, . . . , n kk ukk = a kk − k,k −1 uk−1,k ; ak+1,k ; ukk ak,k +1 = . kk

k+1,k = uk,k +1

4n − 3

2(3n − 2)

L

U

LU U = L T

A A = LL t .

k = 1, . . . , n

   = −  − k −1

kk

akk

2k,r ;

r =1 k−1

ai,k

i,k =

O(

n3

3

)

ir kr

r =1

kk

,

i = k + 1 , . . . , n .

A=

A D−L−U U i>j A D

ij = −aij

D = diag(A) L

uij = − aij

Ax = b

i
(D−L−U)x = b

Dx = (L+U)x+b

x = D −1 (L + U)x + D−1 b.

x(m+1) = B J x(m) + cJ ,

  − )=  −

BJ = D −1 (L + U

0 a21 a22

an1 ann

cJ = D −1 b

bi (m+1)

xi

=

 −

a12 − a11

   =   

···

0 ···

−

an2 ann

b1 a11 b2 a22 bn ann

    

a1n − a a211n − a22

···

aii

(m)

,

,

.

aij x j

j  =i

    0

i = 1, . . . , n .

m

(D−L−U)x = b

Ax = b

(D−L)x = Ux +b

x = (D − L)−1 Ux + (D − L)−1 b.

(D − L)x(m+1) = Ux (m) + b,

m bi − (m+1)

xi

=



(m+1)

aij x j

−

j
 j>i

aii

ω ∈

(m)

aij x j

,

i = 1, . . . , n .

R

A =

1 ω

D−

1 − ω ω

D − L − U.

(D − ω L)x = ((1 − ω )D + ω U)x + ω b.

(D − ωL)x(m+1 = ((1 − ω )D + ω U)x(m) + ω b.

m bi − (m+1)

xi

= ω



(m+1)

aij x j

j
−

 j>i

aii

(m)

aij x j

(m) + (1 − ω )xi ,

i = 1, . . . , n .

A

ω ∈ [0 , 2]

DEPARTAMENTO DE MATEMÁTICA APLICADA www.ugr.es/local/mateapli

INTERPOLACIÓN

2006-2007

José Martínez Aroza

Introducción

Interpolar (D.R.A.E.): Averiguar el valor aproximado de una magnitud en un intervalo cuando se conocen algunos de los valores que toma a uno y otro lado de dicho intervalo, y no se conoce la ley de variación de la magnitud.

Interpolación

2/90

Introducción ventas

800

600

400

200

meses 2

4

6

8

Interpolación

3/90

Introducción ventas

800

600

400

200

meses 2

Interpolación

4

6

8

4/90

Interpolación Lagrangiana

Problema de Interpolación Lagrangiana (versión previa): “Dados los nodos { 0 , x1 ,…, x n } y sus valores respectivos { y 0 , y1 ,…, y n }, obténgase una función p ( x ) “sencilla” que verifique p ( xi ) = yi , i = 0,1,…, n .”

10

5

1

-

2

3

4

5

Interpolación

5/90

Interpolación Lagrangiana p : función interpoladora o interpolante “Sencilla” = hay que fijar un espacio de funciones . Problema bien planteado ⇔ dim F = n + 1 = nº de datos.

Problema de Interpolación Lagrangiana (P.I.L.): “Dados los nodos { 0 , x1 ,…, x n } y sus valores respectivos { y 0 , y1 ,…, y n }, y dado (o fijado) un espacio de funciones con dim F = n + 1, obténgase una función p ∈ F que verifique p ( xi ) = yi , i = 0,1,…, n .”

Problema de Interpolación Polinomial Lagrangiana (P.I.P.L.): , x1 ,…, x n } y sus valores respectivos { y 0 , y1 ,…, y n }, obténgase un polinomio p ∈ Pn tal que ( i ) = y i , i = 0,1, … , n .” “Dados los nodos {

Interpolación

0

6/90


Resolución del P.I.L.: n

Se escoge una base {ϕ 0 , ϕ 1 ,…, ϕ n } de F y se busca p

= ∑ a jϕ j que cumpla las j = 0

condiciones de interpolación

p ( x0 ) = y0 ⇔ a0ϕ 0 ( x0 ) + a1ϕ 1 ( x0 ) +  + anϕ n ( x0 ) = y0 ⎫

⎪

p( x1 ) = y1 ⇔ a0ϕ 0 ( x1 ) + a1ϕ 1 ( x1 ) +  + anϕ n ( x1 ) = y1 ⎬ p( xn ) = y0 ⇔ a0ϕ 0 ( xn ) + a1ϕ 1 ( xn ) +  + anϕ n ( xn ) = yn ⎪ ⎭ (Sistema lineal ( n + 1) × ( n + 1) con incógnitas a j .) El P.I.L. es unisolvente (solución única)

⇔ el sistema es C.D.

Solución = coeficientes de la función interpoladora.

Interpolación

7/90


Ejemplo:

x

−1

y

0

0

1

3

1 −2

4

Polinomio buscado: p ( x ) = a 0 Sistema:

a0

−

a1

, F = P3

= 1, x, x 2 , x 3

+ a1 x + a 2 x 2 + a 3 x 3 . +

a2

−

a3

a0 a0 a0

+ a1 + 3a1

.

+ a2 + 9a 2

+ a3 + 27a 3

= 0⎫ = 1⎪ ⎪ ⎬ = − 2⎪ = 4⎪ ⎭

⇒ p ( x ) = x 3 − 2 x 2 − 2 x + 1.

Interpolación

8/90


5

-

2

-

1

1

-

-

2

4

3

5

10

Interpolación

9/90


Unisolvencia del P.I.P.L.: Usando la base {1, x,…, x n }, el determinante de la matriz de coeficientes es el det. de Vandermonde

⎛ 1 x 0 x 02  x 0n ⎞ ⎜ ⎟ 2 n ⎜1 x1 x1  x1 ⎟ = ∏ ( x i − x j ) ; det ⎜ ⎟      i > j ⎜⎜ ⎟ 2 n ⎟  1 x x x ⎝ n n n det ≠ 0 ⇔ las abscisas x i no se repiten (no hay dos iguales). Demostración: restando a cada columna la anterior × x0 , 0 ⎛ 1 ⎜1 x − x 1 0 = det ⎜   ⎜ ⎝ 1 xn − x0 Interpolación

 0 ⎞ ⎛ 1 x1 n −1  x1 ( x1 − x0 ) ⎟ ⎜  = det ⎟ ⎜⎜   ⎟ n −1 ⎝ 1 xn  xn ( xn − x0 )

n −1

 x1 ⎞ n ⎟   ( xi − x0 ) . n −1 ⎟  xn ⎟ i =1

∏

⎠

10/90

Interpolación Lagrangiana Para no tener que resolver s.e.l. se emplean fórmulas de interpolación, que dan directamente el polinomio de interpolación.

Fórmulas clásicas: la Fórmula de Lagrange y la Fórmula de Newton. Diferentes expresiones, pero el mismo objetivo: el polinomio de interpolación.

Interpolación

11/90

Fórmula de Lagrange

FÓRMULA DE LAGRANGE

Interpolación

12/90

Fórmula de Lagrange Dados los nodos { x 0 , x1 ,…, x n }, los polinomios fundamentales de Lagrange

 0 ,  1 ,…,  n ∈

P n son aquellos que verifican

 0 ( x0 ) = 1,  0 ( x1 ) = 0,  1 ( x0 ) = 0,  1 ( x1 ) = 1,    n ( x0 ) = 0,  n ( x1 ) = 0,

…  0 ( xn ) = 0⎫ …  1 ( xn ) = 0 ⎪ ⎧1 si i = j = = ( x ) δ , es decir,  ⎬ ⎨0 si i ≠ j j i ij   ⎩ ⎪ …  n ( xn ) = 1⎭

(delta de Kronecker). Una vez obtenidos, la solución del P.I.P.L. es la Fórmula de Lagrange

p ( x ) = y 0  0 ( x ) + y1  1 ( x ) +  + y n  n ( x ) =

n

∑ y  j

j

( x) .

j = 0

En efecto, p ( x i )

n

n

j = 0

j = 0

= ∑ y j  j ( x i ) = ∑ y j δ ij = y i .

Interpolación

13/90


Construcción de los polinomios fundamentales de Lagrange:  0 se anula en x1 , … , x n ⇒ es divisible por ( x − x1 ), ( x − x 2 ), … , ( x − x n ) ⇒  0 = K ( x − x1 )( x − x 2 )  ( x − x n ) , con K = cte. Se escoge K para que  0 ( x 0 )

por tanto  0 ( x )

=

x − x1 x − x2

⋅

x0 − x1 x0 − x2 n

En general  j ( x )

= 1 ⇒ K =

=∏ i =0

x − xi x j − xi



x − xn

x0 − xn

1 ( x 0 − x1 )( x 0 − x 2 )  ( x 0 − x n ) n

=∏ i =1

x − xi x0 − xi

.

(hay un error).

Ejercicio: demostrar que { 0 , 1 ,…,  n} es una base de Interpolación

,

n

.

14/90


Ejemplo:

x

−1

y

0 0 =

0

1

3

1 −2

4

x − 0

⋅

x − 1

, unisolvente en P3 .

⋅

x − 3

−1− 0 −1−1 −1− 3 x + 1 x − 1 x − 3 1 = ⋅ ⋅ 0 +1 0 −1 0 − 3 x + 1 x − 0 x − 3 2 = ⋅ ⋅ 1+1 1− 0 1− 3 x + 1 x − 0 x − 1 3 = ⋅ ⋅ 3 +1 3 + 0 3 −1 p = 0 0 + 1 1 − 2 2 + 4 3

=−

1

( x 3 − 4 x 2 + 3 x )

8 1 3 = ( x − 2 x 2 − x + 4) 3 1 3 = − ( x − 2 x 2 − 3 x ) 4 1 = ( x 3 − x ) 24

= x 3 − 2 x 2 − 2 x + 1

Interpolación

15/90


3

2

1

-

Interpolación

2

-

1

1

-

1

-

2

2

3

4

16/90


Ventaja de la fórmula de Lagrange: Los polinomios fundamentales sólo dependen de los nodos { x0 , 1 ,…, x n }; varios P.I.P.L. con el mismo conjunto de abscisas comparten los mismos polinomios fundamentales. 600

Crédito

500 400 300

Contado

200 100

Gastos

1

2

3

4

Inconveniente de la fórmula de Lagrange: Si se añade nueva información (un nuevo punto), hay que reconstruir casi todo.

Interpolación

17/90

Fórmula de Newton

FÓRMULA DE NEWTON

Interpolación

18/90

Fórmula de Newton

Objetivo: construir el pol. de interp. ( x ) gradualmente, partiendo de un solo dato ( x 0 , y 0 ) e ir añadiendo datos progresivamente. p 0 ∈

P 0 : polin. que interpola en

{ x 0 },

p1 ∈


{

p 2 ∈


{ x 0 , x1 , x 2 },

0

, x 1 },



n

= p∈

P n : polin. que interpola en

Se trata de obtener

{ x 0 , x1 , … ,

k conocido

n

}.

pk −1 , k = 1,…, n .

Interpolación

19/90

Fórmula de Newton 4

3

2

1

1

-

2

3

4

5

1

Interpolación

20/90


3

2

1

1

-

2

3

4

5

1

Interpolación

21/90


3

2

1

1

-

2

3

4

5

1

Interpolación

22/90


3

2

1

1

-

2

3

4

5

1

Interpolación

23/90


3

2

1

1

-

2

3

4

5

1

Interpolación

24/90

Fórmula de Newton

Construcción progresiva • p 0 ( x ) = y 0 ∈

p 0 ( x 0 ) = y 0 ). • p1 ( x ) ∈ P1 interpola en { x0 , x1}; q1 = p1 − p0 ∈ P 1 se anula en 0 ⇒ q1 ( x ) = 1 ( − x 0 ) , con 1 adecuado para que y1 − p 0 ( 1 ) p1 ( x1 ) = p 0 ( x1 ) + q1 ( x1 ) = y1 ⇒ A1 = . x1 − x 0 • p 2 = p1 + q 2 ∈ P2 interpola en x0 , x1 , x2 ; q 2 ∈ P2 se anula en 0 y x1 ⇒ q 2 ( x ) = 2 ( x − x 0 )( x − x1 ) con 2 adecuado para que P 0 interpola en

0 (cumple

y 2 − p1 ( x 2 )

p 2 ( x 2 ) = p1 ( x 2 ) + q 2 ( x 2 ) = y 2 ⇒ A2 =

= p k −1 + q k ∈

• en general p k

⇒ q k ( ) =

k

( x 2 − x 0 )( x 2 − x1 ) P k , q k se anula en x 0 , … , x k −1

( x − x 0 )  ( x − x k −1 ) con Ak =

.

y k − p k −1 ( x k )

( x k − x 0 )  ( x k − x k −1 )

Interpolación

.

25/90

Fórmula de Newton Conviniendo

0

= y 0 se tiene

p ( x ) = p n ( x ) =

0

+

1

( x − x 0 ) +

2

( x − x 0 )( x − x1 )

+  + An ( x − x 0 )  ( x − x n −1 ) es decir

p( x ) =

k −1

n

∑ A ∏ ( x − x ) , k

k = 0

i

i =0

que es la fórmula de Newton (preliminar).

Interpolación

26/90

Fórmula de Newton

Ejemplo: A1 =

−1

x y

0

y1 − p0 ( x1 ) x1 − x0

0

1

3

1 −2

=

1− 0 0 +1

y2 − p1 ( x2 )

4

;

0

= y0 = 0 ⇒ p0 ( x ) = 0;

= 1 ⇒ p1 ( x ) = 0 + 1( x + 1) = x + 1;

−2−2 = −2 ⇒ ( x2 − x0 )( x2 − x1 ) (1 + 1)(1 − 0) p2 ( x ) = p1 ( x ) − 2( x + 1)( x − 0) = ( x + 1) − 2( x + 1) x;

A2 =

A3 =

y3 −

2

=

( x3 )

( x3 − x0 )( x3 − x1 )( x3 − x2 )

=

4 + 20 4 ⋅ 3⋅ 2

=1⇒

p3 ( x ) = p2 ( x ) + 1( x + 1) x ( x − 1)

= x 3 − 2

2

− 2 x + 1.

Interpolación

27/90

Diferencias divididas

Notación:

k

= f [ x 0 , … , x k ] diferencia dividida de f en x 0 , … , x k .

Fórmula de Newton: (definitiva)

p( x ) =

∑ f [ x k = 0

Interpolación

k −1

n

0

∏ ( x − x ) .

, … , x k ]

i

i =0

28/90


Propiedad de las diferencias divididas: n

∑

f [ x0 ,…, xn ] =

k = 0

f ( xk ) n

∏ ( x

k

.

− xi )

i =0 i ≠ k

Consecuencia: f [ x0 ,…, xn ] es simétrica respecto de

0

,…, xn .

Demostración: identificando coef. de gr. n en las fórmulas de Lagrange y Newton: p ( x ) =

k −1

n

∑

∏

f [ x0 ,…, xk ]

k = 0

( x − xi ) =

i =0

n

x − xi

n

∑ ∏ x y j

j = 0

i=0 i ≠ j

− xi

j

.

Interpolación

29/90


Propiedad de las diferencias divididas: f [ x0 ,…, xn ] =

f [ x1 ,…, xn ] − f [ x0 ,…, xn −1 ] xn − x0

.

Demostración: identificando coef. de gr. n − 1 en las fórmulas de Newton directa e inversa de gr. n : p ( x ) =

k −1

n

n

∑ f [ x ,…, x ]∏ ( x − x ) = ∑ f [ x ,…, x 0

k = 0

k

i

i =0

k = 0 n −1

∑ x

f [ x0 ,…, xn −1 ] − f [ x0 ,…, xn ]

i

i =0

Interpolación

n

n

n − k

]

∏ ( x − x ) i

i = k +1 n

= f [ xn ,…, x1 ] − f [ x0 ,…, xn ]∑ xi . i =1

30/90

Diferencias divididas Tabla de diferencias divididas

f [ x0 ] = y0

x0

f [ x0 , x1 ]

f [ x1 ] = y1

x1

f [ x0 , x1 , x2 ]

f [ x2 ] = y2

x2



f [ x1 , x2 ] 

f [ x0 ,…, xn ]

 





f [ xn −1 , xn ]

f [ xn ] = yn

xn



Interpolación

31/90


Ejemplo:

x

−1

y

0

0

1

3

1 −2

4

−1

;

0 1

0 1

1

−2

−3

−2 1 2

3 3

4

p ( x ) = 0 + 1( x + 1) − 2( x + 1) x + 1( x + 1) x ( − 1) .

Interpolación

32/90

El error de interpolación

EL ERROR DE INTERPOLACIÓN EN EL P.I.P.L.

Interpolación

33/90

El error de interpolación Sean zi

= f ( xi )

los datos de un P.I.P.L., provenientes de una función f ( x ) .

Definición: Error de interpolación E ( x ) = f ( x ) − p( x ) . Teorema: E ( x ) = f [ x0 , x1, …, xn , x]Π( x) , donde Π ( x ) =

n

∏ ( x − x ). i

i =0

Demostración: Ampliando la fórmula de Newton con xn +1 se tiene f [ x0 ,…, xn +1 ] =

yn+1 − p( xn +1 ) n

∏(

n +1

⇒ f [ x0 ,…, xn+1 ]Π( xn +1 ) = E ( xn +1 ) ,

− xi )

i =0

válido

∀ xn +1 .

Interpolación

34/90

El error de interpolación

Teorema: Sea f ∈ C n [a , b ] y x0 ,…, xn ∈ [ a, b] . Entonces ∃ξ ∈ [a, b] tal que f [ x0 ,…, xn ] =

f ( n ) (ξ ) n!

.

Demostración: E ( x ) ∈ C n [a , b] y se anula en x0 ,… , xn → n + 1 veces; aplicando el T. de Rolle E ′( x ) ∈ C

n −1

[a , b ] se anula n veces;



⇒ E ( n ) ( x ) ∈ C 0 [a , b ]

se anula 1 vez.

Consecuencia: Si f ∈ C n +1 [a , b ] entonces E ( x ) =

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

Π( x ) .

Interpolación

35/90

Interpolación por recurrencia

INTERPOLACIÓN POR RECURRENCIA

Interpolación

36/90


Interpolación por recurrencia:

∀C ⊂ { x0 , x1 ,…, n }, sea pC ( x ) el interpolante en los puntos de C . Teorema (Lema de Aitken): Sean xi , x j ∉ C . Entonces pC ∪{ x i , x j } ( x ) =

( x − x j ) pC ∪{ x i } ( x ) − ( x − xi ) pC ∪{ x j } ( x ) xi − x j

.

Demostración: basta comprobar. Principal aplicación: evaluación del interpolante sin construirlo.

Interpolación

37/90


Método de Neville: x − x0 x − x1 x − x2

{ x 0 }

( x ) = y0

p{ x1 } ( x ) = y1 p{ x 2 } ( x ) = y2

p{ x 0 , x1 } ( x ) p{ x0 , x1 , x 2 } ( x ) 

p{ x1 , x 2 } ( x )  



x − xn

Interpolación



p{ x n } ( x ) = yn

p{ x 0 ,…, x n } ( x ) 



p{ x n −1 , x n } ( x )

38/90


Ejemplo del método de Neville:

3

x

−1

y

0

1

3

1 −2

4

en x

=2

0 3

2

1

1

−2

−1

0

−5

−9 −1

− 3 ⇒ p( 2) = −3

1 4

Interpolación

39/90


Método de Aitken: x − x0 x − x1 x − x2  x − xn

( x ) = y0 p{ x1 } ( x ) = y1 p{ x 2 } ( x ) = y2  p{ x n } ( x ) = yn { x 0 }

p{ x 0 , x1 } ( x ) p{ x 0 , x 2 } ( x ) p{ x 0 , x1 , x 2 } ( x )   p{ x 0 , x n } ( x ) p{ x0 , x1 , x n } ( x )  p{ x 0 ,…, x n } ( x )

Ejemplo del método de Aitken: 3 2 1 −1

Interpolación

0 1 3 ⇒ p ( 2 ) = −3 −2 −3 −9 4 3 3 −3

40/90

Resumen de soluciones para el P.I.P.L.

RESUMEN DE SOLUCIONES PARA EL P.I.P.L.

Interpolación

41/90

Resumen de soluciones para el P.I.P.L. 2

n

Con la base canónica {1, x , x ,…, x } n

Solución: (sistema de ecuaciones lineales) p ( x )

= ∑ a j x j j = 0

Con la base de Lagrange { 0 ,  1 ,…,  n } n

Solución: fórmula de Lagrange p ( x )

= ∑ y j  j ( x ) j = 0

Con base de Newton {1, ( x − x0 ), ( x − x0 )( x − x1 ),…, ( x − x0 )( x − xn −1 )} Solución: fórmula de Newton

p( x ) =

∑ f [ x k = 0

Interpolación

k −1

n

0

∏ ( x − x ) .

, … , x k ]

i

i =0

42/90

Evaluación de polinomios

EVALUACIÓN DE POLINOMIOS

Interpolación

43/90

Evaluación de polinomios Para evaluar p ( x )

= an x n + an −1 x n −1 +  + a1 x + a0 se requieren: n( n − 1) 2

multiplicaciones para potencias

n

multiplicaciones por coeficientes

n

sumas de monomios

n 2 + 3n 2

operaciones

...y los errores de redondeo se amplifican.

Interpolación

44/90

Evaluación de polinomios Expresión equivalente:

p( x ) = a0 + x ( a1 + x ( a2 + x (( an −1 + ( an )))) . Algoritmo de Horner para evaluar p ( x ) en el punto x = r : s = a[[n]]; For[i=n-1, i>=0, i--, s = s*r+a[[i]] ];

Algoritmo de Horner para la fórmula de Newton

p( x ) =

0

+ ( x − x0 )(

1

+ ( x − x1 )((

n −1

+ ( x − xn −1 )( n ))))

s = A[[n]]; For[i=n-1, i>=0, i--, s = s*(r-x[[i]]) + A[[i]] ];

Interpolación

45/90

Otros Problemas de Interpolación

OTROS PROBLEMAS DE INTERPOLACIÓN

Interpolación

46/90

Otros Problemas de Interpolación Problema de Interpolación de Taylor Obténgase

∈ Pn conocidos los valores de p( a ), p′( a ), p′′( a ),…, p ( n ) ( a ) .

Solución 1: Con la base canónica {1, x , x 2 ,…, x n } resulta el sistema

Solución 2: Con la base de Newton {1, ( x − a ), ( x − a ) 2 ,…, ( x − a ) n }

triangular superior:

resulta el sistema diagonal:

⎛ 1 ⎜ ⎜0 ⎜⎜  ⎝ 0

a 1

p( a ) ⎞ ⎟ ′ p ( a ) ⎟

a2  an 2a  na n −1





0

0

 





⎟⎟ p ( a ) (n)

n!

⎛ 0! ⎜0 ⎜ ⎜ ⎝ 0

0  1! 

0 0

p( a ) ⎞ p′( a ) ⎟

⎟



   ⎟ (n) 0  n! p ( a ) ⎠

Solución: el desarrollo polinomial de Taylor.

Interpolación

47/90

Otros Problemas de Interpolación Problema de Interpolación de Hermite clásico Obténgase

∈

P 2 n +1 conocidos los valores de

p ( x0 ), p′( x0 ), p( x1 ), p′( x1 ),…, p ( xn ), p′( xn ) .

Solución: Con la base de Newton {1, ( x − x0 ), ( x − x0 ) 2 , ( x − x0 ) 2 ( x − x1 ), ( x − x0 ) 2 ( x − x1 ) 2 ,… 2 2 …, ( x − x0 ) ( x − x1 ) ( x − xn )} x 1 2 3 2 Ejemplo: y

1

2

3

y ′ −1 −1 −1

1.5

1

unisolvente en P 5.

0.5

0.5

Interpolación

1

1.5

2

48/90


Enfoque de Langrange en el caso polinomial con datos tipo Hermite clásico: Datos: z0 , z1 ,…, zn , z0′ , z1′, …, z n′ ; espacio interpolador: V = Base de Lagrange: {

0

, A1 ,…, An , B0 , B1 , …, Bn} cumpliendo j

( xi ) = δ ij ,

B j ( xi ) = 0, n

Sean  j ( x )

=∏ k =0 k≠ j

2 n +1 .

x − xk x j − xk j

A′j ( xi ) = 0, B ′j ( xi ) = δ ij .

los polinomios fundamentales del problema clásico;

= (1 − 2( x − x j )′j ( x j ) )  2j ( x),

B j = ( x − x j ) 2j ( x ). Fórmula de Lagrange: p ( x )

= ∑ z j A j ( x ) + ∑ z ′j B j ( x) . j

j

Ejercicio: comprobar. Interpolación

49/90


Diferencia dividida generalizada: sup. f ∈ C n [a , b] y x0 ≤ x1 ≤  ≤ xn ; ⎧ f [ x1 ,…, xn ] − f [ x0 ,…, xn−1 ] si x0 < xn ⎪ xn − x0 f [ x0 ,…, xn ] = ⎨

⎪ ⎩

1 n!

f ( n ) ( x0 )

si x0

= xn

Enfoque de Newton: cada serie f ( xi ), f ′( xi ),…, f ( k ) ( xi ) genera el fragmento xi f [ xi ] = f ( xi ) f [ xi , xi ] = f ′( xi ) xi f [ xi ] = f ( xi ) f [ xi , xi , xi ] = f [ xi , xi ] = f ′( xi ) xi f [ xi ] = f ( xi )      f [ xi , xi ] = f ′( xi ) xi f [ xi ] = f ( xi )

Interpolación

1 2!

f ′′( x i )  1 k !



f

( k )

( xi )

50/90

Otros Problemas de Interpolación x

1

2

3

Ejemplo: y

1

2

3 unisolvente en P5 .

y ′ −1 −1 −1 Base: {1,( x − 1),( x − 1) ,( x − 1) ( x − 2),( x − 1) ( x − 2) ,( x − 1) ( x − 2) ( x − 3)} 2

2

2

2

2

2

Tabla de diferencias divididas generalizadas:

1

1

1

1

2

2

2

2

3

3

3

3

−1 1

−1 1

−1

2

−2 2

−2

−4 2

−4

3

−3

−3

Interpolación

51/90

El Problema General de Interpolación

EL PROBLEMA GENERAL DE INTERPOLACIÓN

Interpolación

52/90


Definición: Dado V espacio vectorial sobre R , una forma lineal sobre V es L : V  R verificando L(α f + β g ) = α L( f ) + β L( g ) ∀α , β ∈ R ∀f , g ∈ V .

Ejemplos: (V = espacio de funciones R → R ) L1 ( f ) = f (7), L2 ( f ) = f ( x0 ), L3 ( f ) = f ′( −2), L4 ( f ) =

b

∫ f ( x ) dx, a

L5 ( f ) = f ′′′(0).

Interpolación

53/90


Problema General de Interpolación (P.G.I.): “Sea V un espacio vectorial real de dimensión finita n . Dadas las formas lineales { L1 , L2 , …, Ln} sobre V y los valores reales { z1 , z2 ,…, z n } , obténgase p ∈V que verifique Li ( p ) = zi , i

= 1, 2,…, n .”

Teorema: “Sea {v1 , v2 ,…, vn} una base de V . El P.G.I. es unisolvente si y sólo si el determinante de Gram

det( Li ( v j )) es no nulo.”

Demostración: Si p =

∑a v j

j

j

, entonces

∑ a L (v ) = z , i = 1, 2,…, n es un j

i

j

i

j

sistema de ecuaciones lineales cuyo determinante es el de Gram.

Interpolación

54/90


Casos particulares usuales: V = n → interpolación polinomial V = 1,cos x,sen x,cos 2 x,sen 2 x, … → interpolación trigonométrica Li ( f ) = f ( xi ) → datos lagrangianos L2 i ( f ) = f ( xi ), L2 i−1 ( f ) = f ′( xi ) → datos tipo Hermite clásico Li ( f ) = f ( i ) ( x0 ) → datos tipo Taylor

Características del P.G.I.: • el nº de datos es finito • las condiciones de interpolación son lineales • la dimensión del espacio interpolador es igual al nº de datos.

Interpolación

55/90


Construcción de la solución del P.G.I. • Enfoque directo: con el sistema de Gram

∑ a L (v ) = z , i = 1, 2,…, n . j

i

j

i

j

• Enfoque de Lagrange: sea una base {1 ,  2 ,…,  n} de V tal que Li (  j ) = δ ij ; entonces p

= ∑ z j  j

(fórmula de Lagrange).

j

• Enfoque de Newton: sea una base {w1 , w2 ,… , wn } de V con Li ( w j ) = 0 ∀i < j y L j ( w j )

≠ 0 ∀j ; entonces p = ∑ a j w j (fórmula de Newton), donde j

a1 =

1

L1 ( w1 )

,

⎛ k −1 ⎞ zk − Lk ⎜ ∑ a j w j ⎟ ⎝ j =1 ⎠ , k = 2,…, n . ak = Lk ( wk )

Ejercicio: Comprobar todo. Interpolación

56/90

Interpolación mediante splines

INTERPOLACIÓN MEDIANTE FUNCIONES SPLINE

Interpolación

57/90

Interpolación mediante splines Los polinomios de alto grado interpolación adoptan comportamientos anómalos.

20

15

10

5

5

Interpolación

10

15

20

58/90


20

15

10

5

5

10

15

Interpolación

20

59/90


20

15

10

5

5

Interpolación

10

15

20

60/90

Funciones spline

Definición: Un spline de grado n y clase m con nodos { x0 ,…, xN } es una función s( x ) definida en [a , b] = [ x0 , x ] verificando: 1. si

= s [ x − , x ] es un polinomio de grado ≤ n (i = 1,2,…, N ) i 1

i

2. s ∈ C [ a, b] . m

Convenio: salvo indicación expresa, se sobreentiende m = n − 1. Espacios de splines en los nodos { x0 ,…, x } m • de grado n y clase m : S n ( x0 ,…, x N ) . • de grado n (y clase m = n − 1): S n ( x0 ,…, x N ) . Casos usuales: n = 1: grado 1 (clase 0 = continua) → poligonal spline lineal n = 2 : grado 2 (clase 1 = derivable) spline cuadrático n = 3: grado 3 (clase 2 = derivable 2 veces) spline cúbico n = 3, m = 1: spline cúbico de clase 1 (caso especial). Interpolación

61/90

Funciones spline 5

4

3

2

1

2

Interpolación

4

6

8

62/90

Funciones spline 5

4

3

2

1

2

4

6

8

Interpolación

63/90

Funciones spline 5

4

3

2

1

2

Interpolación

4

6

8

64/90

Funciones spline 4

3

2

1

1

2

3

4

5

Interpolación

6

65/90

Construcción de splines por trozos

CONSTRUCCIÓN DE SPLINES POR TROZOS

Interpolación

66/90


Expresión a trozos:

⎧ s1 ( x ) ∈ Pn si x ∈ [ x 0 , x1 ] ⎪ s 2 ( x ) ∈ Pn si x ∈ [ x1 , x 2 ] s ( x ) = ⎨   ⎪ ⎩ s N ( x ) ∈ Pn si x ∈ [ x N −1 , x N ]

Restricciones de suavidad: si( k ) ( xi ) = si(+k 1) ( xi ) i = 1,…, N − 1, k = 0,…, m por tanto

dim S n m ( x0 ,…, x N ) = ( n + 1) N − ( m + 1)( N − 1)

= N (n − m) + m + 1

Interpolación

67/90


Caso lineal ( n = 1): Cada trozo si ( x ) = ai x + bi aporta dos incógnitas; dos ecuaciones: si ( xi −1 ) = yi −1 (extremo izdo.) y si ( xi ) = yi (extremo dcho.) → subproblema 2 × 2 . Construcción directa:

− xi − l

si ( x ) = yi con hi

= xi − xi −1 , i = 1,…,

Ejemplo:

x

−1

y

0

0

hi

+ yi −1

hi

, i = 1,…, N

.

1

3

1 −2

4

s1 ( −1) = 0, s1 (0) = 1 s2 ( 0) = 1, s2 (1) = −2 s3 ( 1) = −2, s3 (3) = 4 Interpolación

xi − x

⇒ s1 ( x ) = x + 1 ⇒ s2 ( x ) = −3 x + 1 ⇒ s3 ( x ) = 3 x − 5 68/90


Caso cuadrático ( n = 2 ): Cada trozo si ( x ) = ai x 2 + bi x + ci aporta tres incógnitas → total 3 N ; Ecuaciones: • Interpolación (clase 0 ): 2 por trozo (extremos izdo. y dcho.) • Clase 1: si′( xi ) = si′+1 ( xi ) i = 1,…, N − 1 son N − 1 ecuaciones.

→ total 2 + N − 1 = 3N − 1 ecuaciones → ¡falta una ecuación!. → Se añade una ecuación adicional para cuadrar el problema. Construcción directa con s′( x0 ) = d 0 : si ( x ) = yi −1 + d i −1 ( x − xi −1 ) + con wi

=

yi − yi −1 hi

wi − d i −1 hi

( x − xi −1 ) 2 , i = 1,…, N

, d i = 2 wi − d i −1 , i = 1,…, N , d 0 = y0′ .

Interpolación

69/90


Caso cúbico ( n = 3): Cada trozo si ( x ) = ai x 3 + bi x 2 + ci x + d i aporta cuatro incógnitas → total 4 N ; Ecuaciones: • Interpolación (clase 0 ): 2 por trozo (extremos izdo. y dcho.) • Clase 1: si′( xi ) = si′+1 ( xi ) i = 1,…, N − 1 son N − 1 ecuaciones.

• Clase 2 : si′′( xi ) = si′′+1 ( xi ) i

= 1,…, N − 1 son N − 1 ecuaciones.

→ total 2

− 2 ecuaciones → ¡faltan dos ecuaciones!.

+ 2( N − 1) = 4

Ecuaciones adicionales típicas: • s′′( a ) = s′′(b) = 0 • •

spline natural prolongable mediante rectas → de clase 2 en todo R . ′ s′( a ) = y0′ , s′(b) = y N spline de extremo sujeto pendientes fijas de entrada y salida. s′( a ) = s′(b), s′′( a ) = s′′(b) spline periódico prolongable por periodicidad → de clase 2 en todo R .

Interpolación

70/90


Construcción directa del spline cúbico natural o sujeto: si ( x ) = yi −1 + d i −1 ( x − xi −1 ) +

+

wi − d i −1 hi

d i −1 + d i − 2 wi hi2

( x − xi −1 ) 2

( x − xi −1 ) 2 ( x − xi ), i = 1,…, N

siendo d i la solución del s.e.l. tridiagonal simétrico

1 hi

⎛ 1

d i −1 + 2⎜

⎝ hi

+

1 ⎞ hi +1

⎟d i +

1 hi +1

⎛ wi +1

d i +1 = 3⎜

⎝ hi +1

+

wi ⎞ ⎟, i = 1,…, N − 1 hi

con las ecuaciones adicionales

′; • natural: d 0 = y0′ , d N = y N • sujeto: 2d 0 + d 1 = 3w1 , d N −1 + 2d

= 3w

.

Interpolación

71/90


Caso especial: n = 3, m = 1 para datos tipo Hermite clásico (valor y derivada en cada punto). Para cada trozo si ( ) se tienen 4 incógnitas y 4 ecuaciones ( 2 valores y 2 derivadas), formando un subproblema 4 × 4 . Caso general (P.I.L. en S n m ( x0 ,…, x N ) ): Hay N ( n + 1) incógnitas; Ecuaciones: • Interpolación+continuidad: 2 N ( = N + 1 interpolación • Suavidad de clase 1: N − 1 • Suavidad de clase 2: N − 1

 • Suavidad de clase m :

+ N − 1 clase 0).

−1

Total ecuaciones: 2 N + m( N − 1) = ( m + 2) N − m ; si m = n − 1 entonces siempre faltan m ecuaciones adicionales. Se obtiene un sistema N ( n + 1) × N ( n + 1) .

Interpolación

72/90


Teoremas de unisolvencia: • P.I.L. en S 1 ( x0 ,…, x ) • P.I.L. en S 2 (

0

,…, x N ) con s′( xi ) = d i para algún i ∈ {0,1,…, N }

• P.I.L. en S 3 ( x0 ,…, x N ) natural • P.I.L. en S 3 ( x0 ,…, x N ) sujeto • P.I.L. en S 3 ( x0 ,…, x N ) periódico 1

• Hermite clásico en S 3 ( x0 ,…, x )

Interpolación

73/90

Construcción de splines mediante Potencias truncadas

CONSTRUCCIÓN DE SPLINES POR POTENCIAS TRUNCADAS

Interpolación

74/90

Construcción de splines mediante Potencias truncadas Potencia truncada

⎧ x n ( x ) + = ⎨ ⎩0 n

si x si x

≥0 <0

Ejemplos: ( x ) + , ( x − 2) 2+ , ( x 2 − 2) + , ( 2 − x 2 ) + , ( x − 3) 3+ , ( x − 4) 3+ 5

4

3

2

1

-

2

2

4

6

Interpolación

75/90


Ejercicios: x = ( x ) + + ( − x ) + , x = ( x ) + − ( − x ) + . Propiedad: ( x − r ) n+ es un spline de grado n con un nodo en r . Bases de splines:

= 1, x, ( x − x1 ) + , ( x − x2 ) + ,…, ( x − x N −1 ) + 2 2 2 2 … … = − − − ( x , , x ) 1 , x , x , ( x x ) , ( x x ) , , ( x x ) S 2 0 N 1 + 2 + N −1 + S 1 ( x0 ,…, x N )

S 3 ( x0 ,…, x N )

= 1, x, x 2 , x 3 , ( x − x1 ) 3+ , ( x − x2 ) 3+ ,…, ( x − x N −1 ) 3+

S n ( x0 ,…, x N )

= 1, x, x 2 ,…, x n , ( x − x1 ) n+ , ( x − x2 ) n+ ,…, ( x − x N −1 ) n+

1

S 3

( x0 ,…, x N ) = 1, x, x 2 , x 3 , ( x − x1 ) 2+ , ( x − x2 ) 2+ ,…, ( x − x N −1 ) 2+ , ( x − x1 ) 3+ , ( x − x2 ) 3+ ,…, ( x − x N −1 )3+

Interpolación

76/90


Ejemplo: Spline cúbico natural para

x

−1

y

0

0

1

3

1 −2

4

s ( x ) = a + bx + cx 2 + dx 3 + A( x )3+ + B ( x − 1) 3+ sistema lineal

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Solución: s ( x ) =

1 23

1 −1 1 0 1 1 1 3 0 0 0 0

1 −1 0 0 1 1 9 27 2 −6 2 18

(23 − 37 x − 90 x

2

0 0 0 0 1 0 27 8 0 0 18 12

0 ⎞ 1⎟ − 2⎟ 4⎟ 0⎟ 0⎟

− 30 x 3 + 88( x ) 3+ − 72( x − 1) 3+ )

Interpolación

77/90


Ejemplo: Spline cúbico sujeto horizontal para

x

−1

y

0

0

1

3

1 −2

4

s ( x ) = a + bx + cx 2 + dx 3 + A( x )3+ + B ( x − 1) 3+ sistema lineal

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Solución: s ( x ) =

Interpolación

1 22

1 −1 1 −1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 1 3 9 27 27 8 0 1 −2 3 0 0 0 1 6 27 27 12

(22 − 27 x − 120 x

2

0 ⎞ 1⎟ − 2⎟ 4⎟ 0⎟ 0 ⎠⎟

− 71 x 3 + 152( x )3+ − 120( x − 1)3+ ) 78/90

Construcción de splines mediante Potencias truncadas 4

3

2

1

-

1

1

-

1

-

2

2

Interpolación

3

79/90


3

2

1

-

1

Interpolación

1

-

1

-

2

2

3

80/90


3

2

1

-

1

1

-

1

-

2

2

Interpolación

3

81/90


3

2

1

-

1

Interpolación

1

-

1

-

2

2

3

82/90

Interpolación de Hermite con Splines

INTERPOLACIÓN DE HERMITE CON SPLINES

Interpolación

83/90

Interpolación de Hermite con Splines x Ejemplo: y y ′

0 1 0.5 1 0 0

2 2 0

4 1.5 en S 31 (0,1,2,4) 0

2

1.5

1

0.5

1

-

2

3

4

0.5

Interpolación

84/90

Propiedades extremales de los splines cúbicos

PROPIEDADES EXTREMALES DE LOS SPLINES CÚBICOS

Interpolación

85/90


Condiciones comunes: “Sean los nodos a

= x0 < x1 <  < x N = b . Sea s ( x ) el spline cúbico

natural o sujeto que interpola los nodos con datos y 0 , … , y , y adicionales

′. s′′( a ) = s′′(b) = 0 o bien s′( a ) = y0′ , s′(b) = y N 2 Sea f ∈ C [a , b] cualquier otra función que interpole los mismos datos, incluyendo los adicionales. Sea E ( ) = f ( x ) − s ( x ) .” Lema: b

∫ s′′( x) E ′′( x )dx = 0. a

Interpolación

86/90


Demostración:

∫

N

∑∫

b

s′′( x ) E ′′( x )dx =

a

i =1 N

x i

x i −1

s′′( x ) E ′′( x )dx = ⎡

⎤ ⎢⎣u = s′′, dv = E ′′dx ⎥⎦

(

(por partes)

x = ∑ s′′( x ) E ′( x ) x − − ∫ s′′′( x ) E ′( x )dx i =1

x i

i

x i −1

i 1

N

x i

i =1

x i −1

)

= s′′(b) E ′(b) − s′′( a ) E ′( a ) − ∑ K i ∫ E ′( x )dx si es natural: s ′′( a )

= s′′(b) = 0 , si es sujeto: E ′( a ) = ′(b) = 0 , y además

∫

x i

x i −1

E ′( x )dx = E ( x )

x i xi −1

= 0.

Interpolación

87/90


Teorema (propiedad extremal):

∫

b

s ′′( x ) dx ≤ 2

a

∫

b

a

f ′′( x ) 2 dx .

Demostración:

∫

b

a

b

∫ = ∫ E ′′( x ) dx + ∫ s′′( x ) dx + 0 ≥ ∫ s′′( x ) dx.

f ′′( x ) dx = 2

a

( E ′′( x ) + s′′( x ) )2 dx

b

2

a

b

b

2

a

2

a

Interpolación

88/90

Cálculo Numérico-[Miguel Pasadas Fernández]

Recommend Documents