Prévia do material em texto
PONTIFÍCIA UNIVERSIDADE CATÓLICA DO RIO GRANDE DO SUL
FACULDADE DE MATEMÁTICA
CÁLCULO NUMÉRICO
Notas de Aula , Aplicações e Exercícios com Sistema Maple
Eliete Biasotto Hauser
1 –Teoria dos Erros - Sistema de Ponto Flutuante
Um número y na base 2≥β , y = LL 2b1b,0a1a2a3a1nana − pode ser descrito
na forma
y = LL +−+−++++++−
−
+ 22b
1
1b0a1a
2
2a
3
3a
1n
1na
n
na βββββββ .
Por exemplo, 2106110271012105310326,3517 −×+−×++×+×+×=
Em aritmética de ponto flutuante normalizado de t dígitos, y tem a forma:
y = etd2d1d,0 ββ
×
L
i) td2d1d,0 L é a mantissa (uma fração na base β),
ii) 01d,1jd0 ≠−≤≤ β , j=1,2,...,t
iii) “e” é um expoente inteiro que varia no intervalo [m,M].
M e m dependem da máquina utilizada. Em geral, m = -M∈Z.
iv) t define a precisão da máquina, é o número de dígitos da mantissa.
Obs: precisão ≠ exatidão(depende da precisão da máquina e do algorítmo utilizado)
A união de todos os números de ponto flutuante normalizados com o zero:
m
vezest
0000.00 β×= 43421 L
é chamado Sistema de Ponto Flutuante Normalizado e representado por F(β, t, m, M).
Alguns exemplos de máquinas com precisão simples:
a) HP 48 : F(10, 12, -498, 500)
b) IBM 3090 : F(16, 6, -64, 63)
c) Cray1 : F(2, 48, -8192, 8191)
d) Burroughs B6700: F(8, 13, -63, 64)
Em F valem as propriedades:
1) p=0,1 x mβ é o menor número em módulo, não nulo, de F.
2) G=0,
444 3444 21
vezest
)1)...(1)(1( −−− βββ x Mβ é o maior número, em módulo, não nulo, de F.
3) # F = 1)1(1)1(2 +
+−−− mMtββ é o cardinal de F
4) Para qualquer mantissa, vale 1mantissa1 <≤−β .
5) Se Fy ∈ , então Fy ∈− .
6) F1eF0 ∈∈ .
Se o expoente da base não pertencer a [m,M], y não pode ser representado em F.
São os casos de erro de:
- underflow, se e < m (ultrapassa a capacidade mínima)
- overflow, se e > M (ultrapassa a capacidade máxima)
1
E.B.Hauser – Cálculo Numérico
Exemplo: O sistema de ponto flutuante normalizado F(2, 2, -2, 1) é um conjunto discreto
com 17 elementos ( #F=17), β=2, t=2, m=-2, M=1 .
Para obter a representação real dos elementos de F é necessário utilizar desnormalização
seguida do processo de mudança de base 2=β para base decimal.
( ) ( ) 125,03212201200202001,022210,0 =−×+−×+−×+×===−×= ββ
( ) ( ) 25,0221120020201,012210,0 =−×+−×+×===−×= ββ
( ) ( ) 5,012102021,002210,0 =−×+×===×= ββ
( ) ( ) 10212112210,0 =×===×= ββ
( ) ( ) 1875,042132122012002020011,022211,0 =−×+−×+−×+−×+×===−×= ββ
( ) ( ) 375,03212211200202011,012211,0 =−×+−×+−×+×===−×= ββ
( ) ( ) 75,0221122020211,002211,0 =−×+−×+×===×= ββ
( ) ( ) 5,112102121,112211,0 =−×+×===×= ββ
Como p=0,125 e G= 1,5, temos:
underflow: (-0,125 ; 0,125)-{0}
overflow: (-∞ ; -1,5) U (1,5 ; +∞)
Observamos que:
a) 0,75∈F e 1,5∈F mas 0,75+1,5=2,25∉F, isto é, 2,25 está na região de overflow de F;
b) 0,125∈F e 0,25∈F e 0,125+0,25=0,375∈F;
c) 0,25∈F e 1∈F e 0,25+1=1,25∉F mas não está na regiões de overflow e underflow
de F. Neste caso, 1,25 será arredondado para 1 ou 1,5 dependendo to tipo de arredon-
damento definido para F.
2
1 - Teoria dos Erros
Arredondamentos
Se a representação do real y em F não é exata, é necessário utilizar um arredonda-
mento. Os tipos de arredondamento mais conhecidos são:
- Arredondamento para número mais próximo de máquina (Oy).
- Arredondamento por falta, ou truncamento (∇ y).
- Arredondamento por excesso ( ∆ y)
Exemplo: Arredondar y em F(10, 4, -5, 6), se possível. β=10, t=4, m=-5, M=6.
Y Oy ∇ y ∆ y
Π 3,142 3,141 3,142
2000/3 666,7 666,6 666,7
1/3 0,3333 0,3333 0,3334
215=4084101=0, 4084101 x 107∈overflow pois 7>6=M
21-5=1/4084101=0, 2448519270 x 10-6∈underflow pois -6<-5=m
Em F, geralmente, as operações de adição e multiplicação não são comutativas, as-
sociativas e nem distributivas, pois numa série de operações aritméticas, o arredondamento é
feito após cada operação. Ou seja, nem sempre as operações aritméticas válidas para os nú-
meros reais são válidas em F. Isto influi na solução obtida através de um método numérico.
Assim, métodos numéricos matematicamente equivalentes, podem fornecer resultados dife-
rentes.
Exemplo: 7,52x52,142x)x(p ++= tem raízes -7,3471 e -7,1728.
Arredondando 14,52 para 14,5, 7,52x5,142x)x(p ++= tem raízes complexas: i3708,025,7 ±−
Medidas de Exatidão
Quando se aproxima um número real x por x*, o erro que resulta é x-x*. Define-se:
erro absoluto: EA = *xx − e o erro relativo: ER =
x
xx *−
para 0x ≠ .
Exemplo:
Comprimento
real (cm)
Comprimento
medido (cm)
EA ER
Ponte
10000 9999 1 0,0001=0,01%
Rebite
10 9 1 0,1=10%
Na resolução de um problema o valor exato da solução x pode ser desconhecido.
Podemos usar duas aproximações sucessivas de x, definindo:
+
−+++−=+
1
1log3,0)1,(
kx
kxkx
kxkxDIGSE µ
3
E.B.Hauser – Cálculo Numérico
o qual expressa o número de dígitos significativos exatos de kx em relação a 1kx + . Aqui
µ representa a unidade de arredondamento da máquina ( µ = t1
2
1
−β se o arredondamen-
to for Ox ).
No caso de duas aproximações de π, 587,3)1415926535,3;142,3(DIGSE ≅
Algoritmos Numéricos
O Cálculo Numérico tem por objetivo estudar e aplicar algoritmos numéricos para a
solução de problemas, visando o menor "custo" e confiabilidade do resultado (observar
tempo de execução, número de operações aritméticas, quantidade de memória, propagação
do erro, etc).
Um bom algoritmo numérico deve se caracterizado por:
a) Independência da máquina :a execução do programa pode ser realizada em diferentes
máquinas.
b) Inexistência de Erro Aritmético: problemas de overflow e underflow devem ser detecta-
dos a priori
c) Inexistência de Erro Lógico: previsão de tudo o que poderá ocorrer durante o processo
(divisão por zero, por exemplo)
d) Quantidade finita de cálculos.
e) Existência de um critério de exatidão.
f) O erro deve convergir para zero com precisão infinita.
g) Eficiência: produz a resposta correta com economia
Passos para a resolução de um problema: tipos de erros
A resolução de problemas reais envolve diversas etapas:
Estudaremos algoritmos numéricos a fim de obter uma solução numérica do proble-
ma, a qual, geralmente, difere da solução exata. Possíveis fontes de erro que geram essa
diferença são:
a) Simplificação do modelo matemático ( algumas variáveis envolvidas são desprezadas);
b) Erro nos dados de entrada ( com frequência provindos de dados experimentais);
c) Truncamento (por exemplo, substituição de um processo infinito por um finito e lineari-
zações);
d) Erro de arredondamento em aritmética de ponto flutuante.
Curiosidades: Some disasters caused by numerical errors.
http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html
Problema
Real
Modelagem
Matemática Solução
4
1 - Teoria dos Erros
Exercícios
1) No sistema Maple estimar 3.8e− utilizando
∑
∞
=
−
=
−
0i !i
i)x(xe e
∑
∞
=
==
−
0i !i
ix
1
xe
1xe
com 26 termos cada e comparar com 3.8e− ≈ 0.2485168271x10-3.
2) Em F(10 , 3 , -98 , 98) e arredondamento por truncamento estimar p(2.73) se:
a) 55.0x62x53x)x(p ++−=
b) 55.0x)6x)5x(()x(p ++−=
Em ambos os casos estimaro erro absoluto ao comparar com p(2.73) ≈ 0.11917x10-1.
Obs: Estimar p(x) pelo algoritmo usual
01
2
2
3
3...
1
1)( axaxaxa
nx
n
anxnaxp +++++
−
−
+=
exige n adições e (n2+n)/2 multiplicações enquanto que o algoritmo de Horner
{ 0)1)2...)2)1
1
(((...()( axaxaxnaxnaxna
n
xp +++
−
+
−
+
−
=
requer n adições e n multiplicações.
3) Sejam A, B, C e D matrizes genéricas de ordem 10x20, 20x50, 50x1 e 1x100 respecti-
vamente. Utilizando a propriedade associativa, pode-se determinar o produto matricial
AxBxCxD de diversas formas. Qual das duas abaixo é mais eficiente? Porque?
a) Ax(Bx(CxD)) b) (Ax(BxC))xD
4) Representar o número real x na base 2 usando 8 algarismos significativos? Essa represen-
tação é exata? a) x=0.6 b) x=13.25 c) x= 2.47
5) Determinar o cardinal , regiões de underflow e overflow e todos elementos reais de:
a) F(2,3,-1,2) b) F(3,2,-1,2) c) F(2,2,-2,2)
6) Representar, se possível, os números abaixo em utilizando arredondamento por trunca-
mento( x∇ )e arredondamento para número mais próximo de máquina (Ox) em F(10,5,-2,2).
a) 3 b)
3
200
c)
3
2000
d) e e)
3000
1
f) 2
7) Considerando: ∑
=
=
10
12
1
i
i
A
a) Calcular o valor de A utilizando precisão infinita..
b) Utilizando arredondamento por truncamento ( x∇ ) em F(10,3,-98,98), estimar o valor de
A somando da direita para esquerda e após somando da esquerda para a direita. Comparar
os resultados.
5
E.B.Hauser – Cálculo Numérico
Respostas:
1) exp(-8.3);
.0002485168271
> f1:=sum(((-x)^i)/i!, i=0..25):
> f1a:=unapply(f1,x):
> f1a(8.3);
-.001241405905
Obs: Causas desse erro: subtração de grandezas muito próximas e adição de grandezas de
diferentes ordens.
>
> f2:=1/(sum(((x)^i)/i!, i=0..25)):
> f2a:=unapply(f2,x):
> f2a(8.3);
.0002485170000
2) a)p(2.73)= -0.05 ,erro absoluto = 0.061917
b)p(2.73)=0.032 ,erro absoluto = 0.020083
3) (Ax(BxC))xD é mais eficiente pois exige 2200 multiplicações enquanto que para calcular
o produto Ax(Bx(CxD)) são necessárias 125000 multiplicações .
OBS: Se M é de ordem pxq e N de ordem qxr, então MxN, de ordem pxr, é obtida efetuando
pqr operações de multiplicações de elementos de M e N.
4-a) 2)10011001,0(10)6,0( ≈ b) 2)01,1101(10)25,13( ≈ c) 2)011110,10(10)47.2( ≈
5)-a) 0, 1/4, 1/2, 1, 2, 5/16, 5/8, 5/4, 5/2, 3/8, 3/4, 3/2, 3, 7/16, 7/8, 7/4, 7/2 e simétricos.
# F = 33 Região de oferflow: ),2/7()2/7,( +∞∪−−∞
Região de underflow: (-1/4,1/4) - {0}
b) 0, 1/9, 1/3, 1,3, 4/27, 4/9,4/3, 4, 5/27, 5/9, 5/3, 5, 2/9, 2/3, 2, 6, 7/27, 7/9, 7/3, 7, 8/27,
8/9, 8/3, 8 e seus simétricos. # F = 49 Região de oferflow: ),8()8,( +∞∪−−∞
Região de underflow: (-1/9,1/9) - {0}
c) 0, 1/8, 1/4, 1/2, 1, 2, 3/16, 3/8, 3/4, 3/2, 3, e seus simétricos. # F = 21
Região de oferflow: ),3()3,( +∞∪−∞
Região de underflow: (-1/8,1/8) - {0}
6-a) 0,17320*101 e 0,17321 *101 b) 0,66666*102 e 0,66667*102
c) overflow (0,66666*103 ∉ F e 0,66667*103 ∉ F) d) 0.27182*101 e 0.27183*101
e) underflow (0.33333*10 –3 ∉ F) f) 0,14142*101 e 0,14142 *101
7-a) 0,9990234375 b) 0,999 e 0,998
6
2. Resolução de Equações Algébricas e Transcendentes
Objetivo: Resolver equações de forma ( ) 0=xf , isto é, determinar os valores de x para os quais a
igualdade ( ) 0=xf é satisfeita.
Se a função ( )xf só contém operações algébricas repetidas um número finito de vezes, a
equação é dita algébrica.
Exemplos de Equações Algébricas: a) 01x43,43x =−− b) 01xx2 =−−
Exemplos de Equações Transcendentes:
a) 0)x3(senex 5
x
=
−
b) 0)x9,0(sen)x1,0(sen2 =
7
E.B.Hauser – Cálculo Numérico
2.1- Equações Polinomiais
Revisão sobre polinômio:
Seja ( ) 012211 axaxaxaxaxp nnnn +++++= −− K um polinômio de grau n, onde os coefici-
entes , são números reais e 0≠na .
Temos que:
1) ( )xp possui n raízes, contadas as multiplicidades.
2) Se nααα ,,, 21 K forem raízes reais de ( )xp , então ( )xp pode ser fatorado na forma:
( ) ( )( ) ( )nn xxxaxp ααα −−−= K21 .
Ex.: ( ) 2xx2xxp 23 −−+= , tem raízes:
2 ,1 ,1 321 −=−== ααα . Logo,
( ) ( )( )( )2x1x1xxp ++−=
3) Se z= bia + é raiz de ( )xp , então biaz −= também o é.
Ex.: ( ) 10x6xxp 2 +−= tem raízes i±3 .
4) Se bia + é raiz de ( )xp de grau 2≥n , então ( )xp pode ser fatorado:
( ) ( ) ( )xqbaax2xxp 222 ⋅++−= ,onde o grau de ( )xq é 2−n .
Ex.: a) ( ) 2x2xx2xxp 234 −++−= tem raízes i±± 1 ,1 .
( ) ( ) )1x(1x2x22xxp +−⋅
+−=
Ex.: b) ( ) 10167 23 −+−= xxxxp tem raízes i±3 ,1 .
( ) ( ) ( )1x10x6xxp 2 −⋅+−=
8
2 - Resolução de Equações Algébricas e Transcendentes
5)Se ( )xp é de grau ímpar, então ( )xp possui ao menos uma raiz real.
6) Uma raiz α de ( )xp tem multiplicidade m se ( ) ( ) ( ) ( ) 01"' ===== − αααα mpppp K
e ( ) 0≠αmp
Ex.:1) 2=α é raiz de multiplicidade 2 de ( ) 8x6x2xp 23 +−=
2) 2=α é raiz de multiplicidade 3 de ( ) 8465 234 −++−= xxxxxp
9
E.B.Hauser – Cálculo Numérico
7) Valor numérico de um polinômio: para calcular, de forma usual, ( )ixp são necessárias n
adições e ( )
2
1+nn
multiplicações.
O Método de Horner faz esse cálculo com n adições e n multiplicações:
( ) 0121
parênteses 1
)))(((( axaxaxaxaxp nn
n
+++++=
−
−
K
321
KK
Ex.: ( ) 4x)3x)1x)2x)4x3((((4x3xx2x4x3xp 2345 +−−+=−+−−+=
8) Se ( )xp é de grau n , então existe único polinômio de grau ( )xqn ,1− , tal que
( ) ( ) ( ) ( )αα pxqxxp +⋅−= .
Se α é raiz de ( )xp então ( ) 0=αp .
É o algoritmo de Briot-Ruffini utilizado para Deflacionar Raízes.
Ex.
( )
( )( ) ( )
( )( ) ( )
( )( ) ( ) 1483pe14846x10x3x3
22pe26x5x2x2
01pe10x6x1x1
10x16x7xxp
2
3
2
2
2
1
23
−=−−+−+⇒−=
=++−−⇒=
=+−−⇒=
−+−=
α
α
α
2.1.1-Enumeração das Raízes
Enumerar as raízes de ( )xp é dizer quantas são as raízes e se positivas, negativas ou com-
plexas.
Regra de Descartes ou Regras de Sinais
“O número de raízes reais positivas de ( )xp é igual ao número de variações de sinais na
seqüência dos coeficientes ou menor do que este por um número inteiro par, sendo uma raiz de mul-
tiplicidade “m” contada como “m” raízes e não sendo considerados os coeficientes nulos”.
O número de raízes reais negativas é obtido aplicando a regra de Descartes a ( )xp −
Regra de Huat
Se ( ) 00p ≠ e para algum k, 112 +− ×≤ kkk aaa então ( )xp possui raízes complexas.
Regra da Lacuna
Se ( ) 00 ≠p e para algum K, 0 e 0 11 >×= +− kkk aaa , então ( )xp tem raízes complexas.
10
2 - Resolução de Equações Algébricas e Transcendentes
Ex.: ( ) 153 345 −−−+= xxxxxp
Por Descartes, ocorrem as possibilidades.• 1 raiz real positiva, 2 raízes negativas e 2 comple-
xas;
ou
• 1 raiz real positiva, nenhuma negativa, e 4 comple-
xas.
p(x) tem raizes complexas, pois, por Huat: 3a1a
2
2a ×≤ ,
istoé, )5()1(20 −×−< . A Regra da Lacuna também aplica:
03a1a e 02a >×= .
2.1.2-Localização das raízes de p(x)
O objetivo é determinar a região do plano que contém todas as raízes de ( )xp .
Cota de Cauchy: Toda raiz α (real ou complexa) de ( )xp satisfaz βα ≤ , onde
0x,xlim 0k
k
==
∞→
β e n
n
0
k
n
12n
k
n
2n1n
k
n
1n
1k
a
a
x
a
a
x
a
a
x
a
a
x ++++= −−−−+ K .
Ex.: ( ) 1x3136,068,1x37,3x3xxp 234 −+−+−=
4
k
2
k
3
k1k 3136,0x68,1x37,3x3x +++=+
0x0 =
M
M
9575796715,3x
9552764745,3x
59519059204,3x
09469665820,3x
69397307712,3x
40048339019,3x
36165556530,2x
61069395791,2x
64735839615,1x
557483314773,0x
19
18
17
16
15
5
4
3
2
1
=
=
=
=
=
=
=
=
=
=
11
E.B.Hauser – Cálculo Numérico
Então, a cota de Cauchy para p(x) é 4≅β e então todas as raízes α tem 4≤α .
a) Métodos Gráficos
a1) esboçar gráfico da função ( )xf e localizar as abcissas dos pontos onde a curva intercepta o
eixo dos x.
a2) de ( ) 0=xf obter uma equação equivalente ( ) ( )xfxf 21 = . Localizar no mesmo eixo cartesiano
os pontos r onde as duas curvas se interceptem: ( ) ( ) ( ) ( ) ( ) 002121 =⇒=−⇒= rfrfrfrfrf
b) Método Analítico(Teorema de Bolzano)
Seja ( )xf continua no intervalo [ ]ba, . Se ( ) ( ) 0<⋅ bfaf , então existe pelo menos uma raiz de f
em ( )b,a . (Se o sinal de 'f é constante em ( )b,a a raiz é única nesse intervalo).
Ex.: ( ) 393 +−= xxxp
a) Análise gráfica:
2.2- Separação de Raízes Reais de f(x)=0
Logo, existem três raízes reais: ( ) ( ) ( )3,23r,1,02r3,41r ∈∈−−∈
b) analiticamente:
( ) 373923,753113923,13325
323101334
−−−−
−−−−
xp
x
12
2 - Resolução de Equações Algébricas e Transcendentes
Obs: Devemos dar uma atenção especial para os casos de:
� Raízes muito próximas.
� para raízes de multiplicidade par não ocorre troca
de sinal.
Ex: ( ) 3136,068,137,33 234 +−+−= xxxxxp
8,0 e 7,0 21 ≅= rr são raízes de multiplicidade 2.
2.3- Métodos para Resolução de equações algébricas e transcendentes
Qualquer método deve observar um critério de parada, ao qual está associado um estimador
de exatidão. Por exemplo, para onde LC ,,, 21 εε são dados:
• C)x,x(DIGSE 1kk ≥+
• ( ) 1kxf ε≤
• 2
1
1 ε≤
−
+
+
k
kk
x
xx
• Lk > (número máximo de iterações)
2.3.1-Método da Bisseção ou Dicotomia (Algoritmo de quebra)
Seja [ ] ℜ→baf ,: continua e tal que ( ) ( ) 0<⋅ bfaf .
1) Calcula-se o ponto médio
2
ba
xm
+
= , dividindo-se [ ]ba, em dois novos intervalos :
[ ]mxa, , [ ]bxm ,
2) Se ( ) 0≠mxf e:
i) ( ) ( ) 0<⋅ mxfaf então a raiz está em ( )mx,a . Volta-se para (1)
ii) ( ) ( ) 0<⋅ bfxf m então a raiz está em ( )b,mx . Volta-se para (1)
3) Repete-se o processo até obter uma aproximação “razoável” da raiz, isto é, até que
um critério de parada seja satisfeito.
Características: É simples a convergência lenta mas garantida. A velocidade de convergência
é DIGSE⋅3,0 /passo, isto é, a cada 3 ou 4 passos ganha-se um DIGSE .
Ex: ( ) 11205,72 234 −−−+= xxxxxp
a) Enumeração das raízes de ( )xp
Regras de Huat e Lacuna
não aplicam
13
E.B.Hauser – Cálculo Numérico
1r 2r 3r 4r
4211
4031
total⊂ℜℜ
−+
b) Cota de Cauchy:
0,11205,72 04 231 =+++=+ xxxx kkkk
46483007996,4
46481382609,481469532045,4
96478405782,467425618494,3
16472953921,460308011161,3
48211602868,1
17
164
153
142
1
=
==
==
==
=
x
xx
xx
xx
x
M
M
c) Separação de raízes
( ) 5,57617335495,3511515,8775,276
54321012345
−−−−−
−−−−−
xp
x
d) Cálculo da raiz ( )4,34r ∈ utilizando o método da bissecção.
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
M002874148,00352783203,313
20354003906,3;03515625,320152911151,020354003906,312
50356445312,3;03515625,30401045242,050356445312,311
0361328125,3;03515625,30897548754,00361328125,310
037109375,3;03515625,31891500395,0037109375,39
0390625,3;03515625,30095143543,003515625,38
039625,3;03125,33883184832,00390625,37
046875,3;03125,31900454168,1046875,36
0625,3;03125,3405335188,003125,35
0625,3;3871886353,20625,34
125,3;3660400391,9125,33
25,3;300390625,2525,32
35;39375,625,31
kIkxfkxk
−
−
Obtemos 0354039062,312x4r =≈ com ( ) 096,413x,12xDIGSE ≅
2.3.2-Método da Iteração Linear
Consiste em escrever a equação ( ) 0=xf na forma ( )xGx = . Os pontos x* que satis-
fazem a condição ( )** xGx = são ditos pontos fixos de ( )xG e geometricamente repre-
sentam os pontos de intersecção da reta xy = com a curva ( )xGy = .
14
2 - Resolução de Equações Algébricas e Transcendentes
( )xG é dita função de iteração do método. Inicia-se a iteração com um valor 0x
próximo da raiz, e as outras aproximações são dadas
por:
( ) K,2,1,0,1 ==+ ixGx ii
A seqüência de aproximação ix , converge para a
solução x* da equação ( ) 0=xf sob certas condições
.
A construção de G não é única. A escolha de
uma G apropriada é dita “problema do ponto fixo.
Ex. 062 =−+ xx .
( ) 16,
1
6
,6),6)6) 543221 −=+=−−=−=−= xGxGxGcxGbxxGa
Embora não seja preciso usar métodos numéricos para encontrar as duas raízes reais
2,3 21 =−= αα e da presente equação, nota-se que:
i) Tomando 1G e 5,10 =x , a seqüência }{ ix não converge para 2.
( )ii xGx 11 =+
( )
( )
( )
( )
( )
M
8342,12078822
46095276,3475
00390625,59)0625,8(6
0625,8)75,3(6
75,3)5,1(6
415
314
2
213
2
112
2
011
−==
−==
−=−−==
−=−==
=−==
xGx
xGx
xGx
xGx
xGx
ii) Tomando 2G e 5,10 =x , a seqüência }{ ix converge para 2.
( )
( )
( )
( )
( )
( )
( )
M
10000298018,26x2G7x
69998807918,15x2G6x
50004768183,24x2G5x
39980924992,13x2G4x
40076263645,22x2G3x
9694363804,11x2G2x
61213203435,25,160x2G1x
==
==
==
==
==
==
=−==
Teorema da Convergência:
Seja α uma raiz isolada de f em [ ]ba, . Se
i) G e G’são contínuas em [ ]ba, ;
15
E.B.Hauser – Cálculo Numérico
ii) ( )b,ax,1)x('G ∈∀〈 ;
iii) ( ) ,...2,1,0k,b,a)kx(G1kxe0x =∈=+∈ Ι , então a seqüência { kx }, gerada por
)(1 kxGkx =+ ,converge para α .
Ex: Utilizando o método da iteração linear calcule a
raiz de ( ) 3xexf x += , com ( ) 5, 1 ≥+kk xxDIGSE
( )
33
33 00
x
x
xx
eex
exxexf
−=−=⇒
−=⇒=+⇒=
Seja
( )
( ) ( )
( ) 33,00'
24,01';
3
1
'
3
3
≅
≅−−=
−=
G
GexG
exG
x
x
*
G e G’ são continuas em [-1,0] e
( ) ]0,1[1' −∈∀< xxG .
Logo , a seqüência gerada por 31
ix
i ex −=+ converge para α ]0,1[−∈∀x .
Seja 5,00 −=x
M
772882595,0773204044,0
772884374,0771636903,0
772877469,0777723518,0
772904269,0754152577,0
772800243,0846481725,0
105
94
83
72
61
−=−=
−=−=
−=−=
−=−=
−=−=
xx
xx
xx
xx
xx
( )
34,5),(
000003188,0
109
9
≅
=
xxDIGSE
exf
*G não tem Maximo nem mínimo local em [0,1], testa-se entãosó os extremos.
Características do Método da Iteração Linear:
� Não garante a convergência para toda função continua.
� Necessita do calculo de G’(x).
� Pode ocorrer dificuldade para encontrar G(x).
16
2 - Resolução de Equações Algébricas e Transcendentes
� A convergência é linear para raízes simples (a cada passo do método o erro é reduzido por
um fato constante).
� A velocidade de convergência depende de ( )xG' , quanto menor este valor, maior será a
convergência.
Obs.: Dada a equação ( ) 0=xf , existem infinitas funções G(x) que são funções de iteração. A
forma geral destas funções é:
( ) )()( xfxAxxG ⋅+= ,
com a condição que em x*, ponto fixo de G(x), se tenha 0)( * ≠xA . Temos que:
*** )(0)( xxGxf =⇔= Com efeito:
(⇒ ) Seja *x tal que 0)( * =xf .
******* 0)()()()( xxAxxfxAxxG =⋅+=⋅+=
0)(0)(0)()(
)()()()(
****
******
≠=⇒=⋅⇒
=⋅+⇒=⇐
xApoisxfxfxA
xxfxAxxxGSe
2.3.3-Método de Newton-Raphson
Procura garantir e acelerar a convergência do método da Iteração Linear, escolhendo para a
função de iteração a função G(x) tal que G’(x)=0
Dada a função 0)( =xf , tomamos a forma geral para G(x):
),(')(1)(')(')()()('1)('
)(')()()('1)(')()()(
******** xfxAxGxfxAxfxAxG
xfxAxfxAxGxfxAxxG
⋅+=⇒⋅+⋅+=
⇒⋅+⋅+=⇒⋅+=
pois x * é ponto fixo de G )0)()(( *** =⇒= xfxGx .
Agora , )('
1)(0)(')(10)(
*
****
xfxAxfxAxG
−
=⇒=⋅+⇒= .
Tomemos )('
)()()('
1)(
xf
xf
xxGe
xfxA −=−= .
Então dada ( )xf , ( )xG é tal que ( ) 0' * =xG , pois
( ) 22
2
)]('[
)(")('
)]('[
)(")()]('[1'
xf
xfxf
xf
xfxfxf
xG ⋅=⋅−−= e
.0)('0)('0)( *** ≠=⇒= xfsexGxf
Portanto, iniciando-se a iteração com um valor 0x escolhido, a seqüência )( kx é determi-
nada por:
0)(',...2,1,0,)('
)(
1 ≠=−=+ kxfk
kxf
kxfkxkx
Geometricamente , conforme podemos observar na próxima figura, dado kx , 1+kx é abs-
cissa do ponto de intersecção da reta tangente à curva ( )xf em ))(,( kk xfx e o eixo dos “x”. As-
sim,
17
E.B.Hauser – Cálculo Numérico
)('
)()(')( 1
1 k
k
kkk
kk
k
xf
xf
xxxf
xx
xf
tg −=⇒=
−
= +
+
θ
Convergência: (é trabalhoso mostrar que 1)x('G < ).
O método de Newton-Raphson converge se:
2
2 ))('()(")(1))('(
)(")()(' xfxfxf
xf
xfxf
xG <⇒<= .
Para raízes simples a convergência é quadrática e para raízes duplas ou triplas é linear.
Escolha do ponto inicial: Seja )b,a(∈α raiz de f .
Se
bxbfbf
axafaf
=⇒>⋅
=⇒>⋅
0
0
0)(")(
0)(")(
Caso contrário, pode-se considerar
2
)(
0
ba
x
+
= .
Ex.: 1) Estimar o valor da única raiz real de xlnx2)x(f += , utilizando Newton-Raphson.
1kx +
θ
kx 1kx +
)x(f 1k +
)x(f k
xlnx2)x(f +=
x
12
xlnx2
x,x
−
+
−
18
2 - Resolução de Equações Algébricas e Transcendentes
k
kk
k1k
x
12
xlnx2
xx
+
+
−=+
34
3
2
1
0
xx
4263027510,0x
4262969599,0x
4232867952,0x
5,0x
=
=
=
=
=
Logo a aproximação para a raiz é 426302751,0=α com 9 dígitos significativos corretos.
2) Calcular a raiz ]4,3[4 ∈r do polinômio dado anteriormente: 11205,72)( 234 +−+= xxxxxp .
40)4(")4(0)3(")3( 0 =⇒>⋅<⋅ xppepp
201564
11205,72
23
234
1
−−+
−−−+
−=+
kkk
kkkk
kk
xxx
xxxx
xx
56
5
4
3
2
1
0
03524990,3
03525211,3
03709653,3
08982331,3
36397059,3
4
xx
x
x
x
x
x
x
=
=
=
=
=
=
=
Obs: O Método de Newton pode divergir devido ao uso da tangente, oscilando indefinidamente.
Isto acontece quando:
Uma aproximação para a raiz é 03524990,34 =r
com 9 dígitos significativos e 5 iterações.
19
E.B.Hauser – Cálculo Numérico
i) Não há raiz real
ii) Ocorre simetria de ( )xf em torno de α
iii) O valor inicial 0x está longe da raiz exata, fazendo que outra parte da função pren-
da a iteração.
2.3.4-Método da Secante
Uma desvantagem do Método de Newton-Raphson é o calculo do valor numérico de
( )xf ' a cada iteração.
O método da secante contorna este problema, substituindo a derivada pelo quociente das
diferenças:
1
1)()()('
−
−
−
−
≅
kk
kk
k
xx
xfxf
xf
onde kx e 1−kx são duas aproximações para a raiz de ( )xf . A formula iterativa é:
)()(
)()(
1
1
1
−
−
+
−
⋅−
=
kk
kkk
k
xfxf
xfxx
x
Geometricamente, a partir das aproximações para a raiz de kx e 1+kx , o ponto 1+kx é dado pela
abscissa do ponto de intersecção do eixo Ox e da reta secante que passa por
))(,())(,( 11 kkkk xfxexfx −− .
Características do método da secante:
� Garante a convergência para toda função continua
� Pode divergir se )()( 1−≅ kk xfxf
� Convergência mais lenta que o Método de Newton e mais rápida que Bissecção e Itera-
ção Linear.
� São necessárias duas aproximações da raiz a cada iteração.
Ex: )0,1(e21x172x53x)x(p −∈++−= α
)211kx17
2
1kx5
3
1kx()21kx17
2
kx5
3
kx(
)21kx172kx5
3
kx()1kxkx(kx1kx
+
−
+
−
−
−
−++−
++−⋅
−
−
−=+
7x8x
629321148566,07x
89321148568,06x
069321123947,05x
869317876515,04x
599601423487,03x
898888888888,02x
11x,10x
=
−=
−=
−=
−=
−=
−=
=−=
20
2 - Resolução de Equações Algébricas e Transcendentes
Exercícios
1) Uma partícula é arremessada verticalmente, a partir do solo, com uma velocidade inicial ov
.Desprezando a resistência do ar, supomos que a posição p partícula é dada por:
2
2
t
g
tov)t(d −= ,
onde g é aceleração da gravidade. Determinar a altura máxima atingida pela partícula e o instan-
te em que ocorreu.
2) Uma corrente oscilante num circuito elétrico é descrita por I = )t2sen(te9 pi− , t em segun-
dos. Determinar todos os valores positivos de t para os quais I = 4.
3) A concentração de uma bactéria poluente num lago é descrita por
C = t075.0e5.2t5.1e70 −+−
Encontrar o tempo para que a concentração seja reduzida para nove.
4) O deslocamento de uma estrutura está definido pela seguinte equação
D = wtcoskte8 −
onde k = 0.5 e w = 3.
a) Determinar graficamente uma estimativa inicial do tempo necessário para o deslocamento
decrescer para 4.
b) Usar o método de Newton-raphson para determinar essa raiz
5) Enumerar, localizar, separar e calcular (via Newton-Raphson e/ou Bissecção ), se possível, todas
as raízes dos polinomios tendo como critério de parada DIGSE (xk , xk+1) ≥ 5. No caso de raízes
múltiplas, determinar a multiplicidade da raiz e calcular as demais utilizando deflação.
a) p(x) = x x x3 22 3 5− + −
b) p(x) = x5 - 15,5x4 +77,5x3 - 155x2 +124x -31
c) p(x) = x x x x4 3 2121 2247 15043 34300− + − +
d) p(x) = x x x x4 3 2115 1575 7625 12500− + − +
e) p(x) = x x x x4 3 23 3 37 168 0 3136− + − +. . .
f) p(x) = x4-11,101x3+11,1111x2-1,011x+0,001
g) p(x) = x9- 4x8 + 3,9x7 +3,02x6 - 5,5295x5 - 0,84732x4 + 2,83536x3 + 0,37152x2 -0,59616x -
0,15552
h) p(x) = x3 - 2081,93x2 +1424,64x- 3,19728
i) p(x) = x4 + 1,98x3 +1,1424x2 +0,162602x - 0,00174225
21
E.B.Hauser – Cálculo Numérico
6) Localizar graficamente e calcular ( via Newton-Raphson e/ou Método da Iteração Linear) todas
as raízes, com DIGSE(xk , xk+1)≥ 5:
a) f(x) = x2 + ln x
b) f(x) = x + 2 cos x
c) f(x) = 2x + ln x
d) f(x) = cos x + ln x + x
e) f(x)= x + e Bx− 2 para B = 1,5,10,25,50
7) Responder resumidamente:
a) Em que consiste o princípio da bissecção para determinar a primeira aproximação de uma raiz
de uma equação f(x)=0?
b) Explicar o método da iteração linear para calcular uma raiz de uma equação f(x)=0, partindo de
uma primeira aproximação x0.
c) Quando não converge a iteração linear?
d) Quando não converge o Método de Newton Raphson?
e) Interpretar geometricamente o Método de Newton-Raphson.
f) Qual a vantagem de se utilizar o Algoritmo de Horner para se avaliar o valor do polinomio num
ponto? Exemplificar.
Respostas:
1) O deslocamento máximo é vo2/2g e ocorreu em t = vo/g.
2) t = 0.06835432097 e t = 0.4013436927
3) t = 1.556787935
22
2 - Resolução de Equações Algébricas e Transcendentes
4) t = 0 .3151660803
5-a)
r+ r- ¢ T
3 0 0 3
1 0 2 3
p tem raízes complexas.
Existe uma raiz real em (1,2)
Raízes:
x=1,84373427779
x= 0.07813286110 ± 1.644926378i
.
b) Raízes: .4555300547, 1.092601944, 1.940878206, 4.011783883, 7.999205912
c) Raízes: R1=100 e R2= 7(multiplicidade 3), não tem raízes complexas.
d) Raízes: R1=10 e R2=5(multiplicidade 3), não tem raízes complexas.
e) Raízes: R1=0.7(multiplicidade 2), R2=0.8(multiplicidade 2)
f) Raízes: R1=0,001 R2=0,1 R3=1 R4=10
g) Raízes: R1=-0,5(multiplicidade 4), R2=1,2(multiplicidade 5)
h) Raízes: R1= 0.002251681490, R2 = 0.6822607762 , R3= 2081.245488
i) Raízes: R1=-1,01 R2=-0,75 R3=-0,23 R4=0,01
6-a) R ≅ .6529186400 b) R ≅ -1.029866529 c) R ≅ 0 .42630275100 d) R ≅ .2875182754
e) Existe única raiz de f em (-1,0)
23
3. Sistemas de Equações Lineares
O sistema com n equações lineares e n variáveis
bxaxaxaxa
bxaxaxaxa
bxaxaxaxa
nnnn33n22n11n
2nn2323222221
1nn1313212111
=++++
=++++
=++++
L
MMMMM
L
L
pode ser representado matricialmente por BAX = , onde
=
aaa
aaa
aaa
A
nn2n1n
n22221
n11211
K
MMM
K
K
,
=
x
x
x
X
n
2
1
M
,
=
b
b
b
B
n
2
1
M
e
A é a matriz dos coeficientes, X é o vetor das incógnitas e B o vetor dos termos independentes.
Classificação quanto à existência de Solução: Regra de Cramer
A solução do sistema linear AX=B é dada pela fórmula: )det(
)det(
A
A
x ii = , onde, i=1,2,..,.n, )det( iA é o
determinante da matriz obtida de A substituindo a i-ésima coluna de A por B.
O sistema linear :
• tem única solução se 0)det( ≠A ;
• não tem solução se )det(A = 0 e, para pelo menos um i , 0)det( ≠iA .
• tem infinitas soluções se 0)det( =A e para todos i , 0)det( =iA .
Exemplos:
(a)
−=−
=+
7y2x
8yx
tem única solução.
(b)
=+
=+
0y6x10
8y3x5
não tem solução.
(c)
=−
−=+−
3y2x
9y6x3
tem infinitas soluções.
24
E.B.Hauser – Cálculo Numérico
3.1- Introdução à problemática de sistemas
Um SEL pode ser mal condicionado, bem condicionado ou sem solução.
Um sistema é “mal condicionado” se uma pequena alteração nos dados pode provocar uma
grande alteração na solução. Por exemplo:
(a)
=+
=+
5yx
95,4y98,0x
tem solução exata
=
5,2
5,2
x
(b)
=+
=+
5yx
95,4y99,0x
tem solução exata
=
0,5
0,0
x
Uma alteração de 1% (0,01 no coeficiente 0,98) ocasionou uma alteração de 100% na solução.
No caso de um sistema linear de ordem 2, cada equação representa uma reta. Resol-
ver o sistema significa determinar a intersecção das duas retas. Três casos são possíveis:
Bem condicionado Não tem solução. Mal condicionado
0det ≠ 0det = 0det ≅
retas concorrentes retas paralelas (perturbação em 2ℜ )
Exemplo2:
(a)
=+
=+
503,25y501,7x5,1
17y5x
tem solução exata:
=
3
2
X
(b)
=+
=+
501,25y501,7x5,1
17y5x
tem solução:
=
1
12
X
3.2-Medidas de Condicionamento
O determinante normalizado da matriz dos coeficientes A é dado por
2
kn
2
2k
2
1kk
n21
aaaondeAdetANORM L
L
++== α
ααα
, para k = 1, 2, ..., n ( )1,1ANorm −∈
e quanto mais afastado de 1± (isto é, quanto mais próximo de zero) mais mal condicionado é a ma-
triz A.
Retomando o Ex2:
R1 R2 R2
R1
R2
R1
25
Sistemas de Equaçõe Lineares
2)
=
501,75,1
51
A
( )
( ) 586495509853,7501,75,1
90990195135,5251
2
1
22
2
2
1
1
=+=
=+=
α
α
38740000256377,0
0050000128,39
001,0001,0Anorm
001,05,15501,7Adet
21
==
⋅
=
=⋅−=
αα
Agora, como pode ser medido o condicionamento de um sistema linear?
Dado um SEL BAX = , o seu número de condicionamento é dado por:
1AA)A(Cond −= .
Quanto maior for )A(Cond , mais sensível será o sistema linear.
Utilizamos ∑==
=≤≤
∞
n
1i
ij
ni1
amaxAA , a norma do máximo das linhas.
Ex:
=
501,75,1
51
A ,
−
−
=
−
10001500
50007501
A 1
51
1
101501,112521AA)A(Cond
12501A,001,9A
⋅≅==
==
−
−
3.3-Método de Resolução de Sistemas
Métodos Diretos: A solução exata é obtida realizando-se um número finito de operações a-
ritméticas em ℜ (com precisão infinita): Eliminação de Gauss e Gauss Jordan.
Métodos Iterativos: A solução x é obtida como limite de uma seqüência de aproximações
sucessivas x1, x2, ... .
Método de Gauss-Jordan (Matriz Inversa)
A solução do SEL BAX = é calculado utilizando BAX 1−= se 0Adet ≠ .
Obs.: Ver exercício 9.
26
E.B.Hauser – Cálculo Numérico
Método de Eliminação de Gauss
Algoritmo básico de Gauss: A solução de BAX = é calculada em duas etapas:
1º- Triangularização : Mediante operações elementares nas linhas, a matriz A é transforma-
da numa matriz triangular superior.
Algoritmo: para 1n,,1k −= K (indica a linha do pivô)
para n,,1ki K+= (indica a linha a transformar de A)
kk
ik
a
a
m
−
=
0aik =
para nkj ,,1K+= (indica a coluna a transformar da linha i)
amaa kjijij ⋅+=
bmbb kii ⋅−=
Obs.: Se 0aii = é necessário trocar de linha, se possível.
2º- Retro-substituição: A triangularização transforma o sistema BAX = , no sistema equi-
valente:
dxc
dxcxc
dxcxcxc
dxcxcxcxc
nnnn
3nn3333
2nn2323322
1nn1313212111
=
=++
=+++
=++++
LLLLLLL
L
L
L
cuja solução é dada por: ( )
( )
11
nn13132121
1
1n,1n
nn,1n1n
1n
nn
n
n
a
xcxcxcd
x
,
a
xcd
x,
c
d
x
+++−
=
−
==
−−
−−
−
L
Teorema: O método de Gauss produz sempre a solução exata do sistema BAX = (utilizando preci-
são infinita) se 0det ≠A e as linhas de Aforem permutadas quando 0aii = .
27
Sistemas de Equaçõe Lineares
Exemplo: Estimar a solução do sistema linear utilizando o método da eliminação de Gauss.
=−++−
=+++
−=−++−
=+−−
174x3x5,152x221x12
164x5,53x25,62x5,71x
5,64x63x5,62x71x6
124x63x32x21x4
−
−−
−−
−−
17
16
5,6
12
15,152212
5,525,65,71
65,676
6324
L2=L2+(6/4)L1; L3=L3-(1/4)L1; L4=L4+3L1
−−
53
13
5,11
12
175,6160
4780
3240
6324
L3=L3-2L2; L4=L4-4L2
−
−
−
−−
7
10
5,11
12
55,100
2300
3240
6324
L4=L4+(1,5/3)L3
−−
−−
2
10
5,11
12
4000
2300
3240
6324
| Obtemos o sistema triangular, equivalente ao sistema dado.
A solução é obtida por substituição.
=
−=−
=++
=+−−
24x4
104x23x3
5,114x33x22x4
124x63x32x21x4
5,04x,33x,42x,21x =−===
28
E.B.Hauser – Cálculo Numérico
Quando é utilizada aritmética do ponto flutuante, erros de arredondamento podem compro-
meter a solução obtida.
Exemplo: Em ( )98,98,3,10F − , com arredondamento para número mais próximo de máquina
“Ox”, a elminação de Gauss aplicada ao sistema
=+
=+
6y2x2
5y2x0002,0
produz 0x = e 5,2y = o que
não satisfaz a segunda equação do sistema.
(Obs.: multiplicador 000.10m −= por ( )mLLL 122 += )
Gauss com Pivotamento Parcial
O método consiste em trocar linhas (ou colunas) de maneira a minimizar a propagação de
erros nas operações.
Escolhas dos pivôs:
1º pivô é o elemento de maior valor absoluto da coluna 1.
2º pivô é o elemento de maior valor absoluto da coluna 2 e da diagonal para baixo.
Procede-se da mesma forma para os demais pivôs:
Aplicando a técnica ao último exemplo
=+
=+
5y2x0002,0
6y2x2
, em ( )98,98,3,10F − , com
arredondamento para número mais próximo de máquina “Ox”, obtemo 5,0x = e 5,2y = .
(Obs.: Neste caso o multiplicador é menor m= - 0,0001)
1º pivô
2º pivô 3º pivô
29
Sistemas de Equaçõe Lineares
Métodos Iterativos
Os sistemas lineares de grande parte são, em geral, esparsos (muito elementos 0aij = ). Os
métodos diretos não preservam a esparsidade, enquanto que os métodos iterativos sim, além de a-
presentarem relativa insensibilidade ao crescimento dos erros de arredondamento.
No sistema original BAX = , supondo n,,1i,0aii K=≠ , o vetor X é isolado mediante a se-
paração de diagonal principal:
( )
( )
( )1n1n,n22n11nn
nn
n
nn23231212
22
2
nn13132121
11
1
xaxaxab
a
1
x
xaxaxab
a
1
x
xaxaxab
a
1
x
−−
−−−−=
−−−−=
−−−−=
K
MM
K
K
Metodo de Gauss-Jacobi
Dada a aproximação inicial X0, o processo iterativo produz aproximações sucessivas
KK ,,,,
21 kXXX , obtidas de:
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )( )
( ) ( ) ( ) ( )
−−−−=
−−−−=
−−−−=
−
−
+
+
+
k
1n1n,n
k
22n
k
11nn
nn
1k
n
k
nn2
k
323
k
1212
22
1k
2
k
nn1
k
313
k
2121
11
1k
1
xaxaxab
a
1
x
xaxaxab
a
1
x
xaxaxab
a
1
x
K
MM
K
K
Método de Gauss-Seidel
Para X0 dado:
( ) ( ) ( ) ( )( )
( ) ( ) ( ) ( )( )
( ) ( ) ( ) ( )( )
( ) ( ) ( ) ( )( )1k 1n1n,n1k22n1k11nn
nn
1k
n
k
nn3
)k(
434
1k
232
1k
1313
33
1k
3
k
nn2
k
323
1k
1212
22
1k
2
k
nn1
k
313
k
2121
11
1k
1
xaxaxab
a
1
x
xaxaxaxab
a
1
x
xaxaxab
a
1
x
xaxaxab
a
1
x
+
−
−
+++
+++
++
+
−−−−=
−−−−−=
−−−−=
−−−−=
K
MM
K
K
K
30
E.B.Hauser – Cálculo Numérico
Critério de Convergência
O teorema abaixo estabelece uma condição suficiente para a convergência dos métodos
de Gauss-Jacobi e de Gauss-Seidel.
Teorema
Dado o sistema linear YAX = , se a matriz A é Diagonalmente Dominante, isto é,
se iaa
ij ijii
∀∑>
≠
, então tanto o método de Jacobi como o de Gauss-Seidel gera uma
seqüência ( )( )kX convergente para a solução do sistema dado, independente da escolha da
aproximação inicial ( )0X .
Exemplo:
Resolver o sistema abaixo por Gauss- Jacobi e Gauss-Seidel.
Critério de Parada: erro absoluto da solução calculada for menor que 10-3.
−
−
=+
=+
=+−
=−+
5,1
1
3
2
xx5,0
xx2
x2x
xx4x
31
41
43
421
Reordenamos a fim de satisfazer ao critério de convergência.
−
−
=+−
=+
=−+
=+
3
5,1
2
1
x2x
xx5,0
xx4x
xx2
43
31
421
41
12
5,01
114
12
−>
>
−+>
>
Como a matriz dos coeficientes , após a reordenação, é diagonalmente dominante, então a
aplicação dos métodos de Gauss-Jacobi e Gauss-Seidel produzirá uma seqüência ( )( )kX convergen-
te para a solução exata.
31
Sistemas de Equaçõe Lineares
Gauss-Jacobi
Fórmulas iterativas:
( ) ( )
( ) ( ) ( )
( ) ( )( )
( ) ( )k
3
)k(
31k
4
k
1
1k
3
k
4
k
1
)k(
4
)k(
11k
2
k
4
)k(
41k
1
x5,05,1
2
)x3(
x
x5,05,1x
x25,0x25,05,0
4
)xx2(
x
x5,05,0
2
x1
x
+−=
+−
=
−=
+−−=
+−−
=
−=
−
=
+
+
+
+
Aproximação inicial: ( ) 0X 0 = .
k x1k x2k x3k x4k
0 0 0 0 0
1 0,5 0,5 1,5 -1,5
2 1,25 -1,0 1,25 -0,75
3 0,875 -1,0 0,875 -0,875
4 0,9375 -0,9375 1,0625 -10625
5 1,03125 -1,0000 1,03125 -0,96875
6 0,984375 -1,0000 0,984375 -0,984375
7 0,9921875 -0,9921875 1,0078125 -1,0078125
8 1,00390625 -1,0000 1,00390625 0,9960375
9 0,998046875 -1,0000 0,998046875 -0,998046875
10 0,9990234375 -0,999023437 1,000976563 -1,000976563
11 1.000488282 -1.0000000 1.000488281 -0,999511718
: : : : :
40 1 -1 1 -1
na 12a. iteração consegue-se ( ) 51112 10xx −<−
32
E.B.Hauser – Cálculo Numérico
Gauss-Seidel
Fórmulas iterativas:
( ) ( )
( ) ( ) ( )
( ) ( )
( ) ( )1k
3x5,05,1
1k
4x
1k
1x5,05,1
1k
3x
k
4x25,0
1k
1x25,05,0
1k
2x
k
4x5,05,0
1k
1x
++−=+
+
−=
+
++−−=+
−=
+
k x1k x2k x3k x4k
0 0 0 0 0
1 0,5 -0,625 1,25 -0,875
2 0,9375 0,953125 1,03125 -0,984375
3 0,9921875 -0,994140625 1,00390625 0,99806875
4 0,9990234375 -0,9992675781 1,000488281 -0,999755895
5 0,999877929 -0,999908447 1,000061035 -0,999969482
: : : : :
12 1 -1 1 -1
( ) ( ) 545 10xx −<−
33
Sistemas de Equaçõe Lineares
Aplicação
A distribuição de temperatura numa placa, quando a temperatura nas bordas é conhecida, é uma
consideração importante no estudo da transferência de calor . Supomos que a placa da figura repre-
sente um seção transversal de uma barra de metal, com fluxo de calor desprezível na direção per-
pendicular à placa. Sejam T1, T2 , T3 e T4 as temperaturasnos quatro vértices interiores do reticu-
lado da figura. A temperatura num vértice é aproximadamente igual à média aritmética dos quatro
vértices vizinhos mais próximos (à esquerda, acima, à direita e abaixo).
Estimar T1, T2 , T3 e T4, utilizando o método de Gauss Seidel.
10 0 10 0
0 0
T1
T2
0 0
0 0
T3
T4
0 0
10 0 10 0
)1k(3T25,0)1k(2T25,05,2)1k(4T
)k(4T25,0)1k(1T25,05,2)1k(3T
)k(4T25,0)1k(1T25,05,2)1k(2T
)k(2T25,0)k(3T25,05,2)1k(1T
++++=+
+++=+
+++=+
++=+
k T1k T2k T3k T2k
0 0 0 0 0
1 2,5 3,125 3,125 4,0625
2 4,0625 4,53125 4,53125 4,765625
3 4,765625 4,8828125 4,8828125 4,94140625
4 4,94140625 4,970703125 4,970703125 4,9853515625
5 4,9853515625 4,99267578125 4,99267578125 4,996337890625
6 4,996337890625 4,9981689453125 4,9981689453125 4,99908447265625
7 4,99908447265625 4,999542236328125 4,999542236328125 4,9997711181640625
8 4,9997711181640625 4,99988555908203125 4,99988555908203125 4,99994277954101562
: : : : :
5 5 5 5
34
E.B.Hauser – Cálculo Numérico
EXERCÍCIOS
1) Uma consideração importante no estudo da transferência de calor é a de se determinar a
distribuição de temperatura numa placa, quando a temperatura nas bordas é conhecida.
Supomos que a placa da figura represente um seção transversal de uma barra de metal, com
fluxo de calor desprezível na direção perpendicular à placa. Sejam T1, T2 , ..., T6 as tem-
peraturas nos seis vértices interiores do reticulado da figura. A temperatura num vértice é aproxi-
madamente igual à média aritmética dos quatro vértices vizinhos mais próximos (à esquerda, aci-
ma, à direita e abaixo). Por exemplo:
T1 = ( 10 + 20 + T2 + T4 )/4 ou 4T1 = 10 + 20 + T2 + T4
a) Escrever o sistema de seis equações, cuja solução fornece estimativas para as temperatu-
ras T1, T2 , ..., T6.
b) Resolver o sistema.
2) Num experimento num túnel de vento, a força sobre um projétil devido à resis-
tência do ar foi medida para velocidades diferentes.
119
10
3.74
8
6.39
6
8.14
4
90.2
2
0
0
força
velocidade
Expressar a força como função da velocidade aproximando-a a um polinômio de quinto
grau:
5
5
4
4
3
3
2
21o vavavavavaa)v(f +++++=
Verificar a validade de )(vf encontrada e obter uma estimativa para a força sobre o projétil
quando ele estiver se deslocando a uma velocidade de 1, 3, 5, 7 e 9 unidades de velocidade.
3) Considere o sistema (#)
=++
=++
=++
78,0x2,0x25,0x33,0
08,1x25,0x33,0x5,0
83,1x33,0x5,0x
321
321
321
.
a) (#) é bem ou mal condicionado? Porque? O que isso significa?
b) Resolver (#) pelo método de Gauss sem pivotamento.
1
20o 20o 20o
20o 20o 20o
40o
40o
10o
10o
1 2 3
4 5 6
35
Sistemas de Equaçõe Lineares
4) Idem ao 3 para
=++++
=++++
=++++
=++++
=++++
2520/1879x9/1x8/1x7/1x6/1x5/1
840/743x8/1x7/1x6/1x5/1x4/1
420/459x7/1x6/1x5/1x4/1x3/1
60/87x6/1x5/1x4/1x3/1x2/1
60/137x5/1x4/1x3/1x2/1x
54321
54321
54321
54321
54321
5) Resolver por Eliminação de Gauss com pivotamento parcial.
a)
=+−++−+
−=−−−+−+
=+++−++
=−−++++−
=++−++−
=+−++++
=+−−++−
34x9xx7x2xx2x7
10xx4xx11x8x9x2
13xxx2x7xx3x4
8x2xx7x2x3xx
24x11x6xx5x7xx2
17x2x9xx15x9xx
47x8x3x2x4x6x5x3
7654321
7654321
7654321
7654321
7654321
7654321
7654321
b)
=−−+
=−−
=++−
255x10x35x
4414x41x72x
8342x2xx
321
321
321
c)
=−++
−=++−
=−++
=+−+
9.9xx5.2x2.0x
9.3x2.5xxx3.0
9.21x5.8x4x5.0x4.0
7.2xx1.0xx2
4321
4321
4321
4321
6) Referente ao sistema linear AX=B, verificar se a afirmação é Verdadeira ou falsa .
i) Se det A = 0 então o sistema não tem solução .Justificar verificando o que acontece em :
a)
=+−
=+
=+−
2xx2
4xx
5xxx
32
21
321
e b)
=+
=−
=+
0xx
0xx
0xx
32
31
21
ii) Se A não é uma matriz Diagonal Dominante então os métodos de Gauss-Jacobi e Gauss-Seidel
não geram uma sequência convergente para a solução . Justificar verificando o que acontece com
a)
−=−
=+
3y3x
3yx
b)
=+
−=−
3yx
3y3x
7) Resolver pelo Método de Gauss-Seidel. Apresentar as fórmulas iterativas e uma garantia de con-
vergência (se possível).
a)
=+−+
=−++
=+++
=++−
82x7x12x66x17
34x18x11x22x56
26x10x17x16x11
18x20x47x5x3
4321
4321
4321
4321
b)
=++
=+−
=++
31x6x3x2
1x2x9x
30x3x2x8
321
321
321
8) Resolver o sistema de equações algébricas não lineares:
=++
=++
=++
15z4yx
18zy4x
11zyx4
2
2
2
36
E.B.Hauser – Cálculo Numérico
9) Resolver
=++
=++
=++
3321
2321
1321
bx4x3x5
bx2x3x
bx7xx2
para:
a) b1 =16 b2 = -5 b3=11 b) b1 =25 b2 = -11 b3 = -5 c) b1 =3 b2 = 5 b3 = -5
10) Utilizando Eliminação Gaussiana calcular detA.
a) A =
−−−
−−−
−−
−−−−
1322
1413
2101
2132
b) A =
−
01301
25.09110
30350
00023
05.0807
11) Qual o Resíduo produzido pela solução aproximada X’= [ -3 4 0]T de
=++
=++
=++
64.0x25.0x21.0x15.0
52.0x24.0x16.0x12.0
84.0x12.0x36.0x24.0
321
321
321
Respostas:
1) Solução obtida utilizando Maple:
>solve({4*T1=10+20+T2+T4,4*T2=T1+20+T3+T5,4*T3=T2+20+40+T6,
4*T4=10+T1+T5+20,4*T5=T4+T2+T6+20, 4*T6=T5+T3+40+20}, {T1,T2,T3,T4,T5,T6});
{T6 = 190/7, T5 = 150/7, T4 = 120/7, T3 = 190/7, T1 = 120/7, T2 = 150/7}
> evalf(%);
{T6 = 27.14285714, T5 = 21.42857143, T4 = 17.14285714,
T3 = 27.14285714, T1 = 17.14285714, T2 = 21.42857143}
2) Solução utilizando o sistema Maple:
>solve({a0=0,2.90=a0+a1*2+a2*4+a3*8+a4*16+a5*32,
14.8=a0+a1*4+a2*(4^2)+a3*(4^3)+a4*(4^4)+a5*(4^5),
39.6=a0+a1*6+a2*(6^2)+a3*(6^3)+a4*(6^4)+a5*(6^5),
74.3=a0+a1*8+a2*(8^2)+a3*(8^3)+a4*(8^4)+a5*(8^5),
119=a0+a1*10+a2*(10^2)+a3*(10^3)+a4*(10^4)+a5*(10^5)}, {a0,a1,a2,a3,a4,a5});
{a0 = 0, a2 = -1.194791667, a5 = .002604166667, a4 = -.07005208333,
a3 = .6614583333, a1 = 1.712500000}
>f:=x->1.712500000*x-1.194791667*x^2+.6614583333*x^3-.7005208333e-
1*x^4+.2604166667e-2*x^5;
> validade:=[f(0),f(2),f(4),f(6),f(8),f(10)];
validade := [0, 2.899999998, 14.80000000, 39.60000000, 74.29999994, 119.0000000]
> estimativas:=[f(1),f(3),f(5),f(7),f(9)];
estimativas := [1.111718750, 7.202343750, 25.73046873, 55.89609367, 94.9992188]
37
Sistemas de Equaçõe Lineares
3-a) Mal condicionado. NORM A ≅ 0,000181. Uma pequena perturbação nos dados de entrada po-
de causar uma grande alteração na solução.
b) A solução exata é X=[1 1 1]T
4 - NORMA ≅ 0,00257, a solução exata é X=[1 1 1 1 1 1 1]T
5- a)X=[-21,86188 11,46568 2,376447 -8,514801 0,7475478 -15,50981 18,08498]T
b) solução exata X=[1 0 2]T
c) solução exata X=[1 2 3 -1]T
6-i) a) detA=0 e o sistema e incompatível
b) detA = 0 e o sistema tem infinitas soluções dadas por x=z e y=-z
ii) A não é matriz Diagonal Dominante e Gauss-Jacobi e Gauss-Seidel
a) converge (oscilando) para a solução exata [1.5 1.5]T
b) diverge
7 - a) X(10) = [-0.930569 1.901519 1.359500 -1.906078]T
X(35) = [-1.076888 1.990028 1.474477 -1.906078]T
b) solução exata X = [2 1 4]T
8 - A solução exata é x =1 y= 2 z = 3.
9 - soluções exatas: a) X = [3 -4 2]T b) X = [2 -7 4]T c) X = [-3 2 1]T
10 - a) detA = -55 b) det A = 706.875
11 - R = [0.12 0.24 0.25]T
38
4.Interpolação Polinomial
4.1 Objetivo:
Dado um conjunto de n+1 pontos distintos ))x(f,x( ii , n,...,1,0i = , queremos determinar o
polinômio p(x) talque ),x(f)x(p ii = e para os demais pontos do intervalo
)x(f)x(p ≅ , [ ]no x;xx ∈∀ . O polinômio p(x) é dito polinômio aproximante ou interpolador de
f(x) no intervalo [ ]nxox , .
4.2 Aplicações
• Obter uma expressão analítica aproximada de uma função que é conhecida em apenas um
número finito de pontos;
• Avaliar a função num ponto não tabelado [ ]nxoxx ;*∈ ;
• Determinar aproximações para ∫ n
x
ox
dxxf )( , substituindo a f(x) pelo polinômio interpola-
dor;
• Calcular uma aproximação para )(' xf para [ ]nxoxx ;∈ , substituindo f(x) por p(x).
4.2 Existência e Unicidade da Solução
Dados ℜ∈ix e ,)x(f i ℜ∈ n,...,1,0i = , procuramos nP)x(p ∈ tal que
iii x),x(f)x(p ∀= .
Seja ∑
=
=++++=
n
k
kxka
nxnaxaxaoaxp
0
...
2
21)( .
39
E.B.Hauser – Cálculo Numérico
Então de =)x(p i ∑
=
=
n
0k
i
k
ik ),x(fxa obtemos:
=++++
=++++
=++++
n
n
nn
2
n2n10
1
n
1n
2
12110
0
n
0n
2
02010
fxa...xaxaa
..............................................
fxa...xaxaa
fxa...xaxaa
,
o qual representa um sistema linear de ordem n+1 onde as n+1 incógnitas são os ,ka n,0k = e a
matriz dos coeficientes é dada por:
A =
n
n
2
nn
n
2
2
22
n
1
2
11
n
0
2
00
x..xx1
......
......
......
x..xx1
x..xx1
x..xx1
De acordo com Vandermonde, det A =
−∏∏
+=
−
=
n
1ij
ij
1n
0i
)xx( . Como os pontos são
distintos, a diferença ixjx − será sempre diferente de zero, e então detA ≠ 0. Portanto o polinômio
interpolador existe e também é único.
4.3Polinômio Interpolador de Newton Para Diferenças Finitas Ascendentes
Dados ( )ii y,x , )x(fy ii = i=0,1,2,...,n satisfazendo
hxx....xxxx 1nn1201 =−==−=− − ,
isto é , hxx i1i =−+ .
Para k = 1, 2, ...n, e i= 0, 1, 2, ..., n-k , a diferença finita de ordem k é dada por
i
1k
1i
1k
i
k yyy −+
−
−= ∆∆∆ .
E, o polinômio interpolador de Newton para diferenças finitas ascendentes é dado por :
o
n
n
1n1o
o
2
2
1o
0
o
o y
h!n
)xx()xx)(xx(
y
h!2
)xx)(xx(
y
h
)xx(
y)x(p ∆∆∆ −−−−++−−+−+= LL
40
4-Interpolação Polinomial
Exemplo1:
O alongamento de uma mola foi medido em função da carga aplicada. Obteve-se:
3,6
8
0,5
6
5,2
4
0,1
2
)cm(oalongament
)kg(aargc
Estimar o alongamento para o caso de ser aplicada uma carga de 7kg , utilizando o polinômio in-
terpolador de Newton para diferenças finitas.
Solução:
1) Tabela das diferenças finitas:
i ix iy iy∆ i2 y∆ i3 y∆
0 2 1 1.5 1 -2.2
1 4 2.5 2.5 -1.2 --------------------
2 6 5 1.3 -------------------- --------------------
3 8 6.3 -------------------- -------------------- --------------------
2) Polinômio interpolador:
( ) ( )( ) ( )( )( ) )2,2.(
)2(!3
6x4x2x1.
)2(!2
4x2x5,1.
2
2x1)x(p 32 −
−−−
+
−−
+
−
+=
7,2x016666667,2x675,0x30458333333,0)x(p 23 +−+−=
3) Verificação de validade de p(x) : 9999999994,0)2(p = ( 1≅ )
4999999999,2)4(p = ( 5,2≅ )
5)6(p = e p(8)=6.3
4) Estimativa do alongamento ao se aplicar uma carga de 7kg:
O alongamento da mola neste caso é aproximadamente 9375,5)7(p = cm.
5) Análise gráfica do problema:
41
E.B.Hauser – Cálculo Numérico
4.4 Polinômio Interpolador de Newton Para Diferenças Finitas Divididas
Dados ( )iyix , , ),( ixfiy = i= 0, 1, 2, ..., n os pontos ix podem ter um espaçamento
qualquer, não necessariamente eqüidistantes.
O polinômio de Newton para diferenças divididas é dado por:
oy
n
n
xxxxoxxoyxxoxxyoxxoyxp ∆
−
−−−++∆−−+∆−+= )1)...(1)((...
2)1)((0)()( ,
onde, para k = 1, 2, ..., n, e i= 0, 1, 2, ..., n-k a diferença dividida de ordem k é dada por
iki
i
1k
1i
1k
i
k
xx
yy
y
−
−
=
+
−
+
− ∆∆
∆
Por exemplo, para o caso de n = 4:
i ix iy iy∆ i
2 y∆ i
3 y∆ i
4 y∆
0 0x
0y
01
01
xx
yy
−
−
02
01
xx
yy
−
− ∆∆
03
0
2
1
2
xx
yy
−
− ∆∆
04
0
3
1
3
xx
yy
−
− ∆∆
1 1x 1y
12
12
xx
yy
−
−
13
12
xx
yy
−
− ∆∆
14
1
2
2
2
xx
yy
−
− ∆∆
-----------------------
-----------------------
--------------------
2 2x
2y
23
23
xx
yy
−
−
24
23
xx
yy
−
− ∆∆
---------------------
---------------------
------------------
-----------------------
-----------------------
--------------------
3 3x
3y
34
34
xx
yy
−
−
---------------------
---------------------
------------------
---------------------
---------------------
------------------
-----------------------
-----------------------
--------------------
4 4x
4y
--------------- -------------------- -------------------- ----------------------
Observação: Para h constante , a relação entre diferenças divididas e finitas é dada por :
i
k
ki
k y
h!.k
1y ∆∆ = .
Exemplo:
x f(x) iy∆ i
2 y∆ i
3 y∆ i
4 y∆
0 0 4 5 1 0
2 8 19 10 1 -----
3 27 49 14 ----- -----
5 125 91 ----- ----- -----
6 216 ----- ----- ----- -----
42
4-Interpolação Polinomial
Exemplo2:
A velocidade do som na água varia com a temperatura conforme tabela:
532,1
110
538,1
4,104
544,1
9,98
548,1
3,93
552,1
86
)s/m(velocidade
) C(atemperatur o
Solução:
1) Cálculo das diferenças divididas
x y 0y∆ 0
2y∆ 0
3y∆ 0
4y∆
86 1.552 -0.00054794 -0.00001289 -0.114 x 10-5 0.136 x 10-6
93.3 1.548 -0.00071428 -0.00003393 0.213 x 10-5 ---------------
98.9 1.544 -0.00109090 0.175 x 10-5 ------------------------------
104.4 1.538 -0.00107142 --------------- --------------- ---------------
110 1.532 --------------- --------------- --------------- ---------------
2) Construção do polinômio:
)9.98x)(3.93x)(86x()00001289.0)(3.93x)(86x()00054794.0)(86x(552.1)x(p −−−+−−−+−−+=
(-0.114 x 10 5− )+ )4.104x)(9.98x)(3.93x)(86x( −−−− (0.136 x 10 6− )
Simplificando a expressão , encontramos
5036369194.0x830077947032.0x99260000534327.0x10840.13666902p(x) 234-6 −+−×=
3) Verificação da validade de p(x) calculado no item 2:
1.55199999)86(p = ; 548.1)3.93(p = ; 544.1)9.98(p = ;
537999999.1)4.104(p = ; 532.1)110(p =
4) Estimativa da velocidade do som quando a temperatura da água atinge 100 C0 , é
1.54293925)100(p ≅ m/s
6) Análise gráfica do problema:
43
E.B.Hauser – Cálculo Numérico
4.5 Erro de Truncamento
[ ]nxxnf
n
x
xE ,0),()1()!1(
)()( ∈+
+
= ξξϕ com
i
1n
n101n
0
1n
n10 y)*xx)...(xx)(xx(
h)!1n(
f)xx)...(xx)(xx()x( +
+
+
−−−=
+
−−−= ∆∆ϕ ,
pois i
k
ki
k f
h!.k
1y ∆∆ = , para h constante.
4.6 Aplicações utilizando o sistema Maple
APLICAÇAO 1
Cinquenta animais de uma espécie ameaçada de extinção foram colocados numa re-
serva e um controle populacional mostrou os dados:
69
13
73
10
76
7
77
5
73
3
60
1
50
0
animaisdequantidade
)anos(t
a) Determinar a matriz de Vandermonde para o problema e determinar o valor do respectivo
Número de Condicionamento (Cond e Determinante Normalizado). O que podemos conclu-
ir?
b) Determinar o polinômio interpolador utilizando todos os dados tabelados.
c) Verificar a validade do modelo encontrado.
d) Estimar a a população no quarto ano. É possível estimar a população no décimo quinto ano utili-
zando o polinômio determinado no ítem b.
e) Em que ano essa população animal atingiu seu máximo? Qual a população máxima atingida?
f) Plotar num mesmo sistemas de eixos os pontos tabelados e a e o polinômio interpolador determi-
nado no ítem b.
Resposta:
=
65432
65432
65432
65432
65432
1313131313131
1010101010101
7777771
5555551
3333331
1111111
0000001
V
Determinante Normalizado= 0.9072675023 x 10 -11
Cond (V)= 39124291.11
p(t) = .5095984263e-4*t^6-.2488977072e-2*t^5+.4187474562e-1*t^4-.2409954552*t^3-
.6536635972*t^2+10.85522232*t+50.
p(4) ≅ 75.91851648
população máxima ≅ p(5.312779138)= 77.05456458
44
4-Interpolação Polinomial
APLICAÇÃO 2
Para determinar a resistência elétrica de um solo num sistema
de aterramento, enterra-se duas hastes de cobre e aplica-se
uma determinada voltagem, resultando numa corrente elétri-
ca. Numa experiência deste tipo, foram obtidos os seguintes
dados:
5,4
50
3,4
47
5,3
40
8,2
35
2
30
))A(amperecorrente(y
))V(voltsvoltagem(x
−
−
Estimar a corrente se a voltagem aplicada for de 43V usando o Po-
linômio Interpolador de Newton.
Resposta: p(43) ≅ 3,88
APLICAÇÃO 3
“Ao considerar que no Japão a vida média já é superior a 81 anos, a esperança de vida no Brasil de
pouco mais que 71 anos ainda é relativamente baixa. E, de acordo com a projeção mais recente da
mortalidade, somente por volta de 2040 o Brasil estaria alcançando o patamar de 80 anos de espe-
rança de vida ao nascer. (Ver www.ibge.gov.br em População / Tábuas Completas de Mortalidade).
A esperança de vida ao nascer de 71,3 anos coloca o Brasil na 86ª posição no ranking da ONU,
considerando as estimativas para 192 países ou áreas no período 2000-2005 (World Population
Prospects: The 2002 Revision; 2003). Entre 1980 e 2003 a esperança de vida ao nascer, no Brasil,
elevou-se em 8,8 anos: mais 7,9 anos para os homens e mais 9,5 anos para as mulheres. Em 1980,
uma pessoa que completasse 60 anos de idade teria, em média, mais 16,4 anos de vida, perfazendo
76,4 anos. Vinte e três anos mais tarde, um indivíduo na mesma situação alcançaria, em média, os
80,6 anos. Aos 60 anos de idade os diferenciais por sexo já não são tão elevados comparativamente
ao momento do nascimento: em 2003, ao completar tal idade, um homem ainda viveria mais 19,1
anos, enquanto uma mulher teria pela frente mais 22,1 anos de vida”.Na tabela acima obtemos
45
E.B.Hauser – Cálculo Numérico
informações sobre a esperança de vida às idades exatas, especialmente, a esperança de vida ao
nascer que expressa o número médio de anos que se espera viver um recém-nascido (que, ao
longo da vida, estivesse exposto aos riscos de morte da tábua de mortalidade em questão
http://www.ibge.gov.br/home/presidencia/noticias/noticia_visualiza.php?id_noticia=266&id_pagina=1).
Além dos múltiplos usos, não somente pela Demografia e Previdência Provada, mas também pelas
demais Ciências Sociais, a tabela é utilizada para determinar, juntamente com outros parâmetros, o
chamado fator previdenciário para o cálculo das aposentadorias das pessoas regidas pelo Regime
Geral de Previdência Social.
TAREFA: Considerar os resultados de 2003.
0.15
70
4.18
65
1.22
60
0.26
55
1.30
50
3.48
30
0.53
25
8.57
20
6.62
15
5.67
10
2.75
0
Mulheres)2003em,anos(vidadeectativaexp
2003emidade
−
1.13
70
9.15
65
1.19
60
5.22
55
2.26
50
5.42
30
8.46
25
0.51
20
5.55
15
4.60
10
6.67
0
Homens)2003em,anos(vidadeectativaexp
2003emidade
−
A) Determinar o polinômio interpolador utilizando todos os dados tabelados. Sugestão: ?interp
B) Verificar a validade do modelo encontrado.
C) Plotar num mesmo sistema de eixos os pontos tabelados e a e o polinômio interpolador determinado
no item b.
D)Estimar a expectativa de vida para pessoas com idade em 2003, variando de 14 a 30 anos.
> Estimativas_mulheres:= array( [ seq( [i, pexpvidam(i)], i =
14..30)]);
> Estimativas_homens:= array( [ seq( [i, pexpvidah(i)], i =
14..30)]);
:= Estimativas_mulheres
14 63.56097596
15 62.6000001
16 61.6414611
17 60.6831095
18 59.7236187
19 58.7625088
20 57.8000000
21 56.8368317
22 55.8740682
23 54.9129142
24 53.954550
25 53.000000
26 52.050032
27 51.105097
28 50.165308
29 49.230445
30 48.300001
:= Estimativas_homens
14 56.46550812
15 55.5000000
16 54.5559710
17 53.6353083
18 52.7375232
19 51.8602890
20 51.0000001
21 50.1523108
22 49.3126236
23 48.476507
24 47.640028
25 46.800001
26 45.954138
27 45.101122
28 44.240601
29 43.373125
30 42.500001
46
4-Interpolação Polinomial
4.7 Exercícios
(Fonte: Cálculo Numérico Computacional. Claudio, Dalcídio Moraes e Marins, Jussara M. São
Paulo: Ed.Atlas,1994.)
1. A tabela abaixo dá o volume de água num tanque elástico (usado para transpor-te de
óleo, leite, etc. em caminhões)para várias cotas de água. Determinar y(0,12).
x(m) 0,1 0,6 1,1 1,6 2,1
y )3(m 1,1052 1,8221 3,0042 4,9530 8,1662
2. Uma hidroelétrica tem capacidade máxima de 60Mw, a qual é determinada por três geradores de
respectivamente 30Mw, 15Mw e 15Mw.
A demanda de energia varia num ciclo de 24h e é uma função dela que o engenheiro ope-
racional distribui as tarefas dos geradores.
Sabe-se que a demanda mínima ocorre entre 1 e 5h e a demanda máxima entre 13 e 17h .
Pede-se para achar a partir dos dados abaixo essas demandas máximas e mínimas .
H 2 3 4 5 13 14 15 16 17
Demanda
(Mw)
16,4 15,2 14,9 16 28 36,5 43 34 31,2
3. Um paraquedista realizou seis salto; saltando de alturas distintas em cada salto, foi
testada a precisão de seus salto em relação a um alvo de raio de 5m, de acordo com a altura. A dis-
tância apresentada na tabela abaixo é relativa à circunferência.
ALTURA (m) DISTÂNCIA DO ALVO
1O SALTO (1500)
2O SALTO (1250)
3O SALTO (1000)
4O SALTO (750)
5O SALTO (500)
35
25
15
10
7
Levando em consideração os dados acima, a que provável distância do alvo cairia o paraquedista se
ele saltasse de uma altura de 850m ?
4. Um veículo de fabricação nacional, após vários testes, apresentou os resultados abaixo, quando
se analisou o consumo de combustível de acordo com a velocidade média imposta ao veículo. Os
testes foram realizados em rodovia em operação normal de tráfego, numa distância de 72 km. Veri-
ficar o consumo aproximado para o caso de ser desenvolvida a velocidade de 80km/h.
VELOCIDADE (km/h) CONSUMO (km/l)
55 14,08
70 13,56
85 13,28
100 12,27
120 11,3
140 10,4
47
E.B.Hauser – Cálculo Numérico
5. Numa esfera de superfície conhecida, o coeficiente de absorção 0,7 foi mantido à temperatura
de 6000 Ko . Foi calculada a energia irradiada de acordo com o tempo de irradiação, obedecendo à
tabela . Pede-se para obter a possível energia irradiada quando a irradiação atingir o tempo de 25
minutos.
ENERGIA IRRADIADA (Joules) TEMPO DE IRRADIAÇÃO (s)
71,72 . 310 600
94,72 . 310 800
118,4 . 310 1000
142,08 . 310 1200
165,76 . 310 1400
189,44 . 310 1600
6. Uma corda foi tensionada sob a ação de pesos distintos, quando para os respectivos
pesos foram calculadas as devidas velocidades de propagação que estão indicadas abaixo. Pede-se
para calcular a velocidade de propagação quando a corda está tensionada sob a ação de um peso de
7250 gf.
PESO (gf) VELOCIDADE (cm/s)
6000 13728,13
6500 14288,69
7000 14828,07
7500 15348,51
8000 15851,87
Respostas
1) 126904937,1)12,0(p ≅
2) A demanda mínima é 8632,14 Mw. e ocorre entre 3h e 4h da manhã.
A demanda máxima é 101,43 Mw. e ocorre entre 14h e 15h.
3) 4128,11)850(p =
4)
4512685,13)80(p
46701783,13)80(p
=
=
5)
5937,177618)1500(p =
6) p(7250) = 15090,53234
48
5. Ajuste de Funções
Aplicação1: Os dados tabelados descrevem a intensidade da luz como uma função da
distância da fonte, I(d), medida num experimento.
Determinar I (d) ≅ Y (d) =
CBdAd
1
2 ++
.
15.0
75
18.0
70
21.0
65
24.0
60
28.0
55
34.0
50
42.0
45
52.0
40
67.0
35
85.0
30
I
d
Aplicação 2:Segundo a lei de Boyle, o volume de um gás é inversamente proporcional à
pressão exercida, mantendo-se constante a temperatura. Para um certo gás, foram observados
os seguintes valores:
04,1
5,3
15.1
0,3
28,1
5,2
43,1
0,2
60,1
5,1
91,1
0,1
27,2
5,0
85,2
0,0
Volume
essãoPr
Ajustar os dados tabelados a uma hipérbole do tipo: Vol(p) ≅
pB
A)p(Y
+
=
49
E.B.Hauser – Cálculo Numérico
Aplicação 3 : Crescimento Populacional Mundial
908,6
60
2010
089,6
50
2000
276,5
40
1990
454,4
30
1980
708,3
20
1970
04,3
10
1960
555,2
0
1950
)bilhões(MundialPopulação
x
Ano
Ajustar os dados modelo logístico
BxAe1
12)x(Y
+
= .
Obs: 12 bilhões é a capacidade suporte sugerida em :
CHAPRA, Steven C., CANALE, Raymond P. Métodos Numéricos Para Engenhari,. São Paulo: : McGraw-Hill, 2008.
Aplicação 4 : Radiação Solar – Modelo Sinusoidal
A radiação solar, numa determinada cidade, foi medida durante um ano e tabulada como
131
D
159
N
211
0
260
S
287
A
308
J
359
J
351
M
311
A
245
M
188
F
144
J
)2^m/W(solarRadiação
)2015(mes
Supondo que cada mês dure 30 dias, utilizar a técnica dos mínimos quadrados para ajustar os
dados ao modelo sinusoidal e estimar a radiação no mês de agosto do ano subsequente.
)x0(sen.1B)x0cos(.1A0A)x(Y ωω ++= , T
pi
ω
2
0 = , T=360
50
5 - Ajuste de Funções
5.1 Introdução
O ajustamento é uma técnica de aproximação. Conhecendo-se dados experimen-
tais , )y,x( ii , n,...,1,0i = , deseja-se obter a lei )x(fy = relacionando x com y.
Devido aos erros experimentais nos n+1 pares, ))x(f,x( ii , teremos em ge-
ral 11 )x(f ε+ ; 22 )x(f ε+ ; ... ; nn )x(f ε+ , isto é , é impossível calcular exatamente a
função f(x). Por isso , em vez de procurarmos a função f tal que passa por cada um dos pontos
experimentais, determinaremos a função que melhor se ajusta aos pontos dados. O ajusta-
mento traduz um comportamento médio.
Para ajustar uma tabela de dados a uma função devemos:
• Conhecer a natureza física do problema ;
• Determinar o tipo de curva a que se ajustam os valores tabulados (graficamente e/ou
cálculo das diferenças finitas ou divididas) ;
• Calcular os parâmetros da curva.
5.2. Escolha da Função de Ajuste
Plotando os pontos dados ou usando diferenças(finitas ou divididas):
a) Função Linear (regressão linear) : xaa)x(Y 10 += , se iy∆ ≅ cte ou iy∆ ≅ cte .
b) Parábola (ajuste quadrático): 2210 xaxaa)x(Y ++= , se iy2∆ ≅ cte ou .2 cteiy ≅∆
c) Polinômio de grau p: se ip y∆ ≅ cte ou .cteyi
p
≅∆
d) Função Exponencial: xab)x(Y = , se .cte
x
ylog
i
i ≅
∆
∆
d) FunçãoPotência: pax)x(Y = , cte
xlog
ylog
i
i ≅
∆
∆
Outros tipos de Funções Ajuste:
• xaa
y
1
xaa
1)x(Y 10
10
+=⇒
+
= ; ( ) ( ) .cte
x
y/1
y/1
i
i
i ≅= ∆
∆∆
• ;xaa
y
x
xaa
x)x(Y 10
10
+=⇒
+
= ( ) ( )( ) .ctex
y/x
y/x
i
ii
ii ≅= ∆
∆∆
51
E.B.Hauser – Cálculo Numérico
•
2
2102
210
xaxaa
Y
1
xaxaa
1)x(Y ++=⇒
++
= ; ( ) ( ) ( ) .cte
xx
y/1y/1
y/1
i2i
i1i
i
2
≅
−
−
=
+
+ ∆∆∆
•
2cxbxae)x(Y += ; 2cxbxalnyln ++= , .cte
xx
ylnyln
i2i
i1i ≅
−
−
+
+ ∆∆
Exemplo: Considerando a tabela abaixo, como i,1yi
3 ∀=∆ , é adequado o ajuste dos dados
abaixo tabelados a um polinômio de grau 3.
i xi f(xi)=yi iy∆ iy2∆ iy3∆ iy4∆
0 0 0 4 5 1 0
1 2 8 19 10 1 -----
2 3 27 49 14 ----- -----
3 5 125 91 ----- ----- -----
4 6 216----- ----- ----- -----
52
5 - Ajuste de Funções
5.3 Determinação dos Parâmetros da Função de Ajuste: Critério dos Mínimos Quadra-
dos
Para determinar uma função, )x(fY = , que melhor se ajusta ao conjunto de dados
)y,x( ii , n,...,1,0i = , definimos o resíduo(erro) como sendo a diferença entre o valor de
)ix(f da função de ajuste e o valor tabulado de iy :
iy)ix(fi −=δ
Observamos que
∑
=
−=∑
=
n
0i
)iy)ix(f(i
n
0i
δ
pode ser tanto positiva, negativa ou nula e que uma soma de erros nula não significa que o
ajuste realizado seja adequado. Uma maneira de evitar isso é utlizar a soma do quadrado dos
resíduos.
O critério dos mínimos quadrados (CMQ) a soma do quadrado dos resíduos deve ser
mínima, isto é,
∑
=
−=∑
=
n
0i
2)iy)ix(f(
n
0i
2
iδ .
Considerando que o modelo procurado, )x(fY = , depende de parâmetros
m,...,1k,kc = , isto é , )kc,...,2c,1c,x(fy = , temos que
∑
=
−=∑
=
=
n
0i
2)iy)kc,...,2c,1c,ix(f(
n
0i
2
i)kc,...,2c,1c(F δ assume valor mínimos nos
seus pontos críticos. Isso exige que
.0
kc
)ix(fn
0i
]iy)ix(f[
kc
F
=
∂
∂
∑
=
−=
∂
∂
Seja pxpa...2x2ax1a0a)x(fY ++++== a função de ajustamento. Dada uma ta-
bela com n+1 pontos ( )iyix , , chamamos resíduo a diferença entre o valor de iY da equação de
ajustamento e o valor tabulado de iy . n,...,1,0i,yY iii ==− δ .
53
E.B.Hauser – Cálculo Numérico
O critério dos mínimos quadrados estabelece que: ∑
=
n
0i
2
iδ .’
Seja ( )∑
=
−=
n
i
iyiYpaaaF
0
2),...1,0( . Para F ter valor mínimo , é preciso que
0
a
F
0
=
∂
∂
; 0
a
F
1
=
∂
∂
; . . . ; 0=
∂
∂
pa
F
5.31 -Ajuste Linear
A função de ajuste terá a forma xaa)x(Y 10 += . Pelo critério dos Mínimos Quadra-
dos é necessário que :
∑ ∑
= =
−+=−=
n
i
n
i
iyixaaiyiYF
0 0
2)10(2)(
deve ser mínimo .
Sendo F uma função de duas variáveis, 0a e 1a , o menor valor de F será obtido
através de : 0
0
=
a
F
δ
δ
; 0
1
=
a
F
δ
δ
e assim:
0
0
=
a
F
δ
δ
, ( )∑
=
=−+
n
i
iyixaa
0
0102 e
0
1
=
a
F
δ
δ
, ( )∑
=
=−+
n
i
ixiyixaa
0
0102 .
Construímos o sistema de duas equações e duas variáveis:
( )
( )
=
=−+
=
=−+
∑
∑
n
i
ixiyxaa
n
i
iyixaa
0
01102
0
0102
ou
= = =
+=
= = =
+=
∑ ∑ ∑
∑ ∑ ∑
n
i
n
i
n
i
ixaixaiyix
n
i
n
i
n
i
ixaaiy
0 0 0
2
10
0 0 0
10
.
Obtemos:
54
5 - Ajuste de Funções
( )
= = =
=+
= =
=++
∑ ∑ ∑
∑ ∑
n
i
n
i
n
i
iyixixaixa
n
i
n
i
iyixaan
0 0 0
2
10
0 0
101
Resolvendo-se este último sistema linear são obtidos os valores de 0a e 1a e
assim determina-se a função de ajuste : xaaY 10 += .
Exemplo2 : Dada a tabela , achar a equação da reta que se ajusta usando o método dos Míni-
mos Quadrados.
i ix iy iyix 2ix i
Y =Y( ix ) ( )2iyiY −
0 0 2 0 0 2.07142857 0.005102041
1 1 3 3 1 3.14285714 0.02040816
2 2 5 10 4 4.21428571 0.61734694
3 3 5 15 9 5.2857143 0.08163266
4 4 5.5 22 16 6.3571429 0.73469755
5 5 8 40 25 7.4285714 0.32653061
∑ 15 28.5 90 55 --------------- 1.78571440
Seja xaaxY 10)( += a função que ajusta os dados .
Os parâmetros 0a e 1a constituem a solução do sistema :
∑ ∑ ∑
∑ ∑
+=
++=
2
10
10)1(
ixaixaiyix
ixaaniy
⇒
=+
=+
90155015
5.2811506
aa
aa
⇒
=
=
071428571.11
071428572.20
a
a
Logo, xY 071428571.1071428571.2 += .
55
E.B.Hauser – Cálculo Numérico
5.3.2 - Ajuste Quadrático
Seja 2210 xaxaaY ++= a função de ajuste. Assim,
( ) ( )∑ ∑ −++=−=
= =
n
0i
n
0i
2
i
2
i2i10
2
ii210 yxaxaayY)a,a,a(F
assume valor mínimo se 0
a
F
a
F
a
F
210
=== δ
δ
δ
δ
δ
δ
.
Os parâmetros 210 aea,a são obtidos resolvendo:
( )
++=
++=
+++=
∑ ∑ ∑∑
∑ ∑ ∑ ∑
∑ ∑∑
4
ix2a
3
ix1a
2
ix0aiy
2
ix
3
ix2a
2
ix1aix0aiyix
2
ix2aix1a0a1niy
Exemplo3 : Determinar polinômio de 2o grau que se ajusta o sdados tabelados.
i ix iy ii yx 2ix iyix
2
3
ix
4
ix i
y∆
iy
2∆
0 -2 -0.01 0.02 4 -0.04 -8 16 0.52 -0.21
1 -1 0.51 -0.51 1 0.51 -1 1 0.31 -0.25
2 0 0.82 0 0 0 0 0 0.06 -0.13
3 1 0.88 0.88 1 0.88 1 1 -0.07 -0.25
4 2 0.81 1.62 4 3.24 8 16 -0.32 ----------
5 3 0.49 1.47 9 4.41 27 81 ---------- ----------
∑ 3 3.5 3.48 19 9 27 115 ---------- ----------
=++
=++
=++
92115127019
48.322711903
5.32191306
aaa
aaa
aaa
⇒ 806285714.0201.02102142857.0)( ++−= xxxY
56
5 - Ajuste de Funções
5.3.3-Ajuste a polinômio de Grau 3
Considerando como função de ajuste o polinômio 3x3a2x2ax1a0aY +++= , o cri-
tério dos Mínimos Quadrados exige que
( )∑
=
∑
=
−+++=−=
n
0i
n
0i
2
iy
3x3a
2x2ax1a0a2iyiY)3a,2a,1a,0a(F
deve assumir valor mínimo. Isso acontece se 0
3a
F
2a
F
1a
F
0a
F
==== δ
δ
δ
δ
δ
δ
δ
δ
.
Assim, os parâmetros 3ae2a,1a,0a são obtidos resolvendo o sistema linear
( )
∑
=
∑
=
∑
=
=∑
=
+∑
=
++
∑
=
∑
=
∑
=
=∑
=
+∑
=
++
∑
=
∑
=
∑
=
=∑
=
+∑
=
++
∑
=
∑
=
=∑
=
+∑
=
+++
n
0i
n
0i
n
0i
iy3ix
n
0i
6ix3a
n
0i
5ix2a
4
ix1a
3ix0a
n
0i
n
0i
n
0i
iy
2
ix
n
0i
4
ix3a
n
0i
4
ix2a
3
ix1a
2
ix0a
n
0i
n
0i
n
0i
iyix
n
0i
4
ix3a
n
0i
3
ix2a
2
ix1aix0a
n
0i
n
0i
iy
n
0i
3ix3a
n
0i
2ix2aix1a0a1n
Exemplo5 – Utilizando o sistema Maple
> xv:=[0,2,3,5,6,8];
> yv:=[p(0),p(2),p(3),p(5),p(6),p(8)];
>g:=fit[leastsquare[[x,y],y=a*x^3+b*x^2+c*x+d,{a,b,c,d}]]([xv,yv]);
> gll :=unapply(rhs(g),x):
> plots[display]({plots[pointplot]
([0,p(0),2,p(2),3,p(3),5,p(5),6,p(6),8,p(8)]),
plot(gll(x), x=-.5..8.5, y=-50..280),
plot(gll(x), x=-.5..8.5,y=-150..280,thickness=2)});
57
E.B.Hauser – Cálculo Numérico
5.3.4-Ajuste a polinômio de Grau p
O ajuste a um polinômio de grau p. np,pxpa...2x2ax1a0aY <++++= , exige
resolver o sistema linear:
∑ ∑ ∑ ∑
++
∑
+
∑ ∑ ∑ ∑
+
∑
∑ ∑ ∑ ∑ ∑
+
∑ ∑∑∑+
p2
ix
3p
ix
2p
ix
1p
ix
p
ix
2p
ix...
5
ix
4
ix
3
ix
2
ix
1p
ix...
4
ix
3
ix
2
ixix
p
ix...
3
ix
2
ixix)1n(
L
MOMMMM
−
pa
1pa
2a
1a
0a
M
=
∑
∑
−
∑
∑
∑
iy
p
ix
iy
1p
ix
iy
2
ix
iyix
iy
M
5.3.4 -AJUSTE NÃO LINEAR NOS PARÂMETROS: CASOS REDUTÍVEIS AO
LINEAR OU PARABÓLICO POR MUDANÇA DE VARIÁVEIS
5.3.4.1 – Ajuste por Função Exonencial
Seja xabY = a função de ajuste.
Linearizamos Y, aplicando log ou ln (escolher a base adequada)
{ { {
1xa0ay
blnxalnYln +=
Determinamos a e b resolvendo:
=
=
=
+
=
=
=
=
++
∑∑∑
∑∑
iyln
n
0i
ix
n
0i
2
ixbln
n
0i
ixaln
n
0i
iyln
n
0i
ixblnaln)1n(
Então
=→=
=→=
1a1
0a0
ebblna
eaalna
Exemplo 4: Ajustar os dados abaixo a uma função exponencial do tipo xabY = .
58
5 - Ajuste de Funções
i ix iy iyln iyix ln 2ix i
Y 2)iy)ix(Y( −
0 0 3 1.09861 0 0 2.98422 0.0002490
1 0.5 4 1.38629 0.693145 0.25 4.2228031 0.0496412
2 1 6 1.79176 1.79176 1 5.9754291 0.0006025
3 1.5 9 2.19722 3.29583 2.25 8.4555298 0.2964477
4 2 12 2.48491 4.96982 4 11.964948 0.0012286
5 2.5 17 2.83321 7.08303 6.25 16.930930 0.0047706
6 3 24 3.17805 9.53415 9 23.958013 0.0017628
7 3.5 33 3.49651 12.2378 12.25 33.901647 0.8129689
8 4 48 3.87120 15.4848 16 47.972329 0.0007656
∑ 18 ------------- 22.3378 55.0903 51 ------------- 1.1684403
Os gráficos a seguir ilustram o efeito da linearização dos dados.
=+
=+
0903.55151018
3378.2211809
aa
aa
⇒
984218125.20093337778.11
002347015.2169432.00
==→=
==→=
a
eaa
a
eba
⇒
xY )00235.2(98422.2=
-
59
E.B.Hauser – Cálculo Numérico
5.3.4.2 -AJUSTE POR FUNÇÃO POTÊNCIA
Seja baxY = a função de ajuste.
Linearizando a função Y(x) , temos : .xlnbalnyln +=
Resolve-se o sistema de equações lineares e encontra-se a e b:
1
1
0
1
xxln
ab
aaln
yyln
=
=
=
=
⇒
=
=
=
+
=
=
=
=
++
∑∑∑
∑∑
n
0i
iylnixln
n
0i
2)ix(lnb
n
0i
ixlnaln
n
0i
iyln
n
0i
ixlnbaln)1n(
então
10
a
0 abeeaaaln ==⇒= .
Exercício : Os dados abaixo dão a duração de uma broca em função da velocidade de corte.
Pede-se para fazer uma tabela de D=D(v) para 100 (10) 180.
8.2
180
9.7
150
28
120
79
100
.)(
)/(
segD
smV
D(v)=0.813947103e14*1/(x^5.680591334)
-
Efeito da linearização dos dadoss
5.3.4.3 -AJUSTE POR FUNÇÃO HIPERBÓLICA :
xaa
1Y
10 +
= . Linearizando, temos :
∑=∑+∑
∑ ∑=++
⇒+=
===
= =
n
0i
ii
n
0i
2
i1
n
0i
0i
n
0i
n
0i
ii10
10
y/xxaax
y/1xaa)1n(
xaa
Y
1
60
5 - Ajuste de Funções
MODELOS SIGMOIDAIS
5.3.4.4 -AJUSTE POR FUNÇÃO LOGÍSTICA :
xBeA
CY
+
=
1
, onde C, a capacidade de
suporte, é um valor conhecido.
Linearizando, temos :
( )( )
( )( )
−=+
−=++
⇒+=
−
∑∑∑
∑ ∑
===
= =
n
i
ii
n
i
i
n
i
i
n
i
n
i
ii
yCxxBxA
yCxBAn
xBA
Y
C
00
2
0
0 0
1/lnln
1/lnln)1(
ln1ln
Exemplo: Ajustar os dados (0 , 0.0666) e (2, 0.796) a função logística do
tipo
xbea1
1Y
+
= .
Resposta:
−=+
=+
722,2b4aln2
278,1b2aln2
x2e77.141
1)x(Y
−+
≅
5.3.4.4 -AJUSTE POR FUNÇÃO DE GOMPERTZ :
xBeAe
K)x(Y
−
= , onde K, a capaci-
dade de suporte, é um valor conhecido.
Linearizando, temos :
( )( )
( )( )
∑
=
=∑
=
−∑
=
∑
=
∑
=
=−+
⇒−=
n
i
iyKix
n
i
ixB
n
i
ixA
n
i
n
i
iyKixBAn
xBA
Y
K
0
/lnln
0
2
0
ln
0 0
/lnlnln)1(
lnlnln
Exemplo: Ajustar os dados (0.5 , 1.3) e (2 , 3.8) utilizando o modelo de Gompertz .
Resposta:
−=−
−=−
881,525.4ln5.2
853.25.2ln2
ba
ba
xee
xY 058,2145.3
4)(
−
≅
61
E.B.Hauser – Cálculo Numérico
5.3.4.5 – AJUSTE POR FUNÇÃO PERIÓDICA-Modelo Sinusoidal
)x0(sen.1B)x0cos(.1A0A)x(Y ωω ++= , T
pi
ω
2
0 = , T=(n+1)h ix1ixh −+= i=0,1,2..n.
=
+
∑
∑
∑
∑∑∑
∑∑∑
∑∑
)(.
)cos(.
)()(.)cos()(
)(.)cos()(cos)cos(
)()cos(1
0
0
1
1
0
2
000
000
2
0
00
ii
ii
i
iiii
iii
ii
xseny
xy
y
B
A
Ao
xsenxsenxxsen
xsentxx
xsenxn
ω
ω
ωωωω
ωωωω
ωω
+
+
+
2
1n00
0
2
1n0
001n
1
1
0
B
A
A
= ⇒
∑
∑
∑
)(.
)cos(.
0
0
ii
ii
i
xseny
xy
y
ω
ω
∑
∑
∑
+
=
+
=
+
=
)ix0(sen.iy.1n
2
1B
)ix0cos(.iy.1n
2
1A
1n
iy
0A
ω
ω
0
1
)cos( 0
=
+
∑
n
xiω
, 0
1
)( 0
=
+
∑
n
xsen iω
,
2
1
1
)(cos 02
=
+
∑
n
xiω
,
2
1
1
)( 02
=
+
∑
n
xsen iω
,
0
1
)(.)cos( 00
=
+
∑
n
xsenx ii ωω
Exemplo: T=1.5, n+1=10, h=0,15
i ix iy cos(. 0 ii xy ω (. 0 ii xseny ω
0 0 2,200 2,200 0,0000
1 0.15 1,595 1,291 0,938
2 0.30 1,0311 0,319 0,980
3 0.45 0,722 -0,223 0,687
4 0.60 0,786 -0,636 0,462
5 0.75 1,200 -1,200 0,000
6 0.90 1,805 -1,460 -1,061
7 1.05 2,369 -0,732 -2,253
8 1.20 2,678 0,829 -2,547
9 1.35 2,614 2,114 -1,536
∑ 17,00 2,502 -4,330
62
5 - Ajuste de Funções
)x.0(sen.866,0)x.0cos(.500,07,1)x(Y
,Logo
866,0)330,4.(
10
2)ix0(sen.iy.1n
2
1B
500,0502,2.
10
2)ix0cos(.y.1n
2
1A
7,1
10
17
1n
iy
0A
ωω
ω
ω
−+=
−=−=
+
=
==
+
=
==
+
=
∑
∑
∑
63
E.B.Hauser – Cálculo Numérico
5.3.4.6
-OUTROS TIPOS DE FUNÇÕES DE AJUSTE
1)
xaa
xY
10 +
= . Linearizando, temos:
∑=∑+∑
∑=∑++
⇒+=
===
==
n
0i
i
2
i
n
0i
2
i1
n
0i
0i
n
0i
ii
n
0i
i10
10
y/xxaax
y/xxaa)1n(
xaa
Y
x
2) 2
210 xaxaa
1Y
++
= . Linearizando, temos :
∑ ∑ ∑ ∑=++
∑ ∑ ∑=+∑+
∑=∑+∑++⇒++=
i
2
i
4
i2
3
i1
2
i0
ii
3
i2
2
i1i0
i
2
i2i10
2
210
y/xxaxaxa
y/xxaxaxa
y/1xaxaa)1n(
xaxaa
Y
1
3)
{ { 43421
2
x2ax1a
2
0ay
2cxbx
cxbxalnYln
aeY
+
+
++=
=
⇒
∑ ∑ ∑=+∑+
∑ ∑=∑+∑+
∑ ∑ ∑=+++
i
2
i
4
i
3
i
2
i
ii
3
i
2
ii
i
2
ii
ylnxxcxbxaln
ylnxxcxbxaln
ylnxcxbaln)1n(
.
64
5 - Ajuste de Funções
Exercícios:
1) Considerando:
i 0 1 2 3 4 5 6
xi 1 2 3 4 5 6 7
yi 34 45 63 88 120 159 205
a) Mostrar que o ajuste por uma parábola é adequado.
b) Ajustar os dados a uma parábola.
Resposta: a) ∆2 7y ii = ∀,
b) 30x5,02x5,3)x(p ++=
2) Segundo o critério dos Mínimos Quadrados, qual das duas funções x558,0332,2Y1 +=
e x235,02 e037,2Y = melhor ajusta os dados da tabela?
x -2,3 1 3,1
y 1,2 2,5 4,2
Resposta: Y2, pois 2ii1
2
ii2 )yY()yY( −<− ∑∑
2
ii1 )yY( −∑ = 0,194 e 2ii2 )yY( −∑ = 0,0123
3) Um filme vem sendo exibido numa determinada
sala de cinema por cinco semanas e a frequência se-
manal, (aproximada à centena mais próxima) está
dada na tabela abaixo. Utilizar ajuste linear para de-
terminar a frequência esperada na sexta semana.
semana 1 2 3
4 5
frequência 5000 4500 4100 3900 3500
Resposta: Y(x) = 5280 -360x e y(6) = 3120
y = -360x + 5280
R20,9818 =
0
1000
2000
3000
4000
5000
6000
0 1 2 3 4 5 6
65
E.B.Hauser – Cálculo Numérico
4) O número de bactérias, por unidade de volume, existente em uma cultura depois de x ho-
ras é apresentado pela tabela. Ajuste os dados tabulados a uma curva exponencial da forma y
=abx e avaliar y para x=7.
x 0 1 2 3 4 5 6
y 32 47 65 92 132 190 275
Resposta: x42694,11483,32)x(Y ×= e Y (7) = 387, 256
5) Utilizando o critério dos Mínimos Quadrados, ajustar a uma reta os dados tabulados:
xi 3 5 6 8 9 11
yi 2 3 4 6 5 8
Resposta: )x(Y = -0,33328 + 0,714285x
6) A tabela abaixo fornece uma relação entre a temperatura da água e a pressão barométrica.
Ajustar os dados a uma função potência.
p(mm de Hg) 680 690 700 710 720 730 740 780
T( Co ) 96.92 97.32 97.71 98.11 98.49 98.88 99.26 100.73
7) A tabela abaixo fornece uma relação entre a resistência à tração do aço em função da
temperatura.
1500
617
2780
485
4450
412
5260
330
5720
250
)2cm/kg(
)Co(t
σ
Ajustar os dados a um polinômio
de grau 4.
y = 15,476x0,2813
R21 =
96,5
97
97,5
98
98,5
99
99,5
100
100,5
101
660 680 700 720 740 760 780 800
y = 2E-06x40,0033x - 31,8931x + 2466,48x + 47846 -
R21 =
0
1000
2000
3000
4000
5000
6000
7000
0 200 400 600 800
66
5 - Ajuste de Funções
8) Os dados abaixo referem-se a variação do coeficiente de atrito (µ ) ,entre a roda e o tri-
lho seco, com a velocidade.
154.0
70
164.0
60
192.0
40
215.0
30
250.0
20
313.0
10
450.0
0)h/km(V
µ
Ajustar os dados a um polinômio de grau 5.
9) Os dados abaixo relacionam a viscosidade
η em função da temperatura t.
990.0
21
069.1
18
121.1
16
148.1
15
175.1
14
276.1
9.10
409.1
5.7Cot
η
Realizar o ajuste sugerido pelo gráfico ao lado.
Resp: )t(η =-0.030772617131t+1.619873714
y = -9E-10x52E-07x + 42E-05x - 30,0007x + 2 -
0,0196x + 0,45
R21 = 0,000
0,050
0,100
0,150
0,200
0,250
0,300
0,350
0,400
0,450
0,500
0 20 40 60 80
67
E.B.Hauser – Cálculo Numérico
10) Verificar qual das funções
x0001523.0015619.0)x(1Y −=
ou
x)98028.0(017054.0)x(2Y =
melhor se ajusta à tabela dada .
ix 20 40 60 80 100
iy 0.01310 0.006540.004590.00341 0.00263
Respostas:
Aplicação1: Y(d)= 1/(.1329810498e-2*d^2-.2080049421e-1*d+.6161059331)
Aplicação2: Y(p) = 5.779411/(2.043906+p)
Aplicação 3:
Aplicação 4: )x
360
2
cos(22.0)x
360
2
cos(97.10517.246)x(Y pipi −−=
Estimativa para Agosto do ano subsequênte: Y(600) 17,299)600(Y ≅
68
6 . Integração Numérica
Objetivo : Calcular a ∫
b
a
dxxf )( , onde a função integrando )(xf ou é conhecida por sua expressão
analítica ou por uma tabela de valores ( ))x(f,x ii , i = 0, 1,..., n.
Aplicação: Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada, reproduzindo
os dados:
1.75
17
6.78
16
1.81
15
4.86
14
5.86
13
8.84
12
1.83
11
0.77
10
3.75
9
)atemperatur(y
)hora(x
Encontrar a temperatura média da água entre 9h da manhã e 5h da tarde e estimar o erro cometido
nesse cálculo.
]b,a[emcontínuaéf
e0)b(f).a(fse],b,a[em)x(fdemédiovaloroédx)x(f
ab
1fm:OBS
b
a
>∫
−
=
6.1 -Fórmulas de Newton-Cotes
Consideremos ( ))x(f,x ii , i = 0, 1,..., n., hixix =−+1 ,
n
xnxh 0−= , )( ixfiy = .
A integral da função f(x) no intervalo [ ]n0 x;x é dada por :
∫
−−
++
−−
+
−
+=
=∫≅∫
−
nx
0x
0
n
n
1n0
0
2
2
10
0
0
0
nx
0x
n
nx
0x
dxy
h!n
)xx)...(xx(
...y
h!2
)xx)(xx(
y
h
)xx(
y
dx)x(Pdx)x(f
∆∆∆
Se
h
xxR 0−= , então
hdrdxRhxx
0Rxx
0
0
=→+=
=→=
e nRnxx =→= .
Assim:
69
E.B.Hauser – Cálculo Numérico
dR)y
!n
)1nR)...(1R(R
...y
!2
)1R(RyRy(hdR)R(Phdx)x(f 0nn
x
0x
0
2
0
n
0 0
n
0
∆∆∆∫ ∫∫
+−−
++
−
++=≅
=
+−−++−++∫ ∫ ∫∫
n
0
n
0
n
0
0
n
0
2
n
000
dR)1nR)...(1R(R
!n
y
...dR)1R(R
!2
yRdRydRyh ∆∆∆
Na prática não é usual aproximar f(x) por um polinômio de grau n (elevado) devido
ao erro de arredondamento que ocorre no processo.
6.2 – Regra dos Trapézios
Considerando n = 1 na fórmula de Newton-Cotes temos :
∫ ∫ ∫ ∫
+=≅1
x
0x
1
0
1
0
1
000
RdRydRyhdR)R(Phdx)x(f ∆ [ ]1020100 ]2/R[y]R[yh ∆+= =
= [ ]1000 yy2
h
2
y
yh +=
+
∆
, ou para o intervalo [ ]1ii x;x + :
[ ]∫ ∫+ + ++=≅1ixix
1i
i 1ii
yy
2
hdR)R(Phdx)x(f
Generalizamos para n subintervalos:
[ ]n1n21onx
ox
y)yyy(2y
2
hdx)x(f +++++≅∫
−
L
Erro de Truncamento(para n subintervalos):
)(''
],[
max)0(12
2
xf
nxoxx
xnx
h
TE
∈
−≤ ou iy
xnx
TE
2max
12
)0( ∆−≤
Vê-se que a fórmula dos Trapézios é exata para polinômios do 1o grau.
Exemplo1: Determinar h de tal forma que a regra trapezoidal forneça o valor de dxe1
0
2x
∫
− com um
erro de truncamento menor que 410− .
0245.0h10)2(
12
h)x(''fmáx)01(
12
hE 422
T <⇒<=−≤
−
[ ]1;0x∈
h/)xx(n 0n −= , n/)xx(h 0n −= , 0245.0n/)xx( 0n <− , 41n8.40n =→> , 02439.041
1h ==
70
6 -Integração Numérica
6.3 – Fórmula de Simpson
Fazendo n = 2 na fórmula de Newton-Cotes, temos,
[ ]∫ ∫ ∫∫ ++=
−++≅2
x
0x
210
2
0
2
0
0
22
000
yy4y
3
hdR)1R(R
!2
y
RdRydRyhdx)x(f ∆∆
Generalizamos para n subintervalos, n par:
[ ]∫ +++++++++++≅
−−
nx
ox
n2n6421n531o y)yyyy(2)yyyy(4y3
hdx)x(f LL
ERRO DE TRUNCAMENTO PARA A FÓRMULA DE SIMPSON
)x(fmax)xx(
180
hE ''''
]nx,ox[x
0n
4
S
∈
−≤ ou i
40n
S ymax180
)xx(E ∆−≤
Exercícios
1.Aplicação:.Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada, reprodu-
zindo os dados:
1.75
17
6.78
16
1.81
15
4.86
14
5.86
13
8.84
12
1.83
11
0.77
10
3.75
9
)atemperatur(y
)hora(x
Encontrar a temperatura média da água entre 9h da manhã e 5h da tarde e estimar o erro cometido
nesse cálculo.
2. Estimar a área da região hachurada pela regra dos Trapézios e pela de Simpson, usando oito su-
bintervalos.
71
E.B.Hauser – Cálculo Numérico
3. Calcular a integral abaixo pela regra dos Trapézios e pela de Simpson, usando quatro, seis e dez
subintervalos. Comparar os resultados.
xdxsen)xcos2sen(
2/
0
2
∫
pi
4. Calcular por Trapézios:
a) dx
xe
e
2
0
x
∫
−
−
com h = 0.5
b) tgxdxe
1
0
2
x
∫ com h=0.1
5. Calcular por Simpson:
a) ∫
++
2.1
0
x 1xe
dx
com h =0 .3
b) dxe
2
0
2
x
∫
−
com h = 0.25
c) xdxcosx
2/
0
2
∫
pi
com n=0.4
6. O gráfico da figura foi registrado por um instrumento usado para medir uma quantidade física.
Estimar as coordenadas-y dos pontos do gráfico e aproximar a área da região fechada usando seis
subintervalos
72
6 -Integração Numérica
7) A função de Bessel é solução de uma equação que
surge com grande freqüência, em engenharia e/ ou física
matemática, na resolução das equações diferencias parci-
ais pelo método da separação de variáveis. As equações
de Bessel surgem quando aplicamos a técnica de separa-
ção de variáveis a problemas de valor de contorno, por
exemplo, em coordenadas polares, cilíndricas e esféricas.
Constitui exemplos importantes desta modelagem o es-
tudo da evolução da temperatura e reações químicas em
cilindros e esferas. Ao lado a representação gráfica de
)t(0J
Tarefa: Considerar a representação integral da Função de Bessel de primeiro tipo
,...3,2,1,0m,
0
dx)tsenxmxcos(1)t(mJ =−= ∫
pi
pi
Estimar )3(0J com cinco subintervalos e o erro cometido neste cálculo.
Exercício 8: Numa sessão transversal(30pés de comprimento) de um veleiro de competição a for-
ça total do vento que atua mastro é expressa por
dx10
x
e
30
0 x6
x250F
−
∫
+
= .
Estimar o valor da integral F, utilizando o método de
Simpson com quatro subintervalos(n=4)
Exercício 9: Utilizando Trapézios com n=11, estimar a área do Brasil e comparar com o a
área conhecida(8.514.876,599 km²). Criar, com escala de 1cm( Km365≅ ), um sistema de
eixos com origem na fronteira comum do Pará, Tocantins e Mato Grosso.
73
E.B.Hauser – Cálculo Numérico
“O Brasil (oficialmente República Federativa do Brasil) é uma república federativa forma-
da pela união de 26 estados federados e pelo Distrito Federal. O país conta com 5.564 municí-
pios, 189.612.814 habitantes, bem como uma área de 8.514.876,599 km², equivalente a 47%
do território sul-americano. Em comparação com os demais países do globo, dispõe do quinto
maior contingente populacional e da quinta maior área. Nona maior economia do planeta e
maior economia latino-americana, o Brasil tem hoje forte influência internacional, seja em
âmbito regional ou global. O País possui entre 15 e 20% da biodiversidade mundial, sendo
exemplo desta riqueza a Floresta Amazônica, com 3,6 milhões de quilômetros quadrados. Faz
fronteira a norte com a Venezuela, com a Guiana, com o Suriname e com o departamento ul-
tramarino da Guiana Francesa; ao sul com o Uruguai; a sudoeste com a Argentina e com o
Paraguai; a oeste com a Bolívia e com o Peru e, por fim a noroeste com a Colômbia. Os úni-
cos países sul-americanos que não têm uma fronteira comum com o Brasil são o Chile e o
Equador. O país é banhado pelo oceano Atlântico ao longo de toda sua costa norte, nordeste,
sudeste e sul. Além do território continental, o Brasil também possui alguns grandes grupos de
ilhas no oceano Atlântico como exemplo: Penedos de São Pedro e São Paulo, Fernando de
Noronha (território especial do estado de Pernambuco) e Trindade e Martim Vaz no Espírito
Santo. Há também um complexo de pequenas ilhas e corais chamado Atol das Rocas (que
pertence ao estado do Rio Grande do Norte). Apesar de ser o quinto país mais populoso do
mundo, o Brasil apresenta uma das mais baixas densidades populacionais. A maior parte da
população se concentra ao longo do litoral, enquanto o interior do país ainda hoje é marcado
por enormes vazios demográficos.” Fonte: http://pt.wikipedia.org/wiki/Brasil
74
6 -Integração Numérica
Respostas:
1. temperatura média ≅ 81,58 oF ( por Trapézios)
2. Trapézios: 1.761237
Simpson: 1.763624
3.
n 4 6 10
Trapézio 0,481485 0,496396 0,503836
Simpson 0,512682 0,508646 0,508045
4. a) I = 0,636571
b) I = 1,110603
5. a) I = 0,658685
b) I = 0,882065
c) I = 0,466890
7. )3(0J = -0,260052
9. A ≈ 8160031,25 Km2 .( Comparando com a área real, a diferença é de 354845.349 Km2. A área
do estado do Mato Grosso do Sul é ≈ 357125 Km2
Valor da área A obtido considerando os pontos:
A1: (-7,0),(-6,1.5),(-5,3),(-4,4),(-3,4.5),(-2,3.75),(-1,3.5),(0,4),(1,2.5),(2,2),(3,1.5) e (4,1)
A2: (-7,0),(-6,-0.2.5),(-5,-0.5),(-4,-1),(-3,-1.5),(-2,-4),(-1,-7),(0,-5.5),(1,-4.25),(2,-4),(3,-2) e (4,-1)
75
7 - Resolução Numérica de Equações Diferenciais Ordinárias
Objetivo: Resolver numericamente(e generalizar para problemas de ordem mais elevada) o
problema de valor inicial de primeira ordem
(PVI)
=
=
oo y)x(y
)y,x(f'y
Os métodos numéricos são processos que fornecem valores aproximados da solução em
pontos particulares, usando operações algébricas.
Os métodos que estudaremos determinarão estimativas da solução nos pontos
K,x,x,x 21o , onde ihxx oi += , i= 0,1,...,n . A escolha do valor de h é arbitrária e, em geral,
quanto menor h , melhor a estimativa da solução obtida.
Sejam iyixy ≅)( , ihxx oi += , i= 0,1,...,n e n/)xx(h on −= .
Método de Euler(PVI): )y,x(hfyy iii1i +=+ , i= 0,1,...,n-1.
Exemplo- Lei da Radiação de Stefan: A lei de Radiação de Stefan afirma que a taxa de variação
da temperatura de um corpo a T(t) graus, em um meio a M(t) graus, é proporcional a 4T4M − .
Isto é , com K a constante de proporcionalidade,
−= )t(4T)t(4MK
dt
)t(dT
.
Considerar , 440K −= , a temperatura do meio constante o70)t(M = e o100)0(T = , para estimarT(2)
por Euler com h=0,5. Com a notação T=y, t=x, T(t)=y(x), aplicamos o método de Euler, com h=0,5
para estimar y(2), se
=
−
=
100)0(y
440
4y470
y´
.
Solução:
440
4ky470
5,0ky1ky
−
+=+
k xk )kx(yky ≅
0 0 100
1 0,5 85,1582
2 1 79,5761
3 1,5 76,4337
4 2 74,4571
76
E.B.Hauser – Cálculo Numérico
Método de Runge-Kutta de 2 a ordem(PVI) :
( )21i1i kk2
hyy ++=+ , )y,x(fk ii1 = e )hky,hx(fk 1ii2 ++= , i= 0,1,...,n-1.
Exemplo1: Utilizar o método de Runge-Kutta de segunda ordem com h=0,2 para estimar o va-
lor de y(0,8) se
=
−=
2)0(y
yx'y
.
Exemplo2: Utilizar o método de Runge-Kutta de segunda ordem com h=0,1 para estimar o va-
lor de y(0,3) se
=
+=
2)0(y
)yx(sen'y
.
Exemplo3: Utilizar o método de Runge-Kutta de segunda ordem com h=0,2 para estimar o va-
lor de y(1,4) se
=
=
4)1(y
yln2x'y
.
k xk )kx(yky ≅ K1 K2
0 0 2 -2 -1,4
1 0,2 1,66 -1,46 -0.968
2 0,4 1,4172 -1,0172 -0,6137
3 0,6 1,254 -0,6541 -0.3232
4 0,8 1,1563
k xk )kx(yky ≅ K1 K2
0 0 2 0,9092 0,8133
1 0,1 2,0861 0,8165 0,6988
2 0,2 2,1619 0,7030 0,5723
3 0,3 2,2256
k xk )kx(yky ≅ K1 K2
0 1 4 1,3862 2,0927
1 1,2 4,3470 2,1163 3,0626
2 1,4 4,8655
77
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Exemplos Comparando Solução Exata com Euler e Runge-Kutta de Segunda
Método de Euler
Problema: y' = y – x; y(0) = 2
xn
yn
Solução exata
Y(x) = ex + x + 1
h = 0.1 h = 0.05 h = 0.01 h = 0.005
0.0 2.0000 2.0000 2.0000 2.0000 2.0000
0.1 2.2000 2.2025 2.2046 2.2049 2.2052
0.2 2.4100 2.4155 2.4202 2.4208 2.4214
0.3 2.6310 2.6401 2.6478 2.6489 2.6499
0.4 2.8641 2.8775 2.8889 2.8903 2.8918
0.5 3.1105 3.1289 3.1446 3.1467 3.1487
0.6 3.3716 3.3959 3.4167 3.4194 3.4221
0.7 3.6487 3.6799 3.7068 3.7102 3.7138
0.8 3.9436 3.9829 4.0167 4.0211 4.0255
0.9 4.2579 4.3066 4.3486 4.3541 4.3596
1.0 4.5937 4.6533 4.7048 4.7115 4.7183
Método de Runge-Kutta de segunda ordem
Problema: y' = y – x; y(0) = 2
xn
yn Solução exata
Y(x) = ex + x + 1
h = 0.1 h = 0.05 h = 0.01
0.0 2.0000000 2.0000000 2.0000000 2.0000000
0.1 2.2050000 2.2051266 2.2051691 2.2051709
0.2 2.4210250 2.4213047 2.4213987 2.4214028
0.3 2.6492326 2.6496963 2.6498521 2.6498588
0.4 2.8909021 2.8915852 2.8918148 2.8918247
0.5 3.1474468 3.1483904 3.1487076 3.1487213
0.6 3.4204287 3.4216801 3.4221007 3.4221188
0.7 3.7115737 3.7131870 3.7137294 3.7137527
0.8 4.0227889 4.0248265 4.0255115 4.0255409
0.9 4.3561818 4.3587148 4.3595665 4.3596031
1.0 4.7140809 4.7171911 4.7182369 4.7182818
78
E.B.Hauser – Cálculo Numérico
Método de Euler
Problema: y' = y; y(0) = 1
xn
yn
Solução exata
Y(x) = ex
h = 0.1 h = 0.05 h = 0.01 h = 0.005
0.0 1.0000 1.0000 1.0000 1.0000 1.0000
0.1 1.1000 1.1025 1.1046 1.1049 1.1052
0.2 1.2100 1.2155 1.2202 1.2208 1.2214
0.3 1.3310 1.3401 1.3478 1.3489 1.3499
0.4 1.4641 1.4775 1.4889 1.4903 1.4918
0.5 1.6105 1.6289 1.6446 1.6467 1.6487
0.6 1.7716 1.7959 1.8167 1.8194 1.8221
0.7 1.9487 1.9799 2.0068 2.0102 2.0138
0.8 2.1436 2.1829 2.2167 2.2211 2.2255
0.9 2.3579 2.4066 2.4486 2.4541 2.4596
1.0 2.5937 2.6533 2.7048 2.7115 2.7183
Método de Runge-Kutta de segunda ordem
Problema: y' = y; y(0) = 1
xn
yn Solução exata
Y(x) = ex
h = 0.1 h = 0.05 h = 0.01
0.0 1.0000000 1.0000000 1.0000000 1.0000000
0.1 1.1050000 1.1051266 1.1051691 1.1051709
0.2 1.2210250 1.2213047 1.2213987 1.2214028
0.3 1.3492326 1.3496963 1.3498521 1.3498588
0.4 1.4909021 1.4915852 1.4918148 1.4918247
0.5 1.6474468 1.6483904 1.6487076 1.6487213
0.6 1.8204287 1.8216801 1.8221007 1.8221188
0.7 2.0115737 2.0131873 2.0137294 2.0137527
0.8 2.2227889 2.2248265 2.2255115 2.2255409
0.9 2.4561818 2.4587148 2.4595665 2.4596031
1.0 2.7140809 2.7171911 2.7182369 2.7182818
79
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Método de Euler
Problema: y' = 4x3; y(0) = 1
xn
yn
Solução exata
Y(x) = x4
h = 0.1 h = 0.05 h = 0.01
0.1 0.0000 0.0000 0.0000 0.0000
0.1 0.0000 0.0000 0.0001 0.0001
0.2 0.0004 0.0009 0.0014 0.0016
0.3 0.0036 0.0056 0.0076 0.0081
0.4 0.0144 0.0196 0.0243 0.0256
0.5 0.0400 0.0506 0.0600 0.0625
0.6 0.0900 0.1089 0.1253 0.1296
0.7 0.1764 0.2070 0.2333 0.2401
0.8 0.3136 0.3600 0.3994 0.4096
0.9 0.5184 0.5852 0.6416 0.6561
1.0 0.8100 0.9025 0.9801 1.0000
Método de Runge-Kutta de segunda ordem
Problema: y' = 4x3; y(0) = 0
xn yn
Solução exata
Y(x) = x4
h = 0.1 h = 0.05 h = 0.01
0.0 0.0000000 0.0000000 0.0000000 0.0000000
0.1 0.0002000 0.0001250 0.0001010 0.0001000
0.2 0.0020000 0.0017000 0.0016040 0.0016000
0.3 0.0090000 0.0083250 0.0081090 0.0081000
0.4 0.0272000 0.0260000 0.0256160 0.0256000
0.5 0.0650000 0.0631250 0.0625250 0.0625000
0.6 0.1332000 0.1305000 0.1296360 0.1296000
0.7 0.2450000 0.2413250 0.2401490 0.2401000
0.8 0.4160000 0.4112000 0.4096640 0.4096000
0.9 0.6642000 0.6581250 0.6561810 0.6561000
1.0 1.0100000 1.0025000 1.0001000 1.0000000
80
E.B.Hauser – Cálculo Numérico
PVI de ordem maior que um
Sistemas de Equações Diferenciais de primeira ordem com condições iniciais
==
=
=
oooo z)x(z,y)x(y
)z,y,x(g'z
)z,y,x(f'y
Para obtermos sua solução é possível aplicar os métodos de Euler e Runge-Kutta de
segunda ordem. Por exemplo, por Euler as estimativas são obtidas aplicando.
yi+1 = yi + h f (xi, yi, zi ) e zi+1 = zi + h g (xi, yi, zi ).
Mudança de Variável para problemas de valor inicial de segunda ordem:
Para
==
=
oooo z)x('y,y)x(y
)z,y,x(f''y
, fazemos a mudança de varável u'y = , com oo y)x(y = .
Então devemos resolver o sistema de duas equações diferenciais ordinárias, acopladas pelas
condições iniciais:
==
=
=
oooo z)x(u,y)x(y
)z,y,x(f'u
u'y
.
Exemplo1: Considerar um sistema massa-mola-amortecedor descrito pela equação diferencial
ordinária de segunda ordem: m y” + cy’ + ky = L sin x.
Utilizando Euler com h=0,1, estimar y(0,3) ( isto é, o deslocamento para o tempo x=0,3 , para
massa m=1, amortecedor c=0,5, rigidez k=2, amplitude da força L=0,5 , com deslocamento ini-
cial y(0)=1, velocidade inicial y’(0) =0.
Solução:
( )ky2kw5,0)kx(sen5,01,0kw1kw,kw1,0ky1ky,w'y −−+=++=+=
k kx )kx(yky ≅ )kx('ykw ≅
0 0 1 0
1 0,1 1 -0,2
2 0,2 0,98 -0,3805
3 0,3 0,9414 -0,5518
81
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Exemplo2: Utilizar o método de Euler com h=0,1, estimar o valor da soluçãodo problema de
valor inicial
==
=+
0)0´(y,0)0(y
xe2y2''y
, no intervalor [0 ; 0,5].
Solução:
−+=++=+=
2
ky2
kxe1,0kw1kw,kw1,0ky1ky,w'y
Exemplo3: Utilizar o método de Euler com h=0,2, estimar o valor da soluçãodo problema de
valor inicial
===+=++
4)1(''y,2)1('y,3)1(y
yx'y''y2'''y
, no intervalor [ 1,2]
Solução:
)kwkz2kykx(2,0kz1kz,wz2yx'z
kz2,0kw1kw,z'w
,kw2,0ky1ky,w'y
−−++=+−−+=
+=+=
+=+=
k xk )kx(yky ≅ )kx('ykw ≅
0 0 0 0
1 0,1 0 0,1
2 0,2 0,01 0,21051
3 0,3 0,0310 0,3326
4 0,4 0,0643 0,4674
5 0,5 0,1110 0,6157
k xk )kx(yky ≅ )kx('ykw ≅ )kx(''ykz ≅
0 1 3 2 4
1 1.2 3,4 2,8 2,8
2 1,4 3,962 3,36 2,04
3 1,6 4,632 3,768 1,624
4 1,8 5,3856 4,0928 1,4672
5 2 6,2041 4,3862 1,4988
82
E.B.Hauser – Cálculo Numérico
Solução de Problema de Valor de Contorno por Diferenças Finitas
Derivada Diferença Finita
h, k=tamanho do passo na direção x e na direção y (ou t)
)ix´(y atrasada
h
1iyiy
,centrada
h2
1iy1iy
,avançada
h
iy1iy −−−−+−+
)ix´´(y
2h
1iyiy21iy −+−+
centrada
)ix('''y 3h2
2iy1iy21iy22iy −−−++−+
centrada
)ix(IVy 4h
2iy1iy4iy61iy42iy −+−−++−+
centrada
( )jt,ixtu∂∂
k2
1j,iu1j,iu −−+
centrada
( )jt,ix2x
u2
∂
∂
2h
j,1iuj,iu2j,1iu −+−+
centrada
( )jy,ix2y
u2
∂
∂
2k
1j,iuj,iu21j,iu −+−+
centrada
A solução do Problema de valor de contorno de segunda ordem, linear:
(PVC)
==
=++
.y)x(y,y)x(y
)x(fy)x(Q'y)x(P''y
nnoo
por Diferenças Finitas(centradas), para i= 1,...,n-1 , é gerada pelo sistema de n-1 equações li-
neares:
( ) )x(fhy)x(P
2
h1y)x(Qh2y)x(P
2
h1 i21iiii
2
1ii =
−++−+
+
−+ .
83
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Exemplo1: Utilizar diferenças finitas, com h=2 estimar y(3) se
−==
=+−
2)5(y,0)1(y
xlny'xy''y2x
.
Solução:
2
kx
)kxln(4
1ky
kx
11ky22
kx
4
1ky
kx
11 =
−
++
−++
−
Exemplo2: Utilizar diferenças finitas, com n=3, para resolver o PVC
=−=
=+
2)4(y,1)1(y
xy5''y
.
Solução: h=(4-1)/3=1, kx1kyky31ky =−+++
Exemplo3: Utilizar diferenças finitas, com h=1, para resolver o PVC
==
−=+
7)4(y,5,1)1(y
x'y3''y
.
Solução: kx1ky5,0ky21ky5,2 −=−−−+
k kx )kx(yky ≅
0 1 0
1 3 -1,1710
2 5 -2
k kx )kx(yky ≅
0 1 -1
1 2 0,9382
2 3 -0,4022
3 4 2
k kx )kx(yky ≅
0 1 1,5
1 2 10,2380
2 3 7,6904
3 4 7
84
E.B.Hauser – Cálculo Numérico
Exercícios:
1. Utilizando o método de Euler, determinar y(X i ), se :
a)
=
−−=
2)1(y
1xyy`
, h=0,25 , X i =1,75
b)
=
>+=
0)1(y
0y,yxy`
, h=0,1 , X i =1,4
c)
=
+=
2,0)1(y
y2xy` 2
, h=0,1 , X i =1,4
d)
=
−−=
2)0(y
100e101y100y` x
, h=0,1 , X i =0,3
e)
−==
=−−
5,0)1`(y,5.0)1(y
0yx8y```xy 33
, h=0,1 , X i =1,3
f)
==
=+−−
0)0`(y,5,0)0(y
0yy`)y1(``y 2
, h=0,1 , X i =0,3
2.Utilizando o Método de Heun (Runge-Kutta de 2ª ordem), determinar y( X i ), se:
a)
=
+=
1)0(y
y2x3y´
, h=0,2, X i =0,4 b)
=
−=
1)0(y
xyy´
, h=0,1, X i =0,2
c)
=
=
0)0(y
3x4y´
, h=0,1 , X i =0,3
3.Utilizando o Método das diferenças finitas, e o valor indicado de n, resolver o PVC.
a) 4n,
1)2(y,4)0(y
0y9''y
=
==
=+
b) 5n,
0)1(y,4)0(y
x5y'y2''y
=
==
=++
c) 8n,
0)2(y,5)1(y
0y3'xy3''yx2
=
==
=++
d) 10n,
2)1(y,0)0(y
xxy'y)x1(''y
=
==
=+−+
85
7 -Resolução Numérica de Equações Diferenciais Ordinárias
4. Considerar um sistema massa-mola-amortecedor descrito pela equação diferencial ordinária de
segunda ordem:
m y” + cy’ + ky = L sin x.
Utilizando Euler com h=0.25, estimar o deslocamento para o tempo x=0.5, para massa
m=1, amortecedor c=0.5, rigidez k=2, amplitude da força L=0.5 , com deslocamento inicial
y(0)=1, velocidade inicial y’(0) =0.
5. Equação Diferencial Logística na análise do Crescimento Populacional Mundial
Aplicar o método de Euler, com h=10 , na equação diferencial logística
=
−=
555,2)0(y
y
12
y1026,0y´
para simular a população mundial (y bilhões de habitantes) no período de 1950 a 2000.
Comparar com os dados do censo populacional mundial, calculando o erro relativo.
Completar a tabela a seguir:
Ano
k
xk
População estimada por
Euler h=10
yk
Censo
Populacional
Mundial(bilhões)
pk
Desvio Relativo
Percentual
1950 0 0 2,555
2, 555
1960 1 10
3, 040
1970 2 20
3,708
1980 3 30
4, 454
1990 4 40
5, 276
2000 5
50 6, 079
86
E.B.Hauser – Cálculo Numérico
Respostas :
1.a) y(1,75)=1,2018 2.a) y(0,4)=1,1189
b) y(1,4)=0,4558 b) y(0,2)=0,9801
c) y(1,4)=0,7778 c) y(0,3)=0,009
d) y(0,3)=-25,2206
e) y(1,3)= 0, 3647
f) y(0,3)= 0,3991
3.a)
3226,6y
5807,2y
6774,5y
3
2
1
=
−=
−=
b)
1490,04y
5651,03y
3110,12y
4147,,21y
=
=
=
=
c)
2913,0y
6430,0y
0681,1y
5826,1y
2064,2y
9640,2y
8842,3y
7
6
5
4
3
2
1
=
=
=
=
=
=
=
d)
8474,1y
6855,1y
5149,1y
3353,1y
1465,1y
9471,0y
7357,0y
5097,0y
2660,0y
9
8
7
6
5
4
3
2
1
=
=
=
=
=
=
=
=
=
4) y(0.5)=0.875
5)
Ano
k
xk
População estimada
por Euler h=10
yk
Censo
Populacional
Mundial(bilhões)
pk
Desvio Relativo
Percentual
1950 0 0 2,555
2, 555 0
1960 1 10 3,0778
3, 040 1,25%
1970 2 20 3,6728
3,708 0,94%
1980 3 30 4,3355
4, 454 2,65%
1990 4 40 5,0554
5, 276 4,17%
2000
5
50 5,8161 6, 079 4,31%
87
8 – Resolução Numérica de Equações Diferenciais Parciais
Exemplo1:
Consideremos o modelo para determinar a temperatura de estado estacionário para uma pla-
ca quadrada com condições de contorno dadas:
<≤−
<<
==
−==
<<<<=
∂
∂
+
∂
∂
2x1,x2
1x0,x)2,x(u,0)0,x(u
)y2(y)y,2(u,0)y,0(u
2y0,2x0,0)y,x(
y
u)y,x(
x
u
2
2
2
2
Utilizando diferenças finitas centrais com h=k=2/3, para i, j = 0,1,2,3, obtemos a equação
de diferenças finitas:
0u4uuuu j,i1j,ij,1i1j,ij,1i =−+++ −−=+
Para determinar o valor de
22u)3/4,3/4(ue21u)3/2,3/4(u,12u)3/4,3/2(u,11u)3/2,3/2(u≅≅≅≅
Utilizando o sistema Maple para resolver o sistema linear tridiagonal:
−=−+
−=+−
−=+−
=++−
9/1422u412u21u
3/222u12u411u
9/822u21u411u
012u21u11u4
> solve({-4*u11+u21+u12=0, u11-4*u21+u22=-8/9, u11-4*u12+u22=-2/3,
u21+u12-4*u22=-14/9},{u11,u12,u21,u22});
{ }, , , = u12 1336 = u11
7
36 = u22
7
12 = u21
5
12
> evalf(%);
{ }, , , = u12 .3611111111 = u11 .1944444444 = u22 .5833333333 = u21 .4166666667
88
E.B.Hauser – Cálculo Numérico
Exemplo2
Consideremos o problema da determinação do estado estável da distribuição de calor em
uma placa quadrada metálica delgada, com dimensões 0,5 m por 0,5 m. Dois limites adjacentes
são mantidos a 0ºC, e o calor nos outros dois limites aumenta linearmente de 0ºC, em um canto,
para 100ºC no lugar onde ambos os lados se encontram. Se colocarmos os lados com as condi-
ções de limite igual a zero ao longo dos eixos x e y , o problema pode ser expresso como
,0),(),( 2
2
2
2
=
∂
∂
+
∂
∂ yx
y
uyx
x
u
para (x , y) no conjunto },5,0y0,5,0x0/)y,x{(R <<<<= com as condições de fronteira
u (0 , y)=0, u (x , 0)=0, u (x , 0,5)=200x, u (0,5 , y) = 200y
Consideremos n = m = 4. Assim, com h=k=0.125. Construímos a malha da figura 7:
Figura 7
Utilizandos o método das Diferenças Finitas (2.4) obtemos a equação de diferenças
finitas
,0wwwww4 ij,iij,ij,1ij,1ij,i =−−−− +−−+
para todo i= 1,2,3 e j=1,2,3 .
Para expressar isso em função dos interiores da grade, usamos
( )jil y,xP = e ijl ww =
e
)( ii Puw = , e l = i+(m-1-j)(n-1).
89
8 - Resolução Numérica de Equações Diferenciais Parciais
Também, a partir das condições de contorno e usando (2.5)
75ww,e,50ww,25ww
,0wwwwww
3,44,32,44,21,44,1
3,02,01,00,30,20,1
======
======
Assim, para cada ponto interior da grade, iP geramos uma equação linear:
,wwwww4:P
,wwwww4:P
,wwwww4:P
,wwwww4:P
0wwwww4:P
wwwww4:P
wwwww4:P
wwwww4:P
wwwww4:P
1,40,36899
0,257988
0,11,04877
2,493566
824655
2,071544
4,33,46233
4,251322
4,13,04211
+=−−
=−−−
+=−−
=−−−
=−=−−
=−−−
+=−−
=−−−
+=−−
Obtemos um sistema linear cuja a forma matricial é:
Resolvendo o sistema com o sistema de computação algébrica e simbólica Maple , obtemos
as temperaturas aproximadas nos pontos interiores da malha.
90
E.B.Hauser – Cálculo Numérico
> A := matrix( [[4,-1,0,-1,0,0,0,0,0],[-1,4,-1,0,-1,0,0,0,0],[0,-
1,4,0,0,-1,0,0,0],[-1,0,0,4,-1,0,-1,0,0],[0,-1,0,-1,4,-1,0,-
1,0],[0,0,-1,0,-1,4,0,0,-1],[0,0,0,-1,0,0,4,-1,0],[0,0,0,0,-1,0,-
1,4,-1],[0,0,0,0,0,-1,0,-1,4]]):
> b := vector( [25,50,150,0,0,50,0,0,25]):
> linalg[linsolve](A, b);
> W:=evalf(%);
Tabelamos os resultados:
i 1 2 3 4 5 6 7 8 9
iw 18.75 37.50 56.25 12.5 25.00 37.50 6.25 12.50 18.75
Assim, para cada i, )( ii Puw = , isto é, iw é uma estimativa da solução em )y,x(P iii = .
Por exemplo, 5.37)25.0,375.0(u)P(uw 66 ≅==
Exercícios:
1) Com h=k=10, resolver numericamente equação de Laplace
,60y0,80x0,0)t,x(2x
u2)t,x(2y
u2
<<<<=
∂
∂
+
∂
∂
sujeita às condições
100)y,0(ue,0)y,80(u)60,x(u)0,x(u ==== .
j↓ 0 1 2 3 4 5 6 7 8 i←
6 100 0 0 0 0 0 0 0 0 60
5 100 46,993 24,835 13,978 8,068 4,633 2,523 1,110 0 50
4 100 63,138 38,368 23,010 13,660 7,942 4,348 1,917 0 40
3 100 67,192 42,490 26,032 15,620 9,128 5,009 2,211 0 30
2 100 63,138 38,368 23,010 13,660 7,942 4,348 1,917 0 20
1 100 46,993 24,835 13,978 8,068 4,633 2,523 1,110 0 10
0 100 0 0 0 0 0 0 0 0 0
x → 0 10 20 30 40 50 60 70 80 y↓
4
ij,iuij,iuj,1iuj,1iuj,iu +
+
−
+
−
++
=
91
8 - Resolução Numérica de Equações Diferenciais Parciais
2) Resolver numericamente equação de calor
,t0,1x0,0)t,x(
x
u)t,x(
t
u
2
2
≤<<=
∂
∂
−
∂
∂
com as condições de contorno
,t0,0)t,1(u)t,0(u <==
e as condições iniciais
.1x0),x(sen)0,x(u ≤≤= pi
Ajuda:
Diferenças finitas para o problema com h=0,2, k=0.01
−
+
+
+=+ j,1iuj,1iu25,0j,iu5,01j,iu
j↓ 0 1 2 3 4 5 ← i
6 0 0,3219 0,5208 0,5208 0,3219 0 0,06
5 0 0,3559 0,5758 0,5758 0,3559 0 0,05
4 0 0,3934 0,6366 0,6366 0,3934 0 0,04
3 0 0,4350 0,7038 0,7038 0,4350 0 0,03
2 0 0,4809 0,7781 0,7781 0,4809 0 0,02
1 0 0,5317 0,8602 0,8602 0,5317 0 0,01
0 0 0,5878 0,9511 0,9511 0,5878 0 0,00
x → 0 0,2 0,4 0,6 0,8 1 ↑ tempo
Diferenças finitas para o problema com h=0,1, k=0.005
−
+
+
+=+ j,1iuj,1iu1,0j,iu9,01j,iu
j↓
0 1 2 3 4 5 6 7 8 9 10 ← i
6 0 0 0,5189 0,9869 1,3584 1,5969 1,6791 1,5969 1,3584 0,9869 0,5189 0 0,030
5 1 0 0,4759 0,9053 1,2460 1,4647 1,5401 1,4647 1,2460 0,9053 0,4759 0 0,025
4 2 0 0,4365 0,8304 1,1429 1,3435 1,4127 1,3435 1,1429 0,8304 0,4365 0 0,020
3 3 0 0,4004 0,7616 1,0483 1,2324 1,2958 1,2324 1,0483 0,7616 0,4004 0 0,015
2 4 0 0,3673 0,6986 0,9616 1,1304 1,1886 1,1304 0,9616 0,6986 0,3673 0 0,010
1 5 0 0,3369 0,6408 0,8820 1,0369 1,0902 1,0369 0,8820 0,6408 0,3369 0 0,005
0 6 0 0,3090 0,5878 0,8090 0,9511 1,0000 0,9511 0,8090 0,5878 0,3090 0 0,000
x → x 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 ↑ tempo
92
Equação Diferenças Finitas
h, k=tamanho do passo na direção x e na direção y (ou t)
Célula Computacional
Laplace
0)y,x(2y
u2)y,x(2x
u2
=
∂
∂
+
∂
∂
Centradas, com h=k
4
ij,iuij,iuj,1iuj,1iu
j,iu
++−+−++
=
•
••
•
−
+−
+
1j,iu
j,1iuj,iuj,1iu
1j,iu
o
Calor
)t,x(
t
u)t,x(2x
u2
∂
∂
=
∂
∂
α
Material α >0
Prata 1,7
Cobre 1,5
Alumínio 0,85
Ferro 0,15
Concreto 0,005
progressivas
−
+
+
+
−=+ j,1iuj,1iu2h
k2
j,iu2h
k2211j,iu
αα
Converge se
2
1
2h
k
<
•••
+−
+
j,1iuj,iuj,1iu
1j,iu
o
Onda
)t,x(2t
u2)t,x(2x
u22
∂
∂
=
∂
∂β
centradas
1j,iuj,1iuj,1iu2h
2k2
j,iu2h
2k2121j,iu
−
−
−
+
+
+
−=+
ββ
Se 12h
2k
e1 ==β , 1j,iuj,1iuj,1iu1j,iu −−++−=+
•
•••
−
+−
+
1j,iu
j,1iuj,iuj,1iu
1j,iu
o
93
PUCRS – FAMAT - Cálculo Numérico A - Formulário
1) Sistema de Ponto Flutuante ),,,( MmtF β , # F = 1)1mM(1t)1(2 +
+−−− ββ , ),,,( MmtF β .
2) Polinômios: Seja 0ax1a
2
x2a
3
x3a...
1n
x1na
n
xna)x(p +++++−
−
+=
Horner: 0ax)1ax)2a...x)2nax)1naxna
1n
(((...()x(p +++
−
+
−
+
−
=
321
Huat: Se p(0) ≠ 0 e para algum k = 1, 2, ...,n-1, 2ka ≤ 1k1k aa +− , então p tem raízes complexas.
3) Método de Newton-Raphson: 0)kx('fmáx,...,2,1,0k,)kx('f
)kx(fkx1kx ≠=−=+ .
4) Sistema de n equações lineares: AX=B, A =
ija , X =
ix e B =
ib com i,j = 1, 2, ..., n.
2
kna
2
2ka
2
1kakonde
n21
AdetANORM L
L
++== α
ααα
, para k = 1, 2, ..., n.
A matriz A é Diagonal Dominante se .n,..,2,1j,i
n
ij
1j
ijaiia =∀
≠
=
> ∑
5) Interpolação Polinomial: Sejam )iy,ix( , ix1ixh −+= , i= 0,1, 2, ..., n .
Polinômio Interpolador de Newton para Diferenças Finitas Ascendentes:
oy
n
nh!n
)1nxx()1xx)(oxx(
oy
2
2h!2
)1xx)(oxx(
0yh
)oxx(
oy)x(p ∆∆∆ −
−−−
++
−−
+
−
+=
L
L
Se
h
)oxx(z −= ,
!n
oy
n
))1n(z()1z(z
!3
oy
3
)2z)(1z(z
!2
oy
2
)1z(z0yzoy)z(p
∆∆∆
∆ −−−++−−+−++= LL .
Polinômio Interpolador de Newton para Diferenças Divididas
oy
n)1nxx)...(1xx)(oxx(...oy
2)1xx)(oxx(0y)oxx(oy)x(p ∆∆∆ −−−−++−−+−+=
6) Ajuste de Dados: Sejam ,n,...1,0i,)iy,ix( = , ∑ ∑
=
=
n
0i
.
Ajuste Linear: x1aoa)x(Y += ,
+
∑∑
∑
1a
0a
2
ixix
ix1n
=
∑
∑
iyix
iy
,
Ajuste Quadrático: 2x2ax1aoa)x(Y ++= ,
+
∑∑∑
∑∑∑
∑∑
4
ix
3
ix
2
ix
3
ix
2
ixix
2
ixix1n
2a
1a
o
a
=
∑
∑
∑
iy
2
ix
iyix
iy
.
94
7) Integração Numérica: Para
n
oxnxh −= , ihoxix += e )( ixfiy = , i = 0,1,...n,
Trapézios: [ ]n1n21ox
x
y)yyy(2y
2
hn
o
dx)x(f +++++≅
−∫ L , i
20n
T ymax12
)xx(
E ∆−≈ .
Simpson ( n par): ∫
+
−
++++
−
++++≅
nx
ox
ny)2ny4y2y(2)1ny3y1y(4oy3
hdx)x(f LL i40nS ymax180
)xx(
E ∆−≈ .
8) Equações Diferenciais Ordinárias
Derivação Numérica : Três formas de aproximar a função derivada por Diferencas Finitas, ix1ixh −+= :
atrasada
h
1iyiy)iy´(x,centradah2
1iy1iy)iy´(x ;avançadah
iy1iy)iy´(x −
−
≅−
−+≅
−+≅
11) PVI: Sejam
=
=
oy)ox(y
)y,x(f'y
iyixy ≅)( com ihoxix += , i= 0,1,...,n e
n
oxnxh −= .
11 a)Euler: )iy,ix(hfiy1iy +=+ , i= 0,1,...,n-1.
11b) Runge-Kutta de 2
a
Ordem : ( )2k1k2
h
iy1iy ++=+ , )iy,ix(f1k = e )hky,hx(fk 1ii2 ++= , i= 0,1,...,n-1.
12) PVC(Linear): Sejam
==
=++
.y)x(y,y)x(y
)x(fy)x(Q'y)x(P''y
nnoo
, iyixy ≅)( com ihoxix += , i= 0,1,...,n e n/)oxnx(h −= .
Diferenças Finitas: ( ) )x(fhy)x(P
2
h1y)x(Qh2y)x(P
2
h1 i
2
1iiii
2
1ii =
−++−+
+
−+ , i= 1,...,n-1
9) Equações Diferenciais Parciais: Derivadas parciais aproximadas por Diferencas Finitas Centradas , com malha definida por
ix1ixh −+= , it1itk −+= . ( ) k2 1j
,i
u
1j,iu
jt,ixt
u −
−
+
≅
∂
∂
,
( )
2h
j,1iuj,iu2j,1iu
jt,ix2x
u2 −
+−
+
≅
∂
∂
,, ( )jy,ix2y
u2
∂
∂
2k
1j,iuj,iu21j,iu −+−+
≅ .
Se h=k, a equação de Laplace 0)y,x(2y
u2)y,x(2x
u2
=
∂
∂
+
∂
∂
é aproximada por
4
ij,iuij,iuj,1iuj,1iu
j,iu
++−+−++
= .
Equação da Onda )t,x(
2t
u2)t,x(
2x
u22
∂
∂
=
∂
∂β , 1j,iuj,1iuj,1iu2h
2k2
j,iu2h
2k2121j,iu
−
−
−
+
+
+
−=+
ββ
.
95
PUCRS - Faculdade de Matemática - Prof. Eliete Biasotto Hauser
Cálculo Numérico - Laboratório1 - Introdução ao Matlab, Maple, VCN
Notação Comentários Exemplos
1 + Adição 3+5
2 - Subtração 758-195
3 * Multiplicação 7.5*8
4 / Divisão 39/13
5 ^ Potenciação 2^12
6 abs ( x ) Valor Absoluto de x abs(-7)
7 sqrt (x) Raiz Quadrada de x sqrt(7835)
8 x^(1/n) Raiz n-ésima de x 27^(1/3)
9 factorial(n) Fatorial de n factorial(5)
10 pi pi pi
11 sqrt(-1) ou i ou j Unidade Imaginária sqrt (-1)
i
j
12 exp(1) Número de Euler: e exp(1)
13 exp(x) b^x Função Exponencial exp(3.5)
7^10
14 log(x)
logb(x) = log(x)/lob(b).
Logaritmo Natural
Logaritmo de base b
log (7)
log10(3)
15 sin(x) cos(x) tan(x) sec(x) csc(x)
cot(x) sinh(x) cosh(x) tanh(x)
sech(x) csch(x) coth(x) asin(x)
acos(x) atan(x) asec(x) acsc(x)
arccot(x) asinh(x) acosh(x) atanh(x)
asech(x) acsch(x) acoth(x)
Funções trigonométricas, hiperbó-
licas e inversas
(argumento x em radianos)
sin(pi/2)
cos(pi)
acosh(0)
16 help Abre janela com conceitos e
exemplos
help format
help plot
17 realmax
realmin
Maior e menor número positivo em
ponto flutuante
Obs: comparar com calculadora e
Maple_floats(MAX_FLOAT);
0.1 102147483647
Maple_floats(MIN_FLOAT);
0.1 10-2147483645
realmax
realmin
exp(500)
exp(5000)
Para criar um arquivo.m, selecionar no menu principal do Matlab File/New/M-file e após, Save. O nome não
pode iniciar com número, comando do matlab, e conter espaços em branco, caracteres do tipo ( . , % @ ... ) .
Para salvar e executar pressionar F5.
ans
inf
NaN
Resposta mais recente
Infinito
Não-numérico
whos Mostra informações sobre as variáveis armazenadas na memória
clear Apaga todas as varáveis da memória do Matlab
G
r
á
f
i
c
o
s
x = -10:0.1:10;
y=log(abs(-3*x));
plot(x,y)
plot(x,x.^2-4.5*x+9)
plot(x,1-x.*cos(x))
plot(x,1-x.*cos(x),'o')
x = -7*pi:0.3: 7*pi;
plot(x,sin(x),x,sin(x-pi))
plot(x,sin(x),'r',x,sin(x-pi),'b')
x = -2:0.1:2;
plot(x,2*x,x,cos(x))
x = -90:0.1:90;
plot(x,1-x.*cos(x))
plot(x,1-x.*cos(x),'o')
VNC- Utilitários: 1.12Calculadora ; 1.17Formato dos Números;1.19 Salva arquivo
96
96
Bibliografia
ATKINSON, K. E. An Introduction to Numerical Analysis. 2.ed. New York : John Wiley & Sons,
1989.
AYYUB, Bilal M.; McCUEN, Richard H. Numerical Methods For Engineers. New Jersey: Prentice
Hall, 1996.
BARROSO, Leonidas Conceição et al. Cálculo Numérico com Aplicações. 2.ed. São Paulo :
Harbra, 1987.
BURDEN, Richard L., FAIRES, J. Douglas. Análise Numérica. 8.ed. São Paulo : Thomson, 2008.
CHAPRA, Steven C., CANALE, Raymond P. Métodos Numéricos Para Engenharia. 5.ed. São
Paulo: McGraw-Hill, 2008.
CLÁUDIO, Dalcidio Moraes, MARINS, Jussara Maria. Cálculo Numérico Computacional. 2.ed.
São Paulo : Atlas, 2000.
CUNHA, M.Cristina, Métodos Numéricos, Campinas, SP, 2003, Ed. Unicamp.
FRANCO, Neide B.: Cálculo Numérico: São Paulo, Pearson Prentice Hall, 2006.
GANDER, Walter.; HREBICEK, Jiri. Solving Problems in Scientific Computing Using Maple and
Matlab: Berlin, Springer-Verlag, 1995.
GILAT, Amos; SUBRAMANIAM, Vish. Métodos numéricos para engenheiros e cientistas: uma
introdução com aplicações usando o MATLAB. Porto Alegre: Bookman, 2008.
HUMES, Ana Flora P. de Castro et al. Noções de Cálculo Numérico. São Paulo: McGraw-Hill,
1984.
KREYSZIG, Erwin. Advanced Engineering Mathematics. New York, NY : John Wiley & Sons,
1993.
O’NEIL,PeterV. Advanced Engineering Mathematics. 4.ed. Pacific Grove, CA : Brooks/Cole,
1995.
RICE, Richard G.; DO, Duong D. Applied Mathematics and Modelling for Chemical Engineers.
New York : John Wiley & Sons, 1995.
RUGGIERO, Márcia A. Gomes, LOPES, Vera Lúcia da Rocha. Cálculo numérico: aspectos
teóricos e computacionais. São Paulo : McGraw-Hill, 1997.
SCHELEIDER, Maria Amélia N, Cunha, Maria Cristina dC, Métodos Numéricos para Equações
Diferenciais Parciais, Notas em Matemática Aplicada, Volume 4, São Carlos, SP, 2003, SBMAC.
(Disponível em http://www.sbmac.org.br/boletim/pdf_2003/livro_04_2003.pdf)
STROUD, K.A, BOOTH, Dexter J., Advanced Engineering Mathematics, New York, Palgrave
Macmillan, 200
97
Índice
1 Sistema de Ponto Flutuante Normalizado – Teoria dos Erros....................1
2 Resolução de Equações Algébricas e Transcendentes.............................. 7
3 Sistemas de Equações Lineares................................................................24
4 Interpolação Polinomial............................................................................ 39
5 Ajuste de Funções..................................................................................... 49
6 Integração Numérica.................................................................................69
5 Resolução Numérica de Equações Diferenciais Ordinárias..................... 76
8 Resolução Numérica de Equações Diferenciais Parciais.......................... 88
Formulário................................................................................................... 94
Laboratório utilizando Sistema Maple....................................................... 96
Bibliografia ..............................................................................................96