Cálculo II e Cálculo Numérico - Parte IV

Prévia do material em texto

Unidade 04
Aula 01
Ajuste de Curvas (Mínimos Quadrados)
Olá alunos! Sejam bem-vindos a mais uma aula. Continuem com o ritmo! Boa aula.
Considere que se tenha uma tabela de pontos.
ue podem ter sido observados em algum tipo de experimento. Podemos supor que esses pares de pontos 
 quando colocados no plano possuem o comportamento semelhante ao de alguma
função conhecida dada na forma
em que as funções são contínuas e conhecidadas dentro do intervalo . Desse modo, nesta aula
estudaremos um problema do ajuste de curvas que consiste em determinar os valores das constantes 
, de maneira que a função se aproxime o máximo possível da função que modela
os dados observados.
A primeira observação que fazemos é que este modelo de ajuste de curva, dado pela função é linear.
Pois os coe�cientes que devem ser determinados aparecem linearmente na função ,
embora as funções possam ser funções não lineares da variável , como, por
exemplo, Além disso, precisamos de�nir a nossa
noção de proximidade desejada entre as funções e . No caso do método dos mínimos quadrados,
ela consiste em escolher os coe�cientes de tal forma que a soma dos quadrados dos
desvios
x1 x2 x3 ⋯ xk
f (x1) f (x2) f (x3) ⋯ f (xk)
{(xi, f (xi) } ki=1, x − y,
φ (x) = a1g1 (x) + a2g2 (x) + ⋯ + angn (x),
gi (x) [x1,xk]
a1, a2, ⋯ , an φ (x) f (x)
φ (x),
a1,  a2,   … ,  an φ (x)
g1 (x),  g2 (x),   … ,  gn (x) x
g1 (x) =  x4,  g2 (x) = ln (x),   ⋯ , gn (x) = xenx.
φ (x) f (x)
a1,  a2,   … ,  an
h (a1, a2, ⋯ , an) =
k
∑
i=1
[φ (xk) − f (xk)]
2
seja mínima. Para facilitar o entendimento do método, começaremos com um exemplo.
Exemplo 1:
Para exempli�car o método dos mínimos quadrados, vamos, inicialmente, considerar o seguinte conjunto de
pontos que foi gerado a partir da reta adicionando um pouco de ruído nas medidas.
É possível ver como esses pontos se distribuem no plano através do grá�co a seguir. Como o conjunto de
pontos possui ruído, veri�camos que o comportamento deles foge um pouco da reta No
entanto, vamos utilizar o método dos mínimos quadrados para encontrar uma reta que
se aproxima o máximo possível dos pontos dados, segundo o erro quadrático de�nido anteriormente.
Desta forma, temos que encontrar os valores de e que minimizam a função
Para tal, vamos calcular os seus pontos críticos, isto é, encontrar e tais que
Neste caso, teremos
=
k
∑
i=1
[a1g1 (xi) + ⋯ + angn (xi) − f (xi)]
2
y = 3x + 1,
xi 1 2 3 4 5 6 7
f (xi) 3.27 7.58 10.88 15.05 15.76 16.43 21.00
y = 3x + 1.
φ (x) = a1 + a2x
a1 a2
h (a1, a2) =
7
∑
i=1
(a1 + a2xi − f (xi))
2
= (a1 + a2 − 3.27)
2 + (a1 + 2a2 − 7.58)
2 + (a1 + 3a2 − 10.88)
2 + (a1 + 4a2 − 15.05)
2 + (a1 + 5a
+(a1 + 6a2 − 16.43)
2 + (a1 + 7a2 − 21.00)
2.
a1 a2
∂h
∂a1
= 0    e    
∂h
∂a2
= 0.
∂h
∂a1
= 0 ⇒
(a1 + a2 − 3.27) + (a1 + 2a2 − 7.58) + (a1 + 3a2 − 10.88) + (a1 + 4a2 − 15.05) + (a1 + 5a2 − 15.
( ) ( )
e
Portanto, para encontrarmos os valores ótimos para minimizar a distância entre a reta desejada e o
conjunto de pontos dado, precisamos resolver o sistema linear
que pode ser facilmente resolvido e nos dá
Logo, pelo método dos mínimos quadrados, a reta que melhor se ajusta ao conjunto de pontos fornecido é a
reta
que podemos ver a seguir no grá�co.
Agora, se voltamos ao problema mais geral, em que a função que modela o comportamento dos pontos 
 é dada por
+ (a1 + 6a2 − 16.43) + (a1 + 7a2 − 21.00) = 0 ⇒
7a1 + 28a2 = 89.97 ⇒
a1 + 4a2 = 12.85
∂h
∂a2
= 0 ⇒
(a1 + a2 − 3.27) + 2 (a1 + 2a2 − 7.58) + 3 (a1 + 3a2 − 10.88) + 4 (a1 + 4a2 − 15.05) + 5 (a1 + 5a2
+6 (a1 + 6a2 − 16.43) + 7 (a1 + 7a2 − 21.00) = 0 ⇒
28a1 + 140a2 = 435.72 ⇒
a1 + 5a2 = 15.56.
φ (x)
{
a1 + 4a2 = 12.85
a1 + 5a2 = 15.56
a1 = 2.01     e     a2 = 2.71.
φ (x) = 2.71x + 2.01,
{(xi, f (xi) } ki=1
φ (x) = a1g1 (x) + a2g2 (x) + ⋯ + angn (x).
Então, para encontrarmos os coe�cientes que minimizam a distância entre e nos
pontos dados, podemos prosseguir como no exemplo que �zemos anteriormente e buscar pelos pontos
críticos de , para cada calculando
Isso nos levará ao seguinte sistema linear 
Que também pode ser escrito na forma
As equações que formam esse sistema linear são chamadas de equações normais, e a matriz dos coe�cientes
para o caso do método dos mínimos quadrados é uma matriz simétrica. Além disso, é possível mostrar que,
dependendo de como as funções são escolhidas, é sempre possível garantir a existência de uma
solução para o sistema linear e que a solução única do sistema corresponde a um ponto de mínimo da função 
.
Exemplo 2:
Vamos agora seguir a dedução das equações normais realizada anteriormente e utilizar o método dos
mínimos quadrados para ajustar a parábola ao conjunto de pontos tabelado a
seguir. Esses pontos foram gerados a partir da parábola com erro Gaussiano adicionado nas
medidas.
a1, ⋯ ,    an φ (x) f (x),
h j = 1, ⋯ ,n,
∂h
∂aj
= 0 ⇒
k
∑
i=1
[a1g1 (xi) + ⋯ + angn (xi) − f (xi)] ⋅ gj (xi) = 0.
n × n
⎧⎪⎨⎪⎩∑ki=1 [a1g1 (xi) + ⋯ + angn (xi) − f (xi)] ⋅ g1 (xi) = 0∑ki=1 [a1g1 (xi) + ⋯ + angn (xi) − f (xi)] ⋅ g2 (xi) = 0⋮∑ki=1 [a1g1 (xi) + ⋯ + angn (xi) − f (xi)] ⋅ gn (xi) = 0.⎧⎪⎨⎪⎩[ k∑i=1 g1(x1)g1(x1)] ⋅ a1+. . . +[ k∑i=1 g1(x1)gn(x1)]an = k∑i=1 g1(x1)f(x1)[ k∑i=1 g2(x1)g1(x1)] ⋅ a1+. . . +[ k∑i=1 g2(x1)gn(x1)]an = k∑i=1 g2(x1)f(x1)[ k∑i=1 gn(x1)g1(x1)] ⋅ a1+. . . +[ k∑i=1 gn(x1)gn(x1)]an = k∑i=1 gn(x1)f(x1)gi (x)h (a1, a2, ⋯ , an) φ (x) = a1 + a2x + a3x2y = 310 x2 − 12xi −3 −2 −1 0 1 2 3f (xi) 2.47 0.48 −0.18 −0.80 −0.47 0.69 2.58
A seguir, podemos ver na �gura como os pontos se distribuem no plano, e, devido ao erro adicionado, como
esses pontos fogem da parábola 
Para ajustar a parábola aos pontos dados, temos que resolver o sistema linear
em que , e . Assim, os coe�cientes do sistema linear são dados por
y = 310 x
2 − 12 .
⎧⎪⎨⎪⎩[ 7∑i=1 g1 (xi)g1 (xi)] ⋅ a1 + [ 7∑i=1 g1 (xi)g2 (xi)] ⋅ a2 + [ 7∑i=1 g1 (xi)g3 (xi)] ⋅ a3 = 7∑i=1 g1(x1)f(x1)[ 7∑i=1 g2 (xi)g1 (xi)] ⋅ a1 + [ 7∑i=1 g2 (xi)g2 (xi)] ⋅ a2 + [ 7∑i=1 g2 (xi)g3 (xi)] ⋅ a3 = 7∑i=1 g2(x1)f(x1)[ 7∑i=1 g3 (xi)g1 (xi)] ⋅ a1 + [ 7∑i=1 g3 (xi)g2 (xi)] ⋅ a2 + [ 7∑i=1 g3 (xi)g3 (xi)] ⋅ a3 = 7∑i=1 g3(x1)f(x1)g1 (x) = 1 g2 (x) = x g3 (x) = x27∑i=1 g1 (xi)g1 (xi) = 7∑i=1 1 = 77∑i=1 g1 (xi)g2 (xi) = 7∑i=1 xi = (−3) + (−2) + (−1) + 0 + 1 + 2 + 3 = 07∑i=1 g1 (xi)g3 (xi) = 7∑i=1 x2i = (−3)2 + (−2)2 + (−1)2 + 02 + 12 + 22 + 32 = 287∑i=1 g2 (xi)g2 (xi) = 7∑i=1 x2i = (−3)2 + (−2)2 + (−1)2 + 02 + 12 + 22 + 32 = 287∑i=1 g2 (xi)g3 (xi) = 7∑i=1 x3i = (−3)3 + (−2)3 + (−1)3 + 03 + 13 + 23 + 33 = 07∑i=1 g3 (xi)g3 (xi) = 7∑i=1 x4i = (−3)4 + (−2)4 + (−1)4 + 04 + 14 + 24 + 34 = 196.
Enquanto os termos da direita são dados por
Finalmente, podemos reescrever o sistema linear com os valores calculados e temos
Esse sistema pode ser facilmente resolvido e temos que os coe�cientes do polinômio quadrático são dados
por
Logo, e podemos ver no grá�co a seguir a parábola que melhor se
aproxima dos pontos dados.
7
∑
i=1
g1 (xi)f (xi) =
7
∑
i=1
f (xi) = 2.47 + 0.48 + (−0.18) + (−0.80)
+ (−0.47) + 0.69 + 2.58 = 4.77
7
∑
i=1
g2 (xi)f (xi) =
7
∑
i=1
xif (xi) = (−3) × 2.47 + (−2) × 0.48 + (−1) × (−0.18) + (0) × (−0.80)
+1 × (−0.47) + 2 × 0.69 + 3 × 2.58 = 0.46
7
∑
i=1
g3 (xi)f (xi) =
7
∑
i=1
x2i f (xi) = (−3)
2 × 2.47 + (−2)2 × 0.48 + (−1)2 × (−0.18) + (0)2 × (−0.8
+12 × (−0.47) + 22 × 0.69 + 32 × 2.58 = 49.48.
⎧⎪⎨⎪⎩a1 +                   + 4a3 = 0.68              28a2               = 0.46a1 +                   + 7a3 = 1.76.a1 ≈ −0.76a2 ≈ 0.01a3 ≈ 0.36.φ (x) = 0.36x2 + 0.01x − 0.76
Unidade 04
Aula 02
Interpolação Polinomial
Bem-vindos a mais uma aula! Estamos próximos do �m da unidade, continue estudando.
Dada uma tabela de pontos
que foram observados em algum tipo de experimento, uma técnica de ajuste de curvas chamada de
regressão linear. Essa técnica consiste em determinar uma função se aproxime o máximo possível da
função que modela os dados observados, segundo alguma medida de erro prede�nida. Nesta aula,
iremos estudar os conceitosrelacionados à interpolação, que também é uma técnica de ajuste de curvas e
consiste em encontrar uma função que se aproxima da função a qual temos alguns pontos
tabelados. No entanto, a grande diferença da regressão linear para a interpolação é que a função
interpolante deve coincidir com a função nos pontos dados, como podemos ver na �gura a
seguir.
x1 x2 x3 ⋯ xk
f (x1) f (x2) f (x3) ⋯ f (xk)
φ (x)
f (x),
Pn (x) f (x),
Pn (x) f (x)
SAIBA MAIS
Para mais exemplos do método dos mínimos quadrados, clique aqui e acesse. 
Uma vídeoaula com um exemplo do método dos mínimos quadrados: Clique aqui.
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d37214724fd2a4e94a48840_texto03.pdf
https://www.youtube.com/watch?v=hFNajOohP5E
Na �gura, a letra (a) representa o que seria a interpolação, enquanto na letra (b) temos uma representação
da regressão linear.
Nesta aula, iremos nos ater à interpolação polinomial. Isto é, dado um conjunto de pontos 
 obtidos experimentalmente, queremos encontrar um polinômio
de tal forma que o polinômio coincida com os pontos dados satisfazendo . A primeira
observação que fazemos neste é que dados pontos distintos , com então o
grau do polinômio deve ser igual a , isto é, .
Por exemplo, se considerarmos os pontos distintos , e e quisermos encontrar um
polinômio que passa por eles, a nossa melhor opção é a parábola. Pois sabemos que a equação da parábola é
dada na sua forma geral
e possui três constantes a serem determinadas e . Como temos três pontos que devem satisfazer a
parábola, teremos o seguinte sistema linear de três equações e três incógnitas
Em forma matricial, o sistema é escrito na forma
Como os pontos são distintos, é possível provar que a matriz do lado esquerdo tem determinante não
nulo, logo a matriz é inversível e pode-se determinar o valor das constantes , e . Com uma ideia
análoga, podemos mostrar que dados pontos distintos, podemos determinar um único polinômio de
grau que passa por eles. Para um conjunto de pontos, teríamos um sistema linear semelhante ao
sistema dado anteriormente, e a matriz dos coe�cientes que surge é a matriz, conhecida como matriz de
Vandermonde. De modo geral, a Matriz de Vandermonde é dada como
{(x0, f0), (x1, f1), ⋯ , (xn, fn)}
Pn (x) = a0 + a1x + a2x
2 + ⋯ + anx
n,
Pn (xj)  =  fj
n  + 1 xj j =  0,  1,  2,  . . . ,  n
Pn (x) n an ≠ 0
(x0, f0) (x1, f1) (x2, f2)
p (x) = a0 + a1x + a2x
2,
a0, a1 a2
a0 + a1x0 + a2x20 = p (x0) = f0
a0 + a1x1 + a2x21 = p (x1) = f1
a0 + a1x2 + a2x
2
2 = p (x2) = f2.
= .
⎡⎢⎣1 x0 x201 x1 x211 x2 x22⎤⎥⎦⎡⎢⎣a0a1a2⎤⎥⎦ ⎡⎢⎣f0f1f2⎤⎥⎦xi a0 a1 a2n + 1n n + 1
e seu determinante é sempre não nulo e dado por
\[ \det V=\underset{0\le i
Considere agora um conjunto de pontos na forma , o qual queremos
interpolar por uma função polinomial . Como, neste caso, temos pontos, então podemos
encontrar um único polinômio de grau que passa por eles. Nesta aula, iremos nos concentrar em uma
técnica conhecida por interpolação polinomial na forma de Lagrange. Esta talvez seja a maneira mais simples
de se obter um polinômio interpolador. Neste caso, iremos escrever o nosso polinômio no modo da
seguinte combinação linear
O segredo da interpolação de Lagrange está na forma em que as funções são construídas. Elas, na
verdade, são polinômios de grau na variável os quais serão chamados de coe�cientes de interpolação
de Lagrange. Para assegurar que o polinômio interpolador satisfaça para cada ponto o critério 
, dado anteriormente, então os coe�cientes de interpolação de Lagrange deverão ser
construídos de tal forma que eles satisfaçam a seguinte condição
pois
No exemplo a seguir, veremos como construir os coe�cientes de Lagrange para o caso de uma interpolação
de primeiro grau.
Exemplo 1:
Dados dois pontos e , queremos encontrar um polinômio de primeiro grau que passa por
esses pontos e esteja na forma do polinômio interpolador de Lagrange
V =
⎡⎢⎣1 x0 x20 ⋯ xn01 x1 x21 ⋯ xn11 x2 x22 ⋯ xn2⋮ ⋮ ⋮ ⋱ ⋮1 xn x2n ⋯ xnn⎤⎥⎦{(x0, f0), (x1, f1), ⋯ , (xn, fn)}Pn (x) n + 1n Pn (x)Pn (x) =  L0 (x)f0 +  L1 (x)f1 + ⋯   +  Ln (x)fn= n∑i=0 Li (x)fi. Li (x)n x, Pn (x) xjPn (xj) = fj Li (xj) = δi,j = { ,1,   se i = j0,   se i ≠ j
Pn (xj) = L0 (xj)f0 + L1 (xj)f1 + ⋯ + Lj (xj)fj + ⋯ + Ln (xj)fn
= δ0,jf0 + δ1,jf1 + ⋯ + δj,jfj + ⋯ + δn,jfn
= 0 ⋅ f0 + 0 ⋅ f1 + ⋯ + 1 ⋅ fj + ⋯ + 0 ⋅ fn
= fj.
(x0, f0) (x1, f1)
( ) ( ) ( )
Pela de�nição imposta anteriormente, o coe�ciente de Lagrange deve ser do primeiro grau e
satisfazer as seguintes condições e . Neste caso, uma boa opção para esse
coe�ciente seria na forma
É imediato ver que esse coe�ciente satisfaz as condições e . Além disso, como o
coe�ciente também deve satisfazer uma condição semelhante e , então
podemos escrevê-lo de forma similar à descrita anteriormente, sendo
Desse modo, o polinômio interpolador de Lagrange de primeiro grau pode ser escrito como
Vimos como construir o coe�ciente de Lagrange para um polinômio interpolador de primeiro grau. Porém,
não é difícil veri�car que, de forma geral, a seguinte de�nição dos coe�cientes de interpolação de Lagrange
para um conjunto de pontos é dada na forma
e satisfaz a condição imposta anteriormente.
Exemplo 2:
Vamos considerar a seguinte tabela de pontos que pode ser dada por algum experimento. Nosso objetivo é
encontrar um polinômio interpolador de Lagrange que satisfaça os pontos dados.
0 1 2 3
1 2 3 4
1.54 0.58 0.01 0.35
Tendo em vista que temos um conjunto com quatro pontos, então podemos encontrar um polinômio de
terceiro grau que satisfaz esses pontos. Neste caso da interpolação de Lagrange, o polinômio deve ter a
forma
P1 (x) = L0 (x)f0 + L1 (x)f1.
L0 (x)
L0 (x0) = 1 L0 (x1) = 0
L0 (x) =
x − x1
x0 − x1
.
L0 (x0) = 1 L0 (x1) = 0
L1 (x) L1 (x0) = 0 L1 (x1) = 1
L1 (x) =
x − x0
x1 − x0
.
P1 (x) = L0 (x)f0 + L1 (x)f1
=
x − x1
x0 − x1
f0 +
x − x0
x1 − x0
f1.
{(x0, f0), (x1, f1), ⋯ , (xn, fn)}
Li (x)  =
(x − x1) (x − x2) ⋯ (x − xi−1) (x − xi+1) ⋯ (x − xn)
(xi − x1) (xi − x2) ⋯ (xi − xi−1) (xi − xi+1) ⋯ (xi − xn)
Li (xj) = δi,j 
i
xi
fi
P3 (x) = f0L0 (x) + f1L1 (x) + f2L2 (x) + f3L3 (x)
= 1.54L0 (x) + 0.58L1 (x) + 0.01L2 (x) + 0.35L3 (x).
Cada um dos coe�cientes de Lagrange pode ser escrito utilizando a fórmula dada anteriormente, neste caso
os coe�cientes são:
Podemos ver sem grandes di�culdades, por exemplo, que
e que
Finalmente, na sua forma simpli�cada, o polinômio interpolador de Lagrange pode ser escrito como
e, a seguir,temos o grá�co do polinômio calculado realizando a interpolação nos pontos tabelados
dados.
Apesar da sua simplicidade em construir os polinômios interpoladores, uma das desvantagens da
interpolação de Lagrange é que, ao adicionar pontos às medidas e necessitar mudar a interpolação
polinomial de Lagrange de um certo grau para outro, digamos da linear para quadrática, todos cálculos dos
coe�cientes de Lagrange devem ser refeitos. Isso se deve ao fato que os coe�cientes de Lagrange são
altamente dependentes do conjunto de pontos dados. Por exemplo, nenhuma das expressões usadas no caso
da interpolação linear (os coe�cientes lineares de Lagrange e ) é util na construção da
interpolação polinomial quadrática de Lagrange. Em virtude desta observação, muitos esforços foram feitos
na teoria para construir outras formas de interpolações polinomiais que são iterativas, no sentido de que
L0 (x) =
(x − 2) (x − 3) (x − 4)
(1 − 2) (1 − 3) (1 − 4)
= −
1
6
(x − 2) (x − 3) (x − 4)
L1 (x) =
(x − 1) (x − 3) (x − 4)
(2 − 1) (2 − 3) (2 − 4)
=
1
2
(x − 1) (x − 3) (x − 4)
L2 (x) =
(x − 1) (x − 2) (x − 4)
(3 − 1) (3 − 2) (3 − 4)
= −
1
2
(x − 1) (x − 2) (x − 4)
L3 (x) =
(x − 1) (x − 2) (x − 3)
(4 − 1) (4 − 2) (4 − 3)
=
1
6
(x − 1) (x − 2) (x − 3).
L1 (x1) = L1 (2) =
1
2
(2 − 1) (2 − 3) (2 − 4) = 1
L2 (x3) = L2 (4) = −
1
2
(4 − 1) (4 − 2) (4 − 4) = 0.
P3 (x) = 2.37 + 0.591667x − 0.325x2 + 0.0866667x3
P3 (x)
L0 (x) L1 (x)
uma interpolação polinomialde grau superior pode ser obtida simplesmente adicionando termos de maior
grau à interpolação polinomial de grau inferior já existente. Veremos, no �m, que os polinômios de Lagrange
e Newton são os mesmos, no entanto a forma iterativa de se obter o polinômio de Newton é que mostra o
benefício de seu uso.
Continuando estudando a interpolação polinomial, porém iremos estudar uma forma diferente de se determinar o
polinômio interpolador. Estudamos a interpolação polinomial utilizando o método de Lagrange.
Se considerarmos o conjunto de dados formado pelos pontos , então o
polinômio interpolador de Lagrange é dado pela seguinte combinação linear
em que cada um dos coe�cientes de Lagrange é dado na forma
e satisfaz a condição .
Conforme vimos, os coe�cientes de Lagrange são altamente dependentes do conjunto de pontos dado. Por
exemplo, os coe�cientes de Lagrange e , usados no caso da interpolação linear, não servem
para a construção do polinômio interpolador quadrático na forma de Lagrange. Considerando esta situação,
os nossos esforços nesta aula serão para construir um outro método para determinar o polinômio
interpolador de forma iterativa, de modo que uma interpolação polinomial de grau superior possa ser obtida
simplesmente adicionando termos de maior grau à interpolação polinomial de grau inferior já existente. O
método que estudaremos nesta aula é conhecido como método das diferenças divididas de Newton. Na
formulação do método de Newton, não é necessário que os pontos usados sejam igualmente espaçados ou
{(x0, f0), (x1, f1), ⋯ , (xn, fn)}
Pn (x)
Pn (x) =  L0 (x)f0 +  L1 (x)f1 + ⋯   +  Ln (x)fn,
Li (x)    =
(x − x1) (x − x2) ⋯ (x − xi−1) (x − xi+1) ⋯ (x − xn)
(xi − x1) (xi − x2) ⋯ (xi − xi−1) (xi − xi+1) ⋯ (xi − xn)
Li (xj) = δi,j
L0 (x) L1 (x)
SAIBA MAIS
A seguir, é possível assistir a videoaula com exemplos sobre a interpolação de Lagrange: Clique aqui. 
Para mais exemplos, clique aqui e acesse!
https://www.youtube.com/watch?v=khZv01Om6is
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d372d616d9079c25cc3bc96_texto04.pdf
que os valores das abscissas sejam necessariamente em ordem crescente. Assim, considerando o conjunto
de pontos podemos escrever a expressão da forma de Newton que é
dada por:
O desa�o neste momento é como determinar os coe�cientes constantes para que o polinômio
interpolador satisfaça o conjunto de pontos dados, isto é, satisfaça . Assim, considerando a
forma do polinômio dado, se �zermos , podemos encontrar o coe�ciente , pois
Para encontrarmos o coe�ciente substituímos no polinômio para obtermos
Podemos continuar este processo de substituição para determinar os demais. No caso do método das
diferenças divididas de Newton, esses coe�cientes recebem uma notação especial e são usualmente
escritas na forma
Considerando os coe�cientes e calculados anteriormente, podemos, sem maiores di�culdades,
generalizar as funções entre colchetes de�nidas anteriormente, que são chamadas de diferenças divididas e
são de�nidas de acordo com seus índices na forma
{(x0, f0), (x1, f1), ⋯ , (xn, fn)},
Pn (x) = d0 + d1 (x − x0) + d2 (x − x0) (x − x1) + ⋯ + dn (x − x0) (x − x1) ⋯ (x − xn−1).
di
′s
Pn (xj) = fj
x = x0 d0
Pn (x0) = d0 + d1 (x0 − x0) + d2 (x0 − x0) (x0 − x1) + ⋯ + dn (x0 − x0) ⋯ (x0 − xn−1)
⇒ f0 = d0 + d1 ⋅ 0 + d2 ⋅ 0 ⋅ (x0 − x1) + ⋯ + dn ⋅ 0 ⋅ ⋯ (x0 − xn−1)
⇒ f0 = d0.
d1 x = x1
Pn (x1) = d0 + d1 (x1 − x0) + d2 (x1 − x0) (x1 − x1) + ⋯ + dn (x1 − x0) ⋯ (x1 − xn−1)
⇒ Pn (x1) = d0 + d1 (x1 − x0) + d2 (x1 − x0) ⋅ 0 + ⋯ + dn (x1 − x0) ⋅ 0 ⋅ ⋯ (x1 − xn−1)
⇒ f1 = d0 + d1 (x1 − x0)
⇒ d1 =
f1 − d0
x1 − x0
=
f1 − f0
x1 − x0
.
di
d0 = f [x0] = f0
d1 = f [x1,x0]
d2 = f [x2,x1,x0]
⋮
dn = f [xn,  xn−1, ⋯ ,x2,x1,x0].
d0 d1
( ) ( )
O que resulta na expressão do polinômio interpolador de Newton
O interessante do método das diferenças divididas de Newton é que podemos, como desejado, encontrar os
polinômios interpoladores de ordem superior baseados nos coe�cientes calculados anteriormente. Assim,
para facilitar os cálculos, utiliza-se comumente uma tabela recursiva, e, a seguir, para exempli�car, temos a
tabela necessária para encontrar um polinômio interpolador de quarta ordem.
Ordem 0 Ordem 1 Ordem 2 Ordem 3 Ordem 4
       
         
     
       
   
       
     
         
f [xi,xj] =
f (xi) − f (xj)
xi − xj
f [xi,xj,xk] =
f [xi,xj] − f [xj,xk]
xi − xk
f [xi,xj,xk,xl] =
f [xi,xj,xk] − f [xj,xk,xl]
xi − xl
⋮
f [xn,xn−1, ⋯ ,  x1,x0] =
f [xn,xn−1, ⋯ ,x1] − f [xn−1, ⋯ ,x1,x0]
xn − x0
.
Pn (x) = f [x0] + f [x1,x0] (x − x0) + f [x2,x1,x0] (x − x0) (x − x1) + ⋯ + f [xn,  xn−1, ⋯ ,x1,x
x
x0 f [x0]
f [x1,x0]
x1 f [x1] f [x2,x1,x0]
f [x2,x1] f [x3,x2, x1,x0]
x2 f [x2] f [x3,x2,x1] f [x4,x3,x2,x1,x0]
f [x3,x2] f [x4,x3,x2,x1]
x3 f [x3] f [x4,x3,x2]
f [x4,x3]
Ordem 0 Ordem 1 Ordem 2 Ordem 3 Ordem 4
       
Exemplo 1:
Aqui, vamos determinar o polinômio interpolador que interpole os pontos tabelados a seguir
utilizando o polinômio interpolador na forma de Newton.
Como temos um conjunto de três pontos, então devemos encontrar um polinômio de segundo grau. Neste
caso, a forma de Newton para esse polinômio é dada por
Para encontrar os coe�cientes do polinômio, temos que determinar as diferenças divididas e 
. Recursivamente, podemos construir uma tabela de diferenças divididas, que é dada a seguir.
Ordem 0 Ordem 1 Ordem 2
   
     
 
     
   
Observem que, para completar a tabela, foram necessários os cálculos das diferenças divididas 
 e , que têm seus valores calculados a seguir:
Finalmente, temos que o polinômio interpolador é dado por
x
x4 f [x4]
P2 (x)
x −1 0 2
f (x) 4 1 −1
P2 (x) = f [x0] + f [x1,x0] (x − x0) + f [x2,x1,x0] (x − x0) (x − x1).
f [x1,x0]
f [x2,x1,x0]
x
−1 4
−3
0 1 2/3
−1
2 −1
f [x1,x0],
f [x2,x1] f [x2,x1,x0]
f [x1,x0] =
f [x1] − f [x0]
x1 − x0
=
1 − 4
0 + 1
= −3
f [x2,x1] =
f [x2] − f [x1]
x2 − x1
=
−1 − 1
2 − 0
=
−2
2
= −1
f [x2,x1,x0] =
f [x2,x1] − f [x1,x0]
x2 − x0
=
−1 − (−3)
2 − (−1)
=
2
3
.
e seu grá�co pode ser visto na �gura a seguir.
Exemplo 2:
Neste exemplo, vamos estimar o valor de baseado nos valores tabelados a seguir. Para tal, vamos realizar
uma interpolação usando o método de Newton com os nós dados.
Como temos um conjunto de quatro pontos, então devemos encontrar um polinômio de terceiro grau que
interpole esse conjunto de dados. Neste caso, a forma de Newton para o polinômio que interpola os pontos
deve ter a forma
A seguir, temos a tabela formada pelas diferenças divididas que podem ser calculadas usando as fórmulas
P2 (x) = 4 − 3 (x + 1) +
2
3
(x + 1)x
= 4 − 3x − 3 +
2
3
x2 +
2
3
x
= 1 −
7
3
x +
2
3
x2,
e2
x −1 0 1 3
f (x) 0.367879 1 2.718281 20.085536
P3 (x) = f [x0] + f [x1,x0] (x − x0) + f [x2,x1,x0] (x − x0) (x − x1) + f [x3,x2,x1,x0] (x − x0) (x
f [xi,xj] =
f (xi) − f (xj)
xi − xj
f [xi,xj,xk] =
f [xi,xj] − f [xj,xk]
xi − xk
$x$ Ordem 0 Ordem 1 Ordem 2 Ordem 3
     
       
   
     
   
       
     
Portanto, neste caso o polinômio de Newton é dado por
Como queremos um valor aproximado de , então precisamos calcular . Como 
, então o erro de aproximação é de
\(\approx 0.801716.\]
A seguir, podemos ver no grá�co como o polinômio interpolador se aproxima da função exponencial e o erro
envolvido na interpolação entre os pontos e .
A seguir, temos o código para determinar o valor polinômio interpolador através do MATLAB para um
conjunto de pontos dados .
1 % funcao para calcular o polinomio interpolador P(x) considerando 
2 % os dados tabelados (xi,yi) atraves do metodo das diferencas divididas 
f [xi,xj,xk,xl] =
f [xi,xj,xk] − f [xj,xk,xl]
xi − xl
.
−1 0.367879
0.632121
0 1 0.543080
1.718281 0.444675
1 2.718281 2.321782
8.6836275
3 20.085536
P3 (x) = 0.367879 + 0.632121 (x + 1) + 0.543080 (x + 1)x + 0.444675 (x + 1)x (x − 1)
= 1 + 0.730526x + 0.54308x2 + 0.444675x3.
e2 P3 (2) = 8.190770
e2 ≈ 7.389056
erro= P3 (2) − e2∣ ∣x = 1 x = 3 P (x){(xi, fi)}ni=1
3 % de newton 
4 
5 
6 
 
7 function P = interp_newton (xi,yi,x) 
8 
9 
10 
11 % entrada – xi, yi matrizes colunas contendo os nos 
12 % x o ponto no qual o valor do polinomio sera calculado 
13 % 
14 % 
15 % saida – P o valor do polinonio no ponto x 
16 
17 
18 
19 n = length(x); 
20 a(1) = y(1); 
21 for k = 1 : n - 1 
22 d(k, 1) = (yi(k+1) - yi(k))/(xi(k+1) - xi(k));
23 end 
24 for j = 2 : n - 1 
25 for k = 1 : n – j
26 % calcula a tabela de diferencas 
27 d(k, j) = (d(k+1, j - 1) - d(k, j - 1))/(xi(k+j) - xi(k)); 
28 end 
29 end 
30 
31 
32 
33 for j = 2 : n 
34 a(j) = d(1, j-1); 
35 end 
36 Df(1) = 1; 
37 c(1) = a(1); 
38 for j = 2 : n 
39 Df(j)=(x - xi(j-1)) .* Df(j-1); % monta o polinomio P(x) 
40 c(j) = a(j) .* Df(j); 
41 end 
42 P = sum(c);
Unidade 04
Aula 03
Spline Cúbica
Sejam bem-vindos a mais uma unidade! Estamos cada vez mais próximos do �m e, com isso, aumentando o nosso
conhecimento.
Vimos como determinar uma função polinomial de grau que passa por pontos tabelados
usando as formas de Newton e Lagrange. Ocorre que, quando o grau do polinômio interpolador, é muito alto,
o resultado dessa aproximação dada pela interpolação pode ser desastroso, como podemos ver na �gura a
seguir. A curva em vermelho nesta �gura é a chamada função de Runge, e à medida que essa curva é
interpolada por polinômios com graus cada vez maiores, várias oscilações crescentes começam a surgir nas
bordas. Neste caso da �gura, a curva azul é dada por um polinômio de grau e a curva em verde é dada por
um polinômio de grau .
pn (x) n, n + 1
5
9
SAIBA MAIS
Para mais exemplos do método de Newton, clique aqui. Con�ra! 
No seguinte endereço, é possível assistir a videoaula sobre o assunto: Clicando aqui.
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d3730728e4c5b4474706e8b_texto05.pdf
https://www.youtube.com/watch?v=fr1oGgA7QxU
Essas oscilações que ocorrem nas bordas do intervalo são conhecidas como fenômeno de Runge e elas
acontecem, como podemos ver, quando se usa interpolação polinomial com polinômios de ordem muito
elevada. Uma das alternativas que podemos utilizar para mitigar esse problema é interpolar uma função 
tabelada em grupos de poucos pontos, obtendo-se, assim, para cada subintervalo um polinômio de
grau menor, evitando, desse modo, as oscilações. Além disso, impõem-se condições de continuidade e
suavidade sobre a função interpolante por partes, isto é, deseja-se que a aproximação seja contínua e tenha
derivadas contínuas até uma certa ordem.
De forma geral, dada uma tabela de dados
em que \({{x}_{0}}
1. A função spline passa pelos nós dados, isto é, , para todo .
2. Em cada um dos subintervalos a função spline é um polinômio de grau .
3. A função spline tem derivada contínua até a ordem dentro do intervalo na qual a
função é tabelada.
Exemplo 1:
Vamos, nesse exemplo, considerar o caso de interpolação por spline para o caso mais simples, aquele em que
a função de interpolação em cada subintervalo é linear. Sendo assim, chamaremos a função que irá
interpolar os pontos , e de spline linear e ela formada em cada subintervalo por
polinômios de grau . Dessa maneira, para satisfazer as três características de uma spline, como de�nido
anteriormente, a função deve ser uma função por partes dada na forma
f (x) 
x0 x1 x2 ⋯ xn
f0 f1 f2 ⋯ fn
s (x) s (xi) = fi i = 0, 1, ⋯ ,n
[xi−1,xi], s (x) k
s (x) k − 1 [x0,xn],
s (x)
(1, −1) (2, 4) (3, 1)
k = 1
s (x)
s (x) = {
a1x + b1,    x ∈ [x1,x2]
a2x + b2,    x ∈ [x2,x3]
Observe que, neste caso linear, não precisamos impor nenhuma condição sobre a derivada, pois, pela
de�nição da spline, a derivada de ordem já é contínua.
Sendo cada uma das funções que de�ne uma reta, então para escrevê-las respectivamente, são
necessários dois pontos apenas. O modo mais simples de se obter a equação dessas retas é usar para cada
uma delas a forma do polinômio interpolador de Lagrange. Logo, a spline linear para três pontos pode ser
escrita de forma geral como
Portanto, para o conjunto de pontos dados, temos que a spline linear é dada por
Apesar de podermos de�nir uma spline com uma ordem qualquer, na prática o caso mais utilizado da
interpolação por uma spline é o caso cúbico. Isso se deve ao fato da spline cúbica ter um bom grau de
suavidade, o que costuma ser satisfatório nas aplicações práticas. Assim, dada a tabela de nós como a seguir,
uma spline cúbica que interpola esses dados é uma função de�nida por partes como segue
em que cada uma das funções que a de�ne é um polinômio cúbico 
 que satisfaz a condição , para todo . Além disso, para garantir que a
função seja suave nos nós dados (item iii da de�nição), a função deve ter primeira e segunda
derivadas contínuas em todo o intervalo , isto é, as condições e 
 para , devem ser satisfeitas para cada nó .
Exemplo 2:
Tendo em vista que o caso da spline cúbica é o mais utilizado na prática, vamos, nesse exemplo, considerar a
seguinte tabela para encontrar uma spline cúbica que passe pelos pontos dados.
0
s (x)
s (x) =
⎧
⎨⎩
f1
(x−x2)
x1−x2
+ f2
(x−x1)
x2−x1
,    x ∈ [x1,x2]
f2
(x−x3)
x2−x3
+ f3
(x−x2)
x3−x2
,    x ∈ [x2,x3]
s (x) = {
(−1) ⋅ x−21−2 + 4 ⋅
x−1
2−1 ,    x ∈ [1, 2]
4 ⋅ x−32−3 + 1 ⋅
x−2
3−2 ,    x ∈ [2, 3]
= { .
5x − 6,    x ∈ [1, 2]
−3x + 10,    x ∈ [2, 3]
x0 x1 x2 ⋯ xn
f0 f1 f2 ⋯ fn
s (x)
s (x) = ,
⎧⎪⎨⎪⎩ s0 (x),  x ∈ [x0,x1]s1 (x),  x ∈ [x1,x2]⋮sn−1 (x),  x ∈ [xn−1,xn]si (x) = aix3 + bix2 + cix + di(ai ≠ 0) si (xi) = fi i = 0, 1, ⋯ ,ns (x) s (x)[x0,xn] si ′ (xi) = si+1′ (xi)si ′′ (xi) = si+1′′ (xi), i = 0, 1, ⋯ ,n xix 2 3 5
Desta forma, de�nimos um polinômio cúbico em cada um dos subintervalos que se formam com os pontos
dados. Portanto, a função interpoladora deve ser dada por
Assim como nas demais técnicas de interpolação estudadas, é necessário que a spline cumpra a condição de
passar pelos pontos dados na tabela. Logo, temos que
Observe anteriormente que no ponto que se encontra no interior do intervalo , a spline deve
cumprir que cada uma das funções que a de�ne, e , seja igual neste ponto. Caso contrário,
haveria uma perda da continuidade neste ponto. Agora, para garantirmos a suavidade da spline, calculamos a
primeira derivada da função que nos dá
Ao fazermos fomos capazes de garantir a continuidade da spline no ponto Para
garantirmos a suavidade neste ponto, é necessário também que as derivadas neste ponto sejam iguais.
Desse modo, fazemos que fornece a seguinte condição
que é o mesmo que
Analogamente, procedemos com a segunda derivada que nos dá
e pela condição que a derivada também seja suave no ponto , então
y −1 2 −7
s (x) = {
s1 (x),  x ∈ [x1,x2]
s2 (x),  x ∈ [x2,x3]
= { .
a1x
3 + b1x2 + c1x + d1,  x ∈ [2, 3]
a2x
3 + b2x2 + c2x + d2,  x ∈ [3, 5]
s (2) = −1 ⇒ 8a1 + 4b1 + 2c1 + d1 = −1,
s (3) = 2 ⇒ { ,
27a1 + 9b1 + 3c1 + d1 = 2
27a2 + 9b2 + 3c2 + d2 = 2
s (5) = −7 ⇒ 125a5 + 25b2 + 5c2 + d2 = −7.
x = 3, [2, 5]
s1 (x) s2 (x)
s (x),
s′ (x) = { .
3a1x2 + 2b1x + c1,  x ∈ [2, 3]
3a2x2 + 2b2x + c2,  x ∈ [3, 5]
s1 (3) = s2 (3), x = 3.
s1
′ (3) = s2 ′ (3)
3a1(3)
2 + 2b1 (3) + c1 = 3a2(3)
2 + 2b2 (3) + c2
27a1 + 6b1 + c1 = 27a2 + 6b2 + c2.
s (x) = {
6a1x + 2b1,  x ∈ [2, 3]
6a2x + 2b2,  x ∈ [3, 5]
s′ (x) x = 3
6a1 (3) + 2b1 = 6a2 (3) + 2b2 ⇒ 18a1 + 2b1 = 18a2 + 2b2.
Neste ponto, contamos com 6 equações e 8 incógnitas, portanto temos ainda 2 graus de liberdade. Assim,
para termos um sistema com única solução, precisamos de mais duas condições para a spline. De forma geral,
agregam-se ao sistema as duas seguintes condições de fronteira
Essas condições de fronteiran são chamadas de condições de fronteira naturais. Para elas, obtemos as
condições
Finalmente, completamos um sistema quadrado com 8 equações e 8 incógnitas dado pelo seguinte conjunto
de equações
Esse sistema pode serescrito em forma matricial, o que nos dá:
s′′ (x0) = 0   e   s′′ (xn) = 0.
s′′ (2) = 0 ⇒ 12a1 + 2b1 = 0
s′′ (5) = 0 ⇒ 30a2 + 2b2 = 0.
8a1 + 4b1 + 2c1 + d1 = −1,
27a1 + 9b1 + 3c1 + d1 = 2,
27a2 + 9b2 + 3c2 + d2 = 2,
125a5 + 25b2 + 5c2 + d2 = −7,
27a1 + 6b1 + c1 = 27a2 + 6b2 + c2,
18a1 + 2b1 = 18a2 + 2b2,
12a1 + 2b1 = 0,
30a2 + 2b2 = 0.
É possível veri�car, com algum trabalho, que a matriz dos coe�cientes é inversível, então podemos inverter o
sistema para obtermos a seguinte solução
Substituindo esses valores na função inicial, vemos que a spline cúbica que satisfaz os dados da tabela é dada
por
A seguir, temos o grá�co da função e os nós que de�nem essa spline.
Note-se neste exemplo a delicadeza com que os polinômios cúbicos se ligam na spline. Praticamente não é
possível veri�car que se trata de dois polinômios diferentes. Este fato se deve às condições impostas sobre
as derivadas da função . Esta delicadeza quase artística é o que permite aplicar as splines cúbicas em
diversas questões em que a interpolação é necessária, dentro de um problema cuja natureza seja bastante
sensível.
=
⎡⎢⎣ 8 4 2 1 0 0 0 027 9 3 1 0 0 0 00 0 0 0 27 9 3 10 0 0 0 125 25 5 127 6 1 0 27 6 1 018 2 0 0 18 2 0 012 2 0 0 0 0 0 00 0 0 0 30 2 0 0⎤⎥⎦⎡⎢⎣a1b1c1d1a2b2c2d2⎤⎥⎦ ⎡⎢⎣−122−70000 ⎤⎥⎦a1 = −1.25b1 = 7.5c1 = −10.75d1 = 0.5a2 = 0.625b2 = −9.375c2 = 39.875d2 = −50.125.s (x) = { −1.25x3 + 7.5x2 − 10.75x + 0.5,  x ∈ [2, 3]0.625x3 − 9.375x2 + 30.875x − 50.125,  x ∈ [3, 5]s (x)
s (x)
Com o código a seguir, calculam-se os coe�cientes da spline cúbica para um conjunto de pontos 
1 % funcao para calcular os coe�cientes da spline interpolante 
2 % que passa pelos pontos (x,y) usando condicao de contorno natural 
3 % 
4 % a spline e de�nida como um polinomio cubico por partes na forma: 
5 % S(x) = { Sk(x) x(k) <= x <= x(k+1) 
6 % 
7 % cada polinomio cubico Sk(x) e dado na forma: 
8 % Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3 
9 
10 
11 function [s0, s1, s2, s3] = spline_cubica(x,y) 
12 % entrada – x, y matrizes colunas contendo os nos 
13 % 
14 % saida – os coe�cientes sk0, sk1, sk2, e sk3 para cada um dos polinomios 
15 % e retornado atraves dos vetores s0,s1,s2, e s3, respectivamente 
16 
17 
18 n = length(x); 
19 h = x(2:n) - x(1:n-1); 
20 d = (y(2:n) - y(1:n-1))./h; 
21 
22 
23 inferior = h(1:end-1); 
24 principal = 2*(h(1:end-1) + h(2:end)); 
25 superior = h(2:end); 
26 
27 
28 T = spdiags([inferior principal superior], [-1 0 1], n-2, n-2); 
29 ld = 6*(d(2:end)-d(1:end-1)); 
30 
31 
32 m = T\ld; 
33 
34 
35 % uso das condicoes de contorno naturais 
36 m = [ 0; m; 0]; 
37 
38 
39 s0 = y; 
40 s1 = d - h.*(2*m(1:end-1) + m(2:end))/6; 
{(xi, fi } ni=1.
41 s2 = m/2; 
42 s3 =(m(2:end)-m(1:end-1))./(6*h);
Unidade 04
Aula 04
Integração Numérica
Olá alunos! Nesta aula falaremos sobre as fórmulas de Newton-Cotes. Boa aula!
Se em alguma aplicação dos nossos conhecimentos nos depararmos, por exemplo, com alguma das integrais
de�nidas
teremos alguma di�culdade para calcular a área que elas representam. Em ambos os casos temos que o
integrando é contínuo, logo existe uma primitiva. Porém, em nenhum deles seremos capazes de calcular essa
primitiva. Então, para encontrarmos qual o valor aproximado de cada uma dessas integrais será fundamental
construirmos métodos numéricos para isso.
1
∫
−1
e−x
2
dx ou 
1
∫
0
√1 − x3dx,
SAIBA MAIS
Para uma outra aula e mais exemplos resolvidos, clique aqui e acesse. 
Uma videoaula interessante sobre a spline linear: Clique aqui.
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d37332c7822ab09e785270b_texto06.pdf
https://www.youtube.com/watch?v=EN6Dk6p6_ls
Nesta aula serão estudados dois métodos numéricos para determinar a aproximação de uma integral
de�nida, a regra dos trapézios e a regra de Simpson.
Esses métodos se encaixam numa categoria de métodos numéricos conhecidos como fórmulas de Newton-
Cotes. Essas fórmulas fornecem uma aproximação de uma integral de�nida baseadas na avaliação do
integrando em pontos igualmente espaçados dentro do domínio de integração.
Além das duas fórmulas que serão apresentadas nesta aula, existem outras diversas outras fórmulas de
Newton-Cotes e algumas das mais conhecidas serão apresentadas ao �m desta aula. Começaremos aqui com
a regra dos trapézios para calcularmos aproximadamente uma integral de�nida.
Considere a integral de�nida
e suponha que seja contínua em . Para o nosso cálculo aproximado da integral iremos dividir o
intervalo em subintervalos de tamanhos iguais
Essa subdivisão irá nos dar pontos igualmente espaçados
que são as fronteiras dos subintervalos. Em cada um desses pontos iremos calcular o valor da função ,
que nos dá
Nós aproximamos a integral usando trapézios formados pela união dos pontos e por
linhas retas, para , como podemos ver na �gura abaixo.
Como sabemos da geometria plana, a área de um trapézio é dada por
b
∫
a
f (x)dx
f (x) [a, b]
[a, b] n
Δx =
b − a
n
.
n + 1
x0 = a,  x1 = x0 + Δx,  x2 = x0 + 2Δx,   ⋯ ,  xn = x0 + nΔx = b,
f (x)
y0 = f (x0),  y1 = f (x1), ⋯ ,  yn = f (xn).
n (xi−1, yi−1) (xi, yi)
1 ≤ i ≤ n
AT =
1
2
(base maior + base menor) × altura,
e para cada trapézio da �gura, temos que
Finalmente, se adicionarmos as áreas de todos os trapézios, obtemos então a aproximação da integral
de�nida
A seguir temos um exemplo de como utilizar esta fórmula.
Exemplo 1
Neste exemplo vamos usar a regra dos trapézios para aproximar a integral de�nida
com subdivisões do intervalo . Primeiramente, podemos a�rmar que, apesar de muito
trabalhoso, é possível encontrar a primitiva desta integral e encontrar o seu valor exato. Neste caso, a
primitiva da integral acima é dada pela função abaixo
na qual é uma constante de integração e é a função inversa do seno hiperbólico de�nido por
Logo, a integral de�nida pode ser calculada e seu valor exato é
Agora para calcularmos aproximadamente o valor da integral, temos inicialmente que calcular o
comprimento de cada subintervalo que será dado por
ATi =
Δx
2
(yi−1 + yi).
n
b
∫
a
f (x)dx ≈
n
∑
i=1
Δx
2
(yi−1 + yi)
=
Δx
2
(y0 + y1) +
Δx
2
(y1 + y2) + ⋯ +
Δx
2
(yn−1 + yn)
=
Δx
2
(y0 + 2y1 + 2y2 + ⋯ + 2yn−1 + yn).
5
∫
1
√1 + x2dx
n = 8 [1, 5]
∫ √1+x
2
dx =
x
2
√1 + x2 +
1
2
arcsenh (x) + C,
C arcsenh (x)
senh (x) =
ex − e−x
2
 .
5
∫
1
√1 + x2dx =
1
2
[5√26 − √2 + arcsenh (5) − arcsenh (1)] ≈ 12.755974.
Δx =
5 − 1
8
= 0.5.
Pela regra dos trapézios, temos que a aproximação deve ser dada na forma do seguinte somatório
em que os pontos são dados por e Para facilitar os nossos cálculos, tabelamos
os valores da função nos pontos que dividem o intervalo.
Assim, usando os valores da tabela, podemos calcular uma aproximação com duas casas decimais da integral
que será dada por
Observe que o valor aproximado �cou muito perto do valor real da integral. Este comportamento pode ser
justi�cado pelo fato da função ser aproximadamente linear, como podemos ver no grá�co abaixo.
Tendo em vista que a regra dos trapézios aproxima a função linearmente, temos a semelhante entre os
resultados usando poucos trapézios.
Na regra dos trapézio, discretizamos o intervalo de integração em subintervalos e
aproximamos o integrando por segmentos de retas a cada dois pontos. Desta forma, teríamos que calcular
apenas áreas de trapézios para encontramos uma aproximação para a integral. Na regra de Simpson, o
princípio é semelhante: substituir o integrando por outro, de forma que podemos encontrar a área abaixo
desta função aproximada facilmente. Como na regra dos trapézios a aproximação do integrando é por uma
função linear, na regra de Simpson a aproximação do integrando será por uma função quadrática.
5
∫
1
√1 + x2dx ≈
Δx
2
7
∑
i=0
[√1 + x2i +√1 + x2i+1],
xi x0 = 1 xi = x0 + i ⋅ Δx.
xi
x0 x1 x2 x3 x4 x5 x6 x7 x8
x 1 1.5 2 2.5 3 3.5 4 4.5 5
√1 + x2 √2 √3.25 √5 √7.25 √10√13.25 √17 √21.25 √26
5
∫
1
√1 + x2dx ≈
Δx
2
7
∑
i=0
[√1 + x2i +√1 + x2i+1]
=
0.5
2
(√2 + 2√3.25 + 2√5 + 2√7.25 + 2√10 + 2√13.25 + 2√17 + 2√21.25 + √26)
≈ 12.76.
√1 + x2
[a, b] [xi,xi+1]
Assim, para determinarmos a regra de Simpson vamos inicialmente encontrar a área abaixo do grá�co da
parábola de equação que passa pelos três pontos e , como
vemos no grá�co abaixo. Observe que são necessários três pontos para determinar a equação de uma
parábola, diferentemente da regra dos trapézios em que eram necessários apenas dois.
Temos que a área abaixo do grá�co da parábola é dada por
Como os pontos e estão sobre a parábola, eles satisfazem a equação 
 Portanto, podemos encontrar as constantes e baseadas nesses pontos que são
dados por
Agora, observe que a combinação dos pontos e satisfaz
Finalmente, ao compararmos a soma acima, com a área abaixo do grá�co da parábola, temos que esta área
pode ser escrita a partir dos pontos e na forma
Para deduzirmos a regra de Simpson, considere novamente a integral de�nida
y = ax2 + bx + c (−Δx, y0),   (0, y1) (Δx, y2)
A =
Δx
∫
−Δx
(ax2 + bx + c)dx = ( ax
3
3
+
bx2
2
+ cx)
Δx
−Δx
= 2aΔx3 + 2cΔx =
Δx
3
(2aΔx2 + 6c).
(−Δx, y0),   (0, y1) (Δx, y2)
y = ax2 + bx + c. a,  b c
y0 = aΔx
2 − bΔx + c
y1 = c
y2 = aΔx2 + bΔx + c.
y0, y1 y2
y0 + 4y1 + y2 = (aΔx2 − bΔx + c) + 4c + (aΔx2 + bΔx + c) = 2aΔx2 + 6c.
y0, y1 y2
A =
Δx
3
(2aΔx2 + 6c) =
Δx
3
(y0 + 4y1 + y2).
b
∫
a
f (x)dx
e suponha que a função seja contínua em . Para o nosso cálculo aproximado da integral de�nida,
considerando uma aproximação de por uma parábola, devemos dividir o intervalo em um número
par de subintervalos de comprimentos iguais
que nos levará aos pontos igualmente espaçados
que são as fronteiras desses subintervalos. Assim como �zemos na regra dos trapézios, iremos calcular o
valor da função nos pontos discretizados para obtermos
Finalmente, podemos estimar o valor da integral simplesmente adicionando as áreas abaixo dos arcos
parabólicos a cada três pontos sucessivos, como podemos ver na �gura abaixo. Isso nos levará a então
conhecida como regra de Simpson
Exemplo 2
Vamos agora usar a regra de Simpson para aproximar a integral de�nida
com subdivisões do intervalo . Essa integral calcula a área entre o eixo e a elipse 
com , como podemos ver na �gura abaixo.
f (x) [a, b]
f (x) [a, b]
n
Δx =
b − a
n
,
n + 1
x0 = a,  x1 = x0 + Δx,  x2 = x0 + 2Δx,   ⋯ ,  xn = x0 + nΔx = b,
f (x)
y0 = f (x0),  y1 = f (x1), ⋯ ,  yn = f (xn).
b
∫
a
f (x)dx ≈
n
2 −1
∑
i=0
Δx
3
(y2i + 4y2i+1 + y2i+2)
=
Δx
3
(y0 + 4y1 + y2) +
Δx
3
(y2 + 4y3 + y4) + ⋯ +
Δx
3
(yn−2 + 4yn−1 + yn)
=
Δx
3
(y0 + 4y1 + 2y2 + 4y3 + 2y4 + ⋯ + 4yn−1 + yn).
4
∫
1
√4 −
x2
4
dx
n = 6 [1, 4] x x
2
16 +
y2
4 = 1
1 ≤ x ≤ 4
Assim como no primeiro exemplo, é possível determinar a primitiva desta integral e também seu valor exato.
Neste caso, a primitiva da integral acima é dada pela função abaixo
na qual é a função inversa do seno e é uma constante de integração. Logo, a integral de�nida
pode ser calculada e seu valor exato é
Para usarmos a regra de Simpson, começamos calculando o comprimento de cada subintervalo que será
dado por
Assim, a aproximação desejada deve ser dada na forma do seguinte somatório
em que os pontos e são dados por , e . Para facilitar os nossos
cálculos, tabelamos os valores da função nos pontos que dividem o intervalo.
Assim, usando os valores da tabela, temos que a aproximação, com duas casas decimais, da integral é dada
por
∫
√4− x24
dx =
x
2
√4 − x
2
4
+ 4arcsen( x
4
)+ C,
arcsen (x) C
4
∫
1
√4 − x
2
4
dx = 2π −
√15
4
− 4arcsen( 1
4
) ≈ 4.3042.
Δx =
4 − 1
6
= 0.5.
4
∫
1
√4 − x
2
4
dx ≈
Δx
3
2
∑
i=0
[y2i + 4y2i+1 + y2i+2],
xi yi x0 = 1 xi = x0 + i ⋅ Δx  yi = f (xi)
xi
x0 x1 x2 x3 x4 x5 x6
x 1 1.5 2 2.5 3 3.5 4
√4 − x
2
4
√3.75 √3.4375 √3 √2.4375 √1.75 √0.9375 0
4
∫
1
√4 − x
2
4
dx ≈
Δx
3
2
∑
i=0
[y2i + 4y2i+1 + y2i+2]
5
Em parte pelo seu caráter didático, nesta aula deduzimos duas das fórmulas de Newton-Cotes, a regra dos
trapézios e a regra de Simpson. No entanto, existem outras fórmulas de Newton-Cotes conhecidas que
podemos citar e as mais conhecidas são a regra de Simpson
e a regra de Boole
que são construídas também baseadas em polinômios interpoladores de Lagrange, neste caso, polinômios de
grau 3 e 4, respectivamente. Observe que a regra dos trapézios é uma regra que depende de dois pontos, a
regra de Simpson depende de 3 pontos, a regra de Simpson depende de 4 pontos e, �nalmente, a regra
de Boole depende de 5 pontos.
Abaixo temos dois códigos em MATLAB. Um para calcular o valor aproximado de uma integral de�nida
usando a regra dos trapézios e o outro usando a regra de Simpson. Primeiro, segue a regra dos trapézios.
1 % funcao para calcular a numericamente a integral de f(x) no intervalo
2 % [a,b] usando a regra dos trapezio
3
4 function I = trapezios(f, a, b, n) 
5
6 % entrada – f funcao f(x) a ser integrada
7 % a, b intervalo de integracao
8 % n numero de subdivisoes no intervalo
9 %
10 % saida – I valor aproximado da integral
11
12 dx = (b-a)/n;
13 x = a:dx:b;
14 I = 0;
15
16 for i = 1:n
17 I = I + (dx/2)*(f(x(i))+f(x(i+1))); % calcula a integral aproximada
18 end
=
0.5
3
(√3.75 + 4√3.4375 + 2√3 + 4√2.4375 + 2√1.75 + 4√0.9375 + 0)
≈ 4.26.
3/8
x3
∫
x0
f (x)dx ≈
Δx
8
(y0 + 3y1 + 3y2 + y3),
x4
∫
x0
f (x)dx ≈
Δx
90
(7y0 + 32y1 + 12y2 + 32y3 + 7y4),
3/8
 
Em seguida, a regra de Simpson.
1 % funcao para calcular a numericamente a integral de f(x) no intervalo
2 % [a,b] usando a regra de simpson
3
4 function I = rsimpson(f, a, b, n) 
5
6 % entrada – f funcao f(x) a ser integrada
7 % a, b intervalo de integracao
8 % n numero de subdivisoes no intervalo (n deve ser par)
9 %
10 % saida – I valor aproximado da integral
11
12 dx = (b-a)/n;
13 x = a:dx:b;
14 I = 0;
15
16 for i = 1:n/2
17 I = I + (dx/3)*(f(x(2*i-1))+4*f(x(2*i))+f(x(2*i+1))); % calcula a integral aproximada
18 end
 
SAIBA MAIS
Nesse endereço pode ser assistida uma aula sobre a regra dos Trapézios: Clique aqui. 
Nesse endereço pode ser assistida uma aula sobre a regra de Simpson: Clique aqui. 
É possível encontrar nesse endereço uma lista com várias fórmulas de Newton-Cotes: Newton-Cotes
Formulas
https://www.youtube.com/watch?v=HbhO6JCV2Jc
https://www.youtube.com/watch?v=0it9LQkpmtU
http://mathworld.wolfram.com/Newton-CotesFormulas.html
Olá! Iniciamos mais uma aula. Continue seus estudos e se torne um pro�ssional capacitado!
Nesta aula nosso foco ainda é determinar métodos de aproximar uma integral de�nida. Para construir o
nosso novo método vamos, especi�camente, encontrar a área abaixo do grá�co da curva dentro
do intervalo , isto é, calcular a integral de�nida
Nós já vimos a regra dos trapézios, com ela podemos calcular aproximadamente a integral utilizando apenas
duas avaliações da função nos seus pontos extremos, neste caso usamos apenas os pontos extremos, neste
caso usamos apenas os pontos e . Mas, se o grá�co da função for côncavo
para baixo, o erro de aproximação será toda a região que está entre a curva e o segmento de reta unindo os
pontos, o que pode ser um erro muito grande considerando a região. Abaixo no grá�co podemos ver como
seria o erro gerado pela aproximação de uma integral pela regra dos trapézios.
Se diferentemente de utilizarmos os extremos do intervalo, pudéssemos usar dois pontos e que estão
no interior do intervalo , então a reta que passa pelos pontos e irá cruzar a
curva e a área abaixo da reta irá gerar uma melhor aproximação da área abaixo da curva, como podemos ver
na �gura abaixo.
y = f (x)
−1 ≤ x ≤ 1
A =
1
∫
−1
f (x)dx.
(−1, f (−1)) (1, f (1)) y = f (x)
x1 x2
[−1, 1] (x1, f (x1)) (x2, f (x2))
Considerando a reta que passa pelos pontos e , temos que sua equação é dada por
A área do trapézio formado pela reta dada pela equaçãocima é
Perceba que a regra dos trapézios é um caso especial desta última equação, pois se �zermos , 
 e , então a área abaixo do grá�co pela regra dos trapézios é
Podemos usar o método dos coe�cientes a determinar para encontrar as abcissas de coordenadas , , e
os pesos e de tal forma que a integral de�nida seja aproximada pela fórmula
e que essa integral seja exata para polinômios cúbicos (i.e., ). Como
quatro coe�cientes , , e precisam ser determinados na equação acima, podemos escolher quatro
condições para serem satisfeitas. Usando o fato que a integração é aditiva, será su�ciente pedir que a
equação acima seja exata para as quatro funções , , e . As
quatro condições integrais são
(x1, f (x1)) (x2, f (x2))
y = f (x1) +
f (x2) − f (x1)
x2 − x1
(x − x1).
Atrap =
2x2
x2 − x1
f (x1) −
2x2
x2 − x1
f (x2).
x1 = −1
x2 = 1 h = 2
Atrap =
2
2
f (x1) −
−2
2
f (x2) = f (x1) + f (x2).
x1 x2
w1 w2
1
∫
−1
f (x)dx ≈ w1f (x1) + w2f (x2)
f (x) = a3x3 + a2x2 + a1x + a0
w1 w2 x1 x2
f1 (x) = 1 f2 (x) = x f3 (x) = x2 f4 (x) = x3
f1 (x) = 1 ⇒
1
∫
−1
f1 (x)dx =
1
∫
−1
1dx = 2 = w1 + w2
Temos agora agora um sistema não linear de quatro equação e quatro incógnitas a ser resolvido dados por
Como , então substituindo na última equação, temos
Como deve ser diferente de , então temos que Assim, usando esta última condição
podemos encontrar que os pesos devem ser iguais e dados por
 
Substituindo esta última equação na equação , temos
Além disso, temos
 
Portanto, encontramos as coordenadas das abscissas que fornecem uma aproximação cúbica para a função é
 
f2 (x) = x ⇒
1
∫
−1
f2 (x)dx =
1
∫
−1
xdx = 0 = w1x1 + w2x2
f3 (x) = x2 ⇒
1
∫
−1
f3 (x)dx =
1
∫
−1
x2dx =
2
3
= w1x21 + w2x
2
2
f4 (x) = x4 ⇒
1
∫
−1
f4 (x)dx =
1
∫
−1
x3dx = 0 = w1x31 + w2x
3
2.
w1 + w2 = 2
w1x1 + w2x2 = 0
w1x
2
1 + w2x
2
2 =
2
3
w1x
3
1 + w2x
3
2 = 0.
w1x1 = −w2x2
x21 = x
2
2 ⇒ x1 = ±x2.
x1 x2 x1 = −x2.
w1 = w2.
w1 + w2 = 2
w1 = w2 = 1.
w1x
2
1 + w2x
2
2 = x
2
2 + x
2
2 =
2
3
⇒ x22 =
1
3
.
x1 = −
1
√3
    e    x2 =
1
√3
.
Encontramos os nós e os pesos que compõem a quadratura Gaussiana de dois pontos. Uma vez que a
fórmula é exata para equações cúbicas, pode ser mostrado que o termo de erro envolverá uma derivada de
quarta ordem.
Finalmente, a regra da quadratura Gaussiana para dois pontos pode ser então escrita como: dada uma
função é contínua no intervalo , então a integral de�nida é aproximada com erro de quarta
ordem por
Exemplo 1
Neste exemplo, vamos usar a quadratura Gaussiana com dois pontos para calcular uma aproximação para a
interal de�nida
e comparar o resultado com a regra dos trapézios, com e a regra de Simpson, com .
Primeiramente, temos que essa integral pode ser facilmente calculada exatamente e é dada por
Pela quadratura Gaussiana, temos que a integral dada por ser aproximada por
Pela regra dos trapézios com , a integral é aproximada por
Finalmente, pela regra de Simpson com , temos que integral é aproximada por
f (x) [−1, 1]
1
∫
−1
f (x)dx ≈ f (− 1
√3
)+ f ( 1
√3
).
1
∫
−1
1
x + 2
dx
h = 2, h = 1
1
∫
−1
1
x + 2
dx = [ln (x + 2)]1−1 = ln (3) − ln (1) ≈ 1.09861.
1
∫
−1
1
x + 2
dx ≈
1
2 − 1
√3
+
1
2 + 1
√3
≈ 0.70291 + 0.38800 = 1.09091.
h = 2
1
∫
−1
1
x + 2
dx ≈
1
2 − 1
+
1
2 + 1
≈ 1.00000 + 0.33333 = 1.33333.
h = 1
1
∫
−1
1
x + 2
dx ≈
1
3
[ 1
2 − 1
+
4
2 + 0
+
1
2 + 1
] ≈ 1.11111.
Neste caso, os erros de aproximação são dados respectivamente por , e .
Mostrando assim que o método da quadratura Gaussiana fornece uma melhor aproximação para a integral
de�nida quando comparada aos métodos de Newton-Cotes.
Apesar da quadratura Gaussiana para dois pontos ter sido de�nida dentro do intervalo , é possível
utilizar esse método para determinar a aproximação da área de uma função dentro de um intervalo qualquer
 . Para tal, fazemos a mudança de variáveis
 
com isso a integral no intervalo pode ser transformada numa integral no intervalo na forma
Agora, aplicando a quadratura Gaussiana temos que a aproximação de uma integral no intervalo é dada
por
Exemplo 2
Agora vamos usar a quadratura Gaussiana com dois pontos para calcular a área abaixo do grá�co da função 
 no intervalo , como podemos ver no grá�co abaixo. Isto é, queremos calcular a
aproximação para a integral de�nida
Neste caso, a integral de�nida pode ser calculada de forma exata e é dada por
Para usarmos a quadratura Gaussiana, temos que e assim realizamos a mudança de variáveis
de tal forma que
0.00770 0.23472 0.01250
[−1, 1]
[a, b]
t =
a + b
2
+
b − a
2
x
[a, b] [−1, 1]
b
∫
a
f (t)dt =
1
∫
−1
f ( a + b
2
+
b − a
2
x)
(b − a)
2
dx.
[a, b]
b
∫
a
f (t)dt ≈
(b − a)
2
[f ( a + b
2
−
b − a
2
1
√3
)+ f ( a + b
2
+
b − a
2
1
√3
)].
f (x) = √1 + x2 [−1, 3]
3
∫
−1
√1 + x2dx.
3
∫
−1
√1 + x2dx =
1
2
(√1 + x2 + arcsenh (x))
3
−1
≈ 6.8004.
a = −1, b = 3
3
∫
−1
√1 + x2dx =
1
∫
−1
√1 + (1 + 2x)2dx.
Finalmente, temos que a integral de�nida é aproximada por
que fornece um erro de com relação ao resultado exato da integral.
Ao longo desta aula deduzimos a regra da Quadratura Gaussiana para dois pontos, porém é possível
escrever a regra da quadratura para três, ou mais pontos. No caso da regra para três pontos, é possível
mostrar que a integral de�nida de no intervalo é aproximada por
Além disso, no caso de uma integral de�nida dentro do intervalo , realizando a mudança de variáveis
mostrada anteriormente, teríamos a regra da quadratura Gaussiana de três pontos na forma
Exemplo 3
Repetindo o último exemplo para a regra de três pontos teríamos que a integral seria aproximada por
1
∫
−1
√1 + (1 + 2x)2dx ≈ 2 √1 + (1 − 2
√3
)
2
+√1 + (1 + 2
√3
)
2
  ≈ 6.7746,
⎛
⎝
⎞
⎠
0.0258
f (x) [−1, 1]
1
∫
−1
f (x)dx
≈
1
9
[5f(−√ 3
5
)+ 8f (0) + 5f(√ 3
5
)].
[a, b]
b
∫
a
f (t)dt =
1
∫
−1
f ( a + b
2
+
b − a
2
x)
(b − a)
2
dx
≈
(b − a)
18
[5f(
a + b
2
−
b − a
2
√ 3
5
)+ 8f (
a + b
2
)+ 5f(
a + b
2
+
b − a
2
√ 3
5
)]
3
∫
−1
√1 + x2dx =
1
∫
−1
√1 + (1 + 2x)2dx
que fornece um erro de . Apesar da pequena diferença, temos que a quadratura Gaussiana de três
pontos fornece melhor resultado que a regra da quadratura para dois pontos.
Segue o código para a quadratura Gaussiana de dois pontos.
1 % funcao para calcular a numericamente a integral de f(x) no intervalo
2 % [a,b] usando a quadratura gaussiana de dois pontos
3
4 function I = quad_gaussiana(f, a, b, n) 
5
6 % entrada – f funcao f(x) a ser integrada
7 % a, b intervalo de integracao
8 % n numero de subdivisoes no intervalo
9 %
10 % saida – I valor aproximado da integral
11
12 dx = (b-a)/n;
13 x = a:dx:b;
14 I = 0;
15
16 for i = 1:n
17 I = I+((x(i+1)-x(i))/2)*(f((x(i+1)+x(i))/2-sqrt(3)*(x(i+1)-x(i))/6)+...
18 f((x(i+1)+x(i))/2+sqrt(3)*(x(i+1)-x(i))/6));
19 end
 
≈
2
9
5 1 +(1 − 2√ 3
5
)
2
+ 8√2 + 5 1 +(1 + 2√ 3
5
)
2⎛⎜⎝ ⎷ ⎷ ⎞⎟⎠≈ 6.8243,0.0239
SAIBA MAIS
Para mais exemplos sobre a quadratura Gaussiana acessem os links a seguir ou o acervo da disciplina. 
Quadratura de Gauss-Legendre 
Quadratura Gaussian
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d373793b4bc282c5d890ef9_texto07.pdf
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d3737a1c5aacbbc4f673991_texto08.pdf
Unidade 04
Aula 05
Solução de Problemas de Valor Inicial
Até este ponto, praticamente todas as equações diferenciais que lhes foram apresentadas nos cursos de cálculo
puderam ser resolvidas. O grande problema é que essas equações estão mais para exceção do que regra, pois a vasta
maioria das equações de primeira ordem não podem ter suas soluções analíticas encontradas.
Nestes casos, devemos recorrer a métodos numéricos que nos permitam aproximar soluções para essas
equações diferenciais que não possuem soluções analíticas.
Existem muitos métodos diferentesque podem ser usados para aproximar as soluções para uma equação
diferencial e nesta aula vamos olhar para um dos mais antigos e talvez o mais fácil de se usar.
Este método foi originalmente concebido por Leonhard Euler e é chamado, curiosamente, de método de
Euler.
Considerando o problema de valor inicial
\[\frac{{dy}}{{dx}} = f\left( {x,y} \right)\]
com \(y\left( {{x_0}} \right) = {y_0}.\) Começaremos querendo aproximar a solução deste PVI perto do ponto \
(x = {x_0}\). Iremos usar com duas informações que temos sobre a solução. Primeiro, sabemos o valor da
solução em \(x = {x_0}\) da condição inicial. Em segundo lugar, sabemos também o valor da derivada no ponto
\(x = {x_0}.\) Para calcular essa derivada, basta substituir a condição inicial em \(f\left( {x,y} \right)\), isto é, na
equação diferencial propriamente dita. Assim, a derivada neste ponto é
\[{\left. {\frac{{dy}}{{dx}}} \right|_{x0}} = f\left( {{x_0},{y_0}} \right).\]
Lembrem-se do curso de cálculo que essas informações são su�cientes para escrevemos a equação da reta
tangente a solução no ponto \(x = {x_0}\), que é dada por
\[y = {y_0} + f\left( {{x_0},{y_0}} \right)\left( {x - {x_0}} \right).\]
Ao olharmos para a �gura abaixo, se o ponto \({x_1}\) está su�cientemente próximo de \({x_0}\), então o
ponto \({y_1}\) sobre a reta tangente deverá estar próximo do valor real da solução no ponto \({x_1}\), isto é,
próximo de \(y\left( {{x_1}} \right)\). Determinar o ponto \({y_1}\) é relativamente fácil, tudo o que
precisamos é substituir \({x_1}\) na equação da reta tangente para a equação da reta dada acima e temos
\[{y_1} = {y_0} + f\left( {{x_0},{y_0}} \right)\left( {{x_1} - {x_0}} \right).\]
Agora, considerando que \({y_1}\) é uma boa aproximação para a solução e se construirmos uma reta que
passa pelo ponto \(\left( {{x_1},{y_1}} \right)\) e tem inclinação \(f\left( {{x_1},{y_1}} \right)\), então a equação
da reta será dada por
\[y = {y_1} + f\left( {{x_1},{y_1}} \right)\left( {x - {x_1}} \right).\]
Assim, para encontrarmos a aproximação da solução no ponto \(x = {x_2}\), esperamos que essa nova reta
estará su�cientemente próxima a solução real no ponto \({x_2}\) e usar esse valor da reta no ponto \({x_2}\)
como uma aproximação para a solução real. Isso nos dará
\[{y_2} = {y_1} + f\left( {{x_1},{y_1}} \right)\left( {{x_2} - {x_1}} \right).\]
Em geral, se tivermos \({x_n}\), uma aproximação da solução neste ponto, \({y_n},\) e quisermos encontrar
uma aproximação em \({x_{n + 1}}\) o que nós precisamos é repetir o procedimento descrito acima para obter
\[{y_{n + 1}} = {y_n} + f\left( {{x_n},{y_n}} \right)\left( {{x_{n + 1}} - {x_n}} \right).\]
Em muitas vezes, iremos supor que os tamanhos entre os pontos \({t_0}\), \({t_1}\), \({t_2}, \cdots \) são de
uma distância uniforme de tamanho \(h\). Em outras palavras, vamos assumir frequentemente que
\[{x_n} = {x_0} + n \cdot h.\]
Assim, o método de Euler é dado pela sequência
\[{y_{n + 1}} = {y_n} + hf\left( {{x_n},{y_n}} \right),\]
com um passo de tamanho \(h\), \({x_n} = {x_0} + n \cdot h\) com \(n\in \mathbb{N}\), \({x_0} = 0\) e \({y_0} =
y\left( 0 \right)\). Abaixo, veremos exemplos da aplicação do método de Euler.
Exemplo 1
Para este exemplo vamos considerar o problema de valor dado pela equação diferencial
\[\frac{{dy}}{{dx}} = - x{y^2}\]
sujeito a condição inicial \(y\left( 0 \right) = 1.\) A equação diferencial é separável e podemos determinar sem
grandes problemas a sua solução simplesmente fazendo
\[\frac{{dy}}{{{y^2}}} = - xdx \Rightarrow \mathop \int \nolimits^ \frac{{dy}}{{{y^2}}} = - \mathop \int
\nolimits^ xdx \Rightarrow - \frac{1}{y} = - \frac{{{x^2}}}{2} + C \Rightarrow y\left( x \right) = \frac{2}{{{x^2} -
2C}},\]
na qual a constante \(C\) é uma constante de integração. Agora, utilizando a condição inicial, \(y\left( 0 \right)
= 1\), podemos determinar a constante de integração que é \(C = - 1\). Finalmente, a solução do problema de
valor inicial é dada por
\[y\left( x \right) = \frac{2}{{{x^2} + 2}}.\]
Pelo método de Euler, a solução aproximada pode ser construída através da sequência
\[{y_{n + 1}} = {y_n} + hf\left( {{x_n},{y_n}} \right),\]
em que \(h\) é o passo, \({x_n} = {x_0} + n \cdot h\) com \(n\in \mathbb{N}\), \({x_0} = 0,\) \({y_0} = y\left( 0
\right)\) e \(f\left( {x,y} \right) = - x{y^2}.\) Ou seja, a sequência neste caso é dada por
\[{y_{n + 1}} = {y_n} - h{x_n}y_n^2.\]
Agora, escolhendo o passo \(h = 0.1\), podemos então gerar a nossa sequência de pontos. Assim, os três
primeiros pontos da sequência são dados por
\[{y_1} = {y_0} - h{x_0}y_0^2 = 1 - 0.1 \cdot 0 \cdot {1^2} = 1\]
\[{y_2} = {y_1} - h{x_1}y_1^2 = 1 - 0.1 \cdot \left( {0 + 0.1} \right) \cdot {1^2} = 0.99\]
\[{y_3} = {y_2} - h{x_2}y_2^2 = 0.99 - 0.1 \cdot \left( {0 + 2 \cdot 0.1} \right) \cdot {0.99^2} = 0.970398.\]
Continuando a sequência podemos montar a seguinte tabela.
\[{x_n}\] Solução aproximada (\({y_n})\) Solução analítica \(\left( {y\left( {{x_n}} \right)} \right)\) Erro
\[0.0\] \[1.000000\] \[1.000000\] \[0.000000\]
\[0.1\] \[1.000000\] \[0.995025\] \[0.004975\]
\[0.2\] \[0.990000\] \[0.980392\] \[0.009608\]
\[0.3\] \[0.970398\] \[0.956938\] \[0.013460\]
\[0.4\] \[0.942148\] \[0.925926\] \[0.016222\]
\[0.5\] \[0.906642\] \[0.888889\] \[0.017753\]
\[0.6\] \[0.865542\] \[0.847458\] \[0.018085\]
\[0.7\] \[0.820592\] \[0.803213\] \[0.017379\]
\[0.8\] \[0.773456\] \[0.757576\] \[0.015881\]
\[0.9\] \[0.725598\] \[0.711744\] \[0.013854\]
\[1.0\] \[0.678213\] \[0.666667\] \[0.011547\]
Abaixo podemos ver o grá�co da solução analítica do problema de valor inicial dada pela linha contínua em
azul, enquanto a solução aproximada é dada pelos círculos em vermelho.
Apesar do método de Euler ser descrito para uma equação de primeira ordem, ele pode ser facilmente
adaptado para determinar a solução aproximada de um problema de valor inicial envolvendo uma equação
de segunda ordem
\[y'' = f\left( {x,y,y'} \right)\]
com \(y\left( {{x_0}} \right) = {y_0}\) e \(y'\left( {{x_0}} \right) = {y_1}.\) Para utilizarmos o método de Euler é
necessário transformar a equação de segunda ordem em uma equação de primeira ordem. Fazemos isso
através da transformação \(y'\left( x \right) = z\left( x \right).\) Neste caso, a equação de segunda ordem é
transformada em um sistema de equações de primeira ordem o que nos torna possível utilizar o método de
Euler. Teremos então o seguinte sistema de equações de primeira ordem
\[y'\left( x \right) = z\]
\[z'\left( x \right) = f\left( {x,y,z} \right)\]
que pode ser escrito como
\[\frac{d}{{dx}}w = F\left( {x,w} \right),\]
com
\[w=\left( \begin{matrix} y \\ z \\ \end{matrix} \right)~~~\text{e}~~~~F\left( x,y,z \right)=\left( \begin{matrix}
z \\ f\left( x,y,z \right) \\ \end{matrix} \right).\]
Logo o método de Euler pode ser aplicado e as sequências de aproximação são dadas por
\[{w_{n + 1}} = {w_n} + hF\left( {{x_n},{w_n}} \right).\]
Abaixo, veremos um exemplo de como utilizar o método de Euler na solução de um problema de segunda
ordem.
Exemplo 2
Considere a equação
\[\frac{{{d^2}y}}{{d{x^2}}} + \frac{{dy}}{{dx}} + y = 0\]
sujeito as condições iniciais \(y\left( 0 \right) = - 1\) e \(y'\left( 0 \right) = 1\). Podemos determinar a solução
desta equação analiticamente sem di�culdades considerando que a solução deve ser da forma \(y\left( x
\right) = {e^{rx}}\). Substituindo essa última expressão na equação diferencial, chegamos na seguinte
equação característica
\[{r^2} + r + 1 = 0 \Rightarrow {{r}_{1}}=-\frac{1}{2}+i\frac{\sqrt{3}}{2}~~\text{e}~~{{r}_{2}}=-\frac{1}{2}-
i\frac{\sqrt{3}}{2}.\]
Considerando que as soluções da equação característica são complexas, podemos escrever a solução geral
como sendo
\[y\left( x \right) = {c_1}{e^{ - \frac{x}{2}}}\cos \left( {\frac{{\sqrt 3 }}{2}x} \right) + {c_2}{e^{ - \frac{x}{2}}}
{\rm{sen}}\left( {\frac{{\sqrt 3 }}{2}x} \right).\]
Finalmente, utilizando as condições iniciais podemos determinar as constantes \({c_1}\) e \({c_2}\)e concluir
que a solução do problema de valor inicial é dado por
\[y\left( x \right) = \frac{1}{3}\left[ {\sqrt 3 {e^{ - \frac{x}{2}}}{\rm{sen}}\left( {\frac{{\sqrt 3 }}{2}x} \right) - 3{e^{
- \frac{x}{2}}}\cos \left( {\frac{{\sqrt 3 }}{2}x} \right)} \right].\]
Fazendo \(\frac{{dy}}{{dx}} = z\left( x \right)\), isso nos dará a nova equação com duas variáveis dependentes \
(y~\)e \(z\)
\[\frac{{dz}}{{dx}} + z + y = 0.\]
Agora temos duas equações diferenciais de primeira ordem e duas variáveis dependentes a serem
determinadas, podemos então escrevê-las na forma de um sistema
\[\frac{d}{{dx}}\left( {\begin{array}{*{20}{c}} y\\ z \end{array}} \right) = \left( {\begin{array}{*{20}{c}} z\\ { - z - y}
\end{array}} \right),\]
com condição inicial
\[\left( {\begin{array}{*{20}{c}} {y\left( 0 \right)}\\ {z\left( 0 \right)} \end{array}} \right) = \left( {\begin{array}{*
{20}{c}} { - 1}\\ 1 \end{array}} \right).\]
Assim, se considerarmos que \(w = {\left( {y,z} \right)^T}\), então o método de Euler manterá a mesma forma
será dado pela seguinte sequência
\[{w_{n + 1}} = {w_n} + hF\left( {{x_n},{w_n}} \right),\]
em que \(h\) é o passo, \({x_n} = {x_0} + n \cdot h\) com \(n\in \mathbb{N}\), \({x_0} = 0,\) \({w_0} = {\left(
{y\left( 0 \right),z\left( 0 \right)} \right)^T}\) e \(F\left( {x,w} \right) = {\left( {z, - z - y} \right)^T}.\) Ou seja, a
sequência neste caso é dada vetorialmente por
\[\left( {\begin{array}{*{20}{c}} {{y_{n + 1}}}\\ {{z_{n + 1}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{y_n}}\\ {{z_n}} \end{array}} \right) + h\left( {\begin{array}{*{20}{c}} {{z_n}}\\ { - {z_n} - {y_n}} \end{array}} \right)
\Rightarrow \]
\[\left( {\begin{array}{*{20}{c}} {{y_{n + 1}}}\\ {{z_{n + 1}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{y_n} + h{z_n}}\\ {{z_n} - h{z_n} - h{y_n}} \end{array}} \right).\]
Escolhendo para este exemplo o passo \(h = 0.2\), podemos então gerar a nossa sequência de pontos. Assim,
os primeiros termos da sequência são dados por
\[\left( {\begin{array}{*{20}{c}} {{y_1}}\\ {{z_1}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{y_0} +
h{z_0}}\\ {{z_0} - h{z_0} - h{y_0}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 1 + 0.2 \cdot 1}\\ {1 - 0.2
\cdot 1 - 0.2 \cdot \left( { - 1} \right)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 0.8}\\ 1 \end{array}}
\right)\]
\[\left( {\begin{array}{*{20}{c}} {{y_2}}\\ {{z_2}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{y_1} +
h{z_1}}\\ {{z_1} - h{z_1} - h{y_1}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 0.8 + 0.2 \cdot 1}\\ {1 -
0.2 \cdot 1 - 0.2 \cdot \left( { - 0.8} \right)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 0.6}\\ {0.96}
\end{array}} \right)\]
\[\left( {\begin{array}{*{20}{c}} {{y_3}}\\ {{z_3}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{y_2} +
h{z_2}}\\ {{z_2} - h{z_2} - h{y_2}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 0.6 + 0.2 \cdot 0.96}\\
{0.96 - 0.2 \cdot 0.96 - 0.2 \cdot \left( { - 0.6} \right)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { -
0.408}\\ {0.888} \end{array}} \right)\]
Continuando a sequência podemos montar a seguinte tabela.
\[{x_n}\] Solução aproximada (\({y_n})\) Solução analítica \(\left( {y\left( {{x_n}} \right)} \right)\) Erro
\[0.0\] \[ - 1.000000\] \[ - 1.000000\] \[0.000000\]
\[0.2\] \[ - 0.800000\] \[ - 0.802167\] \[0.001267\]
\[0.4\] \[ - 0.600000\] \[ - 0.609605\] \[0.009605\]
\[0.6\] \[ - 0.408000\] \[ - 0.430659\] \[0.022659\]
\[0.8\] \[ - 0.230400\] \[ - 0.268589\] \[0.038189\]
\[1.0\] \[ - 0.072000\] \[ - 0.126193\] \[0.054193\]
Abaixo podemos ver o grá�co da solução analítica do problema de valor inicial dada pela linha contínua em
azul, enquanto a solução aproximada é dada pelos círculos em vermelho.
A seguir temos o código em MATLAB para determinar a solução aproximada de um problema de valor inicial
usando o método de Euler.
1 % calcula a solucao aproximada de um problema de valor inicial
2 % y’=f(x,y), y(x0)=y0
3 % usando o método de euler
4
5 % entrada – funcao f(x,y) a ser integrada
6 % a, b intervalo no qual sera calculada a solucao
7 % y0 condicao inicial
8 %
9 % saida – y solucao do pvi
10
11 clear all, close all, clc
12
13 f = @(x,y) -x.*y.^2; % funcão y’ = f(x,y)
14 h = 1e-1; % passo
15 a = 0;
16 b = 2;
17 xk = a:h:b; 
18 y(1) = 0; % condicao inicial
19
20 for i = 1:length(xk)-1
21 y(i+1) = y(i)+h*f(xk(i),y(i)); % calcula solucao aproximada
22 end
23
24 % imprime a solucao aproximada
25 �gure(1)
26 plot(xk,y);
SAIBA MAIS
Para uma explicação alternativa sobre o método de Euler acesse o vídeo: Clique aqui. 
Outros exemplos da utilização do método de Euler podem ser estudados no documento: 
Solução de Equações Diferenciais Ordinárias Usando Métodos Numéricos
https://www.youtube.com/watch?v=p-pyhXTw8F4
https://iesb.blackboard.com/bbcswebdav/institution/Ead/_disciplinas/EADG568/nova_novo/documents/5d373930c5aacb5f90673da0_texto09.pdf
Bem vindo(a) aluno(a)! Iniciamos aqui mais uma aula. Nela, veremos sobre os métodos de Runge-Kutta.
Esses métodos podem considerados como versões aperfeiçoadas do simplório método de Euler. Dado um
problema de valor inicial
\[y' = f\left( {x,y} \right)\]
com \(y\left( {{x_0}} \right) = {y_0}\), conforme vimos sobre o método de Euler, a aproximação do valor da
função \({y_{n + 1}}\) é calculado a partir do valor de \({y_n}\) e da derivada da função \(y\left( x \right)\) no
ponto \(\left( {{x_n},{y_n}} \right)\), isto é,
\[{y_{n + 1}} = {y_n} + hf\left( {{x_n},{y_n}} \right).\]
Nos métodos de Runge-Kutta, busca-se uma melhor estimativa da derivada da função \(y\left( x \right)\)
através de uma média das avaliações da função \(f\left( {x,y} \right)\) em um ou mais de um ponto dentro do
intervalo \(\left[ {{x_n},{x_{n + 1}}} \right]\).
Começaremos nesta aula discutindo o método de Runge-Kutta de segunda ordem, ou também conhecido
como método de Heun. Neste método, o valor da estimativa de \({y_{n + 1}}\) é encontrada a partir do valor
de \({y_n}\) e também de uma estimativa da derivada no ponto \({x_{n + 1}}.\) Assim, primeiramente se calcula
uma aproximação intermediária \({\tilde y_{n + 1}}\) dada por
\[{\tilde y_{n + 1}} = {y_n} + hf\left( {{x_n},{y_n}} \right)\]
e em seguida calculamos uma aproximação para \({y_{n + 1}}\) utilizando uma média aritmética entre as
inclinações nos pontos \(\left( {{x_n},{y_n}} \right)\) e \(\left( {{x_{n + 1}},{{\tilde y}_{n + 1}}} \right)\) na forma
\[{y_{n + 1}} = {y_n} + \frac{h}{2}\left[ {f\left( {{x_n},{y_n}} \right) + f\left( {{x_{n + 1}},{{\tilde y}_{n + 1}}} \right)}
\right].\]
Na �gura abaixo, podemos ver como a média aritmética das inclinações in�uencia a precisão da solução
aproximada, considerando que o ponto amarelo é a média das inclinações nos pontos \(\left( {{x_n},{y_n}}
\right)\) e \(\left( {{x_{n + 1}},{{\tilde y}_{n + 1}}} \right)\).
SAIBA MAIS
Os métodos de Runge-Kutta formam uma categoria de métodos numéricos para encontrar a solução
aproximada de problemas de valor inicial.
Exemplo 1
Considere o problema de valor inicial dado pela equação diferencial de segunda ordem
\[\frac{{{d^2}y}}{{d{x^2}}} + y = {\rm{sen}}\left( {\frac{{8x}}{{10}}} \right)\]
sujeita as condições iniciais \(y\left( 0 \right) = 0\) e \(y'\left( 0 \right) = 0.\) Assim como os demais exemplos
que resolvemos, para este exemplo também é possível determinar sua solução analítica. Lembre-se que a
solução de uma equação diferencial linear de segunda ordem não-homogênea é escrita como sendo a
combinação da solução homogênea, \({y_h}\left( x \right)\), e a solução particular, \({y_p}\left( x \right)\), isto
é,
\[y\left( x \right) = {y_h}\left( x \right) + {y_p}\left( x \right).\]
A solução homogênea é a solução da equação
\[\frac{{{d^2}{y_h}}}{{d{x^2}}} + {y_h} = 0,\]
enquanto a solução particular é a solução da equação
\[\frac{{{d^2}{y_p}}}{{d{x^2}}} + {y_p} ={\rm{sen}}\left( {\frac{{8x}}{{10}}} \right).\]
Neste caso, é possível mostrar sem grandes esforços que a solução homogênea é dada por
\[{y_h}\left( x \right) = {c_1}{\rm{sen}}\left( x \right) + {c_2}\cos \left( x \right),\]
e a solução particular pode ser encontrada usando o método dos coe�cientes a determinar usando a
expressão
\[{y_p}\left( x \right) = A{\rm{sen}}\left( {\frac{{8x}}{{10}}} \right) + B{\rm{cos}}\left( {\frac{{8x}}{{10}}} \right).\]
Finalmente, considerando as condições iniciais podemos mostrar então que a solução deste problema de
valor inicial é dada por
\[y\left( x \right) = \frac{{100}}{{36}}{\rm{sen}}\left( {\frac{{8x}}{{10}}} \right) - \frac{{20}}{9}{\rm{sen}}\left( x
\right).\]
Conforme �zemos, para podermos utilizar o método de Runge-Kutta de segunda ordem para resolver o
problema de valor inicial de segunda ordem é necessário transformar a equação em uma equação de
primeira ordem. Para isso consideramos a mudança \(\frac{{dy}}{{dx}} = z\left( x \right)\) o que nos dará a
nova equação com duas variáveis dependentes \(y\) e \(z\)
\[\frac{{dz}}{{dx}} + y = {\rm{sen}}\left( {\frac{{8x}}{{10}}} \right).\]
Neste momento, possuímos duas equações diferenciais de primeira ordem e duas variáveis dependentes a
serem determinadas. Escrevemos as equações na forma de um sistema
\[\frac{d}{{dx}}\left( {\begin{array}{*{20}{c}} y\\ z \end{array}} \right) = \left( {\begin{array}{*{20}{c}} z\\
{{\rm{sen}}\left( {\frac{{8x}}{{10}}} \right) - y} \end{array}} \right),\]
com condição inicial
\[\left( {\begin{array}{*{20}{c}} {y\left( 0 \right)}\\ {z\left( 0 \right)} \end{array}} \right) = \left( {\begin{array}{*
{20}{c}} 0\\ 0 \end{array}} \right).\]
Assim, considerando o vetor \(w = {\left( {y,z} \right)^T}\), o método de Runge-Kutta de segunda ordem irá
gerar as suas sequências de aproximações através da fórmula
\[{w_{n + 1}} = {w_n} + \frac{h}{2}\left[ {F\left( {{x_n},{w_n}} \right) + F\left( {{x_n} + h,{w_n} + hF\left( {{x_n},
{w_n}} \right)} \right)} \right],\]
em que \(h\) é o passo, \({x_n} = {x_0} + n \cdot h\) com \(n\in \mathbb{N}\), \({x_0} = 0\), \({w_0} = {\left(
{y\left( 0 \right),z\left( 0 \right)} \right)^T}\) e \(F\left( {x,w} \right) = {\left( {z,{\rm{sen}}\left( {\frac{{8x}}{{10}}}
\right) - y} \right)^T}.\) Vetorialmente, podemos escrever a sequência como
\[\left( {\begin{array}{*{20}{c}} {{y_{n + 1}}}\\ {{z_{n + 1}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{y_n}}\\ {{z_n}} \end{array}} \right) + \frac{h}{2}\left[ {\left( {\begin{array}{*{20}{c}} {{z_n}}\\ {{\rm{sen}}\left(
{\frac{{8{x_n}}}{{10}}} \right) - {y_n}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{{\tilde z}_{n + 1}}}\\
{{\rm{sen}}\left( {\frac{{8{x_{n + 1}}}}{{10}}} \right) - {{\tilde y}_{n + 1}}} \end{array}} \right)} \right]\]
com
\[\left( {\begin{array}{*{20}{c}} {{{\tilde y}_{n + 1}}}\\ {{{\tilde z}_{n + 1}}} \end{array}} \right) = \left(
{\begin{array}{*{20}{c}} {{y_n} + h{z_n}}\\ {{z_n} + h\left( {{\rm{sen}}\left( {\frac{{8{x_n}}}{{10}}} \right) - {y_n}}
\right)} \end{array}} \right).\]
Se considerarmos o passo \(h = 0.1\), podemos então gerar a nossa sequência de pontos. Observe que o
método de Runge-Kutta de segunda ordem necessita calcular duas aproximações. Primeiramente,
calculamos
\[\left( {\begin{array}{*{20}{c}} {{{\tilde y}_1}}\\ {{{\tilde z}_1}} \end{array}} \right) = \left( {\begin{array}{*{20}
{c}} {{y_0} + h{z_0}}\\ {{z_0} + h\left( {{\rm{sen}}\left( {\frac{{8{x_0}}}{{10}}} \right) - {y_0}} \right)} \end{array}}
\right) = \left( {\begin{array}{*{20}{c}} {0 + 0.1 \cdot 0}\\ {0 + 0.1 \cdot \left( {{\rm{sen}}\left( {\frac{{8 \cdot 0}}
{{10}}} \right) - 0} \right)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0\\ 0 \end{array}} \right).\]
Em seguida, determinamos a primeira aproximação que é dada por
\[\left( {\begin{array}{*{20}{c}} {{y_1}}\\ {{z_1}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{y_0}}\\
{{z_0}} \end{array}} \right) + \frac{h}{2}\left[ {\left( {\begin{array}{*{20}{c}} {{z_0}}\\ {{\rm{sen}}\left(
{\frac{{8{x_0}}}{{10}}} \right) - {y_0}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{{\tilde z}_1}}\\
{{\rm{sen}}\left( {\frac{{8{x_1}}}{{10}}} \right) - {{\tilde y}_1}} \end{array}} \right)} \right]\]
\[ = \left( {\begin{array}{*{20}{c}} 0\\ 0 \end{array}} \right) + \frac{{0.1}}{2}\left[ {\left( {\begin{array}{*{20}{c}}
0\\ {{\rm{sen}}\left( {\frac{{8 \cdot 0}}{{10}}} \right) - 0} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0\\
{{\rm{sen}}\left( {\frac{{8 \cdot 0.1}}{{10}}} \right) - 0} \end{array}} \right)} \right]\]
\[ = \left( {\begin{array}{*{20}{c}} {0.000000}\\ {0.003996} \end{array}} \right).\]
Repetindo esse procedimento sucessivamente, podemos tabelarmos alguns pontos da solução e
compararmos o erro, como podemos ver abaixo.
\[{x_n}\] Solução aproximada (\({y_n})\) Solução analítica \(\left( {y\left( {{x_n}} \right)} \right)\) Erro
\[0.0\] \[0.000000\] \[0.000000\] \[0.000000\]
\[0.1\] \[0.000000\] \[0.000134\] \[0.000134\]
\[0.2\] \[0.000799\] \[0.001063\] \[0.000265\]
\[0.3\] \[0.003185\] \[0.003575\] \[0.000389\]
\[0.4\] \[0.007921\] \[0.008423\] \[0.000503\]
\[0.5\] \[0.015729\] \[0.016329\] \[0.000601\]
\[0.6\] \[0.027281\] \[0.027961\] \[0.000680\]
\[0.7\] \[0.043188\] \[0.043925\] \[0.000737\]
\[0.8\] \[0.063985\] \[0.064755\] \[0.000769\]
\[0.9\] \[0.090127\] \[0.090901\] \[0.000774\]
\[1.0\] \[0.121975\] \[0.122724\] \[0.000749\]
Abaixo podemos ver o grá�co da solução analítica do problema de valor inicial dada pela linha contínua em
azul, enquanto a solução aproximada é dada pelos círculos em vermelho. Perceba que o método de Runge-
Kutta de segunda ordem fornece efetivamente uma boa aproximação para a solução do problema de valor
inicial.
Pudemos ver no exemplo anterior que o método de Runge-Kutta de segunda ordem fornece uma ótima
aproximação para a solução do problema de valor inicial. Porém, é conhecido como o método clássico de
Runge-Kutta é o método de quarta ordem. A sua construção é semelhante ao método de segunda ordem, de
forma que a sequência de aproximações é construída com a ajuda de uma média de inclinações. É possível
mostrar que a sequência de inclinações para o método de Runge-Kutta de quarta ordem é dada por
\[{y_{n + 1}} = {y_n} + \frac{h}{6}\left( {{k_1} + 2{k_2} + 2{k_3} + {k_4}} \right),\]
de forma que as inclinações intermediárias são dadas por
\[{k_1} = f\left( {{x_n},{y_n}} \right),\]
\[{k_2} = f\left( {{x_n} + \frac{h}{2},{y_n} + \frac{h}{2}{k_1}} \right),\]
\[{k_3} = f\left( {{x_n} + \frac{h}{2},{y_n} + \frac{h}{2}{k_2}} \right),\]
\[{k_4} = f\left( {{x_n} + h,{y_n} + h{k_3}} \right).\]
Este é o método de Runge-Kutta mais utilizado, inclusive o nome do método de Runge-Kutta é confundido e
difundido na literatura como sendo o método de ordem quatro. Veremos no exemplo abaixo o poder deste
método e como o seu erro de ordem quatro chegam a aproximar quase que perfeitamente a solução.
Exemplo 2
Para exempli�carmos o método de Runge-Kutta de quarta ordem, vamos revisitar um problema de valor
inicial que utilizamos na aula sobre o método de Euler. Assim, considere o problema de valor inicial dado pela
equação diferencial
\[\frac{{dy}}{{dx}} = - x{y^2}\]
sujeito a condição inicial \(y\left( 0 \right) = 1\). Como vimos, essa equação diferencial é separável e sua
solução geral é dada por
\[y\left( x \right) = \frac{2}{{{x^2} - 2C}},\]
na qual a constante \(C\) é uma constante de integração. Utilizando a condição inicial, \(y\left( 0 \right) = 1\),
determina-se a constante de integração e mostra-se que a solução do problema de valor inicial é dada por
\[y\left( x \right) = \frac{2}{{{x^2} + 2}}.\]
Conforme vimos acima, o método de Runge-Kutta de quarta ordem fornece a solução aproximada através da
sequência
\[{y_{n + 1}} = {y_n} + \frac{h}{6}\left( {{k_1} + {k_2} + {k_3} + {k_4}} \right),\]
em que
\[{k_1} = f\left( {{x_n},{y_n}} \right),\]
\[{k_2} = f\left( {{x_n} + \frac{h}{2},{y_n}