ULPGC

Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 1 de 22
Volver
Pantalla completa
Cerrar
Salir
Tutorial de An

alisis Num

erico
Interpolaci on : Splines c ubicos
Jes us Garca Quesada
Departamento de Inform atica y Sistemas
Universidad de Las Palmas de Gran Canaria
35017 Campus de Tara, Espa na
Email : [email protected]
2 de Octubre de 2000, v0.3
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 2 de 22
Volver
Pantalla completa
Cerrar
Salir

Indice General
1 INTERPOLACI

ON POR SPLINES C

UBICOS 3
2 PROBLEMAS 16
Soluciones a los Problemas 19
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 3 de 22
Volver
Pantalla completa
Cerrar
Salir
1. INTERPOLACI

ON POR SPLINES C

UBICOS
Supongamos que tenemos los n + 1 puntos:
P
k
(x
k
, y
k
), donde y
k
= f(x
k
), k = 0, 1, . . . , n
en los cuales se quiere interpolar la funcion f. Las abcisas no es necesario que sean
equidistantes, pero se suponen ordenados, o sea,
x
0
< x
1
< x
2
< < x
n
En esta seccion trataremos de la interpolaci on polinomica a trozos. La idea es encon-
trar polinomios c ubicos q
k
(x) que interpolen la funcion f en el subintervalo [x
k
, x
k+1
], k =
0, 1, . . . , n 1.
Denicion 1. La funcion s(x) se llama c ubica a trozos en [x
0
, x
n
] si existen polinomios
c ubicos q
0
(x), q
1
(x), . . . , q
n1
(x) tales que :
s(x) = q
k
(x) en [x
k
, x
k+1
], para k = 0, 1, . . . , n 1
Para que s(x) interpole en los puntos P
0
, P
1
, . . . , P
n
los q
k
(x) han de vericar :
_
q
k
(x
k
) = y
k
q
k
(x
k+1
) = y
k+1
, k = 0, 1, . . . , n 1
(1)
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 4 de 22
Volver
Pantalla completa
Cerrar
Salir
P
0
P
1
P
2
P
3
P
k
P
k+1
P
n-1
P
n
. . . . . . . . . . . . . . .
q (x)
0
x
n-1
q (x)
x x x x x x x
0 1 3 k
k
q (x)
1
q (x) en [x , x ]
k k+1
solo existe
2 k+1 n-1 n
Figura 1: Interpolaci on polinomica a trozos.
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 5 de 22
Volver
Pantalla completa
Cerrar
Salir
lo cual supone 2n condiciones. Llamaremos a s(x) spline c ubico, o simplemente spline,
si los polinomios q
k
(x) tienen la misma pendiente y la misma concavidad en los nodos
que las unen, o sea :
_
q

k1
(x
k
) = q

k
(x
k
)
q

k1
(x
k
) = q

k
(x
k
), k = 1, 2, . . . , n 1
(2)
lo cual supone 2(n 1) condiciones a cumplir. Al tener que vericar las condiciones
(1) y (2) se asegura que s(x) tiene su primera y segunda derivadas continuas en [x
0
, x
n
].
En este caso se dice que s(x) es un spline interpolador para P
0
, P
1
, . . . , P
n
.
Si s(x) es c ubica a trozos en el intervalo [x
0
, x
n
], su derivada segunda s

(x) es lineal en
el mismo intervalo e interpola en los puntos (x
k
, s

(x
k
)) y (x
k+1
, s

(x
k+1
)) en [x
k
, x
k+1
].
Por tanto, q
k
(x) es un polinomio de grado uno que interpola en los puntos (x
k
, s

(x
k
)) y
(x
k+1
, s

(x
k+1
)):
q

k
(x) = s

(x
k
)
x x
k+1
x
k
x
k+1
+s

(x
k+1
)
x x
k
x
k+1
x
k
, para k = 0, 1, . . . , n 1
Denotando con
h
k
= x
k+1
x
k
, k = 0, 1, . . . , n 1
y

k
= s

(x
k
), k = 0, 1, . . . , n
tenemos :
q

k
(x) =

k
h
k
(x
k+1
x) +

k+1
h
k
(x x
k
), para k = 0, 1, . . . , n 1 (3)
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 6 de 22
Volver
Pantalla completa
Cerrar
Salir
donde h
k
y
k
son constantes (
k
a determinar). Integrando dos veces :
q
k
(x) =

k
h
k
(x
k+1
x)
3
6
+

k+1
h
k
(x x
k
)
3
6
+C
k
+D
k
x (4)
donde el termino lineal lo podemos escribir como :
C
k
+D
k
x = A
k
(x x
k
) +B
k
(x
k+1
x)
siendo A
k
y B
k
constantes arbitrarias, quedando entonces :
q
k
(x) =

k
h
k
(x
k+1
x)
3
6
+

k+1
h
k
(x x
k
)
3
6
+A
k
(x x
k
) +B
k
(x
k+1
x) (5)
Aplicando a (5) las condiciones (1) :
y
k
=

k
h
k
h
3
k
6
+

k+1
h
k
0 +A
k
.0 +B
k
h
k
=

k
6
h
2
k
+B
k
h
k
(6)
y
k+1
=

k+1
h
k
h
3
k
+A
k
h
k
=

k+1
6
h
2
k
+A
k
h
k
(7)
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 7 de 22
Volver
Pantalla completa
Cerrar
Salir
y despejando de aqu A
k
y B
k
y sustituyendo en (5) resulta :
q
k
(x) =

k
6
_
(x
k+1
x)
3
h
k
h
k
(x
k+1
x)
_
+

k+1
6
_
(x x
k
)
3
h
k
h
k
(x x
k
)
_
+y
k
_
x
k+1
x
h
k
_
+y
k+1
_
x x
k
h
k
_
, para k = 0, 1, . . . , n 1
(8)
que es la ecuacion del spline q
k
(x).
Nos falta a un conocer los valores
0
,
1
, . . . ,
n
(n+1 incognitas) para lo cual usamos
(2); derivando en (8) tenemos :
q

k
(x) =

k
6
_
3(x
k+1
x)
2
h
k
+h
k
_
+

k+1
6
_
3(x
k
x)
2
h
k
h
k
_
+
y
k+1
y
k
h
k
Por tanto :
q

k
(x
k
) =

k
6
(2h
k
) +

k+1
6
(h
k
) +
y
k+1
y
k
h
k
(9)
q

k
(x
k+1
) =

k
6
(h
k
) +

k+1
6
(2h
k
) +
y
k+1
y
k
h
k
(10)
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 8 de 22
Volver
Pantalla completa
Cerrar
Salir
Reemplazando k por k 1 en (10) para obtener q

k1
(x
k
) e igualando a (9) nos da :
h
k1

k1
+ 2 (h
k1
+h
k
)
k
+h
k

k+1
= 6
_
y
k+1
y
k
h
k

y
k
y
k1
h
k1
_
, para k = 1, 2, . . . , n 1
(11)
o tambien
h
k1

k1
+ 2 (h
k1
+h
k
)
k
+h
k

k+1
= 6
_
y
k
h
k

y
k1
h
k1
_
, para k = 1, 2, . . . , n 1
(12)
o incluso
h
k1

k1
+ 2 (h
k1
+h
k
)
k
+h
k

k+1
= 6
_
f[x
k
, x
k+1
] f[x
k1
, x
k
]
_
, para k = 1, 2, . . . , n 1
(13)
Como el ndice k vara de 1 a n 1, se producen n 1 ecuaciones lineales con n + 1
incognitas
0
,
1
, . . . ,
n
, lo cual produce un sistema subdeterminado que tiene innitas
soluciones.
Existen varias estrategias para eliminar
0
de la primera ecuacion y
n
de la (n-1)-
esima produciendo un sistema tridiagonal de orden (n-1) en las variables
1
,
2
, . . . ,
n1
.
ALTERNATIVA I Especicar el valor de s

(x) en los puntos extremos :

0
= s

(x
0
)
y
n
= s

(x
n
). Si se pone
0
= 0,
n
= 0 se denomina spline c ubico natural.
ALTERNATIVA II Suponer que s

(x) es constante en los extremos :

0
=
1
y

n
=
n1
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 9 de 22
Volver
Pantalla completa
Cerrar
Salir
ALTERNATIVA III Suponer que s

(x) es lineal cerca de los extremos :

0
=
1
h
1
_
(h
0
+
h
1
)
1
h
0

2
_
y

n
=
1
h
n2
_
h
n1
)
n2
+ (h
n2
+h
n1
)
n1
_
ALTERNATIVA IV Especicar el valor de s

(x) en los puntos extremos :

0
=
3
h
0
_
y
0
s

(x
0
)

1
2

1
y

n
=
3
h
n1
_
s

(x
n
) y
n1

1
2

n1
Si hay que calcular muchas veces s(z) entonces es preferible hacer la sustitucion :
x
k+1
z = (x
k+1
x
k
) (z x
k
) = h
k
(z x
k
)
en q
k
(z) y entonces expresar este en potencias de z x
k
para obtener :
q
k
(z) = y
k
+
1
(z x
k
) +
2
(z x
k
)
2
+
3
(z x
k
)
3
= y
k
+ (z x
k
)(
1
+ (z x
k
)(
2
+ (z x
k
)
3
))
evaluado con solo 4 sumas/restas y 3 productos, donde

1
= f[x
k
, x
k+1
]
h
k
6
(
k+1
+ 2
k
),
2
=

k
2
,
3
=

k+1

k
6h
k
En forma matricial, el sistema tridiagonal que resulta es (caso de spline c ubico natu-
ral):
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 10 de 22
Volver
Pantalla completa
Cerrar
Salir
_

_
2 (h
0
+h
1
) h
1
0
h
1
2 (h
1
+h
2
) 0
0 h
2
0
.
.
.
.
.
. 2 (h
n3
+h
n2
) h
n2
0 0 h
n2
2 (h
n2
+h
n1
)
_

_
_

3
.
.
.

n1
_

_
=
= 6
_

_
f[x
1
, x
2
] f[x
0
, x
1
]
f[x
2
, x
3
] f[x
1
, x
2
]
f[x
3
, x
4
] f[x
2
, x
3
]
.
.
.
f[x
n1
, x
n
] f[x
n1
, x
n2
]
_

_
o tambien
_

_
2 (h
0
+h
1
) h
1
0
h
1
2 (h
1
+h
2
) 0
0 h
2
0
.
.
.
.
.
. 2 (h
n3
+h
n2
) h
n2
0 0 h
n2
2 (h
n2
+h
n1
)
_

_
_

3
.
.
.

n1
_

_
= 6
_

_
y
1
h
1

y
0
h
0
y
2
h
2

y
1
h
1
y
3
h
3

y
2
h
2
.
.
.
y
n1
h
n1

y
n2
h
n2
_

_
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 11 de 22
Volver
Pantalla completa
Cerrar
Salir
Ejemplo.
Interpolar por splines c ubicos la funcion f(x) = 1/x en x = 1.5 tomando los puntos
(0.1, 10.0), (0.2, 5.0), (0.5, 2.0), (1.0, 1.0), (2.0, 0.5), (5.0, 0.2), (10.0, 0.1).
Solucion:
h
0
= 0.2 0.1 = 0.1 h
3
= 2.0 1.0 = 1.0
h
1
= 0.5 0.2 = 0.3 h
4
= 5.0 2.0 = 3.0
h
2
= 1.0 0.5 = 0.5 h
5
= 10.0 5.0 = 5.0
El sistema que resulta es
0.1
0
+ 2 (0.1 + 0.3)
1
+ 0.3
2
= 6
_
2 5
0.5 0.2

5 10
0.2 0.1
_
0.3
1
+ 2 (0.3 + 0.5)
2
+ 0.5
3
= 6
_
1 2
1.0 0.5

2 5
0.5 0.2
_
0.5
2
+ 2 (0.5 + 1.0)
3
+ 1.0
4
= 6
_
0.5 1.0
2.0 1.0

1 2
1.0 0.5
_
1.0
3
+ 2 (1.0 + 3.0)
4
+ 3.0
5
= 6
_
0.2 0.5
5.0 2.0

0.5 1.0
2.0 1.0
_
3.0
4
+ 2 (3.0 + 5.0)
5
+ 5.0
6
= 6
_
0.1 0.2
1.0 5.0

0.2 0.5
5.0 2.0
_
Poniendo
0
=
6
= 0 tenemos
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 12 de 22
Volver
Pantalla completa
Cerrar
Salir
_

_
0.8 0.3 0 0 0
0.3 1.6 0.5 0 0
0 0.5 3.0 1.0 0
0 0 1.0 8.0 3.0
0 0 0 3.0 16.0
_

_
_

5
_

_
= 6
_

_
40
8
1.5
0.4
0.08
_

_
de donde se obtienen

1
= 311.65398570643
2
= 31.077295217152

3
= 8.4549532710280
4
= .82621220450797

5
= 0.18491478834524
Para x = 1.5 habr a que elegir q
3
(x)
q
3
(x) =

3
6
_
(2.0 x)
3
1.0
1.0 (2.0 x)
_
+

4
6
_
(x 1.0)
3
1.0
1.0 (x 1.0)
_
+1.0(2.0x)+0.5(x1.0)
y su valor q
3
(x) = 0.27320367097855 es una mejor estimacion que la obtenida por
interpolaci on polinomica (ver gura 2).
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 13 de 22
Volver
Pantalla completa
Cerrar
Salir
Los diferentes splines que resultan son:
q
0
(x) =519.423309x
3
155.826993x
2
39.611534x + 15
q
1
(x) =190.406267x
3
+ 270.070753x
2
124.791083x + 20.678637
q
2
(x) =13.177416x
3
35.304772x
2
+ 27.896679x 4.769324
q
3
(x) =1.546861x
3
+ 8.868059x
2
16.276152x + 9.954953
q
4
(x) = + 0.0561737x
3
0.750148x
2
+ 2.960264x 2.869324
q
5
(x) =0.00616383x
3
+ 0.1849148x
2
1.715052x + 4.922870
s(x) =
_

_
q
0
(x) , si x [0.1, 0.2], (tambien (, 0.2])
q
1
(x) , si x [0.2, 0.5],
q
2
(x) , si x [0.5, 1.0],
q
3
(x) , si x [1.0, 2.0],
q
4
(x) , si x [2.0, 5.0],
q
5
(x) , si x [5.0, 10.0], (tambien [5.0, +))
En el intervalo [1.0, 2.0], la representacion de ambas funciones es la que aparece en la
gura 2.

Ejemplo. Interpolar por splines c ubicos la funcion f(x) = 1/(x

2
+ 1) en el intervalo
0 x 1 tomando los seis puntos de abcisas x
k
= k/5, k = 0, 1, 2, 3, 4, 5.
Solucion:
Por calculo directo tenemos :
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 14 de 22
Volver
Pantalla completa
Cerrar
Salir
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1 1.2 1.4 1.6 1.8 2
1/x
8.4575/6*((2-x)**3-(2-x))-0.8334/6*((x-1)**3-(x-1))+(2-x)+0.5*(x-1)
Figura 2: La funcion 1/x y q
3
(x) en el intervalo [1.0,2.0].
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 15 de 22
Volver
Pantalla completa
Cerrar
Salir
y
0
= 1.00000000 y
3
= 0.73529412
y
1
= 0.96153846 y
4
= 0.60975610
y
2
= 0.86206896 y
5
= 0.50000000
y es h
0
= h
1
= h
2
= h
3
= h
4
= 1/5 = h y poniendo
0
=
5
= 0 y multiplicando
ambas partes por 6/h :
_

_
4 1 0 0
1 4 1 0
0 1 4 1
0 0 1 4
_

_
_

4
_

_
=
6
h
2
_

_
y
2
2 y
1
+y
0
y
3
2 y
2
+y
1
y
4
2 y
3
+y
2
y
5
2 y
4
+y
3
_

_
=
_

_
9.151194
4.095801
0.185523
2.367288
_

_
y resolviendo
1
= 2.165814,
2
= 0.487920,
3
= 0.022536,
4
= 0.581866 y
tabulando la funcion entre 0 y 1.0 con paso 0.002 la graca de f(x) y el spline c ubico son
indistinguibles (error maximo 0.0040 que se produce entre 0 y 0.2).

ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 16 de 22
Volver
Pantalla completa
Cerrar
Salir
2. PROBLEMAS
Problema 1. Construir el spline c ubico natural que interpola a partir de los datos:
x 0 1 2 2.5 3 4
y 1.4 0.6 1.0 0.65 0.6 1.0
Problema 2. Considerando los datos:
x 0.15 0.76 0.89 1.07 1.73 2.11
y 0.3495 0.2989 0.2685 0.2251 0.0893 0.0431
obtener el spline c ubico natural que interpola en dichos puntos.
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 17 de 22
Volver
Pantalla completa
Cerrar
Salir
Referencias
[Act90] F.S. Acton. Numerical Methods That (Usually) Work. The Mathematical As-
sociation of America, Washington, 1990.
[Atk89] K. E. Atkinson. An Introduction to Numerical Analysis. John Wiley, New York,
2nd. edition, 1989.
[BF80] R.L. Burden and D. Faires. Analisis Numerico. Grupo Editorial Iberoameri-
cana, Mexico, 1980.
[CC89] S.C. Chapra and R.P. Canale. Numerical Methods for Engineers. McGraw-Hill
International, New York, second edition, 1989.
[CdB80] S. D. Conte and C. de Boor. Elementary Numerical Analysis: An Algorithmic
Approach. McGrawHill, New York, third edition, 1980.
[DB74] Germund Dahlquist and

Ake Bjorck. Numerical Methods. Prentice-Hall, En-
glewood Clis, New Jersey, 1974.
[Fad59] V.N. Faddeeva. Computational Methods of Linear Algebra. Dover Publications,
Inc, New York, 1959.
[Fro79] C.-E. Froberg. Introduction to Numerical Analysis. AdisonWesley, Reading,
Massachusetts, 2nd. edition, 1979.
[GW89] C.F. Gerald and P.O. Wheatley. Applied Numerical Analysis. AddisonWesley
Publishing Co., Reading, Massachusets, fourth edition, 1989.
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 18 de 22
Volver
Pantalla completa
Cerrar
Salir
[Hen72] P. Henrici. Elementos de Analisis Numerico. Ed. Trillas, Mexico, 1972.
[Hil74] F. B. Hildebrand. Introduction to Numerical Analysis. McGrawHill, New
York, second edition, 1974.
[KC94] D. Kincaid and W. Cheney. Analisis Numerico : las matematicas del calculo
cientco. Addison-Wesley Iberoamericana, 1994.
[Mar87] M. J. Maron. Numerical Analysis: A Practical Approach. Macmillan Publishing
Co., New York, second edition, 1987.
[ML91] M. J. Maron and R. J. Lopez. Numerical Analysis: A Practical Approach.
Wadsworth, Belmont, California, third edition, 1991.
[RR78] Anthony Ralston and Philip Rabinowitz. A First Course in Numerical Analysis.
McGraw-Hill, New York, 2nd. edition, 1978.
[Sch89] H.R. Schwarz. Numerical Analysis. John Wiley & Sons, Chichester, 1989.
[Wer84] W. Werner. Mathematics of Computation, 43:205217, 1984.
[YG73a] David M. Young and R.T. Gregory. A Survey of Numerical Mathematics, vol-
ume I. Dover Publications, New York, 1973.
[YG73b] David M. Young and R.T. Gregory. A Survey of Numerical Mathematics, vol-
ume II. Dover Publications, New York, 1973.
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 19 de 22
Volver
Pantalla completa
Cerrar
Salir
Soluciones a los Problemas
Problema 1. El sistema es:
_

_
4 1 0 0
1 3 0.5 0
0 0.5 2 0.5
0 0 0.5 3
_

_
_

4
_

_
=
_

_
7.2
6.6
3.6
3
_

_
de donde se obtienen

1
= 2.6788381742739
2
= 3.5153526970954

3
= 2.5344398340249
4
= 0.57759336099585
Los diferentes splines son:
q
0
(x) = + 0.44647303x
3
1.24647303x + 1.4
q
1
(x) =1.03236514x
3
+ 4.436514523x
2
5.68298755x + 2.87883817
q
2
(x) = + 2.01659751x
3
13.85726141x
2
+ 30.90456432x 21.51286307
q
3
(x) =0.65228216x
3
+ 6.15933610x
2
19.13692946x + 20.18838174
q
4
(x) =0.09626556x
3
+ 1.15518672x
2
4.12448133x + 5.17593361
ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 20 de 22
Volver
Pantalla completa
Cerrar
Salir
s(x) =
_

_
q
0
(x) , si x [0, 1], (tambien (, 0])
q
1
(x) , si x [1, 2],
q
2
(x) , si x [2, 2.5],
q
3
(x) , si x [2.5, 3],
q
4
(x) , si x [3, 4], (tambien [3, +))

ULPGC
Inform atica
Pagina Web
Pagina de Inicio
Contenido


Pagina 21 de 22
Volver
Pantalla completa
Cerrar
Salir
Problema 2. El sistema es ahora:
_

_
1.48 0.13 0 0
0.13 0.62 0.18 0
0 0.18 1.68 0.66
0 0 0.66 2.08
_

_
_

4
_

_
=
_

_
0.905372005
0.0435897436
0.212121212
0.50507177
_

_
obteniendose:

1
= 0.61616885710569
2
= 0.050445411325263

3
= 0.029089182290732
4
= 0.23359274520339
Los diferentes splines son:
q
0
(x) =0.16835215x
3
+ 0.075758466x
2
0.0316707558x + 0.35311424
q
1
(x) = + 0.85463368x
3
2.25664921x
2
+ 1.74095908x 0.095951989
q
2
(x) =0.019774286x
3
+ 0.078020050x
2
0.33689656x + 0.52047852
q
3
(x) = + 0.051642314x
3
0.15122724x
2
0.091601967x + 0.43299011
q
4
(x) =0.10245296x
3
+ 0.64852723x
2
1.47517719x + 1.23085182
s(x) =
_

_
q
0
(x) , si x [0.15, 0.76], (tambien (, 0.15])
q
1
(x) , si x [0.76, 0.89],
q
2
(x) , si x [0.89, 1.07],
q
3
(x) , si x [1.07, 1.73],
q
4
(x) , si x [1.73, 2.11], (tambien [1.73, +))