Livro-texto - Unidade II

Prévia do material em texto

207
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Unidade II
5 INTERPOLAÇÃO E AJUSTE DE CURVAS
5.1 Introdução
Muitas observações científicas e cotidianas (indústria, comércio etc.) são baseadas em grandezas 
físicas medidas e tabeladas. Normalmente, tais registros são denominados dados ou pontos experimentais 
e são representados por valores discretos. O conteúdo desse tópico trata das técnicas para ajustar curvas 
a tais dados, com o objetivo de obter estimativas intermediárias. Outro propósito é utilizar os dados 
experimentais para desenvolver ou avaliar funções que os representem a partir de determinação de 
parâmetros relacionados com a curva que mais se aproxime da representação gráfica deles. Ambos os 
procedimentos são conhecidos como ajuste de curvas.
Como observa Chapra (2016), há duas abordagens gerais para o ajuste de curvas, baseadas e distintas 
entre si, com base na grandeza do erro associada com os dados.
A primeira dessas abordagens é utilizada quando os dados exibem um grau significativo de erro ou 
“ruído”; nesse caso, a estratégia adotada é encontrar uma única curva que represente a tendência geral 
dos dados. Como cada ponto pode estar incorreto, a curva não passará obrigatoriamente por todos 
eles; ela deverá seguir o padrão dos pontos considerados como um grupo. Um comportamento dessa 
natureza é chamado de regressão por mínimos quadrados e é utilizado quando há imprecisão nos dados.
O segundo tipo de abordagem é quando se tem a informação de antemão que os dados são muito 
precisos. Nesse caso, busca‑se ajustar uma curva ou uma série de curvas que passem por cada um dos pontos, 
ou com mais frequência utiliza‑se interpolações e extrapolações. A interpolação é um procedimento de 
avaliação dos valores prováveis entre os pontos medidos; a extrapolação estima valores além do intervalo 
no qual foram feitas medições.
Ainda, segundo Chapra (2016), em geral, o procedimento de ajuste de dados experimentais possibilita 
avaliações denominadas análise de tendência e teste de hipótese. A análise de tendência utiliza o padrão 
dos dados para fazer previsões de valores para a variável dependente, tanto por extrapolação (além dos 
limites dos dados observados) ou interpolação (interior do intervalo dos dados). A segunda aplicação 
do ajuste de curvas a dados experimentais é no teste de hipótese. Nesse caso, um exemplo matemático 
já existente é comparado com os dados medidos. Se os coeficientes do modelo forem desconhecidos, 
pode ser necessário determinar valores que melhor ajustem os dados observados. Por outro lado, se 
estimativas para os referidos coeficientes já estiverem disponíveis, pode ser apropriado comparar os 
valores previstos pelo modelo com aqueles observados para testar a adequação do modelo. Em geral, 
são avaliados vários modelos alternativos com base nas observações existentes.
208
Unidade II
5.2 Interpolação
Há situações em que é necessário fazer estimativas de valores intermediários entre dados 
precisos. O método mais comum usado para esse propósito é a interpolação polinomial. Ela consiste 
em, dada uma listagem de pares de dados cuja função f(x) é desconhecida, determinar um valor 
aproximado da função y = f(x) em um dado intervalo [x0, x1], sendo x ∈ [x0, x1]. Essa determinação 
é feita através da aproximação dos dados da tabela por outra função F(x), escolhida entre uma 
classe de funções definidas a priori e que satisfaça algumas propriedades em substituição à função 
original, desconhecida.
 Observação
A determinação da função F(x) é descrita por Chapra (2016) a partir da 
expressão geral para um polinômio de grau n como:
2 3 n
0 1 2 3 ny a a x a x a x ... a x= + + + + +
Para (n + 1) pontos dados, existe um e somente um polinômio de grau n que passa por todos os pontos. 
Por exemplo, uma única reta (isto é, um polinômio de primeiro grau) que liga dois pontos. Analogamente, 
uma única parábola está bem determinada por um conjunto de três pontos. A interpolação polinomial 
consiste em definir o único polinômio de grau n que passa pelos n + 1 pontos dados. Esse polinômio, 
então, fornece uma fórmula para calcular valores intermediários.
No entanto, embora exista um só polinômio de grau n que passa por (n + 1) pontos, há diversas 
formas matemáticas alternativas segundo as quais ele pode ser expresso. O polinômio interpolador por 
diferenças divididas de Newton está entre as fórmulas mais populares e úteis.
5.2.1 Interpolação linear
A forma mais simples de interpolação é ligar dois pontos dados com uma reta. Essa técnica, 
interpolação linear, é descrita graficamente por semelhança de triângulos.
Figura 48
209
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Da figura pode ser retirado o triângulo:
Figura 49
Assim:
0 1 0
0 1 0
F(x) f(x ) f(x ) f(x )
x x x x
− −
=
− −
Reorganizando:
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
− = −
−
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
= + −
−
Essa expressão é a fórmula de interpolação linear. F(x) é um polinômio interpolador de primeiro grau. 
Em geral, quanto menor o intervalo entre os pontos dados, melhor a aproximação, o que se deve ao 
fato de que, conforme o intervalo diminui, uma função contínua será melhor aproximada por uma reta.
Exemplo 28
Observe o gráfico de f(x) ‑ Inx.
a) Faça uma estimativa do logaritmo natural de 2 usando uma interpolação linear, utilizando o 
intervalo [1,6].
b) Repita o procedimento fazendo uso do intervalo [1,4].
210
Unidade II
Figura 50
Resolução:
a) Usando o intervalo [1,6]:
f(x) = Inx
Então:
0x 1
f(1) ln1 0
=
 = =
 e 
1x 6
f(6) ln6 1,7917595
=
 = =
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
= + −
−
( )1,7917595 0F(2) 0 . 2 1
6 1
−
= + −
−
( )1,7917595F(2) . 1
5
=
F(2) 0,3583519=
211
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
a
valor real aproximação atual
100%
valor real
−
ε =
Considerando que o valor correto de ln2 0,6931472= :
a
0,6931472 0,3583519
100% 48,30%
0,6931472
−
ε = =
b) Usando o intervalo [1,4]:
f(x) = Inx
Então:
0x 1
f(1) ln1 0
=
 = = e 
1x 4
f(4) ln4 1,3862944
=
 = =
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
= + −
−
( )1,3862944 0F(2) 0 . 2 1
4 1
−
= + −
−
( )1,3862944F(2) . 1
3
=
F(2) 0,4620981=
a
valor real aproximação atual
100%
valor real
−
ε =
Considerando que o valor correto de ln2 0,6931472= :
a
0,6931472 0,4620981
100% 33,33%
0,6931472
−
ε = =
Logo, o uso do intervalo menor reduz o erro relativo porcentual.
212
Unidade II
5.2.2 Interpolação quadrática
Anteriormente, vimos a interpolação linear; agora, iremos estudar a interpolação quadrática. As 
duas formas de interpolação são casos particulares de uma interpolação denominada interpolação de 
Lagrange, sendo a linear também chamada de interpolação de Lagrange de grau 1 e a interpolação 
quadrática, interpolação de Lagrange de grau 2.
Assim, quando dispomos de um conjunto de dados experimentais que aparece de forma linear, 
conforme vimos anteriormente, utilizamos a interpolação linear; já quando a distribuição dos dados 
(x,y) são formas aproximadas de uma parábola, podemos fazer a interpolação quadrática, ou, como já 
mencionado, interpolação de Lagrange de grau 2.
Para tal, vamos conhecer o caso geral da interpolação de Lagrange. Como caso geral, denominamos 
o interpolador de grau N, sendo que, para tanto, precisamos conhecer (N‑1) pontos experimentais 
(x1, x2, ... xn‑1) e (y1, y2, ... yn‑1) e este polinômio obtido passará pelos (N‑1) pontos.
O polinômio interpolador de Lagrange para um grau n qualquer é apresentado como:
n 0 0 1 1 2 2 n nP (x) f(x ) L (x) f(x ) L (x) f(x ) L (x) ... f(x ) L (x)= ⋅ + ⋅ + ⋅ + + ⋅
Onde:
Lk(x) representam frações com k variando de 0 a n. Essas frações são razões entre os valores possíveis 
de x. Considerando que f(xn) = yn, segue que para um polinômio de grau 1 (n = 1), uma reta, teríamos:
1 0 0 1 1
1
0
0 1
0
1
1 0
P y L y L
x x
L (x)
x x
x x
L (x)
x x
= ⋅ + ⋅
−
=
−
−=
−
Conforme mencionado, só é possível obter um polinômio interpolador de Lagrange quando 
dispomos de um número de pontos maior que o grau desejado. Considerando o de grau 2 anteriormente 
apresentado, teríamos os pontos genéricos:
0 0 1 1 2 2 3 3(x ,y ),(x ,y ),(x ,y ),(x ,y )
A partir do caso de grau 2, nota‑se que as frações Ln são formadas no numerador pelas diferenças 
de (x ‑ xk), com k variando de 0 a n. No entanto, não há este termo quando k é igual a n.
Já o denominador da fração Ln é formado pelas diferenças entre xn e os demais pontos x.
213
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Dessa forma, podemos escrever o caso geral (grau n) do polinômio interpolador de Lagrange:
0 0 1 1 2 2 3 3(x ,y ),(x ,y ),(x ,y ),(x ,y )
n 0 0 1 1 2 2 n nP y L y L y L ... y L= ⋅ + ⋅ + ⋅ + + ⋅
E
0 i 1 i 1 n
i
i 0 i i 1 i i 1 i n
x x x x x x x x
L (x) ... ...
x x x x x x x x
− +
− +
− − − −
= ⋅ ⋅ ⋅ ⋅ ⋅
− − − −
Lembrando que:
i iL (x ) 1= e i j
j i
L (x ) 0
≠
=
Assim, conforme esperado:
n i 0 1 n iP (x ) y 0 ...y 1 y 0 y= ⋅ + ⋅ + ⋅ =
Voltando à interpolação quadrática, ela ocorre quando o gráfico obtido a partir de dados coletados sugere 
que uma estratégia para melhorar a estimativa é introduzir alguma curvatura na curva, ligando os pontos.
Figura 51
Para Chapra (2016), se estiverem disponíveis três pontos, a alternativa é um polinômio de segundo 
grau (cujo gráfico é uma parábola). Considere a expressão:
( ) ( )( )0 1 0 2 0 1F(x) b b x x b x x x x= + − + − −
214
Unidade II
Efetuando os produtos e agrupando os termos:
( )( )0 1 1 0 2 2 0 1F(x) b b x b x b x b x x x= + − + − −
2
0 1 1 0 2 2 1 2 0 2 0 1F(x) b b x b x b x b xx b x x b x x= + − + − − +
2
0 1 0 2 0 1 1 2 0 2 1 2F(x) b b x b x x x(b b x b x ) b x= − + + − − +
2
0 1 0 2 0 1 1 2 0 2 1 2F(x) (b b x b x x ) x(b b x b x ) b x= − + + − − +
Para que a equação anterior tenha a forma de
2
0 1 2F(x) a a x a x= + +
É necessário fazer:
0 0 1 0 2 0 1a (b b x b x x )= − +
1 1 2 0 2 1a (b b x b x )= − −
a2 = b2
Logo, as equações vistas anteriormente são alternativas equivalentes do único polinômio de segundo 
grau ligando os três pontos.
Um procedimento simples pode ser usado para determinar os valores dos coeficientes:
Fazendo:
b0 = f(x0)
Ao obter F(x1) em ( ) ( )( )0 1 0 2 0 1F(x) b b x x b x x x x= + − + − −
( ) ( )( )1 0 1 1 0 2 1 0 1 1F(x ) f(x ) b x x b x x x x= + − + − −
( )1 0 1 1 0F(x ) f(x ) b x x= + −
( )1 0 1 1 0F(x ) f(x ) b x x− = −
( )
1 0
1
1 0
F(x ) f(x )
b
x x
−
=
−
215
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Como F(x1) = f(x1)
Então:
( )
1 0
1
1 0
f(x ) f(x )
b
x x
−
=
−
Substituindo:
( )
0 0
1 0
1
1 0
b f(x )
f(x ) f(x )
b
x x
=
 − = −
Em:
( ) ( )( )2 0 1 2 0 2 2 0 2 1F(x ) b b x x b x x x x= + − + − −
Como F(x2) = f(x2)
( ) ( ) ( )( )
1 0
2 0 2 0 2 2 0 2 1
1 0
f(x ) f(x )
f(x ) f(x ) x x b x x x x
x x
−
= + − + − −
−
Obtém‑se:
( ) ( )
1 02 1
2 1 1 0
2
2 0
f(x ) f(x )f(x ) f(x )
x x x x
b
(x x )
−−
−
− −
=
−
Exemplo 29
Considere a função f(x) = Inx.
Faça uma estimativa do logaritmo natural de 2 utilizando uma interpolação quadrática usando os 
valores a seguir. Ou seja,
0x 1
f(1) ln1 0
=
 = =
 1
x 4
f(4) ln4 1,3862944
=
 = =
 e 2
x 6
f(6) ln6 1,7917595
=
 = =
216
Unidade II
Resolução:
2
0 1 2F(x) a a x a x= + +
0 0 1 0 2 0 1
1 1 2 0 2 1
2 2
a (b b x b x x )
a (b b x b x )
a b
= − +
 = − −
 =
( )
( ) ( )
0 0
1 0
1
1 0
1 02 1
2 1 1 0
2
2 0
b f(x )
f(x ) f(x )
b
x x
f(x ) f(x )f(x ) f(x )
x x x x
b
(x x )



 =

− = −
 −−
−
− − = −
Logo:
( )
( ) ( )
0
1
2
b 0
1,3862944 0
b 0,4620981
4 1
1,7917595 1,3862944 1,3862944 0 0,4054651 1,3862944
6 4 4 1 2 3b
(6 1) 5



=
 − = = −
 − −
− −− − = =
 −
0
1
2
b 0
b 0,4620981
0,2027326 0,4620981
b 0,0518731
5

 =

=
 − = = −

217
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
0 0 1 0 2 0 1
1 1 2 0 2 1
2 2
a (b b x b x x )
a (b b x b x )
a b
= − +
 = − −
 =
0
1
2
a 0 0,4620981.1 ( 0,0518731)1.4 0,4620981 0,2074924
a 0,4620981.1 ( 0,0518731).1 ( 0,0518731).4) 0,4620981 0,0518731 0,2074924
a 0,0518731
= − + − = − −
 = − − − − = + +
 = −
0
1
2
a 0,6695905
a 0,7214636
a 0,0518731
= −
 =
 = −
Assim:
2F(x) 0,6695905 0,7214636x 0,0518731x= − + −
2F(2) 0,6695905 0,7214636.2 0,0518731.2= − + −
F(2) 0,5658443=
a
valor real aproximação atual
100%
valor real
−
ε =
Considerando que o valor correto de ln2 0,6931472= :
a
0,6931472 0,5658443
100% 18,37%
0,6931472
−
ε = =
Logo, a interpolação quadrática reduziu significativamente o erro relativo porcentual.
A partir de manipulação algébrica da teoria anterior é possível estabelecer um dispositivo prático 
para interpolação quadrática.
( ) ( )( )0 1 0 2 0 1F(x) b b x x b x x x x= + − + − −
218
Unidade II
( )
( ) ( )
0 0
1 0
1
1 0
1 02 1
2 1 1 0
2
2 0
b f(x )
f(x ) f(x )
b
x x
f(x ) f(x )f(x ) f(x )
x x x x
b
(x x )



 =

− = −
 −−
−
− − = −
Fazendo:
( ) ( )
( ) ( )
1 02 1
2 1 1 0 1 02 1
2
2 0 2 1 1 0 2 0
f(x ) f(x )f(x ) f(x )
x x x x f(x ) f(x )f(x ) f(x ) 1
b .
(x x ) x x x x (x x )
−−
−
 − − −−
= = − − − − − 
( ) ( )
1 02 1
2
2 1 2 0 1 0 2 0
f(x ) f(x )f(x ) f(x )
b
x x (x x ) x x (x x )
 −−
= − − − − − 
( ) ( )( )0 1 0 2 0 1F(x) b b x x b x x x x= + − + − −
( ) ( ) ( ) ( ) ( )( )
1 0 1 02 1
0 0 0 1
1 0 2 1 2 0 1 0 2 0
f(x ) f(x ) f(x ) f(x )f(x ) f(x )
F(x) f(x ) x x x x x x
x x x x (x x ) x x (x x )
 − −−
= + − + − − − − − − − − 
Efetuando as manipulações algébricas necessárias, obtém‑se:
( )
( )( )
( )
( )( )
0 2 0 11 2
0 1 2
0 1 0 2 1 0 1 2 2 0 2 1
x x (x x ) x x (x x )(x x )(x x )
F(x) f(x ) f(x ) f(x )
(x x )(x x ) x x x x x x x x
− − − −− −
= + +
− − − − − −
Multiplicando e dividindo por termos iguais cada fator:
( )
( )
( ) ( )
( )( )( )
( ) ( )
( )( )( )
1 2 0 0 2 1 0 1 2
0 1 2
0 1 0 2 0 1 0 1 2 1 2 0 2 1 2
(x x )(x x ) x x x x (x x ) x x x x (x x ) x x
F(x) f(x ) f(x ) f(x )
(x x )(x x ) x x x x x x x x x x x x x x
− − − − − − − − −
= + +
− − − − − − − − −
( )
( )
( ) ( )
( )( )( )
( )( )
( )( )( )
1 2 0 1 2 0 1 2 0
0 1 2
0 1 0 2 0 1 0 1 2 1 2 0 2 1 2
(x x )(x x ) x x x x (x x ) x x (x x ) x x x x
F(x) f(x ) f(x ) f(x )
(x x )(x x ) x x x x x x x x x x x x x x
− − − − − − − − −
= + +
− − − − − − − − −
219
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Fazendo:
( )
( )
1 2 0
0 0 1 0 2 0
N (x x )(x x ) x x
D (x x )(x x ) x x
= − − −
= − − −
( )( )( )1 1 0 1 2 1D x x x x x x= − − −
( )( )( )2 2 0 2 1 2D x x x x x x= − − −
0 1 2
0 1 2
N N N
F(x) f(x ) f(x ) f(x )
D D D
= + +
0 1 2
0 1 2
f(x ) f(x ) f(x )
F(x) N
D D D
 
= + +  
Tabela 28
(x ‑ x0) (x0 ‑ x1) (x0 ‑ x2)
(x1 ‑ x0) (x ‑ x1) (x1 ‑ x2)
(x2 ‑ x0) (x2 ‑ x1) (x ‑ x2)
 D0
D1
D2
 N
Exemplo 30
Considere a função f(x) = Inx.
Faça uma estimativa do logaritmo natural de 2 utilizando uma interpolação quadrática usando os 
valores a seguir e o dispositivo prático. Ou seja,
0x 1
f(1) ln1 0
=
 = =
 1
x 4
f(4) ln4 1,3862944
=
 = =
 e 2
x 6
f(6) ln6 1,7917595
=
 = =
( )
( )
1 2 0
0 0 1 0 2 0
N (x x )(x x ) x x
D (x x )(x x ) x x
= − − −
= − − −
( )( )( )1 1 0 1 2 1D x x x x x x= − − −
220
Unidade II
( )( )( )2 2 0 2 1 2D x x x x x x= − − −
( )
( )0
N (2 4)(2 6) 2 1 8
D (1 4)(1 6) 2 1 15
= − − − =
= − − − =
( )( )( )1D 4 1 4 6 2 4 12= − − − =
( )( )( )2D 6 1 6 4 2 6 40= − − − = −
0 1 2
0 1 2
f(x ) f(x ) f(x )
F(x) N
D D D
 
= + +  
0 1,3862944 1,7917595
F(2) 8
15 12 ( 40)
 = + +  −
( )F(2) 8 0,1155245 0,0447940= −
F(2) 0,5658441=
Comparando com o resultado do exemplo anterior, a diferença só acontece na sétima casa decimal.
Como foi visto antes, dado um conjunto de (n+1) pontos distintos (xi, f(xi)), com i = 0, 1, 2, 3, 4, ..., n, 
só é possível determinar um polinômio de grau “n” que satisfaz as condições Pn(xi)= f(xi), ou seja, a 
curva da função polinômio de grau “n” passa por todos os (n+1) pontos.
Generalizando o polinômio (obtido anteriormente para a interpolação quadrática) para um polinômio 
de grau 3, a partir de quatro pontos conhecidos:
( )
( )
( )( )
( )( )( )
( ) ( )
( )( )( )
( ) ( )
( )( )( )
0 2 31 2 3
0 1
0 1 0 2 0 3 1 0 1 2 1 3
0 1 3 0 1 2
2 3
2 0 2 1 2 3 3 0 3 1 3 2
x x x x (x x )x x (x x )(x x )
F(x) f(x ) f(x )
(x x )(x x ) x x x x x x x x
x x (x x ) x x x x (x x ) x x
f(x ) f(x )
x x x x x x x x x x x x
− − −− − −
= + +
− − − − − −
− − − − − −
+ +
− − − − − −
221
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Multiplicando e dividindo por termos iguais cada fator:
( )
( )
( )( )
( )( )( )
( ) ( )
( )( )( )
( ) ( )
( )( )( )
0 2 3 11 2 3 0
0 1
0 1 0 2 0 3 0 1 0 1 2 1 3 1
0 1 3 2 0 1 2 3
2 3
2 0 2 1 2 3 2 3 0 3 1 3 2 3
x x x x (x x )(x x )x x (x x )(x x )(x x )
F(x) f(x ) f(x )
(x x )(x x ) x x (x x ) x x x x x x (x x )
x x (x x ) x x (x x ) x x (x x ) x x (x x )
f(x ) f(x )
x x x x x x (x x ) x x x x x x (x x )
− − − −− − − −
= + +
− − − − − − − −
− − − − − − − −
+ +
− − − − − − − −
( )
( )
( ) ( )
( )( )( )
( ) ( )
( )( )( )
( ) ( )
( )( )( )
0 1 2 30 1 2 3
0 1
0 1 0 2 0 3 0 1 0 1 2 1 3 1
0 1 2 3 0 1 2 3
2 3
2 0 2 1 2 3 2 3 0 3 1 3 2 3
x x (x x ) x x (x x )(x x ) x x (x x )(x x )
F(x) f(x ) f(x )
(x x )(x x ) x x (x x ) x x x x x x (x x )
x x (x x )(x x ) x x x x (x x ) x x (x x )
f(x ) f(x )
x x x x x x (x x ) x x x x x x (x x )
− − − −− − − −
= + +
− − − − − − − −
− − − − − − − −
+ +
− − − − − − − −
Chamando de: ( )( )( )( )0 1 2 3N x x x x x x x x= − − − −
( )( )( )( )0 0 1 0 2 0 3 0D x x x x x x x x= − − − −
( )( )( )( )1 1 0 1 2 1 3 1D x x x x x x x x= − − − −
( )( )( )( )2 2 0 2 1 2 3 2D x x x x x x x x= − − − −
( )( )( )( )3 3 0 3 1 3 2 3D x x x x x x x x= − − − −
0 1 2 3
0 1 2 3
N N N N
F(x) f(x ) f(x ) f(x ) f(x )
D D D D
= + + +
0 31 2
0 1 2 3
f(x ) f(x )f(x ) f(x )
F(x) N
D D D D
 
= + + + 
 
Exemplo 31
Considere a função f(x) = Inx.
Faça uma estimativa do logaritmo natural de 2 usando uma interpolação quadrática usando os valores:
0x 1
f(1) ln1 0
=
 = =
 1
x 4
f(4) ln4 1,3862944
=
 = =
 2
x 5
f(5) ln5 1,6094379
=
 = =
 e 3
x 6
f(6) ln6 1,7917595
=
 = =
222
Unidade II
( )( )( )( )0 1 2 3N x x x x x x x x= − − − −
( )( )( )( )0 0 1 0 2 0 3 0D x x x x x x x x= − − − −
( )( )( )( )1 1 0 1 2 1 3 1D x x x x x x x x= − − − −
( )( )( )( )2 2 0 2 1 2 3 2D x x x x x x x x= − − − −
( )( )( )( )3 3 0 3 1 3 2 3D x x x x x x x x= − − − −
( )( )( )( )N 2 1 2 4 2 5 2 6 1( 2)( 3)( 4) 24= − − − − = − − − = −
( )( )( )( )0D 1 4 1 5 1 6 2 1 ( 3)( 4)( 5).1 60= − − − − = − − − = −
( )( )( )( )1D 4 1 4 5 4 6 2 4 (3)( 1)( 2)( 2) 12= − − − − = − − − = −
( )( )( )( )2D 5 1 5 4 5 6 2 5 (4)(1)( 1)( 3) 12= − − − − = − − =
( )( )( )( )3D 6 1 6 4 6 5 2 6 (5)(2)(1)( 4) 40= − − − − = − = −
0 31 2
0 1 2 3
f(x ) f(x )f(x ) f(x )
F(x) N
D D D D
 
= + + + 
 
0 1,3862944 1,6094379 1,7917595
F(2) ( 24)
( 60) ( 12) 12 ( 40)
 = − + + + − − − 
[ ]F(2) ( 24) 0,1155245 0,1341198 0,0447940= − − + −
[ ]F(2) ( 24) 0,0261987= − −
F(2) 0,6287688=
a
valor real aproximação atual
100%
valor real
−
ε =
223
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Considerando que o valor correto de ln2 0,6931472= :
a
0,6931472 0,6287688
100% 9,29%
0,6931472
−
ε = =
Logo, a interpolação quadrática com quatro pontos reduziu significativamente o erro relativo porcentual.
5.2.2.1 Interpolação polinomial Lagrange – implementação computacional
Esta seção mostrará um exemplo de aplicação do polinômio interpolador de Lagrange. Sabemos que 
esse polinômio é da forma mostrada pela equação a seguir.
( ) ( )i i
i
P x f x . L=∑
Onde 
( )
( )
j
i j i
i j
x x
L
x x
≠
−
=
−
∏ 
Vamos agora ilustrar a aplicação desse método em um problema prático. Imagine que a ocupação 
de um estoque, durante o ano, possa ser representada pelos valores mostrados na tabela a seguir. 
Deseja‑se obter uma função para estimar os valores de estoque para os outros meses do ano, além de 
se determinar qual o mês do ano em que ocorre a maior ocupação do estoque.
Tabela 29
Mês Ocupação do estoque
1 1
6 5
11 3,3
Primeiro devemos encontrar o polinômio de Lagrange para esse problema. Vamos calcular cada 
termo isoladamente.
L =
x x
=
x x +
L =
x x
0
2
1
6 11
1 1 1 11
17 66
1 11
6 1 6
�� � �� �
�� � �� �
�
�� � �� �
�� � �
50
111
12 11
25
1 6
11 1 11 6
7 6
50
2
2
2
� �
�
�
�� � �� �
�� � �� �
�
� �
=
x x +
L =
x x
=
x x +
P x = LL x +L x +L x
P x =
x x +
+
x x +
1 2 3
2 2
1 2 3
17 66
50
1
12 11
� � � � � �
� � �
�
�
��
�
�
��
�
�225
5
7 6
50
3,3
5,7 79,9 24,2
5
2
2
�
�
��
�
�
��
��
�
��
�
�
��
� � � �
× +
x x +
×
P x =
x + x
00
224
Unidade II
L =
x x
=
x x +
L =
x x
0
2
1
6 11
1 1 1 11
17 66
1 11
6 1 6
�� � �� �
�� � �� �
�
�� � �� �
�� � �
50
111
12 11
25
1 6
11 1 11 6
7 6
50
2
2
2
� �
�
�
�� � �� �
�� � �� �
�
� �
=
x x +
L =
x x
=
x x +
P x = LL x +L x +L x
P x =
x x +
+
x x +
1 2 3
2 2
1 2 3
17 66
50
1
12 11
� � � � � �
� � �
�
�
��
�
�
��
�
�225
5
7 6
50
3,3
5,7 79,9 24,2
5
2
2
�
�
��
�
�
��
��
�
��
�
�
��
� � � �
× +
x x +
×
P x =
x + x
00
Com esse polinômio em mãos, é possível calcular o valor da ocupação do estoque para qualquer mês 
do ano. A tabela seguinte mostra esse valor calculado para todos os meses.
Tabela 30
Mês Ocupação do estoque
1 1
2 2.256
3 3.284
4 4.084
5 4.656
6 5
7 5.116
8 5.004
9 4.664
10 4.096
11 3.3
12 2.276
Note que para os meses 1, 6 e 11 os valores de ocupação da última tabela são os mesmos observados 
na tabela anterior a esta.
Para calcular o mês de maior ocupação do estoque, pode‑se calcular a derivada da função e igualá‑la 
a zero. Assim:
�� � � �
�
P x =
x +
= o x = = meses
2 5,7 79,9
50
0
79,9
2 5,7
7,0087log
Ou seja, o estoque terá sua ocupação máxima no mês 7. Essa informação pode ser confirmada 
também na tabela seguinte.
225
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Tabela 31
Mês Ocupação do estoque
1 1
2 2.256
3 3.284
4 4.084
5 4.656
6 5
7 5.116
8 5.004
9 4.664
10 4.096
11 3.3
12 2.276
5.2.2.2 Exemplo de implementação computacional: polinômio interpolador de Lagrange 
no Python
Será visto um algoritmo que calcula o polinômio de Lagrange para um dado conjunto de pontos. 
Esse algoritmo segue a equação ( ) ( )i i
i
P x f x . L=∑ , passando ponto a ponto para criar os termos L1 
 
e finalmente o polinômio P(x). Note que, a cada passo, um valor é calculado e salvo na variável aux. A 
variável aux1 recebe o valor do produto dos polinômios de cada ponto, utilizando a função numpy.polymul.
A variável pol recebe o polinômio final. Na sequência, ele é usado para calcular os valores da função 
para diversos valores de x (salvos na variável XX), armazenando‑os na variável YY. O comando numpy.
polyval é utilizado para realizar essas contas. O código a seguir apresenta essa implementação.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Polinômio interpolador de Lagrange
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
from matplotlib import pyplot as plt
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Dados observados, dispostos nos vetores x e y
x = np.array([‑2, 0, 4])
y = np.array([2, ‑2, 1])
226
Unidade II
#######################################################################
## CÁLCULOS ##
#######################################################################
aux1 = 1;
pol = 0;
# Loop para criação do polinômio
for j in range(np.size(x)):
 aux1 = 1
 for i in range(np.size(x)):
 if i!=j:
 aux = [1, ‑x[i]]/(x[j]‑x[i])
 aux1 = np.polymul(aux1, aux)
 pol = pol + aux1*y[j]
# Aplicação do polinômio a diversos valores de x
XX = np.linspace(np.min(x),np.max(x), 100)
YY = np.polyval(pol, XX)
#######################################################################
## PÓS‑PROCESSAMENTO ##
#######################################################################
plt.clf()
plt.scatter(x,y)
plt.plot(XX,YY)
titulo = ‘Polinômio interpolador de Lagrange’
plt.title(titulo)
A variável pol conterá os coeficientes do polinômio de grau 2 gerado para esse caso. O conteúdo 
dessa variável é mostrado a seguir. Observe que esses números são iguais àqueles mostrados na equação 
211x 13x
P(x) 2
24 12
= − − .
In [1]: pol
Out[1]: array([ 0.45833333, ‑1.08333333, ‑2. ])
Essa execução irá gerar um gráfico, no qual podemos confirmar que o polinômio interpolador passa 
pelos pontos de maneira satisfatória.
Polinômio interpolador de Lagrange
2
1
0
‑1
‑2
‑2 ‑1 0 1 2 3 4
Figura 52 – Polinômio de Lagrange e pontos experimentais para o código programado no Python
227
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Finalmente, criaremos um programa utilizando funções pré‑programadas no Python para realizar o 
polinômio interpolador de Lagrange. Para isso, usaremos a biblioteca scipy, que contém a função scipy.
interpolate.lagrange. O código é apresentado na sequência.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Polinômio interpolador de Lagrange usando o scipy
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
import scipy
from matplotlib import pyplot as plt
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Dados observados, dispostos nos vetores x e y
x = np.array([‑2, 0, 4])
y = np.array([2, ‑2, 1])
#######################################################################
## CÁLCULOS ##
#######################################################################
pol = scipy.interpolate.lagrange(x,y)
XX = np.linspace(np.min(x), np.max(x), 100)
YY = np.polyval(pol, XX)
#######################################################################
## PÓS‑PROCESSAMENTO ##
#######################################################################
plt.clf()
plt.scatter(x,y)
plt.plot(XX,YY)
titulo = ‘Polinômio interpolador de Lagrange’
plt.title(titulo)
Nesse ponto, sugere‑se implementar esse código em seu próprio computador e comparar os 
resultados com os exemplos obtidos anteriormente.
5.2.3 Interpolação de diferenças finitas: Newton‑Gregory
A interpolação polinomial de Newton‑Gregory possui duas particularidades, a primeira é que os 
pontos usados estejam em ordem crescente e de diferença constante (equidistantes), ou seja, podemos 
escrever que (xi+1 ‑ xi = h para i=0,1,2,...,n), sendo o valor de h constante.
A segunda particularidade é a questão de diferença finita (∆), isto é, considerando os pontos:
0 0 1 1 2 2 n n(x ,y ),(x ,y ),(x ,y ), ,(x ,y ) para n = 0, 1, 2, ..., n
Seguem as diferenças para as ordens definidas:
Ordem 0
0
i iy y∆ =
228
Unidade II
Ordem 1
1 0 0
i i 1 i i 1 iy y y y y+ +∆ = ∆ − ∆ = −
Ordem 2
2 1 1
i i 1 iy y y+∆ = ∆ − ∆
Ordem n
n n 1 n 1
i i 1 iy y y
− −
+∆ = ∆ − ∆ ou 
n n 1 n 1f(x) f(x h) f(x)− −∆ = ∆ + − ∆
Segue um exemplo de tabela com o cálculo das diferenças finitas.
Tabela 32 
i xi
yi 
(∆0yi)
(∆1yi) (∆
2yi) (∆
3yi) (∆
4yi)
0 3 6,35
0,77 
(7,12 ‑ 6,35)
1 4 7,12 0,82 
(1,59 – 0,77)
1,59 
(8,71 ‑ 7,12)
‑2,07 
(‑1,25 – 0,82)
2 5 8,71 ‑1,25 
(0,34 – 1,59)
5,07 
(3 ‑ (‑2,07))
0,34 
(9,05 ‑ 8,71)
3 
(1,75 – (‑1,25))
3 6 9,05 1,75 
(2,09 ‑ 0,34)
2,09 
(11,14 ‑ 9,05)
4 7 11,14
Em geral, a forma da interpolação polinomial de Newton‑Gregory:
1 2 n
0 0 0 0
0 0 1 0 1 n 11 2 n
f(x ) f(x ) f(x )
P(x) f(x) (x x ) (x x )(x x ) ... (x x )(x x )...(x x )
1!h 2!h n!h
−
∆ ∆ ∆
= ∆ + − ⋅ + − − ⋅ + + − − − ⋅
Os termos ∆kf(xj) são diferenças ordinárias de ordem k para o ponto xj.
Exemplo 32
Usando a expressão dada para diferenças ordinárias, elabore a tabela de ordem 4 correspondente em 
relação à função f(x) representada a seguir:
229
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Tabela 33 
X ‑1 0 1 2 3
f(x) 2 1 2 5 10
Pode‑se verificar que é respeitado o critério de pontos equidistantes, uma vez que, para os pontos 
da tabela, para n = 4, temos:
4 3 3 2 2 1 1 0x x x x x x x x h 1− = − = − = − = =
Partimos assim para a construção da tabela de diferenças finitas:
Tabela 34 
x F(x) ∆1f(x) ∆2f(x) ∆3f(x) ∆4f(x)
‑1 2
‑1
0 1 2
1 0
1 2 2 0
3 0
2 5 2
5
3 10
A partir da interpolação, o polinômio resultante permite calcular/estimar um ponto que não está 
presente no conjunto de pontos apresentados. Para tanto, utiliza‑se o polinômio obtido através do 
processo de interpolação.
Neste momento será feita a tabela de diferenças finitas para gerar o polinômio de interpolação e 
calcular a aproximação de f(0,5) pelo método de Newton‑Gregory.
Para ordem 4, a forma de Newton‑Gregory define‑se como:
1 2
0 0
4 0 0 0 11 2
3 4
0 0
0 1 2 0 1 2 33 4
f(x ) f(x )
P (x) f(x ) (x x ) (x x )(x x )
1! h 2! h
f(x ) f(x )
(x x )(x x )(x x ) (x x )(x x )(x x )(x x )
3! h 4! h
∆ ∆
= + − + − − +
⋅ ⋅
∆ ∆
+ − − − + − − − −
⋅ ⋅
Substituindo os valores, temos:
230
Unidade II
4 2 3
4
2
4
2
4
1 2 0
P (x) 2 (x 1) (x 1)(x 0) (x 1)(x 0)(x 1)
1! 1 2! 1 3! 1
0
(x 1)(x 0)(x 1)(x 2)
4! 1
P (x) 2 (x 1) (x x)
P (x) x 1
−
= + + + + − + + + − − +
⋅ ⋅ ⋅
+ + − − −
⋅
= − + + +
= +
Assim, P4(x) = x
2 + 1 é polinômio obtido através da interpolação de Newton‑Gregory para o conjunto 
de valores numéricos apresentados na tabela.
Podemos agora fazer um cálculo aproximado da função para um valor que não está presente na 
tabela; ou seja, estimar, a partir do uso adequado do polinômio de interpolação calculado, de forma 
aproximada, o valor da nova função para x = 0,5.
2
4P (0,4) (0,5) 1 0,25 1 1,25= + = + =
Figura 53
Exemplo 33
Usando a expressão dada para diferenças ordinárias, elabore a tabela de ordem 4 correspondente em 
relação à função f(x) representada a seguir:
Tabela 35 
X 3 4 5 6 7 8
f(x) ‑1 5 7 8 6 4
231
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Confirma‑se que:
5 4 4 3 3 2 2 1 1 0x x x x x x x x x x h 1− = − = − = − = − = =
Para cada operador diferença devemos seguir a fórmula apresentada, temos:
Tabela 36 
x F(x) ∆1f(x) ∆2f(x) ∆3f(x) ∆4f(x)
3 ‑1
(5 + 1) = 6
4 5 (2 – 6) = ‑4
(7 ‑ 5) = 2 (‑1 ‑ ‑4) = 3
5 7 (1 ‑2) = ‑1 (‑2 – 3) = ‑5
(8 – 7) = 1 (‑3 ‑ ‑1) = ‑2
6 8 (‑2 – 1) = ‑3 (‑3 ‑ ‑2) = ‑1
(6 – 8) = ‑2 (0 ‑ ‑3) = 3
7 6 (‑2 ‑ ‑2) = 0
(4 – 6) = ‑2
8 4
A interpolação de Gregory‑Newton possui como condição para sua formulação que os pontos usados 
estejam em ordem crescente e de diferença constante, ou seja, podemos escrever que (xi=1 ‑ xi = h, 
i=0,1,2,...,n). O valor de h é constante.
5.2.4 Interpolação de Splines
A origem do nome Spline é referente a uma régua maleável utilizada para desenhos e projetos de 
engenharia. Esta régua pode ser deformada para passar por um conjunto de pontos. Spline pode ser 
entendida também como “curva flexível”.
Figura 54 – Spline
232
Unidade II
A interpolação Spline consiste em dividir um intervalo de interesse em diversos subintervalos e neles 
definir polinômios de grau pequeno, buscando uma aderência aos dados da forma mais suave possível. 
Assim, por retornar polinômios de graus menores – normalmente linear (grau 1), quadrática (grau 2) ou 
cúbica (grau 3), tem‑se preferência em seu uso em detrimento das demais interpolações polinomiais, 
uma vez que os outros métodos podem atingir graus elevados em seus polinômios interpoladores 
quando o conjunto de pontos é muito grande.
Portanto, o conceito da interpolação Spline, ao invés de usar um único polinômio para interpolar 
todos os dados, segmenta o polinômio em diversos intervalos, apresentando uma maior simplicidade, 
acurácia e aproximação aos dados. No entanto, o polinômio só é aplicável ao seu intervalo específico. 
As Splines são as que apresentam melhores resultados nas aplicações computacionais como na área de 
computação gráfica e desenhoauxiliado por computadores (CAD – Computer Aided Design).
Uma visualização da interpolação por Splines para os pontos:
Tabela 37 
x 0 1 2 3 4 5 6 7 8 9
y 1 9 5 10 0 4 1 10 0 8
É conforme segue:
Figura 55 
Por exemplo, para a Spline linear, podemos ver que existem nove equações polinomiais para 
descrevê‑la. Isto é, entre cada dois pontos, há uma equação de reta única para expor este trecho. Já 
para a Spline cúbica busca‑se uma equação polinomial cúbica entre intervalos definidos respeitando‑se 
a transição entre uma equação e outra através da análise de suas derivadas para que a transição entre 
um polinômio e outro (intervalo aplicável) seja a mais suave possível.
233
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
 Lembrete
A spline também pode ser entendida como curva flexível.
Definição matemática
Podemos definir matematicamente a construção da Spline como um conjunto de polinômios S(x) 
que atendem ao seguinte critério:
0 0 1
1 1 2
n 1 n 1 n
S (x), para x x x
S (x), para x x x
S(x)
S (x), para x x x− −
≤ <
 ≤ <= 

 ≤ <

Para Splines de primeira ordem, também chamadas de linear ou grau 1, sendo os nós 1 2 nx ,x , , x , 
em cada subintervalo [xi‑1, xi], i = 1, 2, ..., n:
i i 1
i i 1 i i 1 i
i i 1 i i 1
x x x x
S f(x ) f(x ) , x [x ,x ]
x x x x
−
− −
− −
− −
= + ∀ ∈
− −
Importante mencionar que, no caso da Spline linear, utilizamos sempre um par de pontos para 
determinar a reta Si(x). Significa simplesmente unir cada um dos pontos mediante segmentos de retas 
na forma apresentada a seguir:
... etc
Figura 56 
Exemplo 34
Achar a função Spline linear que interpola a função representada pelo conjunto de pontos a seguir:
234
Unidade II
Tabela 38 
x0 x1 x2 X3
x 1 2 4 8
f(x) 1 2 3 1
Pela definição:
01
1 0 1
1 0 1 0
x xx x 2 x x 1
S f(x ) f(x ) =1 2 (2 x) (2x 2) x 
x x x x 2 1 2 1
−− − −
= + + = − + − =
− − − −
1 S (x) x, x [1,2]• = ∀ ∈
2 1
2 1 2
2 1 2 1
x x x x 4 x x 2 3 x
S f(x ) f(x ) =2 3 (4 x) (x 2) +1 
x x x x 4 2 4 2 2 2
− − − −
= + + = − + − =
− − − −
2
x
 S (x) +1 , x [2,4]
2
• = ∀ ∈
3 2
3 2 3
3 2 3 2
x x x x 8 x x 4 3 1 x
S f(x ) f(x ) =3 1 (8 x) (x 4) 4 
x x x x 8 4 1 3 4 2 4
− − − −
= + + = − + − = −
− − − −
3
x
 S 4 ‑ , x [4,8]
4
• = ∀ ∈
Graficamente temos:
Figura 57 
235
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Spline cúbica
A Spline linear tem como desvantagem a derivada primeira descontínua nos nós de transição entre 
polinômios aplicáveis. No caso da Spline quadrática (grau 2), as derivadas são contínuas apenas até 
o grau 1 e a curvatura nos pontos de intersecção pode apresentar inversão de sentido (crescente/
decrescente) pelo sinal da derivada neste ponto quando aplicada na curva que representa o trecho 
da esquerda e da direita e utilizada no mesmo ponto de interseção. Assim são mais usadas as Splines 
cúbicas, que são polinômios de grau 3 no intervalo nos intervalos [xi‑1, xi], i = 1, 2, ..., n.
Destaca‑se também que a primeira e a segunda derivadas dos polinômios da Spline cúbica são 
contínuas, o que garante que a Spline não tenha picos e não troque de curvatura nos nós de transição.
Sabendo que a Spline cúbica será composta de polinômios de grau 3, teremos quatro constantes a 
serem determinadas nesse processo. Assim, a forma para cada seção da Spline cúbica tem a forma:
2 3
k k k k k k k ks a b (x x ) c (x x ) d (x x ) , k 1,2,...,n= + − + − + − =
Lembre‑se que:
sk(xk) = ak = f(xk) ‑ o polinômio possui o mesmo valor da função quando aplicado no ponto Σ.
sk+1(xk+1) = sk(xk+1) ‑ o valor do polinômio é o mesmo na transição de Sk aplicável.
Considerando ainda a primeira e segunda derivadas de Sk, podemos escrever:
' 2
k k k k k k
''
k k k k
s (x) 3a (x x ) 2b (x x ) c
s (x) 6a (x x ) 2b
= − + − +
= − +
Note que 
''
'' k k
k k k k
s (x )
s (x ) 2b , permitindo escrever b
2
= = . Usando ainda a notação k k k 1h x x −= − , 
tem‑se que:
''
k k 1 k k k ks (x ) 6a h 2b , reescrevemos assim a :− = − +
'' '' ''
k k k 1 k k k k 1
k
k k
2b s (x ) s (x ) s (x )
a 
6h 6h
− −− −= =
Uma vez que dk = f(xk) e já tivermos definido ak e bk, podemos obter ck em função das derivadas 
segunda nos nós:
3 2
2k 1 k k k k k k k 1
k k k k k
k k
'' '' ''
k k 1 k k k k 1 k k
k k
k
f(x ) a h b h d f(x ) f(x )
c (a h b h )
h h
f(x ) f(x ) [s (x ) s (x )] s (x )
h h
h 6 2
− −
− −
− − + + −
= = − −
 − − = − − 
  
236
Unidade II 3 2
2k 1 k k k k k k k 1
k k k k k
k k
'' '' ''
k k 1 k k k k 1 k k
k k
k
f(x ) a h b h d f(x ) f(x )
c (a h b h )
h h
f(x ) f(x ) [s (x ) s (x )] s (x )
h h
h 6 2
− −
− −
− − + + −
= = − −
 − − = − − 
  
Portanto:
'' ''
k k 1 k k k k 1 k 1 k
k
k
f(x ) f(x ) 2s (x )h s (x )h
c
h 6
− − −− − −= −
Com o objetivo de simplificar a apresentação, podemos usar as seguintes definições: 
''
k k k k ks (x ) g e f(x ) y= = que nos permite reescrever:
k k 1
k
k
g g
a
6h
−−=
k
k
g
b
2
=
k k 1 k k k 1 k
k
k
y y 2h g g h
c
h 6
− − − += + 
 
k kd y=
Para k = 1, 2, ..., n, podemos calcular todos os coeficientes de Sk(x) em função de 
''
j j jg s (x ) , j = 0, 1, ... ,n= .
Pode‑se genericamente escrever:
k 1 k k k 1
k k 1 k k 1 k k 1 k 1
k 1 k
y y y y
h g 2(h h )g h g 6
h h
+ −
− + + +
+
 − −
+ + + = −  
Que é um sistema de equações lineares (Ax = b), o qual pode ser reescrito:
1 1 2 3
2 2 3 4
h 2(h h ) h
 h 2(h h ) h
 . . .
A
 . . .
 
+
+
=
n 1 n 1 n n (n 1)x(n 1)
 . . .
 h 2(h +h ) h − − − +
 
 
 
 
 
 
 
  
237
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
e
1 02 1
2 1
3 2 2 1
3 2
n n 1 n 1 n 2
n n 1 (n 1)x1
y yy y
h h
y y y y
h hb
 
y y y y
h h
− − −
− −
−− − 
 
− − − =
 
 
 − −
− 
 
 
Para resolver este sistema de equações será necessária a imposição de duas condições a fim de 
podermos determinar ak, bk, ck e dk a todo sk(x). As condições possíveis são:
Spline natural
Esta escolha é equivalente a dizer que os polinômios em seus extremos são linerares. Assim:
'' ''
3 0 0 3 n ns (x ) g 0 e s (x ) g 0= = = =
Vamos supor que as cúbicas são parábolas nos extremos:
0 1 n n 1g g , g g −= =
Valores determinados de inclinação (derivada primeira) nos extremos ' '3 0 3 ns (x ) A e s (x ) B= = , 
assim obtemos duas novas equações:
' 2
1 0 1 1 1
'
n n n
s (x ) 3a h 2b h c A
s (x ) c B
= − + =
= =
Exemplo 35
Dado o conjunto de pontos a seguir, encontre uma aproximação para f(0,3) por Spline cúbica natural.
Tabela 39 
X 0 0,5 1,0 1,5 2,0
f(x) 3 1,8616 ‑0,5571 ‑4,1987 ‑9,0536
O conjunto de pontos apresentado x varia de 0 a 2,0, temos 4 intervalos, logo, n=4.
238
Unidade II
Determinar 1 2 3 4s (x),s (x),s (x) e s (x) resolvendo para 1 k 3 (n‑1 3)≤ ≤ = , sendo:
k 1 k k k 1
k k 1 k k 1 k k 1 k 1
k 1 k
y y y y
h g 2(h h )g h g 6
h h
+ −
− + + +
+
 − −
+ + + = −  
E dado o conjunto de dados em que kh h 0,5= = , ficando:
k 1 k k 1 k 1 k k 1
0 1 2 2 1 0
1 2 3 3 2 1
2 3 4 4 3 2
6
hg 4hg hg (y 2y y )
h
6
hg 4hg hg (y 2y y )
h
6
hg 4hg hg (y 2y y )
h
6
hg 4hg hg (y 2y y )
h
− + + −+ + = − +
 + + = − +

 + + = − +

 + + = − +
Como pede‑se via Spline cúbica natural, fazemos g0 = g4 = 0, resultando no sistema:
1 2 2 1 0
1 2 3 3 2 1
2 3 4 3 2
2 1 01
2 3 2 1
3 4 3 2
4hg hg (6 / h)(y 2y y )
hg 4hg hg (6 / h)(y 2y y )
 hg 4hg (6 / h)(y 2y y )
y 2y yg4h h 0
6
h 4h h g y 2y y
h
0 h 4h g y 2y y
+ = − +
 + + = − +
 + = − +
− +   
    = − +
   
  − +   
1 2 2 1 0
1 2 3 3 2 1
2 3 4 3 2
2 1 01
2 3 2 1
3 4 3 2
4hg hg (6 / h)(y 2y y )
hg 4hg hg (6 / h)(y 2y y )
 hg 4hg(6 / h)(y 2y y )
y 2y yg4h h 0
6
h 4h h g y 2y y
h
0 h 4h g y 2y y
+ = − +
 + + = − +
 + = − +
− +   
    = − +
   
  − +   
Realizando a substituição de h por yi
1
2
3
g 15,36362 0,5 0
0,5 2 0,5 g 14,6748
0 0,5 2 14,5598g
−    
     = −
    
−    
Resolvendo‑se o sistema de equações, temos:
g1 = ‑6,252
g2 = ‑4,111
g3 = ‑6,654
239
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Com estes valores, podemos obter s1(x), s2(x), s3(x) e s4(x). Como a aproximação é de f(0,3), utiliza‑se s1(x):
3 2
1 1 1 1 1 1 1 1
1 0
1
1
1
1 0 1 0
1
1 1
s (x) a (x x ) b (x x ) c (x x ) d
g g 6,654
a 2,218
6h 3
g
b 3,327
2
y y 2hg g h
c 3,386
h 6
d y 1,8616
= − + − + − +
− −
= = = −
= = −
− +
= + = −
= =
Assim:
3 2
1s (0,3) 2,218( 0,3) 3,327(0,3) 3,386( 0,3) 1,8616 2,777= − − − − − + =
Portanto, por Spline cúbica natural entende‑se:
f(0,3) = s1(0,3)=2,777
 Observação
O estudo de interpolação e regressão possibilita ao estudante perceber 
uma perspectiva nova em relação ao estudo de funções. As funções nem 
sempre são dadas por fórmulas. Muitas delas que ocorrem em aplicações 
não têm fórmulas específicas, alerta Ávila (2003). É o caso da temperatura 
num determinado lugar a cada instante de tempo, a corrente cardíaca 
de certa pessoa, registrada num eletrocardiograma como função de 
tempo etc. Nesses eventos, observa o autor, elas são conhecidas, ao 
menos parcialmente, por tabelas e valores, dados numéricos resultantes 
de registros de observações, os quais levam à construção de gráficos 
representativos dessas funções.
6 REGRESSÃO LINEAR POR MÍNIMOS QUADRADOS
Quando um erro substancial estiver associado aos dados, a estratégia mais adequada é 
determinar uma função que se ajuste à forma ou tendência geral deles sem necessariamente 
passar pelos pontos individuais.
240
Unidade II
Figura 58
A figura anterior ilustra como uma reta pode ser usada para caracterizar a tendência geral dos dados 
sem passar por nenhum dos pontos particulares.
6.1 Regressão linear ou ajuste de curvas com equações lineares
O ajuste de curvas usando uma equação linear (polinômio de primeiro grau) é o processo pelo qual, 
a partir da forma geral,
1 0f(x) a x a= +
obtém‑se uma função que representa o padrão das observações ao se determinar constantes a1 e 
a0, as quais possibilitam o menor erro possível quando os pontos medidos são substituídos na equação 
obtida. Caso só se disponha de duas observações, as constantes podem ser determinadas de forma tal 
que a função obtida forneça os valores exatos nos pontos.
Figura 59
Quando as observações consistirem em mais de dois pontos, as constantes a1 e a0 serão determinadas 
de tal forma que a reta promova o melhor ajuste como um todo.
241
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Figura 60
Em muitas situações cotidianas, uma relação linear entre duas variáveis pode ser identificada por 
meio de gráficos denominados “diagramas de dispersão”. No entanto, essa opção é subjetiva, pode variar 
conforme a forma de avaliação de cada pessoa. Daí que, para se obter a função que desempenha o 
comportamento padrão dos dados, é necessário estabelecer o comportamento do conjunto de duas 
variáveis de forma objetiva, através do que é normalmente dito a “força da relação das variáveis”. Essa 
avaliação de comportamento é denominada “correlação linear entre as variáveis”.
6.2 Diagrama de dispersão
Um diagrama de dispersão é a primeira indicação da possível existência de uma associação entre 
duas variáveis. Smailes e McGrane (2002) simplificam a linguagem quando se referem à forma de análise 
de um diagrama de dispersão, pois ele é uma representação do mundo real e quase sempre sua utilidade 
é auxiliar decisões, mas sem negligenciar o fato de que são análises subjetivas.
O eixo y é utilizado para representar a variável dependente que interessa 
a quem toma decisões, enquanto o eixo x é para representar uma variável 
que pode ser controlada ou medida por quem toma decisões (geralmente 
chamada de variável independente) (SMAILES; MCGRANE, 2002, p. 115).
Nesse momento da análise é importante, observam as autoras, justificar com antecedência na 
sequência de cálculos a possibilidade de que uma correlação entre variáveis ser mera coincidência. Para 
tanto, elas denominam a variável independente de variável causa e a variável dependente de efeito 
resultante das mudanças em x. Daí é necessário identificar se a variável y é de fato o efeito resultante 
das mudanças em x, a variável causa.
242
Unidade II
 Saiba mais
Para saber mais sobre o diagrama da dispersão, veja a obra:
SMAILES, J.; MCGRANE, Â. Estatística aplicada à administração com 
Excel. São Paulo: Atlas, 2002.
A planilha Excel® possibilita a obtenção de diagramas de dispersão:
• Introduzir os valores das duas variáveis em duas células distintas e seus respectivos valores na 
sequência. Selecionar as linhas e colunas utilizadas e a opção “dispersão”. Em seguida, escolher 
“dispersão somente com marcadores”.
Figura 61
• Selecionar “linha de tendência” e “linha de tendência linear”. Há outras linhas de tendência que 
poderão ser selecionadas caso o padrão dos dados no gráfico sugira uma linha de tendência 
exponencial, logarítmica etc.
Figura 62
243
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Gráfico de dispersão final.
Figura 63
Exemplo 36
Um pesquisador inqueriu um grupo de pessoas da mesma faixa etária sobre o nível de escolaridade 
e o número de carros que tiveram até o momento. As perguntas foram:
a) Quantos anos você estudou?
b) Quantos carros você comprou em sua vida?
As respostas foram tabeladas da seguinte forma:
xi representa o número de anos que a pessoa estudou;
yi representa o número de carros que a pessoa teve em sua vida.
244
Unidade II
Tabela 40
Xi 3 5 7 9 10
yi 1 2 3 5 7
Apresentar o diagrama de dispersão dessa pesquisa.
De acordo com o exposto na teoria anteriormente, o gráfico de dispersão pode ser criado pelo Excel® 
e a linha de tendência é linear. Ou seja, há fortes indícios de que o número de anos de escolaridade 
influencie a capacidade de aquisição das pessoas.
Figura 64
O exemplo anterior, acrescido da figura a seguir, poderia sugerir a conclusão: o grau de escolaridade 
das pessoas influencia a situação caótica das grandes cidades, visto que está relacionado com a capacidade 
aquisitiva. Observe que, embora uma inferência sobre o caos das cidades ser decorrente do excesso de carros 
possa ser válida, ela não está relacionada com o modelo proposto na pesquisa citada no exemplo. Não era 
essa a pergunta do pesquisador, as variáveis são diferentes e, portanto, uma afirmativa como essa seria 
arbitrária. Percebe‑se, assim, que o modelo matemático tem relação direta com a motivação de uma pesquisa.
245
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Figura 65
6.3 Coeficiente de correlação
O grau de correlação entre as variáveis é dado pelo coeficiente de Pearson:
( ) ( )
( ) ( )
i i i i
2 22 2
i i i i
n x y x y
r
n x x . n y y
−
=
   − −   
∑ ∑ ∑
∑ ∑ ∑ ∑
Os valores de correlação linear variam de ‑1 a 1. Ou seja: 1 r 1− ≤ ≤
Quando as variáveis do modelo têm variações contrárias, ou seja, se a variável independente é 
crescente em determinado intervalo e a respectiva variável dependente é decrescente (ou vice‑versa), 
o coeficiente de correlação é negativo. Normalmente, os valores a seguir indicam como classificar uma 
correlação negativa.
Quadro 1 
Coeficiente Classificação
r = ‑ 1,00 Correlação negativa perfeita
r = ‑ 0,75 Correlação negativa forte
r = ‑ 0,50 Correlação negativa média
r = ‑ 0,25 Correlação negativa fraca
246
Unidade II
Quando as variáveis do modelo têm variações no mesmo sentido, ou seja, se a variável independente 
é crescente em determinado intervalo e a respectiva variável dependente também é crescente (ou 
se ambas são decrescentes), o coeficiente de correlação é positivo. Assim, na correlaçãopositiva, 
normalmente os valores na sequência indicam como classificar uma correlação.
Quadro 2 
Coeficiente Classificação
r = 1,00 Correlação positiva perfeita
r = 0,75 Correlação positiva forte
r = 0,50 Correlação positiva média
r = 0,25 Correlação positiva fraca
 Lembrete
Quando r = 0,0 a correlação linear não existe. Porém, a relação entre as 
variáveis pode ser não linear.
Exemplo 37
Utilizando os dados a seguir, obter o coeficiente de correlação do problema de pesquisa.
Resolução:
Tabela 41
Xi 3 5 7 9 10
yi 1 2 3 5 7
A construção de uma tabela, como a tabela a seguir, auxilia no cálculo do coeficiente de correlação:
Tabela 42 
Xi yi x
2
i y
2
i xi yi 
3 1 (3)2 = 9 (1)2 = 1 (3 . 1) = 3
5 2 (5)2 = 25 (2)2 = 4 (5 . 2) = 10
7 3 (7)2 = 49 (3)2 = 9 (7 . 3) = 21
9 5 (9)2 = 81 (5)2 = 25 (9 . 5) = 45
10 7 (10)2 = 100 (7)2 = 49 (10 . 7) = 70
Σxi = 34 Σyi = 18 Σx
2
i = 264 Σy
2
i = 88 Σxiyi = 149
247
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
( ) ( )
( ) ( )
i i i i
2 22 2
i i i i
n x y x y
r
n x x . n y y
−
=
   − −   
∑ ∑ ∑
∑ ∑ ∑ ∑
2 2
5.149 34.18
r
(5.264 (34) )(5.88 (18) )
−
=
− −
 Lembrete
Observe que n corresponde ao número de pares dos dados, ou seja 
(pergunta, resposta). Em relação aos dados tabelados, n coincide com o 
número de linhas da tabela.
745 612 133 133 133
r 0,96
137,93(1320 1156)(440 324) 164.116 19024
−
= = = = =
− −
O coeficiente de correlação r = 0,96 indica uma “correlação linear positiva forte”, o que coincide com 
a observação inicial no diagrama de dispersão. Além disso, nesse caso, a correlação é positiva porque 
quando a variável independente xi (número de anos que a pessoa estudou) cresce, a variável dependente 
yi (número de carros que a pessoa teve em sua vida) cresce.
 
Regressão linear é um procedimento muito útil em Estatística. Quando se domina o 
conceito, é possível utilizar a planilha eletrônica Excel® para facilitar o procedimento, uma 
vez que a coleção de dados em pesquisa não é normalmente pequena e os cálculos são 
tediosos. O suporte do Office fornece as instruções:
Sintaxe
PEARSON(matriz1;matriz2)
Matriz1 é um conjunto de valores independentes.
Matriz2 é um conjunto de valores dependentes.
Utilizando o conjunto de dados do exemplo anterior, selecionar uma célula para a variável 
independente e outra para a variável dependente. A matriz coluna da variável independente 
será constituída pelas linhas A2:A6. A matriz coluna da variável dependente será constituída 
pelas linhas B2:B6. Selecionar a opção “mais funções”, “Estatística” e “Pearson”.
248
Unidade II
Figura 66
Ao se inserir os domínios das matrizes colunas das variáveis independente e dependente, 
o resultado da correlação aparece imediatamente no quadro:
r = 0,964274589
Como tal precisão não é necessária: r = 0,96.
 
6.4 Regressão linear
Esse processo está intimamente ligado ao conceito de correlação visto em tópico 
anterior e consiste em estabelecer uma função que relacione variáveis com correlação linear 
considerável. No entanto, como foi comentado anteriormente, o conjunto de dados de um 
problema de pesquisa (seja acadêmica, industrial comercial, jornalística etc.) pode ser grande, 
e inúmeras retas podem se adequar ao padrão estabelecido por ele. Para selecionar a função 
ideal é utilizado o “método dos mínimos quadrados”. A representação gráfica dessa função é 
denominada “reta interpoladora”.
O método dos mínimos quadrados consiste em obter os coeficientes a0 e a1 da função f(x) = a1x + a0.
No entanto, não se pode desprezar o processo de buscar uma aproximação. Logo, seria mais adequado:
f(x) = a1x + a0 + e
O valor de e é o erro ou resíduo, ou seja, a discrepância entre o valor verdadeiro de f(x) e o aproximado, 
a1x + a0, previsto pela equação linear.
249
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Usando a forma:
1 0y a x a e= + +
Então:
1 0e y a x a= − −
Dessa maneira, a soma dos erros residuais para todos os dados disponíveis seria dada pela expressão 
a seguir, a qual representa o total de erros, ou seja, todas as diferenças, ponto a ponto, entre o valor 
verdadeiro e o valor aproximado:
n n
i i 1 i 0
i 1 i 1
e y a x a
= =
= − −∑ ∑
Observação: n é número total de pontos.
Mas as diferenças entre os valores verdadeiros e aproximados podem ser negativas ou positivas, de 
forma que a melhor opção é tomar essas diferenças em valor absoluto:
n n
i i 1 i 0
i 1 i 1
e y a x a
= =
= − −∑ ∑
A partir dessa expressão, o objetivo seria escolher a opção que minimizasse a distância máxima que 
um ponto individual tenha da reta interpoladora. Mas, mesmo assim, haveria a influência indevida de 
pontos com erro grande. Daí a opção por minimizar a soma dos quadrados dos resíduos entre o ymedido e 
o ycalculado computado com o modelo linear:
( )
n n n
22 2
r i i medido imodelo i 1 i 0
i 1 i 1 i 1
S e (y y ) y a x a
= = =
= = − = − −∑ ∑ ∑
Esse critério fornece uma única reta para um dado conjunto de dados.
Derivando expressão:
( )
n
2
r i 1 i 0
i 1
S y a x a
=
= − −∑
( )
n
r
i 1 i 0
0 i 1
S
2( 1) y a x a
a =
∂
= − − −
∂ ∑
250
Unidade II
( )
n
r
i 1 i 0
0 i 1
S
2 y a x a
a =
∂
= − − −
∂ ∑ (1)
( )
n
r
i i 1 i 0
1 i 1
S
2( x ) y a x a
a =
∂
= − − −
∂ ∑
( )
n
r
i i 1 i 0
1 i 1
S
2 x y a x a
a =
∂
= − − −
∂ ∑
( )
n
r
i i 1 i i i 0
1 i 1
S
2 x y a x x x a
a =
∂
= − − −
∂ ∑
( )
n
2r
i i 1 i i 0
1 i 1
S
2 x y a x x a
a =
∂
= − − −
∂ ∑ (2)
Igualando as duas expressões a zero para obter o valor crítico:
Expressão (1):
( )
n
i 1 i 0
i 1
0 2 y a x a
=
= − − −∑
( )
n
i 1 i 0
i 1
y a x a 0
=
− − =∑
Fazendo a distribuição do somatório:
n n n
i 1 i 0
i 1 i 1 i 1
y a x a 0
= = =
− − =∑ ∑ ∑
Como a0 é uma constante da função linear, é possível fazer:
n n
0 0 0
i 1 i 1
a a 1 na
= =
= =∑ ∑
n n
i 1 i 0
i 1 i 1
y a x na 0
= =
− − =∑ ∑
n n
0 i 1 i
i 1 i 1
na y a x
= =
− = − +∑ ∑
251
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
n n
0 i 1 i
i 1 i 1
na y a x
= =
= −∑ ∑
n n
i 1 i
i 1 i 1
0
y a x
a
n
= =
−
=
∑ ∑
n n
i 1 i
i 1 i 1
0
y a x
a
n n
= == −
∑ ∑
Como a1 é uma constante da função linear, é possível tirá‑la do somatório:
n n
i 1 i
i 1 i 1
0
y a x
a
n n
= == −
∑ ∑
Da expressão (2):
( )
n
2
i i 1 i i 0
i 1
0 2 x y a x x a
=
= − − −∑
( )
n
2
i i 1 i i 0
i 1
x y a x x a 0
=
− − =∑
Fazendo a distribuição do somatório:
n n n
2
i i 1 i i 0
i 1 i 1 i 1
x y a x x a 0
= = =
− − =∑ ∑ ∑
n n n
2
1 i i i i 0
i 1 i 1 i 1
a x x y x a
= = =
− = − +∑ ∑ ∑
n n n
2
1 i i i i 0
i 1 i 1 i 1
a x x y x a
= = =
= −∑ ∑ ∑
252
Unidade II
Como a0 e a1 são constantes da função linear, é possível tirá‑las do somatório:
n n n
2
1 i i i 0 i
i 1 i 1 i 1
a x x y a x
= = =
= −∑ ∑ ∑
n n n
2
1 i i i 0 i
i 1 i 1 i 1
a x x y a x
= = =
= −∑ ∑ ∑
Substituindo o valor de 
n n
i 1 i
i 1 i 1
0
y a x
a
n n
= == −
∑ ∑
n n
i 1 in n n
2 i 1 i 1
1 i i i i
i 1 i 1 i 1
y a x
a x x y x
n n
= =
= = =
 
 
 = − −
 
 
 
∑ ∑
∑ ∑ ∑
n n n n
i i i 1 in n
2 i 1 i 1 i 1 i 1
1 i i i
i 1 i 1
x y x a x
a x x y
n n
= = = =
= =
 
 
 = − −
 
 
 
∑ ∑ ∑ ∑
∑ ∑
n n n n n
i i i i i 1 in
2 i 1 i 1 i 1 i 1 i 1
1 i
i 1
n x y x y x a x
a x
n
= = = = =
=
− +
=
∑ ∑ ∑ ∑ ∑
∑
n n n n n n
2
1 i i i i i i 1 i
i 1 i 1 i 1 i 1 i 1 i 1
na x n x y x y x a x
= = = = = =
= − +∑ ∑ ∑ ∑ ∑ ∑
n n n n n n
2
1 i i 1 i i i i i
i 1 i 1 i 1 i 1 i 1 i 1
na x x a x n x y x y
= = = = = =
− = −∑ ∑ ∑ ∑ ∑ ∑
n n n n n n
2
1 i i i i i i i
i 1 i 1 i 1 i 1 i 1 i 1
a n x x x n x y x y
= = = = = =
 
− = − 
 
∑ ∑ ∑ ∑ ∑ ∑
2n n n n n
2
1 i i i i i i
i 1 i 1 i 1 i 1 i 1
a n x x n x y x y
= = = = =
  
 − = −    
∑ ∑ ∑ ∑ ∑
253
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
n n n
i i i i
i 1 i 1 i 1
1 2n n
2
i i
i 1 i 1
n x y x y
a
n x x
= = =
= =
−
=
 
−  
 
∑ ∑ ∑
∑ ∑
 
6.5 Método dos mínimos quadrados – implementação computacional
Seráapresentada uma aplicação do método dos mínimos quadrados, além de sua implementação 
no MS‑Excel® e no Python, tanto programando diretamente as equações como utilizando as fórmulas 
pré‑programadas nesses softwares.
Inicialmente, vamos utilizar o exemplo da quantidade de peças produzidas por uma empresa em 
diferentes intervalos. A tabela a seguir mostra os dados de produção para as seis primeiras horas.
Tabela 43 – Dados de produção para o exemplo
Tempo (h) Quantidade produzida
1 123
2 233
3 358
4 460
5 581
6 697
Deseja‑se uma equação de primeiro grau (y = a . x + b) para estimar a produção em períodos 
diferentes. Para o nosso exemplo, teremos então:
Qtd = a . t + b
Para encontrarem‑se os parâmetros a (inclinação da curva) e b (interceptação da curva no eixo y), 
deve‑se, pelo método dos mínimos quadrados, resolver o sistema de equações dado a seguir.
i i
2
i i i i
b.n a. x y
b. x a. x x y
 + =

+ = +
∑ ∑
∑ ∑ ∑
Desenvolvendo‑se esse sistema linear, é possível calcular os parâmetros a e b com as equações na 
sequência, respectivamente.
254
Unidade II
( )
i i i i
22
i i
i i
n x x y x y
a
n x x x
y a x x
b
n
+ − +
=
−
−
=
∑ ∑
∑ ∑
∑ ∑
Observa‑se, portanto, que para resolver ambas as equações, é necessário calcular o valor dos termos 
de somatórios. Uma forma prática de realizar isso é colocar os dados em uma tabela e criar colunas para 
o cálculo de x2 e x . y. A tabela a seguir mostra essas contas.
Tabela 44 – Tabela auxiliar para o cálculo dos somatórios
x y x2 x . y
1 123 1 123
2 233 4 466
3 358 9 1074
4 460 16 1840
5 581 25 2905
6 697 36 4182
Soma 21 2452 91 10590
Substituindo esses números na equação 
( )
i i i i
22
i i
n x x . y x x y
a
n x x x
−
=
−
∑ ∑
∑ ∑
 e na equação 
 
i iy a x xb
n
−
= ∑ ∑ , chegamos a:
( ) ( )
i i i i
2 22
i i
i i
n x x . y x x y 6 x 10590 21x 2452
a 114,74
6 x 91 21n x x x
y a x x 2452 114,74 x 21
b 7,07
n 6
− −
= = =
−−
− −
= = =
∑ ∑
∑ ∑
∑ ∑
Dessa forma, a equação na sequência mostra a relação entre quantidade de itens produzidos com o 
tempo de produção.
Qtd = 114,74 . t + 7,07
255
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Exemplo de implementação computacional: método dos mínimos quadrados no Excel®
Esta seção mostrará a implementação do método dos mínimos quadrados para equações lineares no 
Excel®, utilizando duas aproximações diferentes: na primeira, as contas serão realizadas programando a 
 
equação 
( )
i i i i
22
i i
n x x . y x x y
a
n x x x
−
=
−
∑ ∑
∑ ∑
 e a equação i i
y a x x
b
n
−
= ∑ ∑ . Na segunda, utilizaremos as 
 
fórmulas pré‑programadas no Excel®.
Programação das equações
Para esse exemplo, vamos tomar uma planilha que contém, inicialmente, apenas os dados constantes 
a seguir.
Figura 67 – Dados de entrada para o problema
Para o próximo passo, é necessário criar colunas para o cálculo de x2 e x . y. Escolhendo‑se a coluna C 
para x2 e a coluna D para x . y, devemos inserir nas células C2 e D2, respectivamente, as fórmulas a seguir:
Célula C2: =A2^2
Célula D2: =A2*B2
Copiando essas fórmulas para as linhas de baixo, chega‑se aos valores mostrados na figura 
na sequência:
Figura 68 – Colunas auxiliares para o cálculo dos somatórios
256
Unidade II
Devemos agora adicionar os valores das variáveis necessárias para resolver o sistema de equações 
 
lineares da equação 
i i
2
i i i i
b.n a. x y
b. x a. x x . y
 + =

+ =
∑ ∑
∑ ∑ ∑
. Esses dados podem ser calculados utilizando a 
 
função SOMA. A quantidade de números de informações pode ser calculada utilizando a função CONT.
NÚM. A figura a seguir mostra essas células com suas respectivas fórmulas.
Figura 69 – Valores dos somatórios e suas respectivas fórmulas
Observe que essa planilha é genérica, uma vez que as somas são realizadas em todos os dados 
numéricos contidos nas colunas A, B, C e D. Para utilizá‑la com qualquer conjunto de dados, basta entrar 
com os valores de x a coluna A, os valores de y a coluna B e realizar os cálculos para x2 e x . y as colunas 
C e D, respectivamente.
Por fim, devemos calcular os valores dos coeficientes a e b, utilizando as equações 
 
( )
i i i i
22
i i
n x .y x y
a
n x x
× ∑ − ∑ × ∑
=
× ∑ − ∑
 e i iy a xb
n
∑ − × ∑
= , respectivamente. Para isso, escolhemos as células 
 
G10 e G11, que deverão receber as seguintes fórmulas:
Célula G10: =(G2*G6‑G3*G4)/(G2*G5‑G3^2)
Célula G11: =(G4‑G10*G3)/G2
A figura a seguir mostra a forma final da planilha. Observe que os valores encontrados para a e b são 
os mesmos encontrados anteriormente, validando assim nossa planilha.
Figura 70 – Valores dos coeficientes a e b encontrados pela planilha
257
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Utilização de funções nativas do Excel®
O Excel® possui nativamente fórmulas que permitem rapidamente que sejam encontrados os valores 
dos coeficientes a inclinação e b interceptação do eixo y para o método dos mínimos quadrados. Vamos 
utilizar os mesmos dados do exemplo anterior, criando uma nova planilha que contenha apenas os 
dados, de acordo com a figura a seguir.
Figura 71 – Dados de entrada para o problema
Podemos adicionar as fórmulas a seguir, respectivamente, nas células E2 e E3:
Célula E2: =INCLINAÇÃO(B:B;A:A)
Célula E3: =INTERCEPÇÃO(B:B;A:A)
Com essas fórmulas, o Excel® calcula diretamente os coeficientes, conforme mostrado na figura na 
sequência. Note que os coeficientes são iguais aos encontrados anteriormente, o que valida a utilização 
desse método.
Figura 72 – Valores dos coeficientes a e b encontrados pelas funções nativas do Excel®
Finalmente, existe outra forma de encontrar esses coeficientes, baseando‑se na inserção de 
um gráfico e adicionando‑se a ele uma linha de tendência. Para isso, seleciona‑se os dados com 
o mouse, clica‑se na aba “Inserir” e procura‑se pelo gráfico do tipo “Dispersão”. A figura a seguir 
mostra esse passo.
258
Unidade II
Figura 73 – Inserção de um gráfico de dispersão a partir dos dados de entrada
O resultado disso será um gráfico de pontos representando cada par (x, y) dos dados observados. A 
figura na sequência mostra esse gráfico.
Figura 74 – Gráfico dos pontos gerados a partir dos dados observados
Ao clicar sobre a figura, um conjunto de abas “Ferramentas de Gráfico” deve aparecer. Entre estas 
encontra‑se a aba “Layout”, que possui o botão “Linha de Tendência”. Selecione a opção “Linha de 
Tendência Linear”, conforme a figura a seguir.
259
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Figura 75 – Inserir a linha de tendência linear aos pontos
Após esse passo, uma linha de tendência deverá aparecer no gráfico. Clicando‑se com o botão direito 
do mouse sobre essa linha, escolha a opção “Formatar Linha de Tendência” e selecione o check‑box 
“Exibir Equação no gráfico”. A figura a seguir mostra esse passo.
Figura 76 – Exibir a equação no gráfico
Após esse passo, a equação com seus coeficientes é apresentada no gráfico, como pode ser 
observado na figura anterior. Mais uma vez, nota‑se que os valores encontrados são os mesmos 
dos casos precedentes.
260
Unidade II
Exemplo de implementação computacional: método dos mínimos quadrados no Python
Essa seção mostra um exemplo de programação do método dos mínimos quadrados no Python. 
Para esse programa, vamos utilizar a forma matricial do método. Segundo essa formulação, temos de 
resolver a equação matricial dada pela equação a seguir.
1T T. . .g
−
 Θ = φ φ φ 
Onde:
Quadro 3
a
b
 
Θ =  
 
coeficientes a serem encontrados
1
2
n
x 1
x 1
x 1
 
 φ =  
  
valores de x dos dados observados
1
2
n
y
g y
y
 
 =  
  
valores de y dos dados observados
Dessa forma, devemos criar um código que disponha os dados observados no formato dessas 
matrizes e realize o cálculo matricial da equação 
1T T. . .g
−
 Θ = φ φ φ  . O código a seguir realiza esse 
procedimento, utilizando a biblioteca numpy.
Em especial, utiliza‑se a função np.transpose() a fim de transpor as matrizes, np.linalg.inv( ) para 
inverter matrizes e np.dot( ) para realizar a multiplicação matricial.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Método dos mínimos quadrados utilizando a formulação matricial
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Dados observados, dispostos nos vetores x e y
x = np.array([[1, 2, 3, 4, 5, 6]])
y = np.array([[123, 233, 358, 460, 581, 697]])
261
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
#######################################################################
## CÁLCULOS ##
#######################################################################
g = np.transpose(y)
phi = np.append(np.transpose(x), np.ones([np.size(x,1),1]), axis=1)
aux1 = np.dot(np.transpose(phi), phi)
aux2 = np.linalg.inv(aux1)
aux3 = np.dot(aux2, np.transpose(phi))
theta = np.dot(aux3, g)
#######################################################################
## PÓS‑PROCESSAMENTO ##
#######################################################################
# Mostra na tela a equação
print(‘y = ‘, theta[0,0], ‘.x + ‘, theta[1,0])
A variável theta conterá os valores de a em sua posição theta [0,0] e de b em theta[1,0]. A execução 
desse programa deverá exibir na tela a equação, conforme mostrado a seguir.
y = 114.742857143 .x + 7.06666666667 
Finalmente, vamos criar um programa utilizando funções pré‑programadas no Python para realizar 
o cálculo dos mínimos quadrados. Para isso, vamos usar a biblioteca scipy, que contém a função stats.
linregress. O código é apresentado a seguir.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Método dos mínimos quadrados utilizando a biblioteca scipy
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
from scipy import stats
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Dados observados, dispostos nos vetores x e y
x = np.array([1, 2, 3, 4, 5, 6])
y = np.array([123, 233, 358, 460, 581, 697])
#######################################################################
## CÁLCULOS ##
#######################################################################
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
#######################################################################
## PÓS‑PROCESSAMENTO ##
#######################################################################
# Mostra na tela a equação
print(‘y = ‘, slope, ‘.x + ‘, intercept)
Os resultados ficarão salvos nas variáveis slope e intercept. A execução desse programa gerará a 
seguinte saída:
y = 114.742857143 .x + 7.06666666667 
262
Unidade II
Neste ponto, sugere‑se implementar esse código em seu próprio computador e comparar os 
resultados com os exemplos obtidos anteriormente.
7 REGRESSÃO POLINOMIAL
Alguns conjuntos de dados não são representados de forma satisfatória por uma reta. Nesses casos, 
uma curva seria mais adequada para ajustar os dados. Uma alternativa é adequar os polinômios aos 
dados usando regressão polinomial.
Observe o gráfico de dispersão a seguir. Obviamente, os dados representados não seriam bem 
representados por uma reta.
Figura 77 
O procedimento dos mínimos quadrados pode ser prontamente estendido para ajustar dados por um 
polinômio de grau mais alto. Por exemplo, suponha que se queira ajustar um polinômio de segundo grau ou 
quadrático. O erro, ou resíduo, em função da adequação de uma curva aos dados isolados será denominado e:
2
0 1 2y a a x a x e= + + +
2
0 1 2e y a a x a x= − − −
Nesse caso, a soma dos quadrados dos resíduos é:
( )
2n
2
r i 0 1 i 2 i
i 1
S y a a x a x
=
= − − −∑
263
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Seguindo o procedimento do tópico anterior, obter a derivada da expressão em relação a cada um 
dos coeficientes desconhecidos do polinômio:
( )
n
2r
i 0 1 i 2 i
0 i 1
S
2( 1) y a a x a x
a =
∂
= − − − −
∂ ∑
( )
n
2r
i 0 1 i 2 i
0 i 1
S
2 y a a x a x
a =
∂
= − − − −
∂ ∑ (1)
( )
n
2r
i i 0 1 i 2 i
1 i 1
S
2( x ) y a a x a x
a =
∂
= − − − −
∂ ∑
( )
n
2r
i i 0 1 i 2 i
1 i 1
S
2 x y a a x a x
a =
∂
= − − − −
∂ ∑ (2)
( )
n
2 2r
i i 0 1 i 2 i
2 i 1
S
2( x ) y a a x a x
a =
∂
= − − − −
∂ ∑
( )
n
2 2r
i i 0 1 i 2 i
2 i 1
S
2 x y a a x a x
a =
∂
= − − − −
∂ ∑ (3)
As três equações são lineares com incógnitas a0, a1 e a2. Igualando as três expressões a zero para 
obter o valor crítico:
Equação (1)
( )
n
2
i 0 1 i 2 i
i 1
0 2 y a a x a x
=
= − − − −∑
( )
n
2
i 0 1 i 2 i
i 1
y a a x a x 0
=
− − − =∑
Distribuindo o somatório:
n n n n
2
i 0 1 i 2 i
i 1 i 1 i 1 i 1
y a a x a x 0
= = = =
− − − =∑ ∑ ∑ ∑
264
Unidade II
Como mostrado anteriormente:
n n
0 0 0
i 1 i 1
a a 1 na
= =
= =∑ ∑
Assim:
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
y na a x a x 0
= = =
− − − =∑ ∑ ∑
Como a0, a1 e a2 são constantes da equação polinomial, então:
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
y na a x a x 0
= = =
− − − =∑ ∑ ∑
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
y na a x a x
= = =
= + +∑ ∑ ∑ (4)
Equação (2)
( )
n
2r
i i 0 1 i 2 i
1 i 1
S
2 x y a a x a x
a =
∂
= − − − −
∂ ∑
( )
n
2
i i 0 1 i 2 i
i 1
x y a a x a x 0
=
− − − =∑
( )
n
2 3
i i 0 i 1 i 2 i
i 1
x y a x a x a x 0
=
− − − =∑
Distribuindo o somatório:
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x 0
= = = =
− − − =∑ ∑ ∑ ∑
Como a0, a1 e a2 são constantes da equação polinomial, então:
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x 0
= = = =
− − − =∑ ∑ ∑ ∑
265
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x
= = = =
= + +∑ ∑ ∑ ∑ (5)
Equação (3)
( )
n
2 2r
i i 0 1 i 2 i
2 i 1
S
2 x y a a x a x
a =
∂
= − − − −
∂ ∑
( )
n
2 2
i i 0 1 i 2 i
i 1
x y a a x a x 0
=
− − − =∑
( )
n
2 2 3 4
i i 0 i 1 i 2 i
i 1
x y a x a x a x 0
=
− − − =∑
Distribuindo o somatório:
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x 0
= = = =
− − − =∑ ∑ ∑ ∑
Como a0, a1 e a2 são constantes da equação polinomial, então:
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x 0
= = = =
− − − =∑ ∑ ∑ ∑
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x
= = = =
= + +∑ ∑ ∑ ∑
 (6)
Portanto, o resultado da manipulação algébrica leva a um sistema de três equações que possibilitam 
o cálculo de a0, a1 e a2 a partir dos dados.
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
y na a x a x
x y a x a x a x
x y a x a x a x
= = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
266
Unidade II
A análise para o caso bidimensional pode ser facilmente estendida para um caso mais geral, e é 
possível determinar os coeficientes de um polinômio de grau n de forma equivalente.
2 n
0 1 2 ny a a x a x ... a x e= + + + + +
Fazendo a mesma manipulação algébrica dos dois casos anteriores é possível chegar a um sistema 
de n equações em que a enésima equação será da forma:
n n n n
n n n 1 n 2
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x ...+ +
= = = =
= + + +∑ ∑ ∑ ∑
 
A planilha Excel® possui uma ferramenta útil na simulação da melhor curva interpoladora. 
Selecionar uma célula para a variável independente e outra para a variável dependente. A 
matriz coluna da variável independente na figura a seguir será constituída pelas linhas 
utilizadas. Obter o gráfico de dispersão como já orientado anteriormente e selecionar a aba 
“Layout” e, em seguida, “Linha de tendência”.
Figura 78 
Para casos em que as opções mostradas não parecerem apropriadas, selecione “Mais 
opções de linha de tendência”.
Figura 79 
 
267
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Exemplo 38
Determine o gráfico de dispersão e a curva interpoladora para o conjuntode dados indicados na 
tabela a seguir.
Tabela 45 
xi 0 1 2 3 4 5
yi 2,1 7,7 13,6 27,2 40,9 61,1
Resolução:
a) Gráfico de dispersão
Figura 80 
Tabela 46 – Tabela auxiliar
xi yi x
2
i x
3
i x
4
i xi yi x
2
i yi
0 2,1 0 0 0 0 0
1 7,7 1 1 1 7,7 7,7
2 13,6 4 8 16 27,2 54,4
3 27,2 9 27 81 81,6 244,8
4 40,9 16 64 256 163,6 654,4
5 61,1 25 125 625 305,5 1527,5
15 152,6 55 225 979 585,6 2488,8
268
Unidade II
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
y na a x a x
x y a x a x a x
x y a x a x a x
= = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
0 1 2
0 1 2
0 1 2
152,6 6a 15a 55a
585,6 15a 55a 225a
2488,8 55a 225a 979a
= + +
 = + +
 = + +
1 2 0
1 2 0
1 2 0
15a 55a 6a 152,6
55a 225a 15a 585,6
225a 979a 55a 2488,8
+ + =
 + + =
 + + =
Resolvendo o sistema de equações pela regra de Cramer:
15 15 6
D 55 225 15 3920
225 979 55
 
 = = 
  
a1
152,6 15 6
D 585,6 225 15 9248,4
2488,8 979 55
 
 = = 
  
a2
15 152,6 6
D 55 585,6 15 7294
225 2488,8 55
 
 = = 
  
a0
15 15 152,6
D 55 225 585,6 9716
225 979 2488,8
 
 = = 
  
Figura 81 
 
269
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
a1
1
D
a
D
=
 
a2
2
D
a
D
=
 
a3
3
D
a
D
=
1
9248,4
a 2,35929
3920
= =
2
7294
a 1,86071
3920
= =
0
9716
a 2,47857
3920
= =
Resolvendo o sistema de equações, obtém‑se a curva interpoladora:
2
1 2 0y a x a x a= + +
2y 1,86071x 2,35929x 2,47857= + +
A planilha do Excel® fornece uma função interpoladora: após a seleção do tipo de curva, selecionar 
também “Exibir equação no gráfico”.
Figura 82 
270
Unidade II
Figura 83 
 
A resolução de um sistema de equações pode ser bastante trabalhosa. A utilização da 
planilha Excel® para obter os determinantes das matrizes, necessários na regra de Cramer, 
facilita muito essa tarefa.
Para tanto, inserir a matriz cujo determinante se deseja obter; selecionar a célula na qual 
se deseja a resposta do determinante procurado (no caso do exemplo anterior, reproduzido 
a seguir, foram utilizadas as células A4, A9, A14, A19 para os valores dos determinantes), 
selecione a aba “Fórmulas”, em seguida “Matemática e trigonometria” e “Matriz.determ”.
Figura 84 
271
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Insira no quadro apresentado o campo da matriz, representado pela célula (primeira linha, primeira 
coluna) e a célula (última linha, última coluna). Selecionar “OK”.
Figura 85 
 
7.1 Regressão polinomial (para grau maior do que 2)
É importante observar que a montagem do novo sistema de equações em uma regressão polinomial 
para grau maior do que 2 não é apenas acréscimo de equações, porque a função interpoladora é um grau 
maior do que aquela (função quadrática) para a qual foi obtido o sistema inicial. É necessário observar 
que os expoentes de xi crescem tanto na horizontal quanto na vertical: os valores das potências iniciam 
com o valor zero e crescem em ambos os sentidos. Reescrevendo o sistema inicial e acrescentando o que 
foi omitido no desenvolvimento teórico, talvez seja possível visualizar melhor.
O sistema obtido anteriormente para função quadrática:
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
y na a x a x
x y a x a x a x
x y a x a x a x
= = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
Ao acrescentar o produto de alguns termos por 1, o sistema não sofrerá alteração (uma vez que 1 é 
o elemento neutro da multiplicação na operação usual em ℜ).
272
Unidade II
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
1y n1a a x a x
x y a x a x a x
x y a x a x a x
= = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
Porém:
0
i
n
0
i
i 1
x 1
x n.1
=
 =


=

∑
Substituindo:
n 1 n n
0 0 2
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
x y a x a x a x
x y a x a x a x
x y a x a x a x
= = = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
Valor da potência de xi 
cresce na descendente da 
vertical (cima para baixo) no 
primeiro termo da equação
Valor da potência de xi cresce 
da esquerda para a direita no 
segundo termo da equação
Assim, para uma função interpoladora de grau m é possível escrever:
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
y na a x a x
x y a x a x a x
x y a x a x a x
.....................................................
...............................
= = =
= = = =
= = = =
= + +
= + +
= + +
∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
n n n n
n n n 1 n 2
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
......................
x y a x a x a x ...+ +
= = = =














 = + + +

∑ ∑ ∑ ∑
273
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Exemplo 39
Ajuste uma equação cúbica aos dados da tabela. Utilizar coeficientes com cinco casas decimais.
Tabela 47 
xi 3 4 5 7 8 9 11 12
yi 1,6 3,6 4,4 3,4 2,2 2,8 3,8 4,6
Resolução:
a) Gráfico de dispersão:
Selecionar “Polinomial, ordem 3” e “Exibir equação no gráfico”.
Figura 86 
Figura 87 
274
Unidade II
b) Equação da curva interpoladora. Como a curva é de grau 3, é necessária mais uma equação, 
partindo da equação genérica:
2 3
0 1 2 3y a a x a x a x= + + +
y na a x a x a x
x y a x
i i
i
n
i
i
n
i
i
n
i
n
i i i
i
n
   
 
  

  

0 1
1
2
2
1
3
3
11
0
1
aa x a x a x
x y a x a x
i
i
n
i i
i
n
i
n
i
n
i i i
i
n
1
2
1
2
3
3
4
111
2
0
2
1
1
 
 
 

 
 ii
i
n
i
i
n
i
i
n
i
n
i i i
i
n
i
a x a x
x y a x a x
3
1
2
4
1
3
5
11
3
0
3
1
1
4
 
  
  

  

ii
n
i
i
n
i
i
n
i
n
a x a x
  
   













 1
2
5
1
3
6
11
Tabela 48 – Tabela auxiliar
xi yi x
2
i x
3
i x
4
i x
5
i x
6
i xi yi x
2
i yi x
3
i yi
3,00 1,60 9,00 27,00 81,00 243,00 729,00 4,80 14,40 43,20
4,00 3,60 16,00 64,00 256,00 1.024,00 4.096,00 14,40 57,60 230,40
5,00 4,40 25,00 125,00 625,00 3.125,00 15.625,00 22,00 110,00 550,00
7,00 3,40 49,00 343,00 2.401,00 16.807,00 117.649,00 23,80 166,60 1.166,20
8,00 2,20 64,00 512,00 4.096,00 32.768,00 262.144,00 17,60 140,80 1.126,40
9,00 2,80 81,00 729,00 6.561,00 59.049,00 531.441,00 25,20 226,80 2.041,20
11,00 3,80 121,00 1.331,00 14.641,00 161.051,00 1.771.561,00 41,80 459,80 5.057,80
12,00 4,60 144,00 1.728,00 20.736,00 248.832,00 2.985.984,00 55,20 662,40 7.948,80
59,00 26,40 509,00 4.859,00 49.397,00 522.899,00 5.689.229,00 204,80 1.838,40 18.164,00
Sistema de equações:
0 1 2 3
0 1 2 3
0 1 2 3
0 1 2 3
26,4 8a 59a 509a 4859a
204,8 59a 509a 4859a 49397a
1838,4 509a 4859a 49397a 522899a
18164 4859a 49397a 522899a 5689229a
= + + +
 = + + +
 = + + +
 = + + +
0 1 2 3
0 1 2 3
0 1 2 3
0 1 2 3
8a 59a 509a 4859a 26,4
59a 509a 4859a 49397a 204,8
509a 4859a 49397a 522899a 1838,4
4859a 49397a 522899a 5689229a 18164
+ + + =
 + + + =
 + + + =
 + + + =
275
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
8 59 509 4859
59 509 4859 49397
D 624678912,00156
509 4859 49397 522899
4859 49397 522899 5689229
 
 
 = =
 
 
 
ao
26,4 59 509 4859
204,8 509 4859 49397
D ‑7176753100,79993
1838,4 4859 49397 522899
18164 49397 522899 5689229
 
 
 = =
 
 
 
a1
8 26,4 509 4859
59 204,8 4859 49397
D 4462591967,99962
509 1838,4 49397 522899
4859 18164 522899 5689229
 
 
 = =
 
 
 
a2
8 59 26,4 4859
59 509 204,8 49397
D ‑650420006,399979
509 4859 1838,4 522899
4859 49397 18164 5689229
 
 
 = =
 
 
 
a38 59 509 26,4
59 509 4859 204,8
D 29157523,20000
509 4859 49397 1838,4
4859 49397 522899 18164
 
 
 = =
 
 
 
a0
0
D ‑7176753100,79993
a ‑11,48871
D 624678912,00156
= = =
a1
1
D 4462591967,99962
a 7,14382
D 624678912,00156
= = =
a2
2
D ‑650420006,399979
a ‑1,04121
D 624678912,00156
= = =
a3
3
D 29157523,20000
a 0,046676
D 624678912,00156
= = =
276
Unidade II
Retomando a equação genérica e substituindo os valores obtidos:
2 3
0 1 2 3y a a x a x a x= + + +
2 3y ‑11,48871 7,14382x‑1,04121x 0,046676x= + +
Logo, a equação da curva interpoladora de grau 3 é:
3 2y 0,046676x ‑1,04121x 7,14382x‑11,48871= +
7.2 Regressão não linear
Há muitos casos nos quais modelos não lineares devem ser ajustados aos dados. Esses modelos são 
definidos como aqueles que têm uma dependência não linear com seus parâmetros. Por exemplo, uma 
função do tipo:
a x1
0y a e=
Cuja representação gráfica é semelhante ao que mostra a figura a seguir:
Figura 88 
277
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Notas históricas 1
Ávila (2003) considera as funções exponenciais como as mais importantes da matemática. 
O rápido crescimento da função exponencial no infinito corresponde ao lento crescimento 
do logaritmo, fatos de grande importância em aplicações no mundo real e em domínios 
científicos. As duas funções estão ligadas à descrição de inúmeros fenômenos naturais, 
como juros compostos, crescimento populacional e desintegração radioativa. Estudos 
envolvendo as funções em crescimento populacional não são exclusivos de populações ou 
grupos humanos, mas estão intimamente ligados a pesquisas de muita importância na área 
da saúde, como é o caso do acompanhamento do crescimento (e decrescimento) de colônias 
de bactérias. A radioatividade foi descoberta no final do século XIX, mas só devidamente 
compreendida a partir do trabalho de Ernest Rutherford (1871‑1937) e seus colaboradores, 
que estudaram quantitativamente a radioatividade através de uma série de experiências 
que permitiram estabelecer a ligação desse fenômeno com a desintegração atômica.
Idades geológicas
Dos estudos de Ernest Rutherford e seus colaboradores surgiu uma constante 
de decaimento radioativo λ, denominada meia‑vida da substância radioativa. Ela é 
definida como o tempo T durante o qual a substância fica reduzida à metade de sua 
massa original M0:
T0
0
M
M e
2
−λ=
Aplicando logaritmo:
T0
0
M
ln lnM e
2
−λ=
T
0 0lnM ln2 lnM lne
−λ− = +
ln2 T− = −λ
ln2 T= λ
1
T ln2=
λ
278
Unidade II
Essa constante de decaimento radioativo λ, que envolve conhecimento de funções 
exponenciais e logarítmicas, possibilitou um método mais preciso de determinação de 
idades geológicas, que tem sido empregado com muito êxito desde o início do século XX, 
observa Ávila (2003). As substâncias radioativas utilizadas devem ter meia‑vida longa, da 
ordem de milhões ou bilhões de anos, como é o caso do urânio‑235 e do urânio‑238 para 
determinações de idades relacionadas com a vida do planeta. Estudos aplicados a rochas 
terrestres e minerais obtidos de meteoritos levam à conclusão de que a idade provável da 
crosta terrestre é de 4,5 bilhões de anos.
A imagem a seguir mostra crânios dos nossos antepassados: do mais antigo (o Australopithecus 
africanus, primeiro à esquerda) até o mais recente (o Homo sapiens, à frente de todos).
Figura 89 
Não só as ciências exatas e da saúde se beneficiam de conhecimentos relacionados 
à determinação de idades de acontecimentos, mas também as artes. A figura a seguir 
mostra desenhos pré‑históricos encontrados na caverna de Lascaux, na França, datados 
aproximadamente de 15000 a.C.
Figura 90 
279
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
A partir da década de 1940, foi desenvolvido um método apropriado à determinação 
de idades de importância histórica e arqueológica. William F. Libby (1908‑1980), na 
Universidade de Chicago, utilizou o isótopo de carbono, o carbono‑14, com vida média de 
aproximadamente 5.730 anos.
Pressão atmosférica
A pressão atmosférica num dado ponto da Terra é igual ao peso da coluna de ar disposta 
verticalmente acima de um elemento de área unitária, no ponto em que se deseja a pressão. 
A partir dessa definição, e de alguns conceitos da Química, como a lei de Boyle‑Mariotte, 
obteve‑se uma fórmula que permite determinar a altura de um lugar medindo‑se a pressão 
barométrica P, quando se conhece a aceleração da gravidade g e a pressão atmosférica da 
Terra em sua superfície, p0.
0pkh ln
g p
=
Também nesse caso, em estudos sobre pressão atmosférica, observa‑se a utilização de 
conhecimento sobre logaritmos, e no desenvolvimento da fórmula sobre exponencial.
Notas históricas 2
Os logaritmos foram inventados por John Napier perto do início do século XVII. 
Boyer (2001) considera que a invenção dos logaritmos deve datar aproximadamente de 
1594. Essa estimativa se deve ao fato de Napier ter afirmado que trabalhou na invenção 
dos logaritmos durante 20 anos antes de publicar seus resultados, em 1614, o Mirifi 
Logarithmorum Canonis Descriptio (uma descrição da maravilhosa regra dos logaritmos).
 
 Saiba mais
Para se aprofundar na história da Matemática, leia a obra a seguir:
BOYER, C. B. História da matemática. São Paulo: Edgard Blucher, 2001.
Como em casos dos mínimos quadrados em funções lineares, a regressão não linear é baseada 
na determinação dos valores dos parâmetros que minimizem a soma dos quadrados dos resíduos. 
A possibilidade em um caso assim é a linearização desta função:
a x1
0y a e=
280
Unidade II
Isso ocorre a partir do pressuposto de que uma função exponencial representa uma boa aproximação 
para os dados disponíveis.
Aplicando uma transformação logarítmica na função anterior:
a x1
0lny ln(a e )=
Usando as propriedades de logaritmo:
a x1
0lny lna lne= +
0 1lny lna a x= +
Chamando:
0
1
Y lny
B lna
A a
X x
=
 =
 =
 =
Obtém‑se a linearização da função anterior, a qual será representada por:
Y B AX= +
Ou:
Y AX B= +
O gráfico anterior seria assumido da seguinte forma:
Figura 91 
281
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Portanto, a partir dos dados da tabela, (xi, yi) e yi = f(x), obtêm‑se novos valores (Xi, Yi), tais que os 
valores de A e B da função Y = AX +B representem uma boa aproximação, ou seja, a verdadeira função 
descrita pelos novos pontos obtidos.
Logo:
a) Os valores A e B serão obtidos pela solução do sistema linear de equações:
n n n
i i
i 1 i 1 i 1
n n n
2
i i i i
i 1 i 1 i 1
Y A X B 1
X Y A X B X
= = =
= = =

= +


 = +

∑ ∑ ∑
∑ ∑ ∑
Ou:
n n
i i
i 1 i 1
n n n
2
i i i i
i 1 i 1 i 1
Y A X nB
X Y A X B X
= =
= = =

= +


 = +

∑ ∑
∑ ∑ ∑
b) A partir dos valores de A e B, obtêm‑se os valores a0 e a1 da função exponencial original 
a x1
0y a e= :
0B lna=
lnaB 0e e=
B
0e a=
Ou: 
B
0a e=
Como A = ai, então:
B
0
1
a e
a A
 =

=
282
Unidade II
Exemplo 40
O número de bactérias, por unidade de volume, existente em uma cultura após x horas é apresentado 
a seguir:
Tabela 49 
xi 0 1 2 3 4 5 6
yi 32 47 65 92 132 190 275
Sendo:
•	 xi número de horas;
•	 yi número de bactérias.
a) Obter o gráfico de dispersão.
b) Ajustar os dados da tabela anterior por uma equação exponencial que descreva o comportamento 
da cultura de bactérias por volume, isto é:
i iy f(x )=
c) Estimar o número de bactérias para 10 horas.
Resolução:
a) Gráfico de dispersão:
Figura 92 
283
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
b) Equação exponencial.
Assumiu‑se que:
0
1
Y lny
B lna
A a
X x
=
 =
 =
 =
Logo, a nova tabela será:
Tabela 50 
Xi = xi 0 1 2 3 4 5 6
Yi = Inyi ln32 ln47 ln65 ln92 ln132 ln190 ln275
Calculando‑se os valores, obtém‑se:
Tabela 51
Xi yi Yi xi yi x
2
i
0 32 3,465736 0 0
1 47 3,850148 3,850148 1
2 65 4,174387 8,348775 4
3 92 4,521789 13,56537 9
4 132 4,882802 19,53121 16
5 190 5,247024 26,23512 25
6 275 5,616771 33,70063 36
21 833 31,75866105,2312 91
n n
i i
i 1 i 1
n n n
2
i i i i
i 1 i 1 i 1
Y A X nB
X Y A X B X
= =
= = =

= +


 = +

∑ ∑
∑ ∑ ∑
31,75866 21A 7B
105,2312 A.91 21B
= +
 = +
284
Unidade II
Resolvendo o sistema linear:
21A 7B 31,75866
91A 21B 105,2312
+ =
 + =
Multiplica‑se a primeira equação por (‑3):
63A 21B 95,27598
91A 21B 105,2312
− − = −
 + =
Somando termo a termo:
28A 9,95522=
A 0,35554=
Substituindo:
21A 7B 31,75866+ =
21.0,35554 7B 31,75866+ =
B 3,47033=
Mas 
B
0
1
a e
A a
 =

=
Então: 
3,47033
0
1
a e 32,14735
a 0,35554
 = =

=
Como:
a x1
0y a e=
Então:
0,35554xy 32,14735e=
0,35554xf(x) 32,14735e=
285
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
c) O número de bactérias para 10 horas será:
0,35554xf(x) 32,14735e=
0,35554.10f(10) 32,14735e 32,14735.35,00182= =
f(10) 1.125,22=
Ou seja, após 10 horas há 1.125 bactérias.
 
A planilha Excel® fornece uma expressão aproximada para a função da curva exponencial 
interpoladora. Para tanto, repetem‑se os procedimentos descritos durante o texto para 
inserção dos dados e obtenção da curva de dispersão.
A seguir, seleciona‑se a aba “Layout” e, em seguida, “Linha de tendência”. Marcar “Exponencial” 
e “Exibir equação no gráfico”.
Figura 93 
286
Unidade II
Figura 94 
Inserindo‑se o número de bactérias em 10 horas, o gráfico que conserva a curva e 
a equação apresentará uma modificação na terceira casa de um dos coeficientes por 
arredondamento, em decorrência de novos cálculos.
Figura 95 
 
7.3 Ajuste polinomial: implementação computacional
Será apresentada uma aplicação do ajuste polinomial para dados que não se encaixem em 
funções de primeiro grau. No exemplo a seguir, veremos que o método dos mínimos quadrados 
pode ser estendido para polinômios de mais alta ordem. Nesse exemplo, vamos considerar os dados 
apresentados na tabela na sequência.
287
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Tabela 52 – Dados para o exemplo de ajuste polinomial
X Y
‑4,90 48,73
‑3,21 21,89
‑1,93 10,45
‑0,53 4,40
1,33 0,55
2,44 3,03
3,62 11,92
5,18 12,44
6,64 25,22
7,92 46,99
A figura na sequência mostra esses pontos em um plano cartesiano. Dessa figura podemos concluir 
que um ajuste linear não representará bem a relação entre x e y.
Figura 96 
Podemos tentar realizar um ajuste polinomial de grau 2, ou seja, ajustar esses pontos para uma 
parábola ( 2y a.x b.x c= + + ). Nesse caso, o sistema linear a ser resolvido, de acordo com o método dos 
mínimos quadrados, é dado pela equação a seguir.
2
i i i
2 3
i i i i i
2 3 4 2
i i i i i
c.n b. x a. x y
c. x b. x a. x x .y
c. x b. x a. x x .y
 + ∑ + ∑ = ∑
 ∑ + ∑ + ∑ = ∑

∑ + ∑ + ∑ = ∑
Nota‑se que a construção do sistema de equações é semelhante ao caso do ajuste linear. No entanto, 
termos adicionais aparecem. Será preciso realizar o somatório para mais termos. A tabela a seguir mostra 
as contas necessárias para encontrar esses somatórios.
288
Unidade II
Tabela 53 – Tabela auxiliar para o cálculo dos somatórios
x y x2 x3 x4 x . y x2 . y
‑4,90 48,73 24,04 ‑117,9 577,96 ‑238,9 1.171,56
‑3,21 21,89 10,29 ‑33,00 105,84 ‑70,22 225,24
‑1,93 10,45 3,71 ‑7,15 13,78 ‑20,14 38,80
‑0,53 4,40 0,29 ‑0,15 0,08 ‑2,35 1,26
1,33 0,55 1,78 2,38 3,17 0,73 0,97
2,44 3,03 5,95 14,52 35,42 7,40 18,05
3,62 11,92 13,13 47,56 172,30 43,18 156,43
5,18 12,44 26,86 139,18 721,30 64,49 334,23
6,64 25,22 44,04 292,30 1939,90 167,38 1110,80
Soma 16,56 185,63 192,74 833,72 7495,58 323,48 6001,63
Dessa forma, o sistema linear a ser resolvido será dado pela equação a seguir:
10.c 16,56.b 192,74.a 185,63
16,56.c 192,74.b 833,72.a 323,48
192,74.c 833,72.b 7495,58.a 6001,63
+ + =
 + + =
 + + =
Ou, na forma matricial, conforme equação na sequência.
10 16,56 192,74 c 185,63
A.x b 16,56 192,74 833,72 . b 323,48
192,74 833,72 7495,58 a 6001,63
     
     = → =     
          


Esse sistema deve ser resolvido de modo a encontrar os parâmetros a, b e c da equação de segundo 
grau que representa a relação entre os pontos. Podemos utilizar qualquer dos métodos estudados 
anteriormente para solucionar esse sistema. Vamos utilizar o método de eliminação de Gauss, por ser 
um método direto, e apresentar a solução exata do problema. A matriz estendida, para esse caso, é 
representada pela equação a seguir.
10 16,56 192,74 185,63
16,56 192,74 833,72 323,48
192,74 833,72 7495,58 6001,63
 
 
 
  
O processo segue através de transformações lineares nas linhas, de modo a transformar a parte 
da esquerda da matriz estendida na matriz identidade. Os passos desse processo estão descritos na 
sequência. As transformações aplicadas estão descritas no lado esquerdo de cada linha.
289
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
L1
L2 x 10
L1
16,56
L3 x 10
L1
192,74
−
−
10 16,56 192,74 185,63
0 99,83 310,72 9,71
0 26,70 196,15 125,75
 
 
 
  
L1
L2
L3 x 99,83
L2
26,70
−
10 16,56 192,74 185,63
0 99,83 310,72 9,71
0 0 422,80 460,54
 
 
 
  
L1
L2
L3
422,80
10 16,56 192,74 185,63
0 99,83 310,72 9,71
0 0 1 1,09
 
 
 
  
L1
L3
192,74
L2
L3
310,72
L3
−
−
0,05 0,09 0 0,13
0 0,32 0 1,06
0 0 1 1,09
 −
 − 
  
L1
L2
0,32
L3
0,05 0,09 0 0,13
0 1 0 3,29
0 0 1 1,09
 −
 − 
  
L1
L2
0,09
L2
L3
− 0,60 0 0 1,82
0 1 0 3,29
0 0 1 1,09
 
 − 
  
L1
0,6
L2
L3
1 0 0 3,02
0 1 0 3,29
0 0 1 1,09
 
 − 
  
 
290
Unidade II
8 INTEGRAÇÃO NUMÉRICA
8.1 Introdução
O problema geral da integração numérica é encontrar um valor aproximado para uma integral definida:
b
a
f(x)dx∫
Mas por que se preocupar em encontrar um método numérico para aproximar o valor da integral se 
já dispomos de métodos analíticos, como o Teorema Fundamental do Cálculo, que, em princípio, permite 
calcular o valor exato da integral?
O cálculo de uma integral por métodos analíticos, usando o Teorema Fundamental do Cálculo, pode 
ser uma tarefa extremamente complicada e, em muitas situações, ser mesmo impossível encontrar uma 
expressão que possa ser expressa como uma combinação finita de funções algébricas, logarítmicas e 
exponenciais. A integral a seguir, que desempenha um importante papel em Estatística, é um exemplo 
bem conhecido de uma integral desse tipo:
21 x
0
e dx−∫
Outra situação, bastante comum, ocorre quando a função f que deve ser integrada representa um 
conjunto de dados obtidos a partir de uma amostra. Nessa situação se conhece apenas um número 
finito de valores da função. Como acontece geralmente, não temos uma expressão analítica para f, e 
a única possibilidade é o cálculo da integral por métodos numéricos, pois esses métodos se utilizam 
apenas de um número finito de valores de f. Antes de estudar os métodos numéricos para o cálculo das 
integrais, será útil recordar a definição de integral, pois os métodos numéricos que serão desenvolvidos 
aqui partem dessa definição, que, por sua vez, tem como base o cálculo de áreas por aproximações. 
Por esse motivo, as fórmulas de integração numérica são chamadas também, por razões históricas, de 
quadratura numérica, pois foi Arquimedes (287‑212 a.C.) que, na tentativa de resolver o problema da 
quadratura do círculo, usou um método semelhante aos de integração numérica.
8.2 Integrais definidas e áreas: a base dos métodos de integração numérica
Se f(x) é uma função definida em um intervalo [a, b], então a integral definida
( )b
a
f x dx∫
calcula a área da região entre o eixo x e o gráfico de f no intervalo [a, b]. Na parte positiva de f a área 
é contada como área positiva e, na parte negativa de f, a área é contada como área negativa:
291
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Figura 97 
O cálculo dessa área é feito por aproximações usando retângulos, e o valor exato é definido por meio 
de um limite. Mais precisamente, escolhemos inicialmente uma partição do intervalo [a, b],ou seja, um 
número finito de pontos 0 1 nx , x , , x… entre a e b:
0 1 2 n 1 na x x x x x b−= < < <…< < =
Os pontos escolhidos podem não ser igualmente espaçados, e inclusive há métodos numéricos que se 
utilizam de maneira vantajosa de pontos que não são igualmente espaçados. Sem perda de generalidade, 
vamos escolher pontos igualmente espaçados. Nesta situação, a partição divide o intervalo [a, b] em n 
subintervalos de mesmo comprimento:
[ ] [ ] [ ]0 1 1 2 n 1 nx ,x , x ,x , , x ,x−…
Uma tal partição é chamada de partição uniforme. Como a largura do intervalo é b ‑ a, a largura de 
cada um dos n subintervalos é:
b a
x
n
−
∆ =
Se mi é o valor mínimo da função f no intervalo [xi‑1, xi], e Mi é o valor máximo da função f no intervalo 
[xi‑1, xi], então a área pode ser aproximada tanto pela soma das áreas dos n retângulos inferiores quanto 
pela soma das áreas dos n retângulos superiores:
n n
i i
i 0 i 0
 I m x ou I M x
= =
= ∆ = ∆∑ ∑
A aproximação da área melhora se tomamos partições com um número cada vez maior de pontos:
292
Unidade II
Figura 98 
Formalmente, fazendo n → ∞ em ambas as somas, se os limites existirem e forem iguais, definiremos 
a integral I como o valor desses limites:
n n
I lim I lim I
→∞ →∞
= =
Podemos, a partir dessa definição de integral, já construir um método numérico simples para o 
cálculo aproximado de integrais. O único inconveniente é que precisamos encontrar os valores máximos 
e mínimos da função f em cada intervalo, o que pode ser um problema. Mas, felizmente, isso pode ser 
evitado se usarmos uma outra propriedade de integrais, demonstrada em cursos avançados de cálculo, 
que diz o seguinte: se a integral existe, então, para qualquer ponto ti no intervalo [xi‑1, xi], tem‑se:
( )
n
i
n i 1
I lim f t x
→∞ =
= ∆∑
Podemos escolher ti como o ponto médio do intervalo [xi‑1, xi], ou seja,
i i 1
i
x x
t
2
−−=
Temos então a seguinte aproximação, chamada de regra do ponto médio:
n
i i 1
i 1
x x
I f x
2
−
=
− ≈ ∆  ∑
Ainda é cedo para tentarmos implementar um algoritmo razoável para o método do ponto médio, 
pois ainda não sabemos qual a precisão desse método, que evidentemente depende do número n de 
subintervalos escolhidos e, de forma mais sutil, de como a função se afasta do valor médio i i 1
x x
f
2
−− 
  em cada intervalo.
293
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Vale a pena nos determos um pouco mais sobre esta última observação: o método do ponto médio 
aproxima a área sob a curva do gráfico da função em cada subintervalo [xi‑1, xi] por meio de um retângulo 
de altura i i 1
x x
f
2
−− 
   :
Figura 99 
Em outras palavras, para cada subintervalo [xi‑1, xi], estamos interpolando a curva do gráfico de f em 
[xi‑1, xi] pela reta horizontal 
i i 1x xy f
2
−− =   
.
Se interpolarmos cada parte do gráfico, em cada subintervalo, por outra curva que se aproxima 
melhor do gráfico da função e, além disso, delimita uma área fácil de calcular, será razoável esperar que 
possamos melhorar nossa aproximação.
Este é o fundamento dos outros dois métodos que estudaremos a seguir:
O método dos trapézios, que interpola cada trecho do gráfico definido no subintervalo [xi‑1, xi] por 
uma reta secante (um polinômio de primeiro grau), definida a partir dos valores dos extremos do intervalo.
O método de Simpson, que interpola cada trecho do gráfico definido no subintervalo [xi‑1, xi] por 
uma parábola (um polinômio de segundo grau).
8.3 Regra do ponto médio
A partir da discussão precedente, deduzimos um método para aproximar o valor de uma integral 
definida, que pode ser enunciado da seguinte maneira:
Regra do ponto médio
Seja f uma função integrável no intervalo [a, b] e 0 1 nx x ,L,x, uma partição uniforme desse intervalo. 
Então, a integral de f pode ser aproximada por:
294
Unidade II
( )
nb
M ia
i 1
f x dx I x f(x )
=
≈ = ∆ ∑∫
Em que 
b a
x
n
−
∆ = e i i 1i
x x
x
2
−−= = ponto médio de [xi‑1, xi].
Em termos práticos:
( )i
x 1
x a i 1 x a i x
2 2
∆   = + + − ∆ = + − ∆   
   
Vamos usar a regra do ponto médio para calcular uma integral conhecida, pois, se compararmos com 
o valor exato da integral, poderemos ter uma ideia da precisão do método.
Exemplo 41
Calcule, usando a regra do ponto médio, com n = 5, a integral 
1
0
senxdx∫ .
Nesse caso, 
1
x
5
∆ = e a = 0.
Resolução:
1
0
1 1 3 5 7 9
senxdx sen sen sen sen sen
5 10 10 10 10 10
 ≈ + + + +  ∫ ≈ 0,460465 
O valor exato dessa integral é cos0 cos1 0,459698− ≈ , o que mostra uma discrepância de 0,000767, ou 
seja, um erro menor do que 0,2%. Vamos, então, aumentar gradativamente o valor de n e observar os resultados:
Tabela 54 
n IM Erro = |I - IM|
5 0,460464751755509 7,671E‑04
10 0,459889290718518 1,916E‑04
20 0,459745582800190 4,789E‑05
40 0,459709665644209 1,197E‑05
80 0,459700686969029 2,993E‑06
160 0,459698442338595 7,482E‑07
Podemos observar que, a cada vez que o valor de n é duplicado, ou seja, o tamanho de cada intervalo 
[xi‑1, xi] é reduzido pela metade, então o erro, que é a diferença em módulo entre o valor exato da integral I e 
o valor aproximado pela regra do ponto médio IM, é reduzido aproximadamente por um fator de 4. Isto não 
é uma coincidência, nem se restringe a esta integral específica. O seguinte resultado, que é demonstrado em 
livros mais avançados de cálculo numérico, mostra que esse é de fato o comportamento do erro na regra do 
ponto médio.
295
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Observe que se n aumenta por um fator de 2, ou seja, n é duplicado, então o erro se reduz por um 
fator de 22, confirmando o que pode ser observado na tabela a seguir.
Tabela 55 
n IM Erro = |I - IM| (Erro anterior)/4
5 0,460464751755509 7,671E‑04 –
10 0,459889290718518 1,916E‑04 1,92E‑04
20 0,459745582800190 4,789E‑05 4,79E‑05
40 0,459709665644209 1,197E‑05 1,20E‑05
80 0,459700686969029 2,993E‑06 2,99E‑06
160 0,459698442338595 7,482E‑07 7,48E‑07
Teorema do erro para a regra do ponto médio
Seja f uma função com primeira e segunda derivadas contínuas num intervalo [a, b]. Suponha que 
( )f '' x D≤ para a x b. ≤ ≤ Se IM é a aproximação da integral ( )
b
a
I f x dx = ∫ pela regra do ponto médio 
com subintervalos, então o erro dessa aproximação é limitado por:
( )3
M 2
b a7
erro= I I D
24 n
−
− ≤
Em geral, não é muito simples estimar o erro, pois precisamos estimar D, ou seja, o valor máximo da 
segunda derivada de f no intervalo [a, b]. Nos casos em que isso é possível, podemos encontrar o número 
n de subintervalos necessários para garantir um erro máximo predeterminado. Vejamos o exemplo.
Exemplo 42
Qual deve ser o valor de n para garantir que o erro no cálculo da integral 
3
21
1
dx
x∫ pela regra do ponto médio seja menor do que 0,0001?
Resolução:
Neste caso, queremos mI I 0,0001− < . Para garantir isso, basta tomar:
( )3
2
3 17
D 0,0001
24 n
−
<
Resolvendo a desigualdade para n, temos:
( )32 3 17n D
24 0,0001
−
>
296
Unidade II
37 2
n D
24 0,0001
>
Precisamos estimar M0, o que não é muito difícil, pois a segunda derivada de ( ) 2
1
f x
x
= é igual a 
( ) 4
6
f " x
x
= , que por sua vez é uma função decrescente. Sendo assim, no intervalo [1, 3], o valor máximo 
de f”(x) ocorre quando x = 1, ou seja, 4
6
D 6
1
= = . Temos então:
37 2
n 6 374,17
24 0,0001
> ≈
Então, n = 375 irá garantir a precisão desejada.
Como implementar um algoritmo para a regra do ponto médio? O problema desse algoritmo é que 
ele depende da escolha de n, ou seja, do número de pontos médios que serão tomados no intervalo de 
integração. Se fixarmos arbitrariamente um valor para n, não teremos certeza da precisão dos resultados. 
O exemplo anterior mostra que podemos usar o teorema sobre o erro para escolher n a partir de uma 
precisão desejada. Poderíamos usar esta ideia para construir um algoritmo que calculasse as integrais 
pela regra do ponto médio com uma precisão preestabelecidapelo usuário. Mas isso não pode ser 
implementado tão facilmente, pois o valor de n depende de uma estimativa de D
Há uma outra maneira de controlarmos a precisão dos cálculos: a ideia é começar com um valor 
inicial para n e escolher valores cada vez maiores e, para cada escolha de n, calcular IM. O algoritmo para 
quando a diferença dos valores obtidos para IM entre duas escolhas sucessivas de n é menor do que uma 
tolerância preestabelecida, que pode ser o erro relativo máximo permitido de uma aproximação IM para 
outra. Este será o padrão geral dos algoritmos de integração numérica.
Mais precisamente, se IM(k) é o valor obtido pela regra do ponto médio quando n = k, n0 é o número 
de subintervalos inicial e ∈ é o erro relativo máximo permitido. Então, podemos tomar as aproximações:
( ) ( ) ( )m 0 m 1 m kI n , I n , ,I n ,… …
Fazendo nk = 2
kn0. Então, para cada k > 1, testamos se:
( ) ( )
( )
m k m k 1
m k 1
I n I n
I n
−
−
−
<∈
Ou, melhor ainda (evitando eventuais divisões por zero), testamos se:
( ) ( ) ( )m k m k 1 m k 1I n I n I n− −− <∈
297
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
O algoritmo para quando a desigualdade anterior é satisfeita. Evidentemente, também é necessário 
limitar k.
8.4 Regra do trapézio
Vimos que, na regra do ponto médio, em cada um dos subintervalos [ ] [ ] [ ]0 1 1 2 n 1 nx ,x , x ,x , , x ,x−… , 
aproximamos a curva do gráfico da função por uma reta horizontal, ou seja, pela função constante 
i i 1x xf
2
−− 
  
, o que faz com que a área sob a curva em cada subintervalo [xi‑1, xi] seja aproximada 
por um retângulo.
Uma outra maneira é aproximar cada trecho da curva do gráfico, em cada subintervalo, por uma 
reta secante que interpola o gráfico nos pontos extremos do subintervalo. A figura a seguir mostra esta 
situação para um caso com quatro subintervalos:
Figura 100 
O nome regra do trapézio vem do fato de que a área é aproximada por trapézios.
Sendo assim, para cada subintervalo [xi‑1, xi], a área do trapézio correspondente é:
( ) ( )i 1 if x f x x
2
− + ∆  
Em que ∆x = |xi ‑ xi‑1|. Tomando n subintervalos de igual comprimento, teremos, como no caso da regra 
 
do ponto médio, 
b a
x
n
−
∆ = . A integral é aproximada pela soma das áreas de todos os trapézios:
( ) ( ) ( )
nb i 1 i
a
i 1
f x f x
f x dx x
2
−
=
 +
≈ ∆  
∑∫
298
Unidade II
( ) ( )
n n
i 1 i
i 1 i 1
x
f x f x
2 += =
 ∆
= + 
 
∑ ∑
( ) ( ) ( ) ( ) ( ) ( )( )
( ) ( ) ( ) ( ) ( )( )
( ) ( ) ( )
( ) ( ) ( )
0 1 n 1 1 2 n
0 1 2 n 1 n
n 1
0 n
i
i 1
n 1
i
i 1
x
f x f x ... f x f x f x ... f x
2
x
f x 2f x 2f x ... 2f x f x
2
f x f x
f x x
2
f a f b
f x x
2
−
−
−
=
−
=
∆
= + + + + + + +
∆
= + + + + +
 +
= + ∆ 
 
 +
= + ∆ 
 
∑
∑
Temos então, a partir dessa discussão, o seguinte:
Regra do trapézio
Seja f uma função integrável no intervalo [a, b] e 0 1 nx x ,L,x, uma partição uniforme desse intervalo. 
Então a integral de f pode ser aproximada por:
( ) ( ) ( ) ( )
n 1b
T ia
i 1
f a f b
f x dx I f x x
2
−
=
 +
≈ = + ∆  
 
∑∫
Em que 
b a
x
n
−
∆ = e xi = a + i∆x.
Antes de analisarmos o erro envolvido nessa aproximação, vamos fazer um exemplo usando a regra 
do trapézio.
Exemplo 43
Aproxime a integral 
x2
1
e dx∫ com n = 6.
Resolução:
Com 
2 1 1
a 1, b 2, x
6 6
−
= = ∆ = = , i
i
x 1
6
= + , a regra do trapézio resulta em:
299
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
( ) ( ) ( )
5x2
T i1
i 1
f 1 f 2 1
e dx I f x
2 6=
 +
≈ = + 
 
∑∫
( ) ( )
7 8 9 10 111 2
6 6 6 6 6
f 1 f 2 7 8 9 10 11 1
f f f f f
2 6 6 6 6 6 6
e e 1
e e e e e
2 6
4,68158125
 +          = + + + + +                    
 +
= + + + + +   
≈
O valor exato dessa integral é igual a 2e e 4,67077427− ≈ . O erro obtido pela regra do trapézio, em 
módulo, é de aproximadamente 0,0108038.
Se calcularmos esta mesma integral usando a regra do ponto médio com n = 6, vamos obter 
MI 4,66537266≈ . O erro, neste caso, é de aproximadamente 0,00540161, ou seja, cerca de metade do 
erro obtido na regra do trapézio.
Podemos esperar mais precisão na regra do ponto médio do que na regra do trapézio, como mostra 
a figura a seguir, que compara as áreas obtidas nos dois métodos.
Figura 101 
Na figura, a área do retângulo AB’C’D, obtido pela regra do ponto médio, é equivalente à área do 
trapézio ABCD, cujo lado superior é tangente ao gráfico em P. Por outro lado, o trapézio AQRD aproxima 
a área pela regra do trapézio. Quando comparamos a área desses dois trapézios com a área sob o gráfico, 
podemos notar que a diferença entre a área de ABCD e a área sob o gráfico (em vermelho) é menor 
do que a diferença entre a área de AQRD e a área sob o gráfico (em azul). A área em vermelho é o erro 
obtido na regra do ponto médio, enquanto a área em azul é o erro obtido na regra do trapézio.
300
Unidade II
O seguinte teorema estima o erro obtido na regra do trapézio e confirma nossas observações quando 
comparamos com o erro na regra do ponto médio.
Teorema do erro para a regra do trapézio
Seja f uma função com primeira e segunda derivadas contínuas num intervalo [a, b]. Suponha que 
( )f '' x D≤ para a x b. ≤ ≤ Se IT é a aproximação da integral ( )
b
a
I f x dx = ∫ pela regra do trapézio com 
n subintervalos, então o erro dessa aproximação é limitado por:
( )3
T 2
b a5
erro= I I D
12 n
−
− ≤
8.5 Regra de Simpson
A regra de Simpson se baseia em aproximar a função f, em cada subintervalo, por um polinômio de 
grau 2, ou seja, em cada subintervalo o gráfico de f é aproximado por uma parábola.
A dedução da fórmula da regra de Simpson é um pouco mais elaborada do que nos métodos anteriores. 
Por isso, vamos começar a análise com uma situação simples e depois generalizar os resultados para a 
situação mais geral.
Vamos considerar inicialmente a seguinte situação: queremos aproximar a curva y = f(x) em 
vermelho na figura) por uma parábola (em azul) que passa por ( ) ( ) ( )0 0 1 1 2 2P h,y , P 0,y e P h,y− , 
conforme a figura a seguir.
Figura 102 
Se y = Ax2 + Bx + C é a equação dessa parábola, então a área sob a curva de f pode ser aproximada 
pela área sob a parábola:
( )h h 2
h h
f x dx Ax Bx C dx
− −
≈ + +∫ ∫
301
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
A integral do lado direito não é difícil de calcular:
h3 2h 2
h
h
Ax Bx
 Ax Bx C dx Cx
3 2−
−
+ + = + +∫
( )
3 2 3 2
3
3
Ah Bh Ah Bx
Ch Ch
3 2 3 2
2Ah h
2Ch 2Ah 6C
3 3
   
= + + − + −      
   
= + = +
Então, a equação ( )h h 2
h h
f x dx Ax Bx C dx
− −
≈ + +∫ ∫ pode ser escrita como:
( ) ( )h 2h
h
f x dx 2Ah 6C
3−
≈ +∫
Agora, precisamos encontrar qual é a relação de A e C com ( )0 0P h,y , − P1 (0,y1) e ( )2 2P h,y . Isso 
pode ser feito lembrando que a parábola passa por P0, P1 e P2, o que significa que:
( ) ( )2 20y A h B h C Ah Bh C= − + − + = − +
( ) ( )21y A 0 B 0 C C= + + =
( ) ( )2 22y A h B h C Ah Bh C= + + = + +
A partir dessas equações, deduzimos que y1 = C e y0 + y2 = 2Ah
2 + 2C. Com isso, podemos reescrever 
nossa aproximação ( ) ( )h 2h
h
f x dx 2Ah 6C
3−
≈ +∫ da seguinte maneira:
( ) ( ) ( )
( )
h 2
0 2 1h
0 1 2
h h
f x dx 2Ah 6C y y 4y
3 3
h
y 4y y
3
−
≈ + = +
= + +
+∫
Observe que f(‑h) = y0, f(0) = y1 e f(h) = y2, pois os pontos P0, P1 e P2 também estão na curva y = f(x). 
Finalmente, podemos escrever a fórmula ( )0 1 2
h
y 4y y 
3
= + + como:
( ) ( ) ( ) ( )h
h
h
f x dx f h 4f 0 f h
3−
 ≈ − + + ∫
302
Unidade II
A fórmula ( ) ( ) ( ) ( )h
h
h
f x dx f h 4f 0 f h
3−
 ≈ − + + ∫ nos diz como aproximar a área sob o gráfico de 
f no intervalo [‑h, h] pela área sob a parábola que passa por P0, P1 e P2. Vamos agora adaptá‑la para a 
situação mais geral. Para fazer isso, observamos o seguinte: se deslocarmos a parábola na horizontal, a 
área não se alterará:
Figura 103 
Além disso, f(‑h) = f(xi), f(0) = f(xi+1) e f(h) = f(xi+2). Combinandoestes dois fatos, temos o seguinte:
( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
x hi 2
x hi
i i 1 i 2
h
f x dx f x dx f h 4f 0 f h
3
x
f x 4f x f x
3
+
−
+ +
 = ≈ − + + = 
∆
 = + + 
∫ ∫
Em que ∆x = xi+1 ‑ xi.
Finalmente, a aproximação para a integral sobre o intervalo total [a, b] pela regra de Simpson é 
construída de maneira análoga aos outros dois métodos que estudamos: o intervalo [a, b] é particionado 
em subintervalos [xi, xi+1] de mesmo comprimento ∆x e, para cada par de subintervalos consecutivos, 
[xi, xi+1] e [xi+1, xi+2], computamos a integral ( ) ( ) ( )i i 1 i 2
x
f x 4f x f x
3 + +
∆
 = + +  , pois cada parábola é 
definida a partir dos três pontos xi, xi+1 e xi+2. Se adicionarmos a soma de todas estas áreas, obteremos 
a aproximação pela regra de Simpson para a integral definida sobre o intervalo total [a, b]. Note que o 
número n de subintervalos deve ser par.
( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
b x x x2 4 n
a x x x0 2 n 2
0 1 2 2 3 4
n 2 n 1 n
f x dx f x dx f x dx f x dx
x x
f x 4f x f x f x 4f x f x
3 3
x
... f x 4f x f x
3
−
− −
≈ + +…+
∆ ∆   = + + + + +   
∆
 + + + 
∫ ∫ ∫ ∫
303
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Podemos melhorar a fórmula anterior observando que os termos com índices pares, f(x2), f(x4), f(x6)..., 
aparecem duas vezes cada um, enquanto os termos com índices ímpares, f(x1), f(x3), f(x7)..., aparecem 
multiplicados por 4:
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )b 0 1 2 3 4 n 2 n 1 na
x
f x dx f x 4f x 2f x 4f x 2f x 2f x 4f x f x 
3 − −
∆  ≈ + + + + +…+ + + ∫
A fórmula anterior é conhecida como regra de Simpson. Podemos reescrevê‑la de uma maneira 
mais compacta:
Regra de Simpson
Seja f uma função integrável no intervalo [a, b] e 0 1 nx x ,L,x, uma partição uniforme desse intervalo, 
com n par. Então, a integral de f pode ser aproximada por:
( ) ( ) ( ) ( ) ( )
n n
1
2 2b
S 2i 1 2ia
i 1 i 1
x
f x dx I f a f b 4 f x 2 f x
3
−
−
= =
 
∆  ≈ = + + + 
  
∑ ∑∫
Em que 
b a
x
n
−
∆ = e xi = a + i∆x.
Vamos usar a regra de Simpson para calcular 
x2
1
e dx∫ . Já havíamos realizado o cálculo dessa mesma 
integral usando a regra dos trapézios. Nosso objetivo ao usar o mesmo exemplo é comparar os resultados 
e/ou desempenhos entre os métodos.
Exemplo 44
Calcule a integral 
x2
1
e dx∫ pela regra de Simpson com n = 6. 
Resolução:
Tomando f(x) = ex, n = 6 e ∆x = 1/6 na fórmula:
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )b 0 1 2 3 4 n 2 n 1 na
x
f x dx f x 4f x 2f x 4f x 2f x 2f x 4f x f x 
3 − −
∆  ≈ + + + + +…+ + + ∫
304
Unidade II
Teremos:
( ) ( )
x2
S1
7 8 9 10 11
2 26 6 6 6 6
1 7 8 9 10 11
e dx I f 1 4f 2f 4f 2f 4f f 2
18 6 6 6 6 6
1
e 4e 2e 4e 2e 4e e e e
18
4,67079423
          ≈ = + + + + + +                    
 
 = + + + + + + − +
  
≈
∫
Esta aproximação é muito melhor do que as obtidas pelas regras do ponto médio e do trapézio com 
o mesmo n = 6.
Tabela 56 
Método Valor obtido ≈ |Erro|
IM 4,66537266 0,0054016
IT 4,68158125 0,0108038
IS 4,67079423 0,00002
Isto se confirma pelo seguinte resultado:
Teorema do erro para a regra de Simpson
Seja f uma função com as primeiras quatro derivadas contínuas num intervalo [a, b]. Suponha que 
( ) ( )4f x D≤ para a x b. ≤ ≤ Se IS é a aproximação da integral ( )
b
a
I f x dx = ∫ pela regra de Simpson com 
n subintervalos, n par, então o erro dessa aproximação é limitado por:
( )
S
5
4
b a49
erro= I I D
2880 n
−
− ≤
Para ilustrar a aplicação da técnica de integração numérica, vamos nos apoiar no seguinte problema:
Um designer em uma fábrica de porcelanas deseja criar um copo com o formato mostrado na 
figura seguinte.
305
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Figura 104 – Design do copo
O formato desse copo é gerado a partir da revolução da linha, mostrada na figura, ao redor do 
eixo y = 0.
5
4
3
2
1
0
Altura (cm)
0 2 4 6 8 10 12
Ra
io
 (c
m
)
Figura 105 – Perfil do design do copo
É possível modelar matematicamente o design desse copo a partir da relação mostrada na 
equação y x = x + x� � � �4 10 0,5sin 2/
Deseja‑se determinar (a) o volume de líquido que cabe dentro do copo e (b) a quantidade de vidro 
necessária para produzir cada copo, considerando uma espessura de 3mm.
Resolução:
(a) Para calcular o volume do sólido gerado pela revolução da função mostrada anteriormente, 
deve‑se calcular a área de cada disco de altura infinitesimal (dV) e raio y(x) e somar desde o início até o 
final do intervalo.
Temos que a área de cada disco (dV) será igual a: dV = y x dx�. 2� �� �
Substituindo y(x) e integrando a equação dV, temos
306
Unidade II
V = x + x dx� 4 10 0,5sin 2
2
0
10
/ � �� ��
Essa integral pode ser calculada, por exemplo, utilizando o método de Simpson. Vamos tomar, por 
exemplo, n = 9. Dessa forma, teremos (2n+1) = 19 pontos. É possível calcular o valor de h = 0,55555 
(refazer as contas para verificar).
A tabela seguinte mostra os cálculos para encontrar os termos auxiliares.
Tabela 57 – Calculando os termos auxiliares da integral
x f(x) Fator Auxiliar
x1 0 0 1 0
x2 0,555556 1,934617 4 7,738468
x3 1,111111 2,996165 2 5,99233
x4 1,666667 2,36455 4 9,458198
x5 2,222222 1,969699 2 3,939397
x6 2,777778 3,152877 4 12,61151
x7 3,333333 6,232396 2 12,46479
x8 3,888889 8,957972 4 35,83189
x9 4,444444 8,537906 2 17,07581
x10 5 6,535266 4 26,14106
x11 5,555556 6,17402 2 12,34804
x12 6,111111 8,751221 4 35,00488
x13 6,666667 13,0535 2 26,10699
x14 7,222222 15,02281 4 60,09123
x15 7,777778 12,98577 2 25,97154
x16 8,333333 10,51225 4 42,049
x17 8,888889 11,10362 2 22,20724
x18 9,444444 15,26436 4 61,05743
x19 10 19,86015 1 19,86015
O valor da integral será dado pela equação
I
h
f x + f x + f x + + f x + f x =n n� � � � � � � � � � ��� ���3 4 2 4 80,731481 2 3 1
Esse valor deve ser multiplicado por n para encontrar o volume, que nesse caso vale: V = 253,63 cm3. 
Esse é o volume interno do copo.
307
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
(b) Para calcular o volume de vidro necessário para a fabricação desse copo, deve‑se encontrar a 
função que represente a parte de fora do copo. O volume do vidro será dado pela subtração entre os 
volumes externos e internos do copo.
Como foi dado no problema que a espessura do copo é constante e igual a 1 mm (0,1 cm), basta 
adicionar esse valor à equação do copo.
5
4
3
2
1
0
Altura (cm)
0 2 4 6 8 10 12
Ra
io
 (c
m
)
Interno
Externo
A modelagem matemática do design interno desse copo foi dado pela equação y x = x + x� � � �4 10 0,5sin 2/ . 
Logo, a modelagem matemática para o design externo será y x = x + x +� � � �4 10 0,5sin 2 0,1/ . 
Sabemos que a área de cada disco (dV) será igual a: dV = y x dx�. 2� �� � . Substituindo y(x) e integrando 
a equação dV, temos V = x + x + dx� 4 10 0,5sin 2 0,1
2
0
10
/ � �� ��
A tabela a seguir mostra os cálculos necessários para calcular essa integral.
Tabela 58 – Calculando os termos auxiliares da integral
x f(x) Fator Auxiliar
x1 0 0,01 1 0,01
x2 0,555556 2,222798 4 8,891193
x3 1,111111 3,352354 2 6,704707
x4 1,666667 2,682091 4 10,72837
x5 2,222222 2,260391 2 4,520781
x6 2,777778 3,518004 4 14,07202
x7 3,333333 6,741691 2 13,48338
x8 3,888889 9,56657 4 38,26628
x9 4,444444 9,1323 2 18,2646
x10 5 7,056549 4 28,2262
x11 5,555556 6,680971 2 13,36194
x12 6,111111 9,35287 4 37,41148
x13 6,666667 13,78609 2 27,57217
308
Unidade II
x14 7,222222 15,80799 4 63,23197
x15 7,777778 13,71649 2 27,43297
x16 8,333333 11,1707 4 44,68281
x17 8,888889 11,78006 2 23,56013
x18 9,444444 16,05575 4 64,223
x19 10 20,76144 1 20,76144
O valor da integral dado pela equação é:
 
I
h
f x + f x + f x + + f x + f x =n n� � � � � � � � � � ��� ���3 4 2 4 86,186191 2 3 1
Esse valor deve ser multiplicado por n para encontrarmos o volume externo, que nesse caso é: 
V = cmexterno 270,7619
3
Esse valor V = 270,7619 cm3 é o volume externodo copo. Subtraindo‑se o volume interno 
V = 253,6254 obtido no item (a) desta aplicação, teremos o volume de vidro necessário para essa 
fabricação, que será:
V = = cmvidro 270,7619 253,6254 17,1365
3−
8.6 Integração numérica: implementação computacional
Esta seção mostrará um exemplo de aplicação da fórmula de Newton‑Cotes para integração 
numérica, tanto usando o método do trapézio quanto o método de Simpson. Para isso, vamos calcular 
como exemplo a integral mostrada na equação a seguir.
2x3
1
I e dx
−
= ∫
Essa integral possui solução analítica conhecida, que é mostrada pela equação na sequência.
( ) ( )
2x3
1
1 1
e dx erf 3 erf 1 0,13938
2 2
−
= π − π = …∫
Para a solução, utilizando o método dos trapézios, devemos dividir o domínio em n pontos e utilizar 
a equação a seguir para aproximar o valor da integral.
( ) ( ) ( ) ( ) ( )1 2 3 n 1 n
h
I f x 2f x 2f x 2f x f x
2 −
 ≅ + + +…+ + 
309
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Vamos utilizar n = 2 para nosso problema, ou seja, vamos utilizar 3 pontos. Dessa forma, teremos 
que h vale:
n 1x x 3 1h 1
n 2
− −
= = =
A tabela a seguir mostra o cálculo dos valores de f(xi).
Tabela 59 – Cálculo dos valores de xi e f(xi) com n= 2
xi f(xi) = e
-x
2
i
1 0,367879
2 0,018316
3 0,000123
Com esses dados, podemos substituir os valores na equação 
( ) ( ) ( ) ( ) ( )1 2 3 n 1 n
h
I f x 2f x 2f x 2f x f x
2 −
 ≅ + + +…+ +  , obtendo a equação a seguir:
[ ]1I 0,367879 2 0,018316 0,000123 0,202317
2
≅ + × + =
Esse valor é distante do valor exato da integral. Caso se deseje melhorar essa aproximação, será 
necessário aumentar o número de pontos. Tomando, por exemplo, n = 10 (11 pontos), temos a próxima 
tabela. O valor de h será dado pela equação a seguir.
n 1x x 3 1h 0,2
n 10
− −
= = =
Tabela 60 – Cálculo dos valores de xi e f(xi) com n= 10
xi f(xi) = e
-x
2
i
1 0,367879
 1,2 0,236928
1,4 0,140858
1,6 0,077305
1,8 0,039164
2 0,018316
2,2 0,007907
310
Unidade II
2,4 0,003151
2,6 0,001159
2,8 0,000394
3 0,000123
Com esses valores, temos que o valor da integral será aproximado pela equação na sequência.
[ ]0,2I 0,36788 2 0,23693 2 0,000394 0,000123 0,14182
2
≅ + × +…+ × + =
Observa‑se que esse valor é mais próximo do valor anterior, com n = 2. Vamos agora resolver 
o mesmo problema utilizando o método de Simpson. Para esse método, temos que a distância 
entre os pontos será dada pela equação a seguir. O número de pontos para esse método será 
igual a (2n + 1).
2n 1 1x xh
2n
+ −=
O valor da integral será aproximado pela equação a seguir. Note que todos os coeficientes 
intermediários (f(xn)) com valores pares de n são multiplicados por 4. Os valores intermediários 
com valores ímpares de n são multiplicados por 2. Os valores extremos, por sua vez, são 
multiplicados por 1.
( ) ( ) ( ) ( ) ( )1 2 3 n 1 n
h
I f x 4f x 2f x 4f x f x
3 −
 ≅ + + +…+ + 
Vamos tomar o exemplo de n = 2, ou seja, o número de pontos sendo 5. Para esse caso, temos h 
calculado pela equação a seguir e os demais valores na tabela na sequência.
2n 1 1x x 3 1h 0,5
2n 2.2
+ − −= = =
Tabela 61 – Cálculo dos valores de xi e f(xi) com n = 2 no método de Simpson
xi f(xi) = e
-x
2
i
1 0,367879
1,5 0,105399
2 0,018316
2,5 0,00193
3 0,000123
311
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
O valor da integral será, portanto, aproximado pela equação a seguir:
[ ]0,5I 0,3679 4 0,1054 2 0,018 4 0,002 0,0001 0,13899
3
≅ + × + × + × + =
Observa‑se que, para n = 2 (5 pontos), o método de Simpson se aproxima muito melhor do valor 
exato do que o método dos trapézios com n = 10 (11 pontos). A tabela a seguir mostra os valores 
encontrados para diferentes números de pontos.
Tabela 62 – Comparação das aproximações conseguidas 
com o método dos trapézios e de Simpson
Método dos trapézios Simpson
n Pontos Integral Erro n Pontos Integral Erro
1 2 0,368002851 0,22862
2 3 0,202317064 0,062934
3 4 0,166998793 0,027616 1 3 0,147088 0,007705
4 5 0,154823372 0,01544
5 6 0,14923601 0,009853 2 5 0,138992 ‑0,00039
6 7 0,146214377 0,006831
7 8 0,144397098 0,005014 3 7 0,139286 ‑9,7E‑05
8 9 0,143219511 0,003836
9 10 0,142413027 0,00303
10 11 0,141836589 0,002453
20 21 0,139995938 0,000613
30 31 0,139655485 0,000272
50 51 0,139481223 9,8E‑05
100 101 0,139407716 2,45E‑05
Nessa tabela é possível observar que o método de Simpson com apenas 7 pontos consegue resultados 
da mesma ordem de erro do método dos trapézios para 51 pontos.
Exemplo de implementação computacional: integração numérica no Python
Essa seção mostra um algoritmo que calcula integrais numéricas no Python. Esse programa recebe 
como entrada os valores limites superior e inferior da integração, além do número de blocos em que se 
deseja dividir o domínio n. Com esses dados, o programa calcula o tamanho de h e cria um vetor x com 
os valores de xi. A variável y armazena os valores de f(xi).
O código a seguir mostra a implementação do método dos trapézios. Note que a variável y2 é criada 
recebendo os mesmos valores da variável y. Na sequência, o primeiro e o último valor de y2 são zerados, 
e y somado com y2. Essa estratégia permite a soma dos coeficientes, de acordo com a equação:
312
Unidade II
( ) ( ) ( ) ( ) ( )1 2 3 n 1 n
h
I f x 2f x 2f x 2f x f x
2 −
 ≅ + + +…+ + 
Para copiar os valores de y para y2, utilizou‑se a função copy. Caso tivéssemos só igualado os vetores, 
o Python entenderia que desejávamos apenas copiar todos os atributos do vetor y, e, dessa forma, ao 
mudar os valores em y2, os valores em y seriam alterados também.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Integração numérica pelo método dos trapézios
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Limites de Integração
limites = np.array([1, 3])
n = 100 # Número de divisões
#######################################################################
## CÁLCULOS ##
#######################################################################
h = (np.max(limites)‑np.min(limites))/n
x = np.linspace(np.min(limites), np.max(limites), n+1)
y = np.exp(‑np.power(x,2))
y2 = y.copy()
y2[0]=0
y2[np.size(y)‑1]=0
Int = h*np.sum(y+y2)/2 # De acordo com a Equação
#######################################################################
## PÓS‑PROCESSAMENTO ##
#######################################################################
# Mostra o resultado na tela
print(‘Valor da integral: ‘, Int)
Para o método de Simpson, uma aproximação parecida pode ser aplicada. Utilizaremos o comando 
y[1::2] = y[1::2]*3 para multiplicar todos os índices pares por 3. Ao somar com a variável y2, de maneira 
semelhante ao caso anterior, teremos a equação ( ) ( ) ( ) ( ) ( )1 2 3 n 1 n
h
I f x 4f x 2f x 4f x f x
3 −
 ≅ + + +…+ +  
satisfeita. O trecho de código a seguir exemplifica essa implementação.
# ‑*‑ coding: utf‑8 ‑*‑
“””
Integração numérica pelo método dos Simpson
@author: Luís Lamas ‑ feb/2018
“””
import numpy as np
#######################################################################
## DADOS DE ENTRADA ##
#######################################################################
# Limites de Integração
limites = np.array([1, 3])
n = 100 # Número de divisões
313
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
#######################################################################
## CÁLCULOS ##
#######################################################################
h = (np.max(limites)‑np.min(limites))/(2*n)
x = np.linspace(np.min(limites), np.max(limites), 2*n+1)
y = np.exp(‑np.power(x,2))
# Cria uma cópia dos valores de y, para zerar o primeiro e o último 
elemento
y2 = y.copy()
y2[0]=0
y2[np.size(y)‑1]=0
# Multiplica os valores de índices pares por 3
y[1::2] = y[1::2]*3
Int = h*np.sum(y+y2)/3
#######################################################################
## PÓS‑PROCESSAMENTO #########################################################################
# Mostra o resultado na tela
print(‘Valor da integral: ‘, Int)
Sugere‑se, nesse ponto, implementar esse código em seu próprio computador e comparar os 
resultados com os exemplos obtidos anteriormente.
 Saiba mais
Recomendamos que você consulte novos exemplos e aplicações, em 
especial as seguintes obras: 
BRASIL, M. L. R. F.; BALTHAZAR, J. M.; GÓIS, W. Métodos numéricos e 
computacionais na prática de engenharia e ciências. São Paulo: Blucher, 2015.
VARGAS, J. V. C.; ARAKI, L. K. Cálculo numérico aplicado. Barueri: 
Manole, 2017.
 Resumo
No ajuste de curvas, dependendo das incertezas dos dados, utiliza‑se 
regressão ou interpolação. Para medidas precisas, a interpolação é usada 
para determinar uma curva que passe diretamente por cada um dos pontos. 
Para medidas imprecisas, a regressão é usada para determinar uma curva 
interpoladora que ajuste a tendência geral dos dados sem necessariamente 
passar por nenhum ponto individual. Todos os métodos de regressão são 
desenvolvidos para ajustar funções que minimizem a soma dos quadrados 
314
Unidade II
dos resíduos entre os dados e a função. Esses métodos são chamados de 
regressão por mínimos quadrados.
A regressão por mínimos quadrados linear é usada nos casos em que 
uma variável dependente e uma independente estão relacionadas entre 
si de uma forma linear. Para casos em que uma variável dependente e 
uma independente exibem uma relação curvilínea no gráfico de dispersão, 
utiliza‑se regressão polinomial. Há casos em que o diagrama de dispersão 
sugere uma curva interpoladora, daí uma função exponencial poder ser 
obtida através da regressão não linear.
O polinômio interpolador de Lagrange para um grau qualquer n é 
apresentado como:
n 0 0 1 1 2 2 n nP (x) f(x ) L (x) f(x ) L (x) f(x ) L (x) ... f(x ) L (x)= ⋅ + ⋅ + ⋅ + + ⋅
Em que: Lk(x) apresentam frações com k variando de 0 a n. Tais frações 
são razões entre os valores possíveis de x.
A determinação do erro é dada por:
a
valor real aproximação atual
100%
valor real
−
ε =
A interpolação linear pode ser obtida pela fórmula:
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
= + −
−
A interpolação quadrática pode ser obtida pela fórmula:
( ) ( ) ( )( )
1 0
0 2 0 2 2 0 2 1
1 0
f(x ) f(x )
F(x) f(x ) x x b x x x x
x x
−
= + − + − −
−
Em que:
( ) ( )
1 02 1
2 1 1 0
2
2 0
f(x ) f(x )f(x ) f(x )
x x x x
b
(x x )
−−
−
− −
=
−
315
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Interpolação de diferenças finitas: Newton-Gregory
As diferenças para as ordens definidas:
Ordem 0
0
i iy y∆ =
Ordem 1
1 0 0
i i 1 i i 1 iy y y y y+ +∆ = ∆ − ∆ = −
Ordem 2
2 1 1
i i 1 iy y y+∆ = ∆ − ∆
Ordem n
n n 1 n 1
i i 1 iy y y
− −
+∆ = ∆ − ∆ ou 
n n 1 n 1f(x) f(x h) f(x)− −∆ = ∆ + − ∆
A formula da interpolação polinomial de Newton‑Gregory é:
1 2 n
0 0 0 0
0 0 1 0 1 n 11 2 n
f(x ) f(x ) f(x )
P(x) f(x) (x x ) (x x )(x x ) ... (x x )(x x )...(x x )
1!h 2!h n!h
−
∆ ∆ ∆
= ∆ + − ⋅ + − − ⋅ + + − − − ⋅
1 2 n
0 0 0 0
0 0 1 0 1 n 11 2 n
f(x ) f(x ) f(x )
P(x) f(x) (x x ) (x x )(x x ) ... (x x )(x x )...(x x )
1!h 2!h n!h
−
∆ ∆ ∆
= ∆ + − ⋅ + − − ⋅ + + − − − ⋅
Interpolação de Splines
Uma Spline é um conjunto de polinômios S(x) que atendem ao critério:
0 0 1
1 1 2
n 1 n 1 n
S (x), para x x x
S (x), para x x x
S(x)
M
S (x), para x x x− −
≤ <
 ≤ <= 

 ≤ <
Uma Spline de primeira ordem, ou Spline de primeira linear: dados os 
nós 1 2 nx ,x ,L,x , para cada subintervalo [xi‑1, xi], i = 1, 2, ..., n:
316
Unidade II
i i 1
i i 1 i i 1 i
i i 1 i i 1
x x x x
S f(x ) f(x ) , x [x ,x ]
x x x x
−
− −
− −
− −
= + ∀ ∈
− −
Regressão linear método dos mínimos quadrados
Coeficiente de correlação de Pearson
( ) ( )
( ) ( )
i i i i
2 22 2
i i i i
n x y x y
r
n x x . n y y
−
=
   − −   
∑ ∑ ∑
∑ ∑ ∑ ∑
Regressão linear
O método dos mínimos quadrados consiste em obter os coeficientes a0 
e a1 da função f(x) = a1x + a0.
Temos que: 
n n
i 1 i
i 1 i 1
0
y a x
a
n n
= == −
∑ ∑
 e 
n n n
i i i i
i 1 i 1 i 1
1 2n n
2
i i
i 1 i 1
n x y x y
a
n x x
= = =
= =
−
=
 
−  
 
∑ ∑ ∑
∑ ∑
No processo de buscar uma aproximação, temos que considerar o erro e:
1 0f(x) a x a e= + +
A soma dos erros residuais é dada por:
n n
i i 1 i 0
i 1 i 1
e y a x a
= =
= − −∑ ∑
Regressão polinomial
É a extensão do procedimento dos mínimos quadrados
2
0 1 2y a a x a x e= + + +
2
0 1 2e y a a x a x= − − −
317
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
A soma dos quadrados dos resíduos é:
( )
2n
2
r i 0 1 i 2 i
i 1
S y a a x a x
=
= − − −∑
Um sistema de três equações que possibilita o cálculo de a0, a1 e a2 a 
partir dos dados.
n n n
2
i 0 1 i 2 i
i 1 i 1 i 1
n n n n
2 3
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
n n n n
2 2 3 4
i i 0 i 1 i 2 i
i 1 i 1 i 1 i 1
y na a x a x
x y a x a x a x
x y a x a x a x
= = =
= = = =
= = = =

= + +


 = + +


 = + +

∑ ∑ ∑
∑ ∑ ∑ ∑
∑ ∑ ∑ ∑
Integração numérica
A base dos métodos de integração numérica apoia‑se na relação entre 
integrais definidas e áreas.
Para todos os casos necessitamos que f seja uma função integrável 
no intervalo [a, b] e 0 1 nx x ,L,x, uma partição uniforme desse intervalo. 
Então, a integral de f em cada caso pode ser aproximada por:
Regra do ponto médio
( )
nb
M ia
i 1
f x dx I x f(x )
=
≈ = ∆ ∑∫
Em que 
b a
x
n
−
∆ = e [ ]i i 1i i 1 i
x x
x ponto médio de x ,x
2
−
−
−
= = .
Regra do trapézio
( ) ( ) ( ) ( )
n 1b
T ia
i 1
f a f b
f x dx I f x x
2
−
=
 +
≈ = + ∆ 
 
∑∫
Em que 
b a
x
n
−
∆ = e xi = a + i∆x.
318
Unidade II
Regra de Simpson
( ) ( ) ( ) ( ) ( )
n n
1
2 2b
S 2i 1 2ia
i 1 i 1
x
f x dx I f a f b 4 f x 2 f x
3
−
−
= =
 
∆  ≈ = + + + 
  
∑ ∑∫
Em que 
b a
x
n
−
∆ = e xi = a + i∆x.
 Exercícios
Questão 1. O gráfico da figura a seguir representa a curva da função y = 5 ln (x).
Figura
A estimativa do valor da função para x = 3 e erro encontrado quando se usa uma interpolação linear, 
com o intervalo [1,4] são, respectivamente:
A) 6,9314718 e 15,88%.
B) 4,6209812 e 31,16%.
C) 5,4930614 e 31,16%.
D) 4,6209812 e 15,88%.
E) 5,4930614 e 15,88%.
Resposta correta: alternativa D.
319
CÁLCULO NUMÉRICO E SIMULAÇÃO PARA ENGENHARIA DE PRODUÇÃO
Análise das alternativas
A) Alternativa incorreta.
Justificativa: o valor apresentado para a função é o valor para x = 4.
B) Alternativa incorreta.
Justificativa: embora o valor da estimativa esteja verdadeiro, o erro é igual ao dobro do correto.
C) Alternativa incorreta.
Justificativa: o valor apresentado para a função é o valor verdadeiro. Não é o da estimativa. Além 
disso, o erro é igual ao dobro do correto.
D) Alternativa correta.
Justificativa: f(x) 5 ln(x)= ⋅ .
Então: 
0x 1
f(1) 5 ln1 0
=

= ⋅ =
 e 1
x 4
f(4) 5 ln4 6,9314718
=
 = ⋅ =
( )1 00 0
1 0
f(x ) f(x )
F(x) f(x ) . x x
x x
−
= + −
−
( )6,9314718 0F(3) 0 3 1
4 1
−
= + −
−
( )6,9314718F(3) 2
3
=
F(3) 4,6209812=
a
valor real aproximação atuall
100%
valor real
−
ε =
Considerando que o valor correto de 5 ln3 5,4930614⋅ = :
a
5,4930614 4,6209812
100% 15,88%
5,4930614
−
ε = =
320
Unidade II
E) Alternativa incorreta.
Justificativa: o valor apresentado para a função é o valor verdadeiro. Não é o da estimativa.
Questão 2. O gráfico da figura a seguir representa a curva da função y = 5 ln (x):
Figura
A estimativa do valor da função para x = 3 e erro encontrado quando se usa uma interpolação 
quadrática, com o intervalo [1,4] são, respectivamente:
A) 6,9314718 e 15,88%.
B) 4,6209812 e 6,43%.
C) 5,1397114 e 31,16%.
D) 4,6209812 e 15,88%.
E) 5,1397114 e 6,43%.
Resolução desta questão na plataforma.
321
FIGURAS E ILUSTRAÇÕES
Figura 2
HOFFMANN, L. D.; BRADLEY, G. L. Cálculo: um curso moderno e suas aplicações. 10. ed. Rio de Janeiro: 
LTC, 2010. p. 15.
Figura 3
CHAPRA, C. S. Métodos numéricos aplicadoscom MATLAB para engenheiros e cientistas. 3. ed. 
Porto Alegre: AMGH, 2013. p. 106.
Figura 16
CHAPRA, C. S. Métodos numéricos para engenharia. 7. ed. São Paulo: AMGH, 2016. p. 1.
Figura 22
CHAPRA, C. S. Métodos numéricos para engenharia. 7. ed. São Paulo: AMGH, 2016. p. 133.
Figura 30
CHAPRA, C. S. Métodos numéricos para engenharia. 7. ed. São Paulo: AMGH, 2016. p. 138.
Figura 35
CHAPRA, C. S. Métodos numéricos para engenharia. 7. ed. São Paulo: AMGH, 2016. p. 143.
Figura 54
800PX‑SPLINE_SMOOTHING.PNG. Disponível em: <https://upload.wikimedia.org/wikipedia/commons/
thumb/f/f7/Spline_Smoothing.png/800px‑Spline_Smoothing.png>. Acesso em: 30 jul. 2018.
Figura 65
Grupo UNIP‑Objetivo.
Figura 95
Grupo UNIP‑Objetivo.
Figura 96
Grupo UNIP‑Objetivo.
322
REFERÊNCIAS
Textuais
ÁVILA, G. Cálculo das funções de uma variável. Rio de Janeiro: LTC, 2003.
BICUDO, M. A. V.; GARNICA, A. V. M. Filosofia da educação matemática. Belo Horizonte: Autêntica, 2001.
BOYER, C. B. História da matemática. São Paulo: Edgard Blucher, 2001.
BRASIL, M. L. R. F.; BALTHAZAR, J. M.; GÓIS, W. Métodos numéricos e computacionais na prática de 
engenharia e ciências. São Paulo: Blucher, 2015.
CHAPRA, C. S. Métodos numéricos aplicados com MATLAB para engenheiros e cientistas. 3. ed. 
Porto Alegre: AMGH, 2013.
___. Métodos numéricos para engenharia. 7. ed. São Paulo: AMGH, 2016.
DOMINGUES, H. H. Fundamentos da matemática. São Paulo: Atual, 1998.
EVES, H. Introdução à história da matemática. Tradução Hygino H. Domingues. Campinas: Editora da 
Unicamp, 2004.
FOUCAULT, M. A verdade e as formas jurídicas. Rio de Janeiro: Nau, 1999.
GILAT, A. Métodos numéricos para engenheiros e cientistas: uma introdução com aplicações usando o 
MATLAB. Porto Alegre: Bookman, 2008.
HOFFMANN, L. D.; BRADLEY, G. L. Cálculo: um curso moderno e suas aplicações. 10. ed. São Paulo: LTC, 2010.
MACHADO, D. I.; NARDI, R. Construção de conceitos de física moderna e sobre a natureza da ciência com o 
suporte da hipermídia. Revista Brasileira de Ensino de Física, v. 28, n. 4, São Paulo, 2006. Disponível em: <http://
www.scielo.br/scielo.php?script=sci_arttext&pid=S1806‑11172006000400010>. Acesso em: 1º ago. 2018.
MAGAGNATO, P. C. Fundamentos teóricos da atividade de estudo como modelo didático para o ensino 
das disciplinas científicas. 2011. Dissertação (Mestrado em Educação para a Ciência). Universidade 
Estadual Paulista, Bauru, 2011.
MARCONDES, G. A. B. Matemática com Pyton: um guia prático. São Paulo: Novatec, 2018.
MEGLIORINI, E. Administração financeira. São Paulo: Prentice Hall, 2012.
RUGGIERO, M. A. G.; ROCHA LOPES, V. L. Cálculo numérico: aspectos teóricos e computacionais. 2. ed. 
São Paulo: Pearson Makron Books, 2014.
323
SISTEMAS lineares. [s.d.]. Disponível em: <http://www.mat.ufmg.br/~rodney/notas_de_aula/sistemas_
lineares.pdf>. Acesso em: 23 jul. 2018.
SMAILES, J.; MCGRANE, A. Estatística aplicada à administração com Excel. São Paulo: Atlas, 2002.
UNIVERSIDADE FEDERAL DO PARÁ. Pesagem da amostra. Belém, [s.d.]. Disponível em <http://www.
ufpa.br/quimicanalitica/sbalancas.htm>. Acesso em: 19 jul. 2018.
VARGAS, J. V. C.; ARAKI, L. K. Cálculo numérico aplicado. Barueri: Manole, 2017.
WINPLOT. Baixaki, 2010. Disponível em: <http://www.baixaki.com.br/download/winplot.htm>. 
Acesso em: 19 jul. 2018.
324
325
326
327
328
Informações:
www.sepi.unip.br ou 0800 010 9000