interpolação Calculo Numérico

Prévia do material em texto

Interpolação e ajuste não-segmentados
1 Introdução
O problema geral da interpolação pode ser de�nido da seguinte forma:
Seja F uma família de funções f : D → E e {(xi, yi)}Ni=1 um conjunto de pares ordenados tais que xi ∈ D
e yi ∈ E, encontrar uma função f da família dada tal que f(xi) = yi para cada 1 ≤ i ≤ N .
Exemplo 1. Encontrar uma função f(x) da forma f(x) = aebx onde a e b são constantes tal que f(1) = 1 e
f(2) = 5. Este problema equivale a resolver o seguinte sistema de equações:
aeb = 1
ae2b = 5
Dividindo a segunda equação pela primeira, temos eb = 5, logo, b = ln(5). Substituindo este valor em qualquer
das equações, temos a = 1
5
. Assim
f(x) =
1
5
eln(5)x =
1
5
5x = 5x−1.
Exemplo 2. Encontrar a função polinomial do tipo f(x) = a+bx+cx2 que passe pelos pontos (−1, 2), (0, 1),
(1, 6). Observamos que podemos encontrar os coe�cientes a, b e c através do seguinte sistema linear:
a− b+ c = 2
a = 1
a+ b+ c = 6
cuja solução é dada por a = 1, b = 2 e c = 3. Portanto
f(x) = 1 + 2x+ 3x2.
2 Interpolação polinomial
Interpolação polinomial é o caso particular do problema geral de interpolação quando a família de funções é
constituída de polinômios.
Teorema 1. Seja {(xi, yi)}ni=0 um conjunto de n+ 1 pares ordenados de números reais tais que
i 6= j =⇒ xi 6= xj (i.e. as abscissas são distintas)
então existe um único polinômio P (x) de grau igual ou inferior a n que passa por todos os pontos dados.
Demonstração. Observamos que o problema de encontrar os coe�cientes a0, a1,. . . , an do polinômio
P (x) = a0 + a1x+ a2x
2 + · · · anxn =
n∑
k=0
akx
k
tal que P (xi) = yi é equilavente ao seguinte sistema linear de n+1 equações e n+1 incógnitas:
a0 + a1x0 + a2x
2
0 + · · ·+ anxn0 = y0
a0 + a1x2 + a2x
2
1 + · · ·+ anxn1 = y1
.
.
.
a0 + a1xn + a2x
2
n + · · ·+ anxnn = yn
1
que pode ser escrito na forma matricial como
1 x0 x
2
0 · · · xn0
1 x1 x
2
1 · · · xn1
1 x2 x
2
2 · · · xn2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 xn x
2
n · · · xnn


a0
a1
a2
.
.
.
an
 =

y0
y1
y2
.
.
.
yn

A matriz envolvida é uma matriz de Vandermonde de ordem n+ 1 cujo determinante é dado por∏
0≤i<j≤n
(xj − xi)
É fácil ver que se as abscissas são diferentes dois a dois, então o determinante é não-nulo. Disto decorre que
o sistema possui um a solução e que esta solução é única.
Exemplo 3. Encontre o polinômio da forma P (x) = a0 + a1x+ a2x
2 + a3x
3
que passa pelos pontos
(0, 1), (1, 2), (2, 4), (3, 8)
Este problema é equivalente ao seguinte sistema linear:
a0 = 1
a0 + a1 + a2 + a3 = 2
a0 + 2a1 + 4a2 + 8a3 = 4
a0 + 3a1 + 9a2 + 27a3 = 8
cuja solução é a0 = 1, a1 =
5
6
, a2 = 0 e a3 =
1
6
. Portanto
P (x) = 1 +
5
6
x+
1
6
x3
Exemplo 4. Encontre o polinômio da forma P (x) = a0 + a1x+ a2x
2 + a3x
3
que passa pelos pontos
(0, 0), (1, 1), (2, 4), (3, 9)
Este problema é equivalente ao seguinte sistema linear:
a0 = 0
a0 + a1 + a2 + a3 = 1
a0 + 2a1 + 4a2 + 8a3 = 4
a0 + 3a1 + 9a2 + 27a3 = 9
cuja solução é a0 = 0, a1 = 0, a2 = 1 e a3 = 0. Portanto
P (x) = x2
Esta abordagem direta que �zemos ao calcular os coe�cientes do polinômio na base canônica se mostra ine-
�ciente quando o número de pontos é grande e quando existe grande discrepância em as abscissas. Neste caso
a matriz de Vandermonde é mal-condicionada (ver [1]), acarretando um aumento dos erros de arredondamento
na solução do sistema.
Uma maneira de resolver este problema é escrever o polinômio em uma base que produza um sistema mais
bem-condicionado.
2
3 Diferenças divididas de Newton
O método das diferenças divididas de Newton consistem em contruir o polinômio interpolador da seguinte
forma:
P (x) = a0 + a1(x− x0) + a2(x− x0)(x− x1) + · · ·+ an(x− x0)(x− x1) · · · (x− xn−1)
Assim, o problema de calcular os coe�cientes a0, a1, . . . , an é equivalente ao seguinte sistema linear:
a0 = y0
a0 + a1(x1 − x0) = y1
a0 + a1(x2 − x0) + a2(x2 − x0)(x2 − x1) = y2
.
.
.
a0 + a1(xn − x0) + a2(xn − x0)(xn − x1) + · · ·+ an(xn − x0)(xn − x1) · · · (xn − xn−1) = yn
Equivalente à sua forma matricial:
1 0 0 · · · 0
1 (x1 − x0) 0 · · · 0
1 (x2 − x0) (x2 − x0)(x2 − x1) · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 (xn − x0) (xn − x0)(xn − x1) · · · (xn − x0)(xn − x1) · · · (xn − xn−1)


a0
a1
a2
.
.
.
an
 =

y0
y1
y2
.
.
.
yn

Este é um sistema triangular inferior que pode ser facilmente resolvido conforme:
a0 = y0
a1 =
y1 − a0
x1 − x0 =
y1 − y0
x1 − x0
a2 =
y2 − a1(x2 − x0)− a0
(x2 − x0)(x2 − x1) =
y2−y1
(x2−x1) −
y1−y0
(x1−x0)
(x2 − x0)
. . .
A solução deste sistema pode ser escrita em termos das Diferenças Divididas de Newton, de�nidas recur-
sivamente conforme:
f [xj] = yj
f [xj, xj+1] =
f [xj+1]− f [xj]
xj+1 − xj
f [xj, xj+1, xj+2] =
f [xj+1, xj+2]− f [xj, xj+1]
xj+2 − xj
.
.
.
Nesta notação, temos ak = f [x0, x1, x2, . . . , xk]
3
Podemos esquematizar o método na seguinte tabela:
j xj f [xj] f [xj−1, xj] f [xj−2, xj−1, xj]
0 x0 f [x0]
f [x0, x1] =
f [x1]−f [x0]
x1−x0
1 x1 f [x1] f [x0, x1, x2] =
f [x1,x2]−f [x0,x1]
x2−x0
f [x1, x2] =
f [x2]−f [x1]
x2−x1
2 x2 f [x2]
Exemplo 5. Encontrar o polinômio que passe pelos seguintes pontos
(−1, 3), (0, 1), (1, 3), (3, 43)
j xj f [xj] f [xj−1, xj] f [xj−2, xj−1, xj] f [xj−3, xj−2, xj−1, xj]
0 −1 3
1−3
0−(−1) = −2
1 0 1 2−(−2)
1−(−1) = 2
3−1
1−0 = 2
6−2
3−(−1) = 1
2 1 3 20−2
3−0 = 6
43−3
3−1 = 20
3 3 43
Portanto
P (x) = 3− 2(x+ 1) + 2(x+ 1)x+ (x+ 1)x(x− 1)
= x3 + 2x2 − x+ 1
Problema 1. Considere o seguinte conjunto de pontos:
(−2,−47), (0,−3), (1, 4)(2, 41)
Encontre o polinômio interpolador usando os métodos vistos. Trace os pontos no Scilab usando o comando
'plot2d' e trace o grá�co do polinômio usando comandos de plotagem e a estrutura de polinômio.
Resp: 5x3 + 2x− 3
4
4 Polinômios de Lagrange
Outra maneira clássica de resolver o problema da interpolação polinomial é através do polinômios de Lagrange.
Dado um conjunto de pontos {xj}nj=1 distintos dois a dois, de�nimos os polinômios de Lagrange como os
polinômios de grau n− 1 que satisfazem as seguintes condições:
Lk(xj) =
{
1, k = j
0, k 6= j
Assim, a solução do problema de encontrar os polinômios de grau n− 1 tais P (xj) = yj, j = 1, · · · , n é dado
por
P (x) = y1L1(x) + y2L2(x) + · · ·+ ynLn(x) =
n∑
j=1
yjLj(x)
Para construir os polinômios de Lagrange, basta olhar para sua forma fatorada, ou seja:
Lk(x) = Ck
∏
1≤j 6=k≤n
(x− xj)
onde o coe�ciente Ck é obtido da condição Lk(xk) = 1:
Lk(xk) = Ck
∏
1≤j 6=k≤n
(xk − xj) =⇒ Ck = 1∏
1≤j 6=k≤n(xk − xj)
Portanto,
Lk(x) =
∏
1≤j 6=k≤n
(x− xj)
(xk − xj)
Observação 1. O problema de interpolação quando escrito usando como base os polinômios de Lagrange
produz um sistema linear diagonal.
Exemplo 6. Encontre o polinômio da forma P (x) = a0 + a1x+ a2x
2 + a3x
3
que passa pelos pontos
(0, 0), (1, 1), (2, 4), (3, 9)
Escrevemos:
L1(x) =
(x− 1)(x− 2)(x− 3)
(0− 1)(0− 2)(0− 3) = −
1
6
x3 + x2 − 11
6
x+ 1
L2(x) =
x(x− 2)(x− 3)
1(1− 2)(1− 3) =
1
2
x3 − 5
2
x2 + 3x
L3(x) =
x(x− 1)(x− 3)
2(2− 1)(2− 3) = −
1
2
x3 + 2x2 − 3
2
x
L4(x) =
x(x− 1)(x− 2)
3(3− 1)(3− 2) =
1
6
x3 − 1
2
x2 +
1
3
x
Assim temos:
P (x) = 0 · L1(x) + 1 · L2(x) + 4 · L3(x) + 9 · L4(x) = x2
Exemplo 7. Encontre o polinômio da forma P (x) = a0 + a1x+ a2x
2 + a3x
3
que passa pelos pontos
(0, 0), (1, 1), (2, 0), (3, 1)
Comos as absissas são as mesmas do exemploanterior, podemos utilizar os mesmos polinômios de Lagrange,
assim temos:
P (x) = 0 · L1(x) + 1 · L2(x) + 0 · L3(x) + 1 · L4(x) = 2
3
x3 − 3x2 + 10
3
x
5
5 Aproximação de funções reais por polinômios interpoladores
Teorema 2. Dados n+1 pontos distintos, x0, x1, · · · , xn, dentro de um intervalo [a, b] e uma função f com
n+1 derivadas contínuas nesse intervalo (f ∈ Cn+1[a, b]), então para cada x em [a, b], existe um número ξ(x)
em (a, b) tal que
f(x) = P (x) +
f (n+1)(ξ(x))
(n+ 1)!
(x− x0)(x− x1) · · · (x− xn),
onde P (x) é o polinômio interpolador. Em especial, pode-se dizer que
|f(x)− P (x)| ≤ M
(n+ 1)!
|(x− x0)(x− x1) · · · (x− xn)| ,
onde
M = max
x∈[a,b]
|f (n+1)(ξ(x))|
Exemplo 8. Considere a função f(x) = cos(x) e o polinômio P (x) de grau 2 tal que P (0) = cos(0) = 1,
P (1
2
) = cos(1
2
) e P (1) = cos(1). Use a fórmula de Lagrange para encontrar P (x). Encontre o erro máximo
que se assume ao aproximar o valor de cos(x) pelo de P (x) no intervalo [0, 1]. Trace os grá�cos de f(x) e P (x)
no intervalo [0, 1] no mesmo plano cartesiano e, depois, trace o grá�co da diferença cos(x)− P (x). Encontre
o erro efetivo máximo | cos(x)− P (x)|.
P (x) = 1
(x− 1
2
)(x− 1)
(0− 1
2
)(0− 1) + cos
(
1
2
)
(x− 0)(x− 1)
(1
2
− 0)(1
2
− 1) + cos(1)
(x− 0)(x− 1
2
)
(1− 0)(1− 1
2
)
≈ 1− 0.0299720583066x− 0.4297256358252x2
L1=poly([.5 1],'x');L1=L1/horner(L1,0)
L2=poly([0 1],'x');L2=L2/horner(L2,0.5)
L3=poly([0 .5],'x');L3=L3/horner(L3,1)
P=L1+cos(.5)*L2+cos(1)*L3
x=[0:.05:1]
plot(x,cos)
plot(x,horner(P,x),'red')
plot(x,horner(P,x)-cos(x))
Para encontrar o erro máximo, precisamos estimar |f ′′′(x)| = | sin(x)| ≤ sin(1) < 0.85 e
max
x∈[0,1]
∣∣∣∣x(x− 12
)
(x− 1)
∣∣∣∣
O polinômio de grau três Q(x) = x
(
x− 1
2
)
(x − 1) tem um mínimo (negativo) em x1 = 3+
√
3
6
e um máximo
(positivo) em x2 =
3−√3
6
. Logo
max
x∈[0,1]
∣∣∣∣x(x− 12
)
(x− 1)
∣∣∣∣ ≤ max{|Q(x1)|, |Q(x2)|} ≈ 0.0481125.
Portanto,
|f(x)− P (x)| < 0.85
3!
0.0481125 ≈ 0.0068159 < 7 · 10−3
Para encontrar o erro efetivo máximo, basta encontrar o máximo de |P (x)− cos(x)|. O mínimo (negativo)
de P (x)− cos(x) acontece em x1 = 4.29 · 10−3 e o máximo (positivo) acontece em x2 = 3.29 · 10−3. Portanto,
o erro máximo efetivo é 4.29 · 10−3.
6
Exemplo 9. Considere o problema de aproximar o valor da integral
∫ 1
0
f(x)dx pelo valor da integral do
polinômio P (x) que coincide com f(x) nos pontos x0 = 0, x1 =
1
2
e x2 = 1. Use a fórmula de Lagrange para
encontrar P (x). Obtenha o valor de
∫ 1
0
f(x)dx e encontre uma expressão para o erro de truncamento.
O polinômio interpolador de f(x) é
P (x) = f(0)
(x− 1
2
)(x− 1)
(0− 1
2
)(0− 1) + f
(
1
2
)
(x− 0)(x− 1)
(1
2
− 0)(1
2
− 1) + f(1)
(x− 0)(x− 1
2
)
(1− 0)(1− 1
2
)
= f(0)(2x2 − 3x+ 1) + f
(
1
2
)
(−4x2 + 4x) + f(1)(2x2 − x)
e a integral de P (x) é∫ 1
0
P (x)dx =
[
f(0)
(
2
3
x3 − 3
2
x2 + x
)
+ f
(
1
2
)(
−4
3
x3 + 2x2
)
+ f(1)
(
2
3
x3 − 1
2
x2
)]1
0
= f(0)
(
2
3
− 3
2
+ 1
)
+ f
(
1
2
)(
−4
3
+ 2
)
+ f(1)
(
2
3
− 1
2
)
=
1
6
f(0) +
2
3
f
(
1
2
)
+
1
6
f(1)
Para fazer a estimativa de erro usando o teorema (2), e temos∣∣∣∣∫ 1
0
f(x)dx−
∫ 1
0
P (x)dx
∣∣∣∣ = ∣∣∣∣∫ 1
0
f(x)− P (x)dx
∣∣∣∣
≤
∫ 1
0
|f(x)− P (x)|dx
≤ M
6
∫ 1
0
∣∣∣∣x(x− 12
)
(x− 1)
∣∣∣∣ dx
=
M
6
[∫ 1/2
0
x
(
x− 1
2
)
(x− 1)dx−
∫ 1
1/2
x
(
x− 1
2
)
(x− 1)dx
]
=
M
6
[
1
64
−
(
− 1
64
)]
=
M
192
.
Lembramos que M = maxx∈[0,1] |f ′′′(x)|.
Observação 2. Existem estimativas melhores para o erro de truncamento para este esquema de integração
numérica. Veremos com mais detalhes tais esquemas na teoria de integração numérica.
Problema 2. Use o resultado do exemplo anterior para aproximar o valor das seguintes integrais:
a)
∫ 1
0
ln(x+ 1)dx
b)
∫ 1
0
e−x
2
dx
Usando a fórmula obtida, temos que ∫ 1
0
ln(x+ 1)dx ≈ 0.39± 1
96∫ 1
0
e−x
2
dx ≈ 0.75± 3.87
192
7
Problema 3. Use as mesmas técnicas usadas o resultado do exemplo (9) para obter uma aproximação do
valor de ∫ 1
0
f(x)dx
através do polinômio interpolador que coincide com f(x) nos pontos x = 0 e x = 1.
Resp:
∫ 1
0
P (x)dx = f(0)+f(1)
2
,
1
12
maxx∈[0,1] |f ′′(x)|
6 Ajuste de curvas pelo método dos mínimos quadrados
No problema de interpolação, desejamos encontrar uma função f(x) tal que
f(xj) = yj
para um conjunto de pontos dados.
Existem diversas situações em que desejamos encontrar uma função que se aproxime desses pontos.
0 10.2 0.4 0.6 0.80.1 0.3 0.5 0.7 0.90.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95
1
0.2
0.4
0.6
0.8
1.2
0.1
0.3
0.5
0.7
0.9
1.1
Figura 1: Conjunto de 15 pontos e a reta que melhor se ajuste a eles pelo critério do mínimos quadrados.
No problema de ajuste de curvas, busca-se a função f(x) de família de funções dadas que melhor se
aproxima de um conjunto de pontos dados. O critério mais usado para o ajuste é critério dos mínimos
quadrados, ou seja, buscamos a função f(x) da família que minimiza a soma dos erros elevados ao quadrado:
Eq = [f(x1)− y1]2 + [f(x2)− y2]2 + · · ·+ [f(xn)− yn]2 =
n∑
j=1
[f(xj)− yj]2
Exemplo 10. Encontre a função do tipo f(x) = ax que melhor se aproxima dos seguintes pontos:
(0,−.1), (1, 2), (2, 3.7) e (3, 7).
De�na
Eq = [f(x1)− y1]2 + [f(x2)− y2]2 + [f(x3)− y3]2 + [f(x4)− y4]2
temos que
Eq = [f(0) + .1]
2 + [f(1)− 2]2 + [f(2)− 3.7]2 + [f(3)− 7]2
= [.1]2 + [a− 2]2 + [2a− 3.7]2 + [3a− 7]2
8
Devemos encontrar o parâmetro a que minima o erro, portanto, calculamos:
∂Eq
∂a
= 2[a− 2] + 4[2a− 3.7] + 6[3a− 7] = 28a− 60.8
Portanto o valor de a que minimiza o erro é a = 60.8
28
.
x=[0 1 2 3]'
y=[-.1 2 3.7 7]'
plot2d(x,y,style=-4)
Problema 4. Encontre a função do tipo f(x) = bx+a que melhor aproxima os pontos do problema anterior.
Resp: f(x) = −.3 + 2.3x
Eq = [f(0) + .1]
2 + [f(1)− 2]2 + [f(2)− 3.7]2 + [f(3)− 7]2
= [a+ .1]2 + [a+ b− 2]2 + [a+ 2b− 3.7]2 + [a+ 3b− 7]2
Devemos encontrar os parâmetros a b que minimizam o erro, por isso, calculamos as derivadas parciais:
∂Eq
∂a
= 2[a+ .1] + 2[a+ b− 2] + 2[a+ 2b− 3.7] + 2[a+ 3b− 7]
∂Eq
∂b
= 2[a+ b− 2] + 4[a+ 2b− 3.7] + 6[a+ 3b− 7]
O erro mínimo acontece quando as derivadas são nulas, ou seja:
8a+ 12b = 25.2
12a+ 28b = 60.8
Cuja solução é dada por a = −0.3 e b = 2.3. Portanto a função que procuramos é f(x) = −.3 + 2.3x.
7 O caso linear
7.1 Revisão de Álgebra Linear - O método dos mínimos quadrados para proble-
mas lineares impossíveis
Considere o sistema linear dado por Ax = b onde A é uma matriz n×m e b é um vetor de n linhas. Assumimos
as seguintes hipóteses:
• n ≥ m. O número de linhas é igual ou superior ao número de colunas. (Mais equações que incógnitas)
• O posto de A é m, i.e., existem m linhas L.I. Isso implica que Av = 0 apenas quando v = 0
Neste caso, não seremos necessariamente capazes de encontrar um vetor x que satisfaça exatamente a
equação Ax = b, pelo que estamos interessamos no problema de encontrar o vetor x (ordem m) que minimiza
o erro quadrático dado por:
E :=
n∑
i=1
[zi − bi]2 (1)
onde z = Ax e zi é linha i do vetor z, dado por:
zi = (Ax)i =
m∑
j=1
aijxj, i = 1, · · · , n (2)
9
onde aij é o elemento de A na linha i e coluna j. Substituindo (2) em (1)
E :=
n∑
i=1
[
m∑
j=1
aijxj − bi
]2
(3)
Esta é uma função diferenciável nos coe�cientes xj e portanto todo ponto de mínimo acontece quando∇E = 0,
ou seja, quando
∂
∂xl
E = 0,∀1 ≤ l ≤ m
O que implica a seguinte condição0 =
∂
∂xl
E =
n∑
i=1
2
[
m∑
j=1
aijxj − bi
]
ail, l = 1, · · · ,m
Equivalente a
n∑
i=1
m∑
j=1
ailxjaij =
n∑
i=1
ailbi, l = 1, · · · ,m
que pode ser reescrito na forma vetorial como:
∑n
i=1
∑m
j=1 ai1xjaij∑n
i=1
∑m
j=1 ai2xjaij
.
.
.∑n
i=1
∑m
j=1 aimxjaij
 =

∑n
i=1 ai1bi∑n
i=1 ai2bi
.
.
.∑n
i=1 aimbi
 (4)
Observamos agora que a expressão (4) é equivalente ao seguinte problema matricial:
ATAx = AT b (5)
Teorema 3. A matriz M = ATA é quadrada de ordem m e é invertível sempre que o posto da matriz A é
igual a número de colunas m.
Demonstração. Para provar que M é invertível precisamos mostrar que Mv = 0 implica v = 0:
Mv = 0 =⇒ ATAv = 0
tomando o produto interno da expressão 0 = ATAv com v, temos:
0 =
〈
ATAv, v
〉
= 〈Av,Av〉 = ‖Av‖2
Então se Mv = 0 Av = 0, como o posto de A é igual ao número de colunas, v = 0.
Outra propriedade importante é que M é simétrica, ou seja, M = MT . Isso é facilmente provado pelo
seguinte argumento:
MT = (ATA)T = (A)T (AT )T = ATA = M
10
7.2 Ajuste linear de curvas pelo método dos mínimos quadrados
Seja f1(x), f2(x), . . . , fm(x) um conjunto dem funções e (x1, y1), (x2, y2), . . . , (xn, yn) um conjunto de n pontos.
Procuram-se os coe�cientes a1, a2, . . . , am tais que a função dada por
f(x) = a1f1(x) + a2f2(x) + . . .+ amfm(x)
minimiza o erro dado por
Eq =
n∑
i=1
[f(xi)− yi]2
como f(x) =
∑m
j=1 ajfj(x), temos
Eq =
n∑
i=1
[
m∑
j=1
ajfj(xi)− yi
]2
Este problema é equivalente a resolver pelo métodos dos mínimos quadrados o seguinte sistema linear:
f1(x1) f2(x1) · · · fm(x1)
f1(x2) f2(x2) · · · fm(x2)
f1(x3) f2(x3) · · · fm(x3)
.
.
.
.
.
.
.
.
.
.
.
.
f1(xn) f2(xn) · · · fm(xn)


a1
a2
.
.
.
am
 =

y1
y2
y3
.
.
.
yn

Exemplo 11. Encontrar a reta que melhor aproxima o seguinte conjunto de dados:
xi yi
0.01 1.99
1.02 4.55
2.04 7.2
2.95 9.51
3.55 10.82
Desejamos então encontrar os valores de a e b tais que a função f(x) = ax+ b melhor se ajusta aos pontos da
tabela. A�m de usar o critério dos mínimos quadrados, escrevemos o problema na forma matricial dada por:
0.01 1
1.02 1
2.04 1
2.95 1
3.55 1

[
a
b
]
=

1.99
4.55
7.2
9.51
10.82

Multiplicamos agora ambos os lados pela transposta
[
0.01 1.02 2.04 2.95 3.55
1 1 1 1 1
]
:
[
0.01 1.02 2.04 2.95 3.55
1 1 1 1 1
]
0.01 1
1.02 1
2.04 1
2.95 1
3.55 1

[
a
b
]
=
[
0.01 1.02 2.04 2.95 3.55
1 1 1 1 1
]
1.99
4.55
7.2
9.51
10.82

11
[
26.5071 9.57
9.57 5
] [
a
b
]
=
[
85.8144
34.07
]
A solução desse sistema é a = 2.5157653 e b = 1.9988251
A tabela abaixo mostra os valores dados e os valores ajustados:
xi yi axi + b axi + b− yi
0.01 1.99 2.0239828 0.0339828
1.02 4.55 4.5649057 0.0149057
2.04 7.2 7.1309863 - 0.0690137
2.95 9.51 9.4203327 - 0.0896673
3.55 10.82 10.929792 0.1097919
Problema 5. Encontrar a parábola y = ax2 + bx+ c que melhor aproxima o seguinte conjunto de dados:
xi yi
0.01 1.99
1.02 4.55
2.04 7.2
2.95 9.51
3.55 10.82
e complete a tabela
xi yi ax
2
i + bxi + c ax
2
i + bxi + c− yi
0.01 1.99
1.02 4.55
2.04 7.2
2.95 9.51
3.55 10.82
Resposta y = −0.0407898x2 + 2.6613293x+ 1.9364598
xi yi ax
2
i + bxi + c ax
2
i + bxi + c− yi
0.01 1.99 1.963069 - 0.0269310
1.02 4.55 4.6085779 0.0585779
2.04 7.2 7.1958206 - 0.0041794
2.95 9.51 9.4324077 - 0.0775923
3.55 10.82 10.870125 0.0501249
Problema 6. Dado o seguinte conjunto de dados
xi yi
0. 31
0.1 35
0.2 37
0.3 33
0.4 28
0.5 20
0.6 16
0.7 15
0.8 18
0.9 23
1. 31
12
• Encontre a função do tipo f(x) = a+ b sin(2pix) + c cos(2pix) que melhor aproxima os valores dados.
• Encontre a função do tipo f(x) = a+ bx+ cx2 + dx3 que melhor aproxima os valores dados.
Resp: a = 25.638625, b = 9.8591874, c = 4.9751219 e a = 31.475524, b = 65.691531, c = −272.84382,
d = 208.23621
8 Problemas não lineares que podem ser aproximados por proble-
mas lineares
Eventualmente, problemas de ajuste de curvas podem recair num sistema não linear. Por exemplo, se dese-
jamos ajustar a função y = Aebx ao conjunto de pontos (x0, y0), (x1, y1) e (x2, y2), temos que minimizar o
funcional
Eq = (Ae
x0b − y0)2 + (Aex1b − y1)2 + (Aex2b − y2)2
ou seja, resolver o sistema
∂Eq
∂A
= 2(Aex0b − y0)ex0b + 2(Aex1b − y1)ex1b + 2(Aex2b − y2)ex2b = 0
∂Eq
∂b
= 2Ax0(Ae
x0b − y0)ex0b + 2Ax1(Aex1b − y1)ex1b + 2x2A(Aex2b − y2)ex2b = 0
que é não linear em A e b. Esse sistema pode ser resolvido pelo método de Newton-Raphson, o que pode se
tornar custoso, ou mesmo inviável quando não dispomos de uma boa aproximação da solução para inicializar
o método.
Felizmente, algumas famílias de curvas admitem uma transformação que nos leva a um problema linear.
No caso da curva y = Aebx, observe que ln y = lnA+ bx. Assim, em vez de ajustar a curva original y = Aebx
a tabela de pontos, ajustamos a curva submetida a transformação logaritmica
z = lnA+ bx := B + bx.
Usamos os três pontos (x0, ln y0) := (x0, y˜0), (x1, ln y1) := (x1, y˜1) e (x2, ln y2) := (x2, y˜2) e resolvemos o
sistema linear
ATA
[
B
b
]
= AT
 y˜0y˜1
y˜2
 ,
onde
A =
 1 x01 x1
1 x2

Exemplo 12. Encontre uma curva da forma y = Aex que melhor ajusta os pontos (1, 2), (2, 3) e (3, 5).
Temos
A =
 1 11 2
1 3

e a solução do sistema leva em B = 0.217442 e b = 0.458145. Portanto, A = e0.217442 = 1.24289.
Observação 3. Os coe�cientes obtidos a partir dessa linearização são aproximados, ou seja, são diferen-
tes daqueles obtidos quando aplicamos mínimos quadrados não linear. Observe que estamos minimizando∑
i
[ln yi − ln(f(xi))]2 em vez de
∑
i
[yi − f(xi)]2. No exemplo resolvido, a solução do sistema não linear
original seria A = 1.19789 e B = 0.474348
13
Observação 4. Mesmo quando se deseja resolver o sistema não linear, a solução do problema linearizado
pode ser usada para construir condições iniciais.
A próxima tabela apresenta algumas curvas e transformações que linearizam o problema de ajuste.
curva transformação problema linearizado
y = aebx Y = ln y Y = ln a+ bx
y = axb Y = ln y Y = ln a+ b lnx
y = axbecx Y = ln y Y = ln a+ b lnx+ cx
y = ae(b+cx)
2
Y = ln y Y = ln a+ b2 + bcx+ c2x2
y = a
b+x
Y = 1
y
Y = b
a
+ 1
a
x
y = A cos(ωx+ φ)
ω conhecido
− y = a cos(ωx)− b sin(ωx),
a = A cos(φ), b = A sin(φ)
Exemplo 13. Encontre a função f da forma y = f(x) = A cos(2pix+ φ) que ajusta a tabela de pontos
xi yi
0 9.12
0.1 1.42
0.2 - 7.76
0.3 - 11.13
0.4 - 11.6
0.5 - 6.44
0.6 1.41
0.7 11.01
0.8 14.73
0.9 13.22
1.0 9.93
Usando o fato que y = A cos(2pix + φ) = a cos(2pix) − b sin(2pix), onde a = A cos(φ) e b = A sin(φ),
z = [ a b ]T é solução do problema
BTBz = BTy,
onde
B =

cos(2pix0) − sin(2pix0)
cos(2pix1) − sin(2pix1)
.
.
.
cos(2pix10) − sin(2pix10)
 =

1. 0.
0.8090170 −0.5877853
0.3090170 −0.9510565
−0.3090170 −0.9510565
−0.8090170 −0.5877853
−1. 0
−0.8090170 0.5877853
−0.3090170 0.9510565
0.3090170 0.9510565
0.8090170 0.5877853
1. 0

.
Assim, a = 7.9614704 e b = 11.405721 e obtemos o seguinte sistema:{
A cos(φ) = 7.9614704
A sin(φ) = 11.405721
.
Observe que
A2 = 7.96147042 + 11.4057212
14
e, escolhendo A > 0, A = 13.909546 e
sin(φ) =
11.405721
13.909546
= 0.8199923
Assim, como cosφ também é positivo, φ é um ângulo do primeiro quadrante:
φ = 0.9613976
Portanto f(x) = 13.909546 cos(2pix+ 0.9613976). Observe que nesse exemplo a soluçãodo problema linear é
a mesma do problema não linear.
Problema 7. Encontre a função f da forma y = f(x) = a
b+x
que ajusta a tabela de pontos
xi yi
0 101
0.2 85
0.4 75
0.6 66
0.8 60
1.0 55
usando uma das transformações tabeladas.
Usando o fato que Y = 1
y
= b
a
+ 1
a
x, z = [ b
a
1
a
]T é solução do problema
ATAz = ATY,
onde
A =

1 x1
1 x2
1 x3
1 x4
1 x5
1 x6
 =

1 0
1 0.2
1 0.4
1 0.6
1 0.8
1 1.0

e
Y =

1/y1
1/y2
1/y3
1/y4
1/y5
1/y6
 =

0.0099010
0.0117647
0.0133333
0.0151515
0.0166667
0.0181818

Assim,
1
a
= 0.0082755 e b
a
= 0.0100288 e, então, a = 120.83924 e b = 1.2118696, ou seja, f(x) = 120.83924
1.2118696+x
.
9 Interpolação linear segmentada
Considere o conjunto (xi, yi)
n
j=1 de n pontos. Assumiremos que xi+1 > xi, ou seja, as abscissas são distintas
e estão em ordem crescente. A função linear que interpola os pontos xi e xi+1 no intervalo i é dada por
Pi(x) = yi
(xi+1 − x)
(xi+1 − xi) + yi+1
(x− xi)
(xi+1 − xi)
O resultado da interpolação linear segmentada é a seguinte função contínua de�nida por partes no intervalo
[x1, xn]:
f(x) = Pi(x), x ∈ [xi, xi+1]
15
Figura 2: Interpolação linear segmentada dos pontos (0,0), (1,4), (2,3), (3,0), (4,2), (5,0)
10 Interpolação cúbica segmentada - spline
Dado um conjunto de n pontos (xj, yj)
n
j=1 tais que xj+1 > xj, ou seja, as abscissas são distintas e estão em
ordem crescente; um spline cúbico que interpola estes pontos é uma função s(x) com as seguintes propriedades:
i Em cada segmento [xj, xj+1], j = 1, 2, . . . n− 1 s(x) é um polinômio cúbico.
ii para cada ponto, s(xj) = yj, i.e., o spline interpola os pontos dados.
iii s(x) ∈ C2, i.e., é função duas vezes continuamente diferenciável.
Da primeira hipótese, escrevemos
s(x) = sj(x), x ∈ [xj, xj+1], j = 1, . . . , n− 1
com
sj(x) = aj + bj(x− xj) + cj(x− xj)2 + dj(x− xj)3
O problema agora consiste em obter os 4 coe�cientes de cada um desses n− 1 polinômios cúbicos.
Veremos que a simples de�nição de spline produz 4n− 6 equações linearmente independentes:
sj(xj) = yj, j = 1, . . . , n− 1
sj(xj+1) = yj+1, j = 1, . . . , n− 1
s′j(xj+1) = s
′
j+1(xj+1), j = 1, . . . , n− 1
s′′j (xj+1) = s
′′
j+1(xj+1), j = 1, . . . , n− 1
Como
s′j(x) = bj + 2cj(x− xj) + 3dj(x− xj)2 (6)
e
s′′j (x) = 2cj + 6dj(x− xj), (7)
16
temos, para j = 1, . . . , n− 1, as seguintes equações
aj = yj,
aj + bj(xj+1 − xj) + cj(xj+1 − xj)2 + dj(xj+1 − xj)3 = yj+1,
bj + 2cj(xj+1 − xj) + 3dj(xj+1 − xj)2 = bj+1,
cj + 3dj(xj+1 − xj) = cj+1,
Por simplicidade, de�nimos
hj = xj+1 − xj
e temos
aj = yj,
aj + bjhj + cjh
2
j + djh
3
j = yj+1,
bj + 2cjhj + 3djh
2
j = bj+1,
cj + 3djhj = cj+1,
que podem ser escrita da seguinte maneira
aj = yj, (8)
dj =
cj+1 − cj
3hj
, (9)
bj =
yj+1 − yj − cjh2j − cj+1−cj3hj h3j
hj
,
=
3yj+1 − 3yj − 3cjh2j − cj+1h2j + cjh2j
3hj
=
3yj+1 − 3yj − 2cjh2j − cj+1h2j
3hj
(10)
Trocando o índice j por j − 1 na terceira equação (8), j = 2, . . . , n− 1
bj−1 + 2cj−1hj−1 + 3dj−1h2j−1 = bj
e, portanto,
3yj − 3yj−1 − 2cj−1h2j−1 − cjh2j−1
3hj−1
+ 2cj−1hj−1 + cjhj−1 − cj−1hj−1 =
3yj+1 − 3yj − 2cjh2j − cj+1h2j
3hj
.
Fazendo as simpli�cações, obtemos:
cj−1hj−1 + cj(2hj + 2hj−1) + cj+1hj = 3
yj+1 − yj
hj
− 3yj − yj−1
hj−1
. (11)
É costumeiro acrescentar a incógnita cn ao sistema. A incógnita cn não está relacionada a nenhum dos
polinômios interpoladores. Ela é uma construção arti�cial que facilita o cálculo dos coe�cientes do spline.
Portanto, a equação acima pode ser resolvida para j = 2, . . . , n− 1.
Para determinar unicamente os n coe�cientes cn precisamos acrescentar duas equações linearmente in-
dependentes às n − 2 equações dadas por (11). Essas duas equações adicionais de�nem o tipo de spline
usado.
17
10.1 Spline natural
Uma forma de de�nir as duas equações adicionais para completar o sistema (11) é impor condições de fronteira
livres (ou naturais), ou seja,
S ′′(x1) = S ′′(xn) = 0. (12)
Substituindo na equação (7)
s′′1(x1) = 2c1 + 6d1(x1 − x1) = 0 =⇒ c1 = 0.
e
s′′n−1(xn) = 2cn−1 + 6dn−1(xn − xn−1) = 0 =⇒ c1 = 0.
Usando o fato que
cn−1 + 3dn−1hn−1 = cn
temos que
cn = −3dn−1(xn − xn−1) + 3dn−1hn−1 = 0.
Essa duas equações para c1 e cn juntamente com as equações (11) formam um sistema de n equações AX = B,
onde
A =

1 0 0 0 · · · 0 0
h1 2h2 + 2h1 h2 0 · · · 0 0
0 h2 2h3 + 2h2 h3 · · · 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · hn−2 2hn−2 + 2hn−1 hn−1
0 0 0 · · · 0 0 1

X =

c1
c2
.
.
.
cn
 e B =

0
3y3−y2
h2
− 3y2−y1
h1
3y4−y3
h3
− 3y3−y2
h2
.
.
.
3yn−1−yn−2
hn−2
− 3yn−2−yn−3
hn−3
0

Observe que a matriz A é diagonal dominante estrita e, portanto, o sistema AX = B possui solução única.
Os valores dos an, bn e dn são obtidos diretamente pelas expressões (8), (10) e (9), respectivamente.
Exemplo 14. Construa um spline cúbico natural que passe pelos pontos (2, 4.5), (5,−1.9), (9, 0.5) e (12,−0.5).
O spline desejado é uma função de�nida por partes da forma:
f(x) =

a1 + b1(x− 2) + c1(x− 2)2 + d1(x− 2)3 2 ≤ x < 5
a2 + b2(x− 5) + c2(x− 5)2 + d2(x− 5)3 5 ≤ x < 9
a3 + b3(x− 9) + c3(x− 9)2 + d3(x− 9)3 9 ≤ x < 12
.
Os coe�cientes c1, c2 e c3 resolvem o sistema AX = B, onde
A =

1 0 0 0
3 2 · 3 + 2 · 4 4 0
0 4 2 · 4 + 2 · 3 3
0 0 0 1
 =

1 0 0 0
3 14 4 0
0 4 14 3
0 0 0 1

18
X =

c1
c2
c3
c4
 e B =

0
30.5−(−1.9)
4
− 3 (−1.9)−4.5
3
3−0.5−0.5
3
− 30.5−(−1.9)
4
0
 =

0
8.2
−2.8
0

Observe que c4 é um coe�ciente arti�cial para o problema. A solução é c1 = 0, c2 = 0.6847826, c3 =
−0.3467391 e c4 = 0. Calculamos os demais coe�cientes usando as expressões (8), (10) e (9):
a1 = y1 = 4.5
a2 = y2 = −1.9
a3 = y3 = 0.5
d1 =
c2 − c1
3h1
=
0.6847826− 0
3 · 3 = 0.0760870
d2 =
c3 − c2
3h2
=
−0.3467391− 0.6847826
3 · 4 = −0.0859601
d3 =
c4 − c3
3h3
=
0 + 0.3467391
3 · 3 = 0.0385266
b1 =
y2 − y1
h1
− h1
3
(2c1 + c2) =
−1.9− 4.5
3
− 3
3
(2 · 0− 0.6847826) = −2.8181159
b2 =
y3 − y2
h2
− h2
3
(2c2 + c3) =
0.5− (−1.9)
4
− 4
3
(2 · 0.6847826− 0.3467391) = −0.7637681
b3 =
y4 − y3
h3
− h3
3
(2c3 + c4) =
−0.5− 0.5
3
− 3
3
(2 · (−0.3467391) + 0) = 0.3601449
Portanto,
f(x) =

4.5− 2.8181159(x− 2) + 0.0760870(x− 2)3 2 ≤ x < 5
−1.9− 0.7637681(x− 5) + 0.6847826(x− 5)2 − 0.0859601(x− 5)3 5 ≤ x < 9
0.5 + 0.3601449(x− 9)− 0.3467391(x− 9)2 + 0.0385266(x− 9)3 9 ≤ x < 12
.
Comandos no scilab:
x=[2 5 9 12]'
y=[4.5 -1.9 0.5 -0.5]
h=x(2:4)-x(1:3)
A=[1 0 0 0;h(1) 2*h(1)+2*h(2) h(2) 0;0 h(1) 2*h(1)+2*h(2) h(2);0 0 0 1 ]
B=[0 3*(y(3)-y(2))/h(2)-3*(y(2)-y(1))/h(1) 3*(y(4)-y(3))/h(3)-3*(y(3)-y(2))/h(2) 0]'
c=A\B
for i=1:3
a(i)=y(i)
d(i)=(c(i+1)-c(i))/(3*h(i))
b(i)=(y(i+1)-y(i))/h(i)-h(i)/3*(2*c(i)+c(i+1))
end
z=zeros(3,1000)
for i=1:3
P(i)=poly([a(i) b(i) c(i) d(i)],'x','coeff')
z(i,:)=linspace(x(i),x(i+1),1000)
plot(z(i,:), horner(P(i),z(i,:) -x(i)))
end
19
10.2 Spline com condições de contorno �xadas
Alternativamente, para completar o sistema (11), podemos impor condições de contorno �xadas, ou seja,
S ′(x1) = f ′(x1)
S ′(xn) = f ′(xn).
Substituindo na equação (6)
s′1(x1) = b1 + 2c1(x1 − x1) + 3dj(x1 − x1)2 = f ′(x1) =⇒ b1 = f ′(x1)
e
s′n−1(xn) = bn−1 + 2cn−1(xn − xn−1) + 3dj(xn − xn−1)2 = bn−1 + 2cn−1hn−1 + 3dn−1h2n−1 = f ′(xn)
Usando as equações (9) e (10) para j = 1 e j = n− 1, temos:
2c1h1 + c2h1 = 3
y2 − y1h1
− 3f ′(x1)
cn−1hn−1 + cnhn−1 = 3f ′(xn)− 3yn − yn−1
hn−1
Essas duas equações juntamente com as equações (11) formam um sistema de n equações AX = B, onde
A =

2h1 h1 0 0 · · · 0 0
h1 2h2 + 2h1 h2 0 · · · 0 0
0 h2 2h3 + 2h2 h3 · · · 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · hn−2 2hn−2 + 2hn−1 hn−1
0 0 0 · · · 0 hn−1 2hn−1

X =

c1
c2
.
.
.
cn
 e B =

3y2−y1
h1
− 3f ′(x1)
3y3−y2
h2
− 3y2−y1
h1
3y4−y3
h3
− 3y3−y2
h2
.
.
.
3yn−1−yn−2
hn−2
− 3yn−2−yn−3
hn−3
3f ′(xn)− 3yn−yn−1hn−1

Observe que a matriz A é diagonal dominante estrita e, portanto, o sistema AX = Y possui solução única.
Os valores dos an, bn e dn são obtidos diretamente pelas expressões (8), (10) e (9), respectivamente.
Exemplo 15. Construa um spline cúbico com fronteira �xada que interpola a função y = sin(x) nos x = 0,
x = pi
2
, x = pi, x = 3pi
2
e x = 2pi.
O spline desejado passa pelos pontos (0, 0), (pi/2, 1), (pi, 0), (3pi/2,−1) e (2pi, 0) e tem a forma:
f(x) =

a1 + b1x+ c1x
2 + d1x
3 0 ≤ x < pi
2
a2 + b2(x− pi2 ) + c2(x− pi2 )2 + d2(x− pi2 )3 pi2 ≤ x < pi
a3 + b3(x− pi) + c3(x− pi)2 + d3(x− pi)3 pi ≤ x < 3pi2
a4 + b4(x− 3pi2 ) + c4(x− 3pi2 )2 + d4(x− 3pi2 )3 3pi2 ≤ x < 2pi
.
Observe que ele satisfaz as condição de contorno f ′(0) = cos(0) = 1 e f ′(2pi) = cos(2pi) = 1.
20
Os coe�cientes c1, c2, c3 e c4 resolvem o sistema AX = B, onde
A =

pi pi/2 0 0 0
pi/2 2pi pi/2 0 0
0 pi/2 2pi pi/2 0
0 0 pi/2 2pi pi/2
0 0 0 pi/2 pi

X =

c1
c2
c3
c4
c5
 e B =

31−0
pi/2
− 3 · 1
30−1
pi/2
− 31−0
pi/2
3−1−0
pi/2
− 30−1
pi/2
30−(−1)
pi/2
− 3 (−1)−0
pi/2
3 · 1− 30−(−1)
pi/2
 =

6/pi − 3
−12/pi
0
12/pi
3− 6/pi

Aqui c5 é um coe�ciente arti�cial para o problema. A solução é c1 = −0.0491874, c2 = −0.5956302, c3 = 0,
c4 = 0.5956302 e c5 = 0.0491874. Calculamos os demais coe�cientes usando as expressões (8), (10) e (9):
a1 = y1 = 0
a2 = y2 = 1
a3 = y3 = 0
a4 = y3 = −1
d1 =
c2 − c1
3h1
=
−0.5956302− (−0.0491874)
3 · pi/2 = −0.1159588
d2 =
c3 − c2
3h2
=
0− (−0.5956302)
3 · pi/2 = 0.1263967
d3 =
c4 − c3
3h3
=
0.5956302− 0
3 · pi/2 = 0.1263967
d4 =
c5 − c4
3h4
=
0.0491874− 0.5956302
3 · pi/2 = −0.1159588
b1 =
y2 − y1
h1
− h1
3
(2c1 + c2) =
1− 0
pi/2
− pi/2
3
(2 · (−0.0491874)− 0.5956302) = 1
b2 =
y3 − y2
h2
− h2
3
(2c2 + c3) =
0− 1
pi/2
− pi/2
3
(2 · (−0.5956302) + 0) = −0.0128772
b3 =
y4 − y3
h3
− h3
3
(2c3 + c4) =
−1− 0
pi/2
− pi/2
3
(2 · 0 + 0.5956302) = −0.9484910
b4 =
y5 − y4
h4
− h4
3
(2c4 + c5) =
0− (−1)
pi/2
− pi/2
3
(2 · 0.5956302 + 0.0491874) = −0.0128772
Portanto,
f(x) =

x− 0.0491874x2 − 0.1159588x3 0 ≤ x < pi
2
1 +−0.0128772(x− pi
2
)− 0.5956302(x− pi
2
)2 + 0.1263967(x− pi
2
)3 pi
2
≤ x < pi
−0.9484910(x− pi) + 0.1263967(x− pi)3 pi ≤ x < 3pi
2−1− 0.0128772(x− 3pi
2
) + 0.5956302(x− 3pi
2
)2 − 0.1159588(x− 3pi
2
)3 3pi
2
≤ x < 2pi
.
21
Referências
[1] Gautschi, W., and Inglese, G. Lower bounds for the condition number of vandermonde matrices.
Numerische Mathematik 52, 3 (1987/1988), 241�250.
22