PONTIFÍCIA UNIVERSIDADE CATÓLICA DO RIO GRANDE DO SUL
FACULDADE DE MATEMÁTICA
CÁLCULO NUMÉRICO
Notas de Aula – Aplicações – Exercícios
Eliete Biasotto Hauser
Índice
1 Sistema de Ponto Flutuante Normalizado – Teoria dos Erros....................4
2 Resolução de Equações Algébricas e Transcendentes.............................. 9
3 Sistemas de Equações Lineares................................................................ 25
4 Interpolação Polinomial............................................................................ 36
5 Ajuste de Funções..................................................................................... 46
6 Integração Numérica.................................................................................60
5 Resolução Numérica de Equações Diferenciais Ordinárias..................... 65
8 Resolução Numérica de Equações Diferenciais Parciais.......................... 72
9 Bibliografia ..............................................................................................85
Formulário................................................................................................... 86
Laboratórios utilizando Sistema Maple....................................................... 90
1 –Teoria dos Erros - Sistema de Ponto Flutuante
Um número y na base β ≥ 2 , y = a n a
L a a a a , b b2 L pode ser descrito
n −1
3 2 1 0 1
na forma
β n − 1 + L + a 3 β 3 + a 2 β 2 + a1 β + a0 + b1 β − 1 + b2 β − 2 + L .
y = an β n + a
n −1
Por exemplo, 3517 ,26 = 3 × 10 3 + 5 × 10 2 + 1 × 10 + 7 + 2 × 10 − 1 + 6 × 10 − 2
Em aritmética de ponto flutuante normalizado de t dígitos, y tem a forma:
y = ⎛⎜ 0 , d d L d t ⎞⎟ × β e
1 2
⎝
⎠β
i)
ii)
0 , d d L d t é a mantissa (uma fração na base β),
1 2
0 ≤ d j ≤ β − 1, d 1 ≠ 0 , 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
0 = 0 .000
L
1
42
4
30 × β
t
vezes
é 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) 0,1 x β m é o menor número em módulo, não nulo, de F.
2) 0, ( β − 1)( β − 1)...( β − 1) x β M é o maior número, em módulo, não nulo, de F.
144424443
t vezes
3) # F = ⎡⎢2( β − 1) β t −1 ( M − m + 1)⎤⎥ + 1 é o cardinal de F
⎦
⎣
−
4) Para qualquer mantissa, vale β 1 ≤ mantissa < 1 .
5) Se y ∈ F , então − y ∈ F .
6) 0 ∈ F e 1 ∈ F .
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)
4
E.B.Hauser – Cálculo Numérico
Se a representação do real y em F não é exata, é necessário utilizar um
arredondamento. 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)
Em F, geralmente, as operações de adição e multiplicação não são comutativas,
associativas 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 diferentes.
Medidas de Exatidão
Quando se aproxima um número real x por x*, o erro que resulta é x-x*. Define-se:
x − x*
para x ≠ 0 .
erro absoluto: EA = x − x * e o erro relativo: ER =
x
A fim de ver o tipo de situação que pode ocorrer um erro relativo de grande
magnitude, vamos considerar a diferença entre os números a seguir, por exemplo:
x = 0,3721478693
y = 0,3720230572
x − y = 0,0001248121=0,1248121 x 10-3
Se os cálculos forem feitos em F(10, 5, -499, 499) com arredondamento Ox:
x*= 0,37215, y*= 0,37202 e x*-y*= 0,00013 = 0,13000 x 10-3
Assim, o erro relativo entre os dois resultados é grande:
x − y − ( x * -y*)
0,1248121 x 10 - 3 − 0 ,13000 x 10 - 3
=
≈ 4%
x− y
0,1248121 x 10 - 3
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:
⎡
⎛
x k +1 − x k ⎞⎤
⎟⎥
DIGSE ( x k , x k +1 ) = − ⎢0,3 + log⎜⎜ µ +
⎟
x
k +1 ⎠⎦⎥
⎝
⎣⎢
o qual expressa o número de dígitos significativos exatos de x k em relação a x k + 1 . Aqui
µ representa a unidade de arredondamento da máquina ( µ =
to for Ox ).
5
1 1− t
β
se o arredondamen2
1 - Teoria dos Erros
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 detectados 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:
Problema
Real
Modelagem
Matemática
Solução
Estudaremos algoritmos numéricos a fim de obter uma solução numérica do
problema, 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
linearizaçõ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
Aplicação
O volume de uma esfera de raio r pode ser obtido de
4π r 3
V=
.
3
Estimar o volume do hemisfério de raio “e” representado
graficamente ao lado. Utilizar arredondamento para número
mais próximo de máquina (0x) em F( 10 , 4 , -98 , 98).
6
E.B.Hauser – Cálculo Numérico
Exercícios
1) No sistema MapleV estimar e − 8.3 utilizando
∞
( − x )i
e
e− x =
i!
i =0
∑
e− x =
1
1
=
x
∞
e
xi
i!
i =0
∑
com 26 termos cada e comparar com e − 8.3 ≈ 0.2485168271x10 .
-3
2) Em F(10 , 3 , -98 , 98) e arredondamento por truncamento estimar p(2.73) se:
a)
p( x ) = x 3 − 5 x 2 + 6 x + 0.55
b)
p( x ) = (( x − 5 )x + 6 )x + 0.55
Em ambos os casos estimar o erro absoluto ao comparar com p(2.73) ≈ 0.11917x10
Obs: Estimar p(x) pelo algoritmo usual
p ( x) = a n x n + a
x n −1 + ... + a x 3 + a x 2 + a x + a
n −1
3
2
1
0
-1
.
2
exige n adições e (n +n)/2 multiplicações enquanto que o algoritmo de Horner
p ( x) = (((...(
{ a n x + a n −1 ) x + a n − 2 ) x + ...a 2 ) x + a1 ) x + a 0
n −1
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
respectivamente. Utilizando a propriedade associativa, pode-se determinar o produto
matricial AxBxCxD de diversas formas. Qual das duas abaixo é mais eficiente? Porque?
b) (Ax(BxC))xD
a) Ax(Bx(CxD))
4) Representar o número real x na base 2 usando 8 algarismos significativos? Essa
representaçã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
truncamento( ∇x )e arredondamento para número mais próximo de máquina (Ox) em
F(10,5,-2,2).
200
2000
1
b)
e)
c)
d) e
f) 2
a) 3
3
3
3000
10 1
∑ i
i =1 2
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.
7) Considerando: A =
7
1 - Teoria dos Erros
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) (0,6)10 ≈ (0,10011001)2
b) (13,25)10 ≈ (1101,01) 2 c) (2.47)10 ≈ (10,011110) 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: ( −∞,−7 / 2) ∪ (7 / 2,+∞)
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, 5/16, 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
d) 0.27182*101 e 0.27183*102
c) overflow (0,66666*103 ∉ e 0,66667*103 ∉ F)
–3
1
e) underflow (0.33333*10 ∉ F) f) 0,14142*10 e 0,14142 *101
7-a) 0,9990234375
b) 0,999 e 0,998
8
2. Resolução de Equações Algébricas e Transcendentes
Objetivo: Resolver equações de forma f ( x ) = 0 , isto é, determinar os valores de x para os quais a
igualdade f ( x ) = 0 é satisfeita.
Se a função f ( x ) só contém operações algébricas repetidas um número finito de vezes, a
equação é dita algébrica.
1 5
Ex.: polinômios, x 2 − x − 1 = 0 ,
− x =0
x3
Equações Transcendentes são aquelas em que a variável independente x está submetida a
uma operação não algébrica um número finito de vezes.
Ex.: e x + x 2 + 1 = 0 ,
ln x + tgx = 0 ,
senx − e x = 0
2.1- Equações Polinomiais
Revisão sobre polinômio:
Seja p( x ) = a n x n + a n −1 x n −1 + K + a 2 x 2 + a1 x + a 0 um polinômio de grau n, onde os
coeficientes , são números reais e an ≠ 0 .
Temos que:
1) p( x ) possui n raízes, contadas as multiplicidades.
2) Se x1 , x2 ,K, xn forem raízes reais de p( x ) , então p( x ) pode ser fatorado na forma:
p( x ) = an ( x − x1 )( x − x2 )K( x − xn ) .
Ex1: p( x ) = x 3 + 2 x 2 − x − 2 , tem raízes:
x1 = 1, x2 = −1, x3 = −2 . Logo,
p ( x ) = ( x − 1)(x + 1)( x + 2 )
3) Se a + bi é raiz de p( x ) então z = a − bi também o é.
Ex.: p( x ) = x 2 − 6 x + 10 tem raízes 3 ± i .
9
E.B.Hauser – Cálculo Numérico
4) Se a + bi é raiz de p( x ) de grau n ≥ 2 , então p( x ) pode ser fatorado:
(
)
p( x ) = x 2 − 2ax + a 2 + b 2 ⋅ q( x )
onde o grau de q( x ) é n − 2 .
Ex.: a) p( x ) = x 4 − 2 x 3 + x 2 + 2 x − 2 tem raízes ± 1, 1 ± i .
(
)(
(
)
)
p( x ) = x 2 − 2 x + 2 ⋅ x 2 − 1
Ex.: b) p( x ) = x 3 − 7 x 2 + 16 x − 10 tem raízes 1, 3 ± i .
p( x ) = x 2 − 6 x + 10 ⋅ ( x − 1)
5)Se p ( x ) é de grau ímpar, então p ( x ) possui ao menos uma raiz real.
6) Uma raiz x0 de p( x ) tem multiplicidade m se p( x0 ) = p ' ( x0 ) = p " ( x0 ) = K = p m−1 ( x0 ) = 0
e p m ( x0 ) ≠ 0
Ex.:1) x0 = 2 é raiz de multiplicidade 2 p( x ) = 2 x 3 − 6 x 2 + 8
2) x0 = 2 é raiz de multiplicidade 3 de p( x ) = x 4 − 5 x 3 + 6 x 2 + 4 x − 8
10
2 - Resolução de Equações Algébricas e Transcendentes
7) Valor numérico de um polinômio: para calcular, de forma usual, p(xi ) são necessárias n
n(n + 1)
adições e
multiplicações.
2
O Método de Horner faz esse cálculo com n adições e n multiplicações:
p( x ) = (((KK (a n x + a n −1 ) x + K + a 2 ) x + a1 ) x + a 0
123
n −1 parênteses
Ex.: p( x ) = 3 x 5 + 4 x 4 − 2 x 3 − x 2 + 3 x − 4 = (((( 3 x + 4 )x − 2 )x − 1 )x + 3 )x4
8) Se p( x ) é de grau n , então existe único polinômio de grau n − 1, q( x ) , tal que
p ( x ) = ( x − α ) ⋅ q( x ) + p (α ) .
Se α é raiz de p( x ) então p (α ) = 0 .
É o algoritmo de Briot-Ruffini utilizado para Deflacionar Raízes.
Ex.
p ( x ) = x 3 − 7 x 2 + 16 x − 10
(
) e p(1) = 0
α 2 = 2 ⇒ ( x − 2 )(x 2 − 5 x + 6 )+ 2 e
p (2 ) = 2
α 3 = −3 ⇒ ( x + 3 )(x 2 − 10 x + 46 )− 148 e
p(− 3) = −148
α 1 = 1 ⇒ ( x − 1) x 2 − 6 x + 10
2.1.1-Enumeração das Raízes
Enumerar as raízes de p( x ) é dizer quantas são as raízes e se positivas, negativas ou
complexas.
Regra de Descartes ou Regras de Sinais
“O número de raízes reais positivas de p( x ) é 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
multiplicidade “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 p (− x )
Regra de Huat
Se p (0 ) ≠ 0 e para algum k, a k2 ≤ a k −1 × ak +1 então p( x ) possui raízes complexas.
Regra da Lacuna
Se p (0 ) ≠ 0 e para algum K, a k = 0 e a k −1 × ak +1 > 0 , então p ( x ) tem raízes complexas.
11
E.B.Hauser – Cálculo Numérico
Ex.: p( x ) = 3 x 5 + x 4 − 5 x 3 − x − 1
•
p( x ) pode ter: 1 raiz real
positiva, 2 raízes negativas e 2 complexas;
ou
• 1 raiz real positiva, nenhuma negativa, e 4
complexas.
2.1.2-Localização das raízes de p(x)
Localizar as raízes de p( x ) é determinar a região do plano que contém todas as raízes.
Cota de Cauchy:
Toda raiz α (real ou complexa) de p( x ) satisfaz α ≤ β .
Onde β = lim x k , x0 = 0 e
k →∞
x k +1 =
a n −1 n −1 a n − 2 n − 2
a
a
xk + K + 1 xk + 0
xk +
an
an
an
an
Ex.: p( x ) = x 4 − 3x 3 + 3,37 x 2 − 1,68 + 0,3136 x − 1
x k +1 = (3 x + 3,37 x + 1,68 x k + 0,3136)
x0 = 0
e
3
k
2
k
x1 = 0,748331477355
1
4
x15 = 3,93973077126
x 2 = 1,47358396156
x16 = 3,94696658200
x3 = 2,10693957916
x17 = 3,95190592045
x 4 = 2,61655565303
x18 = 3,9552764745
x5 = 3,00483390194
x19 = 3,9575796715
M
M
∴α ≤ 4
12
2 - Resolução de Equações Algébricas e Transcendentes
2.2- Separação de Raízes Reais de f(x)=0
a) Métodos Gráficos: Utiliza-se um dos seguintes processos:
esboçar gráfico da função f (x ) e localizar as abcissas dos pontos onde a curva
i)
intercepta o eixo dos x.
ii)
de f ( x ) = 0 obter uma equação equivalente f1 ( x ) = f 2 ( x ) . Localizar no mesmo eixo
cartesiano
os
pontos
r
onde
as
duas
curvas
se
interceptem:
f1 (r ) = f 2 (r ) ⇒ f1 (r ) − f 2 (r ) = 0 ⇒ f (r ) = 0
b) Método Analítico: Seja f ( x ) continua no intervalo [a, b] . Se f (a ) ⋅ f (b ) < 0 , então
existe pelo menos uma raiz de f em ( a , b ) . (Se o sinal de f ' é constante em ( a , b ) a
raiz é única nesse intervalo).
Ex.: p( x ) = x 3 − 9 x + 3
a) Análise gráfica:
Logo, existem três raízes reais: r1 ∈ (− 4 , − 3) r2 ∈ (0 , 1) , r3 ∈ (2 , 3)
b) analiticamente:
x
−4 −3 − 3
−1 0 1
3
2 3
p ( x ) − 25 3 13,3923 11 3 − 5 − 7,3923 − 7 3
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.
13
E.B.Hauser – Cálculo Numérico
Ex1:
p( x ) = x 4 − 3 x 3 + 3,37 x 2 − 1,68 x + 0,3136
r1 = 0,7 e r2 ≅ 0,8 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 C , ε 1 , ε 2 , L são dados:
•
•
•
•
DIGSE( x k , x k + 1 ) ≥ C
f (xk ) ≤ ε 1
x k +1 − x k
≤ ε2
x k +1
k > L (número máximo de iterações)
2.3.1-Método da Bisseção ou Dicotomia (Algoritmo de quebra)
Seja f : [a, b] → ℜ continua e tal que f (a ) ⋅ f (b ) < 0 .
1) Calcula-se o ponto médio x m =
[a,
a+b
, dividindo-se [a, b ] em dois novos intervalos :
2
x m ] , [x m , b]
2) Se f ( x m ) ≠ 0 e:
i)
f (a ) ⋅ f ( x m ) < 0 então a raiz está em ( a , xm ) . Volta-se para (1)
ii)
f ( x m ) ⋅ f (b ) < 0 então a raiz está em ( xm ,b ) . 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 é 0,3 ⋅ DIGSE /passo, isto é, a cada 3 ou 4 passos ganha-se um DIGSE .
Ex: p( x ) = x 4 + 2 x 3 − 7,5 x 2 − 20 x − 11
a) Enumeração das raízes de p (x )
ℜ+
Regras de Huat e Lacuna
não aplicam
1
1
14
ℜ−
3
1
⊂ total
0
4
2
4
2 - Resolução de Equações Algébricas e Transcendentes
b) Cota
de
Cauchy:
x k +1 = 4 2 x k3 + 7,5 2k + 20 x k + 11, x 0 = 0
x1 = 1,82116028684
M
x 2 = 3,03080111616
x14 = 4,64729539211
x3 = 3,74256184946
x15 = 4,64784057829
x 4 = 4,14695320458
x16 = 4,64813826094
M
x17 = 4,64830079964
c) Separação de raízes
x
−5
p( x ) 276,5
− 4 − 3 − 2 −1
77
−1
8,5
r1
0
2
3
4
5
− 11 − 35,5 − 49 − 35 173 576,5
5
r2
1
r4
r3
d) Cálculo da raiz r4 ∈ (3, 4 ) utilizando o método da bissecção.
k
xk
f ( xk )
Ik
(3; 35 )
1
3,5
62 ,9375
(3; 3,25 )
2
3,25
25 ,00390625
(3; 3,125 )
3
3,125
9 ,660400391
(3; 3,0625 )
4
3,0625
2 ,871886353
(3,03125; 3,0625 )
− 0 ,405335188
5
3,03125
(3,03125; 3,046875 )
6
3,046875
1,1900454168
(3,03125; 3,039625 )
7
3,0390625
0 ,3883184832
(3,03515625; 3,0390625 )
8
3,03515625
− 0 ,0095143543
(3,03515625; 3,037109375 )
9
3,037109375
0 ,1891500395
(3,03515625; 3,0361328125 )
10 3,0361328125 0 ,0897548754
11 3,03564453125 0 ,0401045242 (3 ,03515625; 3 ,03564453125 )
12 3,03540039062 0 ,01529111512 (3 ,03515625; 3 ,03540039062 )
13 3 ,0352783203
0 ,002874148
M
Obtemos r4 ≈ x12 = 3,0354039062 com DIGSE ( x12 , x13 ) ≅ 4 ,096
2.3.2-Método da Iteração Linear
Consiste em escrever a equação f ( x ) = 0 na forma x = G ( x ) . Os pontos x* que
satisfazem a condição x * = G x * são ditos pontos fixos de G (x ) e geometricamente
representam os pontos de intersecção da reta y = x com a curva y = G ( x ) .
G ( x ) é dita função de iteração do método. Inicia-se a iteração com um valor x0
próximo da raiz, e as outras aproximações são dadas por:
xi +1 = G ( xi ), i = 0,1,2,K
( )
15
E.B.Hauser – Cálculo Numérico
A seqüência de aproximação xi , converge para
a solução x* da equação f ( x ) = 0 sob certas
condições .
A construção de G não é única. A escolha de
uma G apropriada é dita “problema do ponto fixo.
Ex. x 2 + x − 6 = 0 .
a) G1 ( x ) = 6 − x 2
b) G 2 = 6 − x ,
c ) G3 = − 6 − x ,
G4 =
6
,
x +1
G5 =
6
−1
x
Embora não seja preciso usar métodos numéricos para encontrar as duas raízes reais
α 1 = −3, e α 2 = 2 da presente equação, nota-se que:
i)
Tomando G1 e x0 = 1,5 , a seqüência {xi } não converge para 2.
xi +1 = G1 ( xi )
x1 = G1 ( x0 ) = 6 − (1,5) 2 = 3,75
x 2 = G1 ( x1 ) = 6 − (3,75) 2 = −8,0625
x3 = G1 ( x 2 ) = 6 − (−8,0625) 2 = −59,00390625
x 4 = G1 ( x3 ) = −3475,46095276
x5 = G1 (x 4 ) = −12078822,8342
M
ii)
Tomando G2 e x0 = 1,5 , a seqüência {xi } converge para 2.
x1 = G2 (x0 ) = 6 − 1,5 = 2,12132034356
x 2 = G2 ( x1 ) = 1,9694363804
x 3 = G2 ( x2 ) = 2,00762636454
x 4 = G2 ( x3 ) = 1,99809249923
x 5 = G2 ( x4 ) = 2,00047681835
x 6 = G2 ( x5 ) = 1,99988079186
x 7 = G2 ( x6 ) = 2,00002980181
M
Teorema da Convergência:
Seja α uma raiz isolada de f em [a, b] . Se
i) G e G’são contínuas em [a, b] ;
ii) G' ( x ) 〈1,∀x ∈ (a ,b ) ;
iii) x0 ∈ Ι e xk + 1 = G( xk ) ∈ (a ,b ), k = 0 ,1,2 ,... , então a seqüência { xk }, gerada por
xk +1 = G ( xk ) ,converge para α .
16
2 - Resolução de Equações Algébricas e Transcendentes
Ex: Utilizando o método da iteração linear calcule a
raiz de f ( x ) = e x + x 3 , com DIGSE (x k , x k +1 ) ≥ 5
f ( x ) = 0 ⇒ e x + x 3 = 0 ⇒ x 3 = −e x
⇒ x = − e = −e
Seja
3
G ( x ) = −e
x
x
3
x
3
x
1
G ' ( x ) = − e 3 ; G ' (− 1) ≅ 0,24 *
3
G ' (0 ) ≅ 0,33
G e G’ são continuas em [-1,0] e
G ' ( x ) < 1 ∀x ∈ [−1,0] .
xi
Logo , a seqüência gerada por xi +1 = −e 3 converge para α ∀x ∈ [−1,0] .
Seja x0 = −0,5
x1 = −0,846481725
x 6 = −0,772800243
x 2 = −0,754152577
x7 = −0,772904269
x3 = −0,777723518
x8 = − 0,772877469
f ( x9 ) = 0,000003188 e
x 4 = −0,771636903
x9 = −0,772884374
DIGSE ( x9 , x10 ) ≅ 5,34
x5 = −0,773204044
x10 = − 0,772882595
M
*G não tem Maximo nem mínimo local em [0,1], testa-se então só 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).
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 G ' (x ) , quanto menor este valor, maior será a
convergência.
17
E.B.Hauser – Cálculo Numérico
Obs.: Dada a equação f ( x ) = 0 , existem infinitas funções G(x) que são funções de iteração. A
forma geral destas funções é:
G ( x ) = x + A( x) ⋅ f ( x) ,
com a condição que em x*, ponto fixo de G(x), se tenha A( x * ) ≠ 0 . Temos que:
f ( x * ) = 0 ⇔ G ( x * ) = x * Com efeito:
( ⇒ ) Seja x * tal que f ( x * ) = 0 .
G ( x * ) = x * + A( x * ) ⋅ f ( x * ) = x * + A( x * ) ⋅ 0 = x *
(⇐) Se G ( x * ) = x * ⇒ x * + A( x * ) ⋅ f ( x * ) = x *
⇒ A( x * ) ⋅ f ( x * ) = 0 ⇒ f ( x * ) = 0 pois A( x * ) ≠ 0
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 f ( x) = 0 , tomamos a forma geral para G(x):
G ( x) = x + A( x) ⋅ f ( x) ⇒ G ' ( x) = 1 + A' ( x) ⋅ f ( x) + A( x) ⋅ f ' ( x) ⇒
G ' ( x * ) = 1 + A' ( x * ) ⋅ f ( x * ) + A( x * ) ⋅ f ' ( x * ) ⇒ G ' ( x * ) = 1 + A( x * ) ⋅ f ' ( x * ),
pois x * é ponto fixo de G ( x * = G ( x * ) ⇒ f ( x * ) = 0) .
Agora , G ( x * ) = 0 ⇒ 1 + A( x * ) ⋅ f ' ( x * ) = 0 ⇒ A( x * ) =
−1
.
f ' (x* )
1
f ( x)
e G ( x) = x −
.
f ' ( x)
f ' ( x)
Então dada f ( x ) , G ( x ) é tal que G ' x * = 0 , pois
Tomemos A( x) = −
( )
G ' (x ) = 1 −
[ f ' ( x)] − f ( x) ⋅ f " ( x) f ' ( x) ⋅ f " ( x)
=
e
[ f ' ( x)]2
[ f ' ( x)]2
2
f ( x * ) = 0 ⇒ G ' ( x * ) = 0 se f ' ( x * ) ≠ 0.
Portanto, iniciando-se a iteração com um valor x0 escolhido, a seqüência ( x k ) é
determinada por:
f ( xk )
xk +1 = xk −
, k = 0,1,2,... f ' ( xk ) ≠ 0
f ' ( xk )
Geometricamente , conforme podemos observar na próxima figura, dado x k , x k +1 é
abscissa do ponto de intersecção da reta tangente à curva f (x ) em ( x k , f ( x k )) e o eixo dos “x”.
Assim,
f ( xk )
f ( xk )
tgθ =
= f ' ( x k ) ⇒ x k +1 = x k −
x k − x k +1
f ' ( xk )
18
2 - Resolução de Equações Algébricas e Transcendentes
f ( xk )
f ( xk + 1 )
θ
xk + 1
xk
x k +1
Convergência: (é trabalhoso mostrar que G' ( x ) < 1 ).
O método de Newton-Raphson converge se:
f ( x) f " ( x)
G ' ( x) =
< 1 ⇒ f ( x) f " ( x) < ( f ' ( x)) 2 .
2
( f ' ( x))
Para raízes simples a convergência é quadrática e para raízes duplas ou triplas é linear.
Escolha do ponto inicial: Seja α ∈ ( a ,b ) raiz de f .
f (a) ⋅ f " (a) > 0 ⇒ x0 = a
Se
f (b) ⋅ f " (b) > 0 ⇒ x0 = b
(a + b)
.
Caso contrário, pode-se considerar x0 =
2
Ex.: 1) Estimar o valor da única raiz real de f ( x ) = 2 x + ln x , utilizando Newton-Raphson.
f ( x ) = 2 x + ln x
x,
19
x−
2 x + ln x
1
2−
x
E.B.Hauser – Cálculo Numérico
xk + 1 = x k −
x0 = 0 ,5
2 xk + ln xk
1
2+
xk
x1 = 0 ,42
x2 = 0 ,42699599
x3 = 0 ,426302751
x4 = x 3
Logo a aproximação para a raiz é α = 0 ,426302751 com 9 dígitos significativos corretos.
2) Calcular a raiz r4 ∈ [3,4] do polinômio dado anteriormente: p( x) = x 4 + 2 x 3 − 7,5 x 2 20 x + 11 .
p(3) ⋅ p" (3) < 0 e p (4) ⋅ p" (4) > 0 ⇒ x 0 = 4
x k +1
x k4 + 2 x k3 − 7,5 x k2 − 20 x k − 11
= xk −
4 x k3 + 6 x k2 − 15 x k − 20
x0 = 4
x1 = 3,36397059
x 2 = 3,08982331
Uma aproximação para a raiz é r4 = 3,03524990
com 9 dígitos significativos e 5 iterações.
x3 = 3,03709653
x 4 = 3,03525211
x5 = 3,03524990
x 6 = x5
Obs:
O Método de Newton pode divergir devido ao uso da tangente, oscilando indefinidamente.
20
2 - Resolução de Equações Algébricas e Transcendentes
Isto acontece quando:
i)
Não há raiz real
Ocorre simetria de f ( x ) em torno de α
ii)
O valor inicial x0 está longe da raiz exata, fazendo que outra parte da função
iii)
prenda 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
f ' ( x ) a cada iteração.
O método da secante contorna este problema, substituindo a derivada pelo quociente das
diferenças:
f ( x k ) − f ( x k −1 )
f ' ( xk ) ≅
x k − x k −1
onde xk e x k −1 são duas aproximações para a raiz de f ( x ) . A formula iterativa é:
( x k − x k −1 ) ⋅ f ( x k )
f ( x k ) − f ( x k −1 )
Geometricamente, a partir das aproximações para a raiz de xk e x k +1 , o ponto x k +1 é dado pela
abscissa do ponto de intersecção do eixo Ox e da reta secante que passa
por
( x k −1 , f ( x k −1 )) e ( x k , f ( x k )) .
x k +1 =
Características do método da secante:
9 Garante a convergência para toda função continua
9 Pode divergir se f ( x k ) ≅ f ( x k −1 )
9 Convergência mais lenta que o Método de Newton e mais rápida que Bissecção e
Iteração Linear.
9 São necessárias duas aproximações da raiz a cada iteração.
Ex: p( x ) = x 3 − 5 x 2 + 17 x + 21 e α ∈ ( −1,0 )
x k +1
( x k − x k −1 ) ⋅ ( x k3 − 5 x k2 + 17 x k + 21)
= xk − 3
( x k − 5 x k2 + 17 x k + 21) ⋅ ( x k3−1 − 5 x k2−1 + 17 x k −1 + 21)
x0 = −1, x1 = 1
x1 = −0,888888888889
x 2 = −0,960142348759
x3 = −0,931787651586
x 4 = −0,932112394706
x5 = −0,93211485688
x6 = −0,932114856662
x7 = x6
21
E.B.Hauser – Cálculo Numérico
Exercícios
1) Uma partícula é arremessada verticalmente, a partir do solo, com uma velocidade inicial vo
.Desprezando a resistência do ar, supomos que a posição p partícula é dada por:
g
d ( t ) = vo t − t 2 ,
2
onde g é aceleração da gravidade. Determinar a altura máxima atingida pela partícula e o
instante em que ocorreu.
2) Uma corrente oscilante num circuito elétrico é descrita por I = 9e − t sen( 2 π t ) , t em
segundos. Determinar todos os valores positivos de t para os quais I = 3.5.
3) A concentração de uma bactéria poluente num lago é descrita por
C = 70 e − 1.5t + 2.5e − 0.075t
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 = 8 e − kt cos wt
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 3 − 2 x 2 + 3x − 5
b) p(x) = x5 - 15,5x4 +77,5x3 - 155x2 +124x -31
c) p(x) = x 4 − 121x 3 + 2247 x 2 − 15043x + 34300
d) p(x) = x 4 − 115x 3 + 1575x 2 − 7625x + 12500
e) p(x) = x 4 − 3x 3 + 3.37 x 2 − 168
. x + 0.3136
4
3
f) p(x) = x -11,101x +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
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:
22
2 - Resolução de Equações Algébricas e Transcendentes
a)
b)
c)
d)
f(x) = x2 + ln x
f(x) = x + 2 cos x
f(x) = 2x + ln x
f(x) = cos x + ln x + x
2
e) f(x)= x + e − Bx 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
23
E.B.Hauser – Cálculo Numérico
4) t = 0 .3151660803
5-a)
r+ r3
0
1
0
¢
0
2
T
3
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)
f)
g)
h)
i)
Raízes: R1=0.7(multiplicidade 2), R2=0.8(multiplicidade 2)
Raízes: R1=0,001 R2=0,1 R3=1 R4=10
Raízes: R1=-0,5(multiplicidade 4), R2=1,2(multiplicidade 5)
Raízes: R1= 0.002251681490, R2 = 0.6822607762 , R3= 2081.245488
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
e) Existe única raiz de f em (-1,0)
24
d) R ≅ .2875182754
3. Sistemas de Equações Lineares
O sistema com n equações lineares e n variáveis
a11 x1 + a12 x 2 + a13 x 3 + L + a1n x n = b1
a 21 x 2 + a 22 x 2 + a 23 x 3 + L + a 2 n x n = b 2
M
M
M
M
M
a n1 x1 + a n 2 x 2 + a n3 x 3 + L + a nn x n = b n
pode ser representado matricialmente por AX = B , onde
⎡ a 11 a 12
⎢
a 21 a 22
A=⎢
⎢ M
M
⎢
⎣ a n1 a n 2
K
K
K
a 1n ⎤
⎥
a 2n⎥
,
M ⎥
⎥
a nn ⎦
⎡ x1⎤
⎢ ⎥
x2
X =⎢ ⎥ ,
⎢M⎥
⎢ ⎥
⎣ x n⎦
⎡ b1 ⎤
⎢ ⎥
b2
B=⎢ ⎥
⎢M⎥
⎢ ⎥
⎣b n ⎦
e
A é a matriz dos coeficientes, X é o vetor das incógnitas e B o vetor dos termos independentes.
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:
⎧ x + 0 ,98 y = 4 ,95
(a) ⎨
⎩x + y = 5
⎛ 2 ,5 ⎞
tem solução exata x = ⎜⎜ ⎟⎟
⎝ 2 ,5 ⎠
⎧ x + 0 ,99 y = 4 ,95
(b) ⎨
⎩x + y = 5
⎛ 0 ,0 ⎞
tem solução exata x = ⎜⎜ ⎟⎟
⎝ 5 ,0 ⎠
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.
Resolver o sistema significa determinar a intersecção das duas retas. Três casos são possíveis:
R1
R1
R2
R2
R1
Bem condicionado
det ≠ 0
retas concorrentes
Não tem solução.
det = 0
retas paralelas
Mal condicionado
det ≅ 0
(perturbação em ℜ 2 )
25
R2
E.B.Hauser – Cálculo Numérico
Exemplo2:
⎧ x + 5 y = 17
(a) ⎨
⎩1,5 x + 7 ,501 y = 25 ,503
⎛2⎞
tem solução exata: x = ⎜⎜ ⎟⎟
⎝3⎠
⎧ x + 5 y = 17
(b) ⎨
⎩1,5 x + 7 ,501 y = 25 ,501
⎛ 12 ⎞
tem solução: x = ⎜⎜ ⎟⎟
⎝1⎠
3.2-Medidas de Condicionamento
O determinante normalizado da matriz dos coeficientes A é dado por
det A
2
NORM A =
onde α k = ak21 + ak22 + L akn
, para k = 1, 2, ..., n
α 1α 2 Lα n
Norm A ∈ (− 1 , 1) e quanto mais afastado de ± 1 (isto é, quanto mais próximo de zero) mais mal
condicionado é a matriz A.
Retomando o Ex2:
1
5 ⎞ α 1 = (1 + 25 ) 2 = 5 ,09901951359
⎛ 1
⎟⎟
2) A = ⎜⎜
⎝ 1,5 7 ,501⎠
(
α 2 = 1,5
2
)
1
2 2
+ 7 ,501
det A = 7 ,501 − 5 ⋅ 1,5 = 0 ,001
norm A =
= 7 ,649550985358
0 ,001
0 ,001
=
= 0 ,00002563773874
α 1 ⋅ α 2 39 ,0050000128
Agora, como pode ser medido o condicionamento de um sistema linear?
Dado um SEL AX = B , o seu número de condicionamento é dado por:
Cond ( A ) = A
A−1 .
Quanto maior for Cond ( A ) , mais sensível será o sistema linear.
n
Utilizamos A = A ∞ = max ∑ aij , a norma do máximo das linhas.
1≤ i ≤ n i = 1
Ex:
5 ⎞
⎛ 1
⎛ 7501 − 5000 ⎞
⎟⎟ , A −1 = ⎜⎜
⎟⎟
A = ⎜⎜
⎝ 1,5 7 ,501 ⎠
⎝ − 1500 1000 ⎠
A = 9 ,001 , A − 1 = 12501
Cond ( A ) = A
A − 1 = 112521,501 ≅ 1 ⋅ 10 5
26
Sistemas de Equaçõe Lineares
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
aritmé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 Eliminação de Gauss
Algoritmo básico de Gauss: A solução de AX = B é calculada em duas etapas:
1º- Triangularização : Mediante operações elementares nas linhas, a matriz A é
transformada numa matriz triangular superior.
para k = 1,K ,n − 1
(indica a linha do pivô)
(indica a linha a transformar de A)
para i = k + 1,K ,n
− aik
m=
akk
aik = 0
(indica a coluna a transformar da linha i)
para j = k + 1,K , n
a ij = a ij + m ⋅ a kj
bi = bi − m ⋅ b k
Obs.: Se a ii = 0 é necessário trocar de linha, se possível.
Algoritmo:
2º- Retro-substituição: A triangularização transforma o sistema AX = B , no sistema
equivalente:
c 11 x 1 + c 12 x 2 + c 13 x 3 + L + c 1n x n = d 1
c 22 x 3 + c 23 x 3 + L + c 2 n x n = d 2
c 33 x 3 + L + c 3n x n = d 3
LLLLLLL
c nn x n = d n
cuja solução é dada por:
xn =
dn
,
cnn
xn − 1 =
(d n −1 − cn −1,n xn ) ,
an −1,n −1
d − (c12 x2 + c13 x3 + L + c1n xn )
x1 = 1
a11
Teorema: O método de Gauss produz sempre a solução exata do sistema AX = B (utilizando
precisão infinita) se det A ≠ 0 e as linhas de A forem permutadas quando aii = 0 .
27
E.B.Hauser – Cálculo Numérico
Quando é utilizada aritmética do ponto flutuante, erros de arredondamento podem
comprometer a solução obtida.
Ex.: Em F (10 ,3 ,−98 ,98 ) , com arredondamento para número mais próximo de máquina “ox”, a
⎧0 ,0002 x + 2 y = 5
elminação de Gauss aplicada ao sistema ⎨
produz x = 0 e y = 2 ,5 o que não
2x + 2 y = 6
⎩
satisfaz a segunda equação do sistema.
(Obs.: multiplicador m = −10.000 por L2 = L2 + L1 (m ) )
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:
3º pivô
1º pivô
2º pivô
Aplicando a técnica ao último exemplo
2x + 2 y = 6
⎧
, em F (10 ,3 ,−98 ,98 ) , com
⎨
⎩0 ,0002 x + 2 y = 5
arredondamento para número mais próximo de máquina “ox”, obtemo x = 0 ,5 e y = 2 ,5 .
(Obs.: Neste caso o multiplicador é menor m= - 0,0001)
Método de Gauss-Jordan (Matriz Inversa)
A solução do SEL AX = B é calculado utilizando X = A −1 B se det A ≠ 0 .
Obs.: Ver exercício 9.
Métodos Iterativos
Os sistemas lineares de grande parte são, em geral, esparsos (muito elementos aij = 0 ). Os
métodos diretos não preservam a esparsidade, enquanto que os métodos iterativos sim, além de
apresentarem relativa insensibilidade ao crescimento dos erros de arredondamento.
No sistema original AX = B , supondo aii ≠ 0 ,i = 1,K ,n , o vetor X é isolado mediante a
separação de diagonal principal:
28
Sistemas de Equaçõe Lineares
1
(b1 − a12 x2 − a13 x3 − K − a1n xn )
a11
1
(b2 − a21 x1 − a23 x3 − K − a2 n xn )
x2 =
a22
M
M
1
(bn − an1 x1 − an2 x2 − K − an ,n −1 xn −1 )
xn =
ann
x1 =
Metodo de Gauss-Jacobi
Dada a aproximação inicial X0, o processo iterativo produz aproximações sucessivas
X 1 , X 2 , K, X k ,K , obtidas de:
x1(k +1) =
1 ⎛
(k )
⎜ b1 − a12 x 2(k ) − a13x 3 − K − a1n x n (k ) ⎞⎟
⎠
⎝
a11
1
b 2 − a 21x1(k ) − a 23x 3(k ) − K − a 2n x n (k )
x (2k +1) =
a 22
M
M
1
(k ) ⎞
⎛⎜ b − a x (k ) − a x (k ) − K − a
x (nk +1) =
n ,n −1x n −1 ⎟⎠
n2 2
n
n1 1
⎝
a nn
)
(
Método de Gauss-Seidel
Para X0 dado:
x1(k + 1) =
1
a11
1
x2(k + 1) =
a22
1
x3(k + 1) =
a33
M
1
xn(k + 1) =
ann
(b − a
(b − a
(b − a
(k ) − a x (k ) − K − a x (k )
1n n
13 3
1
12 x2
2
21 x1
3
31 x1
(b
)
(k + 1) − a x (k ) − K − a x (k )
23 3
2n n
)
(k + 1) − a x (k + 1) − a x( k ) − K − a x (k )
3n n
32 2
34 4
M
(k + 1) − a x (k + 1) − K − a
(k + 1)
n ,n −1 xn − 1
n2 2
n − an1 x1
)
)
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 AX = Y , se a matriz A é Diagonalmente Dominante, isto é,
se a ii > ∑ a ij
j ≠i
∀ i , então tanto o método de Jacobi como o de Gauss-Seidel gera uma
seqüência (X (k ) ) convergente para a solução do sistema dado, independente da escolha da
aproximação inicial X (0 ) .
29
E.B.Hauser – Cálculo Numérico
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.
⎧ x1
⎪
⎪
⎨
⎪ 2 x1
⎪⎩0 ,5 x1
+ 4 x2
− x3
− x4
=− 2
+ 2 x4
+ x4
=−3
= 1
+ x3
= 1,5
Reordenamos a fim de satisfazer ao critério de convergência.
⎧ 2 x1
⎪ x
⎪
1
⎨
0
,
5
x
1
⎪
⎪⎩
+ x4
+ 4 x2
= 1
=− 2
− x4
+ x3
− x3
+ 2 x4
2
>
1
= 1,5
4
1
>
>
1 + −1
0 ,5
=−3
2
>
−1
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 (X (k ) )
convergente para a solução exata.
Gauss-Jacobi
Fórmulas iterativas:
(k ) ⎞
⎛
⎟
⎜1 − x
4 ⎠
(
k + 1) ⎝
= 0 ,5 − 0 ,5 x4(k )
x1
=
2
x2(k + 1) =
(k )
(k )
( −2 − x1
+ x4
4
(
x3(k + 1) = 1,5 − 0 ,5 x1(k )
(k )
( −3 + x3
x (k + 1) =
4
2
)
)
= −0 ,5 − 0 ,25 x1(k ) + 0 ,25 x4(k )
)
= −1,5 + 0 ,5 x3(k )
Aproximação inicial: X (0 ) = 0 .
30
Sistemas de Equaçõe Lineares
k
0
1
2
3
4
5
6
7
8
9
10
11
:
40
x1
x2
x3
x4
0
0,5
1,25
0,875
0,9375
1,03125
0,984375
0,9921875
1,00390625
0,998046875
0,9990234375
1.000488282
:
1
0
0,5
-1,0
-1,0
-0,9375
-1,0000
-1,0000
-0,9921875
-1,0000
-1,0000
-0,999023437
-1.0000000
:
-1
0
1,5
1,25
0,875
1,0625
1,03125
0,984375
1,0078125
1,00390625
0,998046875
1,000976563
1.000488281
:
1
0
-1,5
-0,75
-0,875
-10625
-0,96875
-0,984375
-1,0078125
0,9960375
-0,998046875
-1,000976563
-0,999511718
:
-1
na 12a. iteração consegue-se x (12 ) − x 11
< 10 −5
Gauss-Seidel
Fórmulas iterativas:
x1(k + 1) = 0 ,5 − 0 ,5 x4(k )
x2(k + 1) = −0 ,5 − 0 ,25 x1(k + 1) + 0 ,25 x4(k )
x3(k + 1) = 1,5 − 0 ,5 x1(k + 1)
x4(k + 1) = −1,5 + 0 ,5 x3(k + 1)
k
0
1
2
3
4
5
:
12
x1
0
0,5
0,9375
0,9921875
0,9990234375
0,999877929
:
1
x2
0
-0,625
0,953125
-0,994140625
-0,9992675781
-0,999908447
:
-1
x3
0
1,25
1,03125
1,00390625
1,000488281
1,000061035
:
1
x (5 ) − x (4 ) < 10 − 5
31
x4
0
-0,875
-0,984375
0,99806875
-0,999755895
-0,999969482
:
-1
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 temperaturas nos seis
vértices interiores do reticulado da figura. A temperatura num vértice é aproximadamente igual à
média dos quatro vértices vizinhos mais próximos (à esquerda, acima, à 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
temperaturas T1, T2 , ..., T6.
b) Resolver o sistema utilizando o sistema MapleV.
20o
10o
1
20o
20o
1
2
3
4
5
6
10o
40o
40o
20o
20o
20o
2) Num experimento num túnel de vento, a força sobre um projétil devido à resistência do ar
foi medida para velocidades diferentes.
velocidade 0
2
4
6
8 10
força
0 2.90 14.8 39.6 74.3 119
Expressar a força como função da velocidade aproximando-a a um polinômio de quinto
grau:
f ( v ) = a o + a 1 v + a 2 v2 + a 3 v3 + a 4 v4 + a 5 v5
Verificar a validade de f (v) 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.
⎧ x1 + 0 ,5 x2 + 0 ,33 x3 = 1,83
⎪
3) Considere o sistema (#) ⎨0 ,5 x1 + 0 ,33 x2 + 0 ,25 x3 = 1,08 .
⎪0 ,33 x + 0 ,25 x + 0 ,2 x = 0 ,78
1
2
3
⎩
a) (#) é bem ou mal condicionado? Porque? O que isso significa?
b) Resolver (#) pelo método de Gauss sem pivotamento.
32
Sistemas de Equaçõe Lineares
⎧ x1 + 1 / 2 x2 + 1 / 3 x3 + 1 / 4 x4 + 1 / 5 x5 = 137 / 60
⎪
⎪⎪1 / 2 x1 + 1 / 3 x2 + 1 / 4 x3 + 1 / 5 x4 + 1 / 6 x5 = 87 / 60
4) Idem ao 3 para ⎨1 / 3 x1 + 1 / 4 x2 + 1 / 5 x3 + 1 / 6 x4 + 1 / 7 x5 = 459 / 420
⎪1 / 4 x + 1 / 5 x + 1 / 6 x + 1 / 7 x + 1 / 8 x = 743 / 840
1
2
3
4
5
⎪
⎪⎩1 / 5 x1 + 1 / 6 x2 + 1 / 7 x3 + 1 / 8 x4 + 1 / 9 x5 = 1879 / 2520
5) Resolver por Eliminação de Gauss com pivotamento parcial.
⎧3 x1 − 5 x2 + 6 x3 + 4 x4 − 2 x5 − 3 x6 + 8 x7 = 47
⎪ x + x + 9 x + 15 x + x − 9 x + 2 x = 17
2
3
4
5
6
7
⎪ 1
⎪2 x1 − x2 + 7 x3 + 5 x4 − x5 + 6 x6 + 11x7 = 24
⎪
a) ⎨− x1 + x2 + 3 x3 + 2 x4 + 7 x5 − x6 − 2 x7 = 8
⎪4 x + 3 x + x − 7 x + 2 x + x + x = 13
2
3
4
5
6
7
⎪ 1
⎪2 x1 + 9 x2 − 8 x3 + 11x4 − x5 − 4 x6 − x7 = −10
⎪
⎩7 x1 + 2 x2 − x3 + 2 x4 + 7 x5 − x6 + 9 x7 = 34
⎧− x1 + 2x 2 + 42x 3 = 83
⎪
b) ⎨72x1 − 41x 2 − 14x 3 = 44
⎪35x + 10x − 5x − = 25
2
3
⎩ 1
⎧2 x1 + x2 − 0.1x3 + x4 = 2.7
⎪0.4 x + 0.5 x + 4 x − 8.5 x = 21.9
⎪
1
2
3
4
c) ⎨
⎪0.3 x1 − x2 + x3 + 5.2 x4 = −3.9
⎪⎩ x1 + 0.2 x2 + 2.5 x3 − x4 = 9.9
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 :
⎧ x1 − x2 + x3 = 5
⎧ x1 + x2 = 0
⎪
⎪
a) ⎨ x1 + x2 = 4
e
b) ⎨ x1 − x3 = 0
⎪− 2 x + x = 2
⎪x + x = 0
2
3
3
⎩
⎩ 2
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
⎧x + y = 3
⎧ x − 3 y = −3
a) ⎨
b) ⎨
⎩ x − 3 y = −3
⎩x + y = 3
7) Resolver pelo Método de Gauss-Seidel. Apresentar as fórmulas iterativas e uma garantia de
convergência (se possível).
⎧3 x1 − 5 x2 + 47 x3 + 20 x4 = 18
⎧8 x1 + 2 x2 + 3 x3 = 30
⎪11x + 16 x + 17 x + 10 x = 26
⎪
⎪ 1
2
3
4
a) ⎨
b) ⎨ x1 − 9 x2 + 2 x3 = 1
⎪2 x + 3 x + 6 x = 31
⎪56 x1 + 22 x2 + 11x3 − 18 x4 = 34
2
3
⎩ 1
⎪⎩17 x1 + 66 x2 − 12 x3 + 7 x4 = 82
33
E.B.Hauser – Cálculo Numérico
⎧4 x + y 2 + z = 11
⎪
⎪
8) Resolver o sistema de equações algébricas não lineares: ⎨ x + 4 y + z 2 = 18
⎪ 2
⎪⎩ x + y + 4 z = 15
⎧2 x1 + x2 + 7 x3 = b1
9) Resolver ⎪⎨ x1 + 3 x2 + 2 x3 = b2 para:
⎪5 x + 3 x + 4 x = b
2
3
3
⎩ 1
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.
⎡7 0
⎡− 2 − 3 − 1 − 2⎤
⎢3 − 2
⎥
⎢− 1 0
⎢
1
2
−
⎥
a) A = ⎢
b) A = ⎢0 5
⎢− 3 − 1 − 4 1 ⎥
⎢
⎥
⎢
⎢0 1
⎣ − 2 2 − 3 − 1⎦
⎢⎣1 0
8 0.5
0 ⎤
0 0
0 ⎥⎥
3 0
3 ⎥
⎥
1 9 0.25⎥
3 1
0 ⎥⎦
11) Qual o Resíduo produzido pela solução aproximada X’= [ -3 4 0]T de
⎧0.24 x1 + 0.36 x2 + 0.12 x3 = 0.84
⎪
⎨0.12 x1 + 0.16 x2 + 0.24 x3 = 0.52
⎪0.15 x + 0.21x + 0.25 x = 0.64
1
2
3
⎩
Respostas:
1) Solução obtida utilizando o MapleV:
>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 MapleV:
>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-.7005208333e1*x^4+.2604166667e-2*x^5;
> validade:=[f(0),f(2),f(4),f(6),f(8),f(10)];
34
Sistemas de Equaçõe Lineares
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]
3-a) Mal condicionado. NORM A ≅ 0,000181. Uma pequena perturbação nos dados de entrada
pode causar uma grande alteração na solução.
b) A solução exata é X=[1 1 1]T
4 - NORM A ≅ 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
35
4.Interpolação Polinomial
4.1 Objetivo:
Dado um conjunto de n+1 pontos distintos ( xi , f ( xi )) , i = 0 ,1,...,n , queremos determinar o
polinômio p(x) talque
p( xi ) = f ( xi ), e para os demais pontos do intervalo
p( x ) ≅ f ( x ) , ∀x ∈ [xo ; xn ] . O polinômio p(x) é dito polinômio aproximante ou interpolador de
f(x) no intervalo [xo , xn ] .
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 x* ∈ [xo ; xn ] ;
xn
• Determinar aproximações para
∫xo f ( x)dx , substituindo a f(x) pelo polinômio
interpolador;
• Calcular uma aproximação para f ' ( x) para x ∈ [xo ; xn ] , substituindo f(x) por p(x).
4.2 Existência e Unicidade da Solução
Dados xi ∈ ℜ
p( xi ) = f ( xi ), ∀xi .
e
f ( xi ) ∈ ℜ ,
i = 0 ,1,...,n , procuramos
Seja p ( x) = ao + a1x + a2 x 2 + ... + an x n =
36
n
∑ ak x k
k =0
.
p( x ) ∈ Pn tal que
E.B.Hauser – Cálculo Numérico
Então de p( xi ) =
n
∑ ak xik = f ( xi ),
obtemos:
k =0
⎧a0 + a1 x0 + a2 x02 + ... + an x0n =
⎪
⎪a0 + a1 x1 + a2 x12 + ... + an x1n =
⎪
⎨
⎪..............................................
⎪
⎪⎩a0 + a1 xn + a2 xn2 + ... + an xnn =
f0
f1 ,
fn
o qual representa um sistema linear de ordem n+1 onde as n+1 incógnitas são os ak , k = 0 , n e a
matriz dos coeficientes é dada por:
⎡1 x0
x 02
. . x 0n ⎤
⎥
⎢
x 12
. . x 1n ⎥
⎢1 x1
⎢1 x
x 22
. . x 2n ⎥
2
⎥
⎢
A=
.
.
. . . ⎥
⎢.
⎥
⎢
.
.
. . . ⎥
⎢.
⎢.
.
.
. . . ⎥
⎥
⎢
2
xn
. . x nn ⎦
⎣1 x n
n −1 ⎡ n
⎤
det A = ∏ ⎢ ∏ ( x j − xi )⎥ . Como os pontos são
⎥
i = 0 ⎢⎣ j = i + 1
⎦
distintos, a diferença x j − xi será sempre diferente de zero, e então detA ≠ 0. Portanto o polinômio
interpolador existe e também é único.
De acordo com Vandermonde,
4.3Polinômio Interpolador de Newton Para Diferenças Finitas Ascendentes
Dados (xi , yi ) , yi = f ( xi ) i=0,1,2,...,n satisfazendo
x1 − x0 = x2 − x1 = .... = xn − xn −1 = h ,
isto é , xi + 1 − xi = h .
Para k = 1, 2, ...n, e i= 0, 1, 2, ..., n-k , a diferença finita de ordem k é dada por
∆k yi = ∆k − 1 yi + 1 − ∆k − 1 yi .
E, o polinômio interpolador de Newton para diferenças finitas ascendentes é dado por :
p( x ) = y o +
( x − xo )( x − x1 )L( x − x n − 1 ) n
( x − xo )( x − x1 ) 2
( x − xo )
∆y 0 +
∆ yo + L +
∆ yo
2
h
2! h
n! h n
37
4-Interpolação Polinomial
Exemplo1:
O alongamento de uma mola foi medido em função da carga aplicada. Obteve-se:
c arg a( kg )
2 4 6 8
alongamento( cm ) 1,0 2 ,5 5 ,0 6 ,3
Estimar o alongamento para o caso de ser aplicada uma carga de 7kg , utilizando o polinômio
interpolador de Newton para diferenças finitas.
Solução:
1) Tabela das diferenças finitas:
i
yi
xi
0
1
2
3
2
4
6
8
1
2.5
5
6.3
∆yi
∆2 yi
∆3 yi
1.5
1
-2.2
2.5
-1.2
-------------------1.3
-------------------- --------------------------------------- -------------------- --------------------
2) Polinômio interpolador:
(x − 2) .1,5 + (x − 2 )(x − 4 ) .1 + (x − 2 )(x − 4 )(x − 6 ) .( −2 ,2 )
p( x ) = 1 +
2
3! ( 2 )3
2! ( 2 )2
p( x ) = −0 ,04583333333 x 3 + 0 ,675 x 2 − 2 ,016666667 x + 2 ,7
3) Verificação de validade de p(x) : p( 2 ) = 0 ,9999999994 ( ≅ 1 )
p( 4 ) = 2 ,4999999999 ( ≅ 2 ,5 )
p( 6 ) = 5 e p(8)=6.3
4) Estimativa do alongamento ao se aplicar uma carga de 7kg:
O alongamento da mola neste caso é aproximadamente p( 7 ) = 5 ,9375 cm.
5) Análise gráfica do problema:
38
E.B.Hauser – Cálculo Numérico
4.4 Polinômio Interpolador de Newton Para Diferenças Finitas Divididas
Dados (xi , yi ) , yi = f ( xi ), i= 0, 1, 2, ..., n
qualquer, não necessariamente eqüidistantes.
os pontos xi podem ter um espaçamento
O polinômio de Newton para diferenças divididas é dado por:
2
p( x) = yo + ( x − xo )∆ y0 + ( x − xo )( x − x1 ) ∆ yo + ... + ( x − xo )( x − x1 )...( x − xn −1)∆ n yo ,
onde, para k = 1, 2, ..., n, e i= 0, 1, 2, ..., n-k a diferença dividida de ordem k é dada por
∆ yi =
k
∆
k −1
yi + 1 − ∆ k − 1 yi
xi + k − xi
Por exemplo, para o caso de n = 4:
i
xi
yi
∆ yi
∆ yi
2
∆ yi
0
x0
y0
y1 − y0
x1 − x0
∆ y1 − ∆ y0
∆ y1 − ∆ y0
y2 − y1
x2 − x1
∆ y 2 − ∆ y1
y3 − y 2
x3 − x 2
∆ y3 − ∆y2
-----------------------------------------------------------------------------
1
2
x1
x2
y1
y2
3
2
x2 − x0
3
x3
y3
y4 − y3
x4 − x3
4
x4
y4
---------------
2
x3 − x0
2
2
∆ y2 − ∆ y1
x3 − x1
x4 − x 2
4
∆ yi
x4 − x1
--------------------------------------------------------------------------------------------------------------------------------------
3
3
∆ y1 − ∆ y0
x4 − x0
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Observação: Para h constante , a relação entre diferenças divididas e finitas é dada por :
1
∆ k yi =
∆k yi .
k
k !.h
Exemplo:
x
f(x)
∆ yi
∆ yi
2
∆ yi
3
∆ yi
0
2
3
5
6
0
8
27
125
216
4
19
49
91
-----
5
10
14
---------
1
1
-------------
0
-----------------
39
4
4-Interpolação Polinomial
Exemplo2:
A velocidade do som na água varia com a temperatura conforme tabela:
temperatura( oC ) 86 93,3 98 ,9 104 ,4 110
velocidade( m / s ) 1,552 1,548 1,544 1,538 1,532
Solução:
1) Cálculo das diferenças divididas
x
y
∆ y0
86
93.3
98.9
104.4
110
1.552
1.548
1.544
1.538
1.532
-0.00054794
-0.00071428
-0.00109090
-0.00107142
---------------
2
∆ y0
-0.00001289
-0.00003393
0.175 x 10-5
-----------------------------
3
∆ y0
-0.114 x 10-5
0.213 x 10-5
-------------------------------------------
4
∆ y0
0.136 x 10-6
---------------------------------------------------------
2) Construção do polinômio:
p( x ) = 1.552 + ( x − 86 )( −0.00054794 ) + ( x − 86 )( x − 93.3 )( −0.00001289 ) + ( x − 86 )( x − 93.3 )( x − 98.9 )
(-0.114 x 10 −5 )+ ( x − 86 )( x − 93.3 )( x − 98.9 )( x − 104.4 ) (0.136 x 10 −6 )
Simplificando a expressão , encontramos
p(x) = 0.13666902 84 × 10 -6 x 4 − 0 . 0000534327 9926 x 3 + 0 . 0077947032 83 x 2 − 0 . 5036369194
3) Verificação da validade de p(x) calculado no item 2:
p( 86 ) = 1.55199999 ; p( 93.3 ) = 1.548 ;
p( 98.9 ) = 1.544 ;
p( 104.4 ) = 1.537999999 ; p( 110 ) = 1.532
4) Estimativa da velocidade do som quando a temperatura da água atinge 100 0 C , é
p( 100 ) ≅ 1.54293925 m/s
6) Análise gráfica do problema:
40
E.B.Hauser – Cálculo Numérico
4.5 Erro de Truncamento
E ( x) =
ϕ ( x)
(n + 1)!
f (n +1) (ξ ),ξ ∈ [x0 , xn ] com
ϕ ( x ) = ( x − x0 )( x − x1 )...( x − xn )
pois ∆ k yi =
1
k !.h
k
∆n + 1 f 0
( n + 1 )! h
n +1
= ( x − x0 )( x − x1 )...( x − xn )* ∆ n + 1 yi ,
∆k f i , 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
reserva e um controle populacional mostrou os dados:
t( anos )
0 1 3 5 7 10 13
quantidade de animais 50 60 73 77 76 73 69
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
concluir?
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
utilizando 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
determinado no ítem b.
Resposta:
0
0
0
0
0 ⎤
⎡1 0
⎢1 1
1
1
1
1
1 ⎥⎥
⎢
⎢1 3 3 2
33
34
35
36 ⎥
⎢
⎥
V = ⎢1 5 5 2
53
54
55
56 ⎥
⎢1 7 7 2 7 3 7 4 7 5 7 6 ⎥
⎢
⎥
2
10 3 10 4 10 5 10 6 ⎥
⎢1 10 10
⎢1 13 13 2 13 3 13 4 13 5 136 ⎥
⎣
⎦
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
41
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étrica. Numa experiência deste tipo, foram obtidos
os seguintes dados:
x ( voltagem − volts( V )) 30 35 40 47 50
y ( corrente − ampere( A )) 2 2 ,8 3 ,5 4 ,3 4 ,5
Estimar a corrente se a voltagem aplicada for de 43A usando o
Polinô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
esperanç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 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
42
E.B.Hauser – Cálculo Numérico
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.
idade em 2003
0 10 15 20 25 30 50 55 60 65 70
exp ectativa de vida( anos ,em2003 ) − Mulheres 75.2 67.5 62.6 57.8 53.0 48.3 30.1 26.0 22.1 18.4 15.0
idade em 2003
exp ectativa de vida ( anos , em 2003 ) − Homens
0
10
15
20
25
30
50
55 60
65 70
67 . 6 60 . 4 55 . 5 51 . 0 46 . 8 42 . 5 26 . 2 22 . 5 19 . 1 15 . 9 13 . 1
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)],
14..30)]);
> Estimativas_homens:= array( [ seq( [i, pexpvidah(i)], i
14..30)]);
⎡14
⎡14 63.56097596⎤
⎢
⎢
⎥
⎢15
⎢15 62.6000001 ⎥
⎢
⎢
⎥
⎢16
⎢16 61.6414611 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢17
⎢17 60.6831095 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢18
⎢18 59.7236187 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢19
⎢19 58.7625088 ⎥
⎢
⎢
⎥
⎢20
⎢20 57.8000000 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢21
⎢21 56.8368317 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
Estimativas_homens := ⎢⎢22
Estimativas_mulheres := ⎢⎢22 55.8740682 ⎥⎥
⎢
⎢
⎥
⎢23
⎢23 54.9129142 ⎥
⎢
⎢
⎥
⎢24
⎢24 53.954550 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢25
⎢25 53.000000 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢26
⎢26 52.050032 ⎥
⎢
⎢
⎥
⎢27
⎢27 51.105097 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢28
⎢28 50.165308 ⎥
⎢
⎢
⎥
⎢
⎢
⎥
⎢29
⎢29 49.230445 ⎥
⎢
⎢
⎥
⎢⎢
⎢⎢
⎥⎥
⎣30
⎣30 48.300001 ⎦
43
i =
=
56.46550812⎤
⎥
55.5000000 ⎥⎥
54.5559710 ⎥⎥
⎥
53.6353083 ⎥⎥
⎥
52.7375232 ⎥⎥
⎥
51.8602890 ⎥⎥
51.0000001 ⎥⎥
⎥
50.1523108 ⎥⎥
⎥
49.3126236 ⎥⎥
⎥
48.476507 ⎥
⎥
47.640028 ⎥⎥
⎥
46.800001 ⎥⎥
⎥
45.954138 ⎥⎥
45.101122 ⎥⎥
⎥
44.240601 ⎥⎥
⎥
43.373125 ⎥⎥
⎥
42.500001 ⎥⎦
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
1,1052
1,8221
3,0042
4,9530
8,1662
y (m3 )
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 operacional 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
distância apresentada na tabela abaixo é relativa à circunferência.
ALTURA (m)
DISTÂNCIA DO ALVO
O
1 SALTO (1500)
35
2O SALTO (1250)
25
O
3 SALTO (1000)
15
4O SALTO (750)
10
O
5 SALTO (500)
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.
Verificar 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
44
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 o K . 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)
3
600
71,72 . 10
94,72 . 10 3
800
118,4 . 10 3
1000
142,08 . 10 3
1200
165,76 . 10 3
1400
189,44 . 10 3
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) p( 0 ,12 ) ≅ 1,126904937
2)
A demanda mínima é 14 ,8632 Mw. e ocorre entre 3h e 4h da manhã.
A demanda máxima é 43,101 Mw. e ocorre entre 14h e 15h.
3) p( 850 ) = 11,4128
4)
5)
p( 80 ) = 13,46701783
p( 80 ) = 13,4512685
p( 1500 ) = 177618 ,5937
6) p(7250) = 15090,53234
45
5. Ajuste de Funções
Aplicação1: Os dados acima tabelados descrevem a intensidade da luz como uma função
da distância da fonte, I(d), medida num experimento.
1
Determinar I (d) ≅ Y (d) =
.
Ad 2 + Bd + C
d 30 35 40 45 50 55 60 65 70 75
I 0.85 0.67 0.52 0.42 0.34 0.28 0.24 0.21 0.18 0.15
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:
Pr essão 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5
Volume 2,85 2,27 1,91 1,60 1,43 1,28 1.15 1,04
Ajustar os dados tabelados a uma hipérbole do tipo: Vol(p) ≅ Y ( p ) =
46
A
B+ p
E.B.Hauser – Cálculo Numérico
5.1 Introdução
O ajustamento é uma técnica de aproximação. Conhecendo-se dados
experimentais , , ( xi , yi ) , i = 0 ,1,...,n , deseja-se obter a lei y = f ( x ) relacionando x com y.
Devido aos erros experimentais nos n+1 pares, ( xi , f ( xi )) , teremos em
geral f ( x1 ) + ε 1 ; f ( x2 ) + ε 2 ; ... ; f ( xn ) + ε n , 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
ajustamento 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
a) Função Linear (regressão linear) : Y ( x ) = a0 + a1 x , se ∆yi ≅ cte ou ∆ yi ≅ cte .
2
b) Parábola (ajuste quadrático): Y ( x ) = a0 + a1 x + a2 x 2 , se ∆2 yi ≅ cte ou ∆ yi ≅ cte.
p
c) Polinômio de grau p: se ∆ p yi ≅ cte ou ∆ yi ≅ cte.
d) Função Exponencial: Y ( x ) = ab x , se
d) FunçãoPotência: Y ( x ) = ax p ,
∆ log yi
≅ cte.
∆xi
∆ log yi
≅ cte
∆ log xi
Outros tipos de Funções Ajuste:
•
Y( x ) =
1
1
⇒ = a0 + a1 x ;
a0 + a1 x
y
∆ (1 / yi ) =
∆ (1 / yi )
≅ cte.
∆xi
•
Y( x ) =
x
x
⇒ = a0 + a1 x;
a0 + a1 x
y
∆ ( xi / yi ) =
∆ ( xi / yi )
≅ cte.
∆ ( xi )
•
47
5 - Ajuste de Funções
•
•
Y( x ) =
1
a0 + a1 x + a2 x
2
⇒
1
= a0 + a1 x + a2 x 2 ;
Y
2
Y ( x ) = ae bx + cx ; ln y = ln a + bx + cx 2 ,
∆ (1 / y i ) =
2
∆ ln yi + 1 − ∆ ln yi
xi + 2 − xi
∆ (1 / y i +1 ) − ∆ (1 / y i )
xi + 2 − xi
≅ cte.
≅ cte.
3
Exemplo: Considerando a tabela abaixo, como ∆ yi = 1,∀i , é adequado o ajuste dos dados
abaixo tabelados a um polinômio de grau 3.
xi f(xi)=yi
i
0
1
2
3
4
0
2
3
5
6
0
8
27
125
216
∆ yi
4
19
49
91
-----
2
∆ yi
5
10
14
---------
3
∆ yi
1
1
-------------
4
∆ yi
0
-----------------
5.3 - Determinação dos Parâmetros da Função de Ajuste
CRITÉRIO DOS MINÍMOS QUADRADOS
Seja Y = a0 + a1x + a2 x 2 + ... + a p x p a função de ajustamento. Dada uma
tabela com n+1 pontos (xi , yi ) , chamamos resíduo a diferença entre o valor de Yi da equação
de ajustamento e o valor tabulado de yi . Yi − yi = δ i , i = 0 ,1,...,n .
n
O critério dos mínimos quadrados estabelece que:
∑ δ i2 = mínimo .
i =0
Seja F (a0 , a1,...a p ) =
n
∑ (Yi − yi )2 . Para F ter valor mínimo , é preciso que
i =0
∂F
∂F
∂F
=0;
=0; . . . ;
=0
∂a0
∂a1
∂a p
48
E.B.Hauser – Cálculo Numérico
5.31 -Ajuste Linear
A função de ajuste terá a forma Y ( x ) = a0 + a1 x . Pelo critério dos
Quadrados é necessário que :
n
n
F = ∑ (Yi − yi )2 = ∑ (a0 + a1xi − yi ) 2
i =0
i =0
Mínimos
deve ser mínimo .
Sendo F uma função de duas variáveis, a0 e a1 , o menor valor de F será obtido
δF
δF
através de :
= 0;
= 0 e assim:
δa0
δa1
n
δF
= 0 , ∑ 2(a0 + a1xi − yi ) = 0 e
δa0
i =0
n
δF
= 0 , ∑ 2(a0 + a1xi − yi )xi = 0 .
δa1
i =0
Construímos o sistema de duas equações e duas variáveis:
⎧ n
⎪2 ∑ (a0 + a1xi − yi ) = 0
⎪⎪
i =0
⎨
⎪ n
⎪2 ∑ (a0 + a1x1 − yi )xi = 0
⎪⎩ i =0
ou
n
n
⎧ n
⎪ ∑ yi = ∑ a0 + ∑ a1xi
⎪⎪
i =0
i =0
i =0
.
⎨
n
n
⎪ n
2
⎪ ∑ xi yi = ∑ a0 xi + ∑ a1xi
⎪⎩i =0
i =0
i =0
Obtemos:
n
n
⎧
⎪(n + 1)a0 + a1 ∑ xi = ∑ yi
⎪⎪
i =0
i =0
⎨
n
n
n
⎪
2
⎪a0 ∑ xi + a1 ∑ xi = ∑ xi yi
⎪⎩ i =0
i =0
i =0
Resolvendo-se este último sistema linear são obtidos os valores de a0 e a1 e
assim determina-se a função de ajuste : Y = a0 + a1x .
49
5 - Ajuste de Funções
Exemplo2 : Dada a tabela , achar a equação da reta que se ajusta usando o método dos
Mínimos Quadrados.
i
0
1
2
3
4
5
∑
xi
0
1
2
3
4
5
15
yi
xi yi
2
3
5
5
5.5
8
28.5
0
3
10
15
22
40
90
xi2
0
1
4
9
16
25
55
Yi
(Yi − yi )2
2.07142857
3.14285714
4.21428571
5.2857143
6.3571429
7.4285714
---------------
0.005102041
0.02040816
0.61734694
0.08163266
0.73469755
0.32653061
1.78571440
Seja Y ( x) = a0 + a1x a função que ajusta os dados .
Os parâmetros a0 e a1 constituem a solução do sistema :
∑ yi = (n + 1)a0 + a1∑ xi
∑ xi yi = a0 ∑ xi + a1∑ xi2
⎧6a0 + 15a1 = 28.5
⎧a0 = 2.071428572
⇒ ⎨
⇒ ⎨
⎩a1 = 1.071428571
⎩15a0 + 55a1 = 90
Logo, Y = 2.071428571 + 1.071428571x .
50
E.B.Hauser – Cálculo Numérico
5.3.2 - Ajuste Quadrático
Seja Y = a0 + a1 x + a2 x 2 a função de ajuste. Pelo critério dos Mínimos
Quadrados :
(
i =0
n
n
F ( a0 ,a1 ,a2 ) = ∑ (Yi − yi )2 = ∑ a0 + a1 xi + a2 xi2 − yi
i =0
)2
δF δF δF
=
=
= 0.
δa0 δa1 δa2
assume valor mínimo se
Assim, os parâmetros a0 ,a1 e a2 são obtidos resolvendo:
∑ yi = (n + 1) a0 + a1∑ xi + a2 ∑ xi2
∑ xi yi = a0 ∑ xi + a1∑ xi2 + a2 ∑ xi3
∑ xi2 yi = a0 ∑ xi2 + a1∑ xi3 + a2 ∑ xi4
⎧
⎪
⎪⎪
⎨
⎪
⎪
⎪⎩
Exemplo3 : Encontre a expressão do polinômio de 2o grau que se ajusta aos dados da tabela
abaixo.
i
xi
yi
xi yi
0
1
2
3
4
5
-2
-1
0
1
2
3
3
-0.01
0.51
0.82
0.88
0.81
0.49
3.5
0.02
-0.51
0
0.88
1.62
1.47
3.48
∑
xi2
4
1
0
1
4
9
19
xi2 yi
-0.04
0.51
0
0.88
3.24
4.41
9
xi3
-8
-1
0
1
8
27
27
xi4
16
1
0
1
16
81
115
∆yi
∆2 yi
0.52
0.31
0.06
-0.07
-0.32
-------------------
-0.21
-0.25
-0.13
-0.25
----------------------------
⎧6a0 + 3a1 + 19a2 = 3.5
⎪
2
⎨3a0 + 19a1 + 27 a2 = 3.48 ⇒ Y ( x) = −0.102142857 x + 0.201x + 0.806285714
⎪19a + 27 a + 115a = 9
1
2
⎩ 0
51
5 - Ajuste de Funções
5.3.3-Ajuste a polinômio de Grau p
O ajuste a um polinômio de grau p. Y = a0 + a1x + a2 x 2 + ... + a p x p , p < n , exige
resolver o sistema linear:
⎡( n + 1 )∑ x ∑ x 2 ∑ x 3 ...∑ x p ⎤
i
i
i
i
⎥
⎢
+
p
2
3
4
⎢∑ x ∑ x ∑ x ∑ x ...∑ x 1 ⎥
i
i
i
i
⎥
⎢ i
⎢..................................................⎥
⎥
⎢
p
p +1
2p
⎥⎦
⎢⎣∑ xi ∑ xi .................∑ xi
⎡a0 ⎤
⎡∑ yi ⎤
⎢ ⎥
⎥
⎢
⎢a1 ⎥
⎢∑ xi yi ⎥
⎢..... ⎥ = ⎢
⎥
........ ⎥
⎢ ⎥
⎢
⎢..... ⎥
⎢
p ⎥
⎢a ⎥
⎢⎣∑ xi yi ⎥⎦
p
⎣ ⎦
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,y
v]);
> 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=150..280),plot(gll(x), x=-.5..8.5,
y=-150..280,thickness=2)});
52
E.B.Hauser – Cálculo Numérico
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 Y = ab x a função de ajuste.
Linearizamos Y, aplicando log ou ln (escolher a base adequada)
ln
Y = ln
x ln b
{
{a + {
y
a0
xa1
n
n
⎧
⎪( n + 1 ) ln a + ln b
xi =
ln yi
⎪
⎪
i =0
i =0
Determinamos a e b resolvendo: ⎨
n
n
n
⎪
2
⎪ln a
xi + ln b
xi ln yi
xi =
⎪
i =0
i =0
⎩ i =0
⎧⎪a = ln a → a = e a0
0
Então ⎨
⎪⎩a1 = ln b → b = e a1
∑ ∑
∑
∑ ∑
Exemplo 4: Ajustar os dados abaixo a uma função exponencial do tipo Y = ab x .
i
xi
yi
ln yi
0
1
2
3
4
5
6
7
8
0
0.5
1
1.5
2
2.5
3
3.5
4
18
3
4
6
9
12
17
24
33
48
-------------
1.09861
1.38629
1.79176
2.19722
2.48491
2.83321
3.17805
3.49651
3.87120
22.3378
∑
xi2
0
0
0.693145 0.25
1.79176 1
3.29583 2.25
4.96982 4
7.08303 6.25
9.53415 9
12.2378 12.25
15.4848 16
55.0903 51
xi ln yi
Os gráficos a seguir ilustram o efeito da linearização dos dados.
53
Yi
2.98422
4.2228031
5.9754291
8.4555298
11.964948
16.930930
23.958013
33.901647
47.972329
-------------
(Yi − yi ) 2
0.0002490
0.0496412
0.0006025
0.2964477
0.0012286
0.0047706
0.0017628
0.8129689
0.0007656
1.1684403
5 - Ajuste de Funções
⎧9a0 + 18a1 = 22.3378
⎨
⎩18a0 + 51a1 = 55.0903
⇒
a0 = 0.69432 → b = e a1 = 2.002347015
⇒
a
0
a1 = 1.093337778 → a = e = 2.984218125
Y = 2.98422(2.00235) x
-
5.3.4.2 -AJUSTE POR FUNÇÃO POTÊNCIA
Seja Y = axb a função de ajuste.
Linearizando a função Y(x) , temos : ln y = ln a + b ln x.
Resolve-se o sistema de equações lineares e encontra-se a e b:
n
n
⎧
⎪
ln y = y1
( n + 1 ) ln a + b
ln xi =
ln yi
⎪
ln a = a0
⎪
i =0
i =0
⇒
⎨
n
n
n
b = a1
⎪
2
⎪ln a
ln xi + b
(ln xi ) =
ln xi ln yi
ln x = x1
⎪
i =0
i =0
⎩ i =0
∑
∑
∑
∑
então
∑
ln a = a0 ⇒ a = e a0 e b = a1 .
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.
V (m / s ) 100 120 150 180
D(v)=0.813947103e14*1/(x^5.680591334)
D( seg.) 79 28 7.9 2.8
-
54
E.B.Hauser – Cálculo Numérico
Efeito da linearização dos dadoss
5.3.4.3 -AJUSTE POR FUNÇÃO HIPERBÓLICA : Y =
1
. Linearizando, temos :
a0 + a1 x
n
n
⎧
⎪( n + 1 )a0 + a1 ∑ xi = ∑ 1 / yi
1
⎪
i =0
i =0
= a0 + a1 x ⇒ ⎨
n
n
n
Y
⎪∑
xi a0 + a1 ∑ xi2 = ∑ xi / yi
⎪⎩i =0
i =0
i =0
5.3.4.4 -OUTROS TIPOS DE FUNÇÕES DE AJUSTE
x
1) Y =
. Linearizando, temos:
a0 + a1 x
n
n
⎧
⎪( n + 1 )a0 + a1 ∑ xi = ∑ xi / yi
x
⎪
i =0
i =0
= a0 + a1 x ⇒ ⎨
n
n
n
Y
⎪∑
xi a0 + a1 ∑ xi2 = ∑ xi2 / yi
⎪⎩i =0
i =0
i =0
2) Y =
1
a0 + a1 x + a2 x 2
. Linearizando, temos :
⎧( n + 1 )a0 + a1 ∑ xi + a2 ∑ xi2 = ∑ 1 / yi
⎪
1
⎪
= a0 + a1 x + a2 x 2 ⇒ ⎨a0 ∑ xi + a1 ∑ xi2 + a2 ∑ xi3 = ∑ xi / yi
Y
⎪
2
3
4
2
⎪⎩a0 ∑ xi + a1 ∑ xi + a2 ∑ xi = ∑ xi / yi
3)
Y = ae bx + cx
2
2
ln
Y = ln
{
{a + bx
1
4+
2cx
4
3
y
a0
a1 x + a 2 x 2
⇒
⎧( n + 1 ) ln a + b ∑ xi + c ∑ xi2 = ∑ ln yi
⎪
⎪
2
3
⎨ln a ∑ xi + b∑ xi + c ∑ xi = ∑ xi ln yi
⎪
2
3
4
2
⎪⎩ln a ∑ xi + b∑ xi + c ∑ xi = ∑ xi ln yi
⇒ ln a = a0 → ae a0 .
55
⇒
5 - Ajuste de Funções
Exercícios:
1) Considerando:
i
0
1
34
x
y
1
2
45
2
3
63
3
4
88
4
5
6
5
6
7
120 159 205
a) Mostrar que o ajuste por uma parábola é adequado.
b) Ajustar os dados a uma parábola.
Resposta: a) ∆2 y i = 7, ∀i
b) p( x) = 3,51194762 x 2 + 0,4047619048x + 30,14285714
2) Segundo o critério dos Mínimos Quadrados, qual das duas funções Y1 = 2 ,332 + 0 ,558 x
e Y2 = 2 ,037 e0 ,235 x melhor ajusta os dados da tabela?
x
y
-2,3 1
1,2 2,5
3,1
4,2
∑ ( Y2i − yi )2 < ∑ ( Y1i − yi )2
∑ ( Y1i − yi )2 = 0,194 e ∑ ( Y2i − yi )2 = 0,0123
Resposta: Y2, pois
3) Um filme vem sendo exibido numa determinada
sala de cinema por cinco semanas e a frequência
semanal, (aproximada à centena mais próxima) está
dada na tabela abaixo. Utilizar ajuste linear para
determinar a frequência esperada na sexta semana.
6000
5000
4000
y = -360x + 5280
R20,9818 =
3000
semana
2000
1000
1
2
3
4
5
frequência 5000 4500 4100 3900 3500
0
0
1
2
3
4
5
6
Resposta: Y(x) = 5280 -360x e y(6) = 3120
56
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
horas é apresentado pela tabela. Ajuste os dados tabulados a uma curva exponencial da
forma y =abx e avaliar y para x=7.
x
y
0
32
1
47
2
65
3
92
4
5
6
132 190 275
Resposta: y = 32,1483 × 1,42694 x e y (7) = 387, 256
5) Utilizando o critério dos Mínimos Quadrados, ajustar a uma reta os dados tabulados:
xi
yi
3
2
5
3
6
4
8
6
9
5
11
8
Resposta: 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
96.92
T( oC )
690
97.32
101
100,5
100
99,5
99
98,5
98
97,5
97
96,5
660
700
97.71
710
98.11
720
98.49
730
98.88
740
99.26
780
100.73
y = 15,476x 0,2813
R21 =
680
700
720
740
760
780
800
7)
A tabela abaixo fornece uma relação entre a resistência à tração do aço em função da
t (oC )
250 330 412 485 617
temperatura.
Ajustar os dados a um polinômio de
σ (kg / cm2 ) 5720 5260 4450 2780 1500
grau 4.
7000 y = 2E-06x4 0,0033x - 3 1,8931x + 2 466,48x + 47846 2
R 1=
6000
5000
4000
3000
2000
1000
0
0
57
200
400
600
800
5 - Ajuste de Funções
8) Os dados abaixo referem-se a variação do coeficiente de atrito ( µ ) ,entre a roda e o
trilho seco, com a velocidade.
10
20
30
40
60
70
V (km / h) 0
µ
0.450 0.313 0.250 0.215 0.192 0.164 0.154
Ajustar os dados a um polinômio de grau 5.
0,500
0,450
0,400
0,350
0,300
0,250
0,200
0,150
5
4
3
2
0,100 y = -9E-10x 2E-07x + 2E-05x - 0,0007x + 0,0196x + 0,45
0,050
2
R 1=
0,000
0
20
40
60
9) Os dados abaixo relacionam a viscosidade
η em função da temperatura t.
t oC 7.5 10.9 14 15 16 18
21
η 1.409 1.276 1.175 1.148 1.121 1.069 0.990
Realizar o ajuste sugerido pelo gráfico ao lado.
Resp: η( t ) =-0.030772617131t+1.619873714
58
80
E.B.Hauser – Cálculo Numérico
10)
Verificar qual das funções
Y1 = 0.015619 − 0.0001523 x
ou
Y2 = 0.017054(0.98028) x
melhor se ajusta à tabela dada .
20
40
60
80
100
xi
yi 0.0131 0.00654 0.00459 0.0034 0.0026
Respostas:
Aplicação1: Y(d)= 1/(.1329810498e-2*d^2-.2080049421e-1*d+.6161059331)
Aplicação2: Y(p) = 5.779411/(2.043906+p)
59
6 . Integração Numérica
b
Objetivo : Calcular a ∫ f ( x)dx , onde a função integrando f (x) ou é conhecida por sua expressão
a
analítica ou por uma tabela de valores (xi , f ( xi )) , i = 0, 1,..., n.
Aplicação: Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada, reproduzindo
os dados:
x( hora )
9 10 11 12 13 14 15 16 17
y( temperatura ) 75.3 77.0 83.1 84.8 86.5 86.4 81.1 78.6 75.1
Encontrar a temperatura média da água entre 9h da manhã e 5h da tarde e estimar o erro cometido
nesse cálculo.
OBS : fm =
f é
1 b
∫ f ( x )dx
b−aa
é o valor médio de f ( x ) em [ a ,b ], se
f ( a ). f ( b ) > 0 e
contínua em [ a ,b ]
6.1 -Fórmulas de Newton-Cotes
x −x
Consideremos (xi , f ( xi )) , i = 0, 1,..., n., xi +1 − xi = h , h = n 0 , yi = f ( xi ) .
n
A integral da função f(x) no intervalo [x0 ; xn ] é dada por :
xn
xn
x0
x0
∫ f ( x )dx ≅
=
∫ Pn ( x )dx
=
( x − x0 )...( x − xn −1 ) n ⎞
( x − x0 )( x − x1 ) 2
( x − x0 )
⎛
∆y0 +
∆ y0 + ... +
∆ y0 ⎟⎟ dx
∫ ⎜⎜ y0 +
2
h
n! h n
2! h
⎠
x0 ⎝
xn
Se R =
x=
x − x0
, então
x=
h
x0 → R = 0
x0 + Rh → dx = hdr
Assim:
60
e
x = xn → R = n .
E.B.Hauser – Cálculo Numérico
xn
∫x0
n
n
f ( x )dx ≅ h ∫ P( R )dR = h ∫ ( y0 + R∆y0 +
0
0
R( R − 1 ) 2
R( R − 1 )...( R − n + 1 ) n
∆ y0 + ... +
∆ y0 )dR
2!
n!
⎤
⎡
n
n
∆2 y0 n
∆n y0 n
R
(
R
1
)
dR
...
R
(
R
1
)...(
R
n
1
)
dR
= h ⎢ y0 ∫ dR + ∆y0 ∫ RdR +
−
+
+
−
−
+
⎥
0
2! ∫0
n! ∫0
⎥⎦
⎢⎣ 0
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 :
x1
1
1
1
f ( x )dx ≅ h P( R )dR = h ⎡⎢ y0 dR + ∆y0 RdR ⎤⎥
x0
0
0
⎦
⎣ 0
∫
∫
∫
∫
[
]
= h y0 [ R ]01 + ∆y0 [ R 2 / 2 ]01 =
∆y ⎤ h
⎡
= h ⎢ y0 + 0 ⎥ = [ y0 + y1 ], ou para o intervalo [xi ; xi + 1 ] :
2 ⎦ 2
⎣
xi + 1
∫xi
f ( x )dx ≅ h ∫
i +1
i
P( R )dR =
h
[ yi + yi + 1 ]
2
Generalizamos para n subintervalos:
xn
∫ f ( x )dx ≅
xo
[
h
yo + 2( y1 + y 2 + L + y n −1 ) + y n
2
]
Erro de Truncamento(para n subintervalos):
ET ≤
h2
( x n − x0 )
max
f ' ' ( x)
12
x∈[ x o , x n ]
ou
ET ≤
( x n − x0 )
max ∆2 yi
12
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
1 − x2
e
dx com um
0
∫
erro de truncamento menor que 10 −4 .
ET ≤
h2
h2
( 2 ) < 10 − 4 ⇒ h < 0.0245
( 1 − 0 )máx f ' ' ( x ) =
12
12
x ∈ [0 ;1]
n = ( xn − x0 ) / h , h = ( xn − x0 ) / n , ( xn − x0 ) / n < 0.0245 , n > 40.8 → n = 41 , h =
61
1
= 0.02439
41
6 -Integração Numérica
6.3 – Fórmula de Simpson
Fazendo n = 2 na fórmula de Newton-Cotes, temos,
⎡
⎤ h
2
2
x2
∆2 y0 2
≅
+
+
−
∆
R
(
R
1
)
dR
f
(
x
)
dx
h
y
dR
y
RdR
⎢
⎥ = [ y0 + 4 y1 + y 2 ]
0
0
∫x0
∫0
∫0
2! ∫0
⎣⎢
⎦⎥ 3
Generalizamos para n subintervalos, n par:
xn
∫ f ( x )dx ≅
xo
[
h
yo + 4( y1 + y3 + y5 + L + y n −1 ) + 2( y 2 + y4 + y6 + L + y n − 2 ) + y n
3
]
ERRO DE TRUNCAMENTO PARA A FÓRMULA DE SIMPSON
( x − x0 )
h4
( xn − x0 ) max
f ' '' ' ( x ) ou E S ≤ n
max ∆4 yi
ES ≤
180
180
x∈[ xo , x n ]
Exercícios
1.Aplicação:.Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada,
reproduzindo os dados:
x( hora )
9 10 11 12 13 14 15 16 17
y( temperatura ) 75.3 77.0 83.1 84.8 86.5 86.4 81.1 78.6 75.1
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
subintervalos.
62
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.
π /2
∫ sen( 2 cos x ) sen
2
xdx
0
4. Calcular por Trapézios:
e− x
∫ e − x dx
0
2
a)
1
2
b) ∫ e x tgxdx
com
h = 0.5
com
h=0.1
0
5. Calcular por Simpson:
1.2
a)
∫
0
dx
ex + x + 1
2
2
b) ∫ e − x dx
com
h =0 .3
com
h = 0.25
0
π/2
c)
∫x
2
cos xdx
com
n=0.4
0
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
63
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
parciais 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 estudo da evolução da temperatura e
reações químicas em cilindros e esferas. Ao lado a
representação gráfica de J 0 ( t )
Tarefa: Considerar a representação integral da Função de Bessel de primeiro tipo
J m( t ) =
1
π
π
∫ cos( mx − tsenx )dx ,
m = 0 ,1,2 ,3 ,...
0
Estimar J 0 ( 3 ) com cinco subintervalos e o erro cometido neste cálculo.
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. J 0 ( 3 ) = -0,260052
64
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
⎧⎪ y' = f ( x , y )
(PVI) ⎨
,
⎪⎩ y( xo ) = y o
e o problema de valor de contorno de segunda ordem, linear:
⎧⎪ y' ' + P( x ) y' +Q( x ) y = f ( x )
(PVC) ⎨
⎪⎩ y( xo ) = yo , y( xn ) = y n .
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
xo , x1 , x2 , K , onde xi = xo + ih , 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 y ( xi ) ≅ y i , xi = xo + ih , i= 0,1,...,n e h = ( x n − xo ) / n .
Método de Euler(PVI):
yi + 1 = yi + hf ( xi , yi ) , i= 0,1,...,n-1.
Método de Runge-Kutta de 2 a ordem(PVI) :
h
yi + 1 = yi + (k1 + k 2 ) , k1 = f ( xi , yi ) e k 2 = f ( xi + h , yi + hk1 ) , i= 0,1,...,n-1.
2
Diferenças Finitas(p/diferenças centrais)(PVC): Para i= 1,...,n-1 , é gerado um sistema de n-1
equações lineares:
⎛ h
⎞
⎛ h
⎞
2
2
⎜ 1 + P( xi ) ⎟ yi + 1 + − 2 + h Q( xi ) yi + ⎜ 1 − P( xi ) ⎟ yi −1 = h f ( xi )
⎝ 2
⎠
⎝ 2
⎠
(
)
Sistemas de Equações Diferenciais de primeira ordem com condições iniciais
⎧ y' = f ( x , y , z )
⎪⎪
⎨ z' = g ( x , y , z )
⎪
⎪⎩ y( xo ) = yo , z( xo ) = z o
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:
⎧⎪ y' ' = f ( x , y , z )
Para ⎨
, fazemos a mudança de varável y' = u , com y( xo ) = yo .
⎪⎩ y( xo ) = yo , y' ( xo ) = zo
Então devemos resolver o sistema de duas equações diferenciais ordinárias, acopladas
⎧ y' = u
⎪⎪
pelas condições iniciais: ⎨u' = f ( x , y , z )
.
⎪
⎪⎩ y( xo ) = yo , u( xo ) = z o
65
E.B.Hauser – Cálculo Numérico
Derivada
Diferença Finita
h, k=tamanho do passo na direção x e na direção y (ou t)
y
y´( xi )
−y
i +1
2h
yi + 2 − 2 y i + 1 + 2 yi − 1 − yi − 2
y' ' ' ( xi )
2h 3
i −1 centrada ,
centrada
y i + 2 − 4 y i + 1 + 6 y i −4 y i − 1 + y i − 2
centrada
h4
y IV ( xi )
(
y
yi + 1 − 2 y i + y i − 1
centrada
h2
y´´( xi )
∂u
xi , t j
∂t
− yi
i +1
avançada ,
h
)
(
∂ 2u
xi ,t j
∂ x2
(
u ,
−u ,
i j +1
i j −1
centrada
2k
u
)
∂ 2u
xi , y j
∂ y2
)
i +1 , j
− 2u
i,j
h2
+u
i −1 , j
centrada
− 2u , + u ,
u ,
i j +1
i j
i j −1
centrada
k2
Exemplos:
66
yi − y
i − 1 atrasada
h
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Método de Euler
Problema: y' = y – x; y(0) = 2
yn
xn
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
h = 0.1
h = 0.05
h = 0.01
h = 0.005
2.0000
2.2000
2.4100
2.6310
2.8641
3.1105
3.3716
3.6487
3.9436
4.2579
4.5937
2.0000
2.2025
2.4155
2.6401
2.8775
3.1289
3.3959
3.6799
3.9829
4.3066
4.6533
2.0000
2.2046
2.4202
2.6478
2.8889
3.1446
3.4167
3.7068
4.0167
4.3486
4.7048
2.0000
2.2049
2.4208
2.6489
2.8903
3.1467
3.4194
3.7102
4.0211
4.3541
4.7115
xn
h = 0.1
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
2.0000000
2.2050000
2.4210250
2.6492326
2.8909021
3.1474468
3.4204287
3.7115737
4.0227889
4.3561818
4.7140809
Método de Runge-Kutta de segunda ordem
Problema: y' = y – x; y(0) = 2
yn
h = 0.05
h = 0.01
2.0000000
2.2051266
2.4213047
2.6496963
2.8915852
3.1483904
3.4216801
3.7131870
4.0248265
4.3587148
4.7171911
2.0000000
2.2051691
2.4213987
2.6498521
2.8918148
3.1487076
3.4221007
3.7137294
4.0255115
4.3595665
4.7182369
67
Solução exata
Y(x) = ex + x + 1
2.0000
2.2052
2.4214
2.6499
2.8918
3.1487
3.4221
3.7138
4.0255
4.3596
4.7183
Solução exata
Y(x) = ex + x + 1
2.0000000
2.2051709
2.4214028
2.6498588
2.8918247
3.1487213
3.4221188
3.7137527
4.0255409
4.3596031
4.7182818
E.B.Hauser – Cálculo Numérico
Método de Euler
Problema: y' = y; y(0) = 1
yn
xn
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
xn
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
h = 0.1
1.0000
1.1000
1.2100
1.3310
1.4641
1.6105
1.7716
1.9487
2.1436
2.3579
2.5937
h = 0.05
1.0000
1.1025
1.2155
1.3401
1.4775
1.6289
1.7959
1.9799
2.1829
2.4066
2.6533
h = 0.01
1.0000
1.1046
1.2202
1.3478
1.4889
1.6446
1.8167
2.0068
2.2167
2.4486
2.7048
h = 0.005
1.0000
1.1049
1.2208
1.3489
1.4903
1.6467
1.8194
2.0102
2.2211
2.4541
2.7115
Método de Runge-Kutta de segunda ordem
Problema: y' = y; y(0) = 1
yn
h = 0.1
h = 0.05
h = 0.01
1.0000000
1.1050000
1.2210250
1.3492326
1.4909021
1.6474468
1.8204287
2.0115737
2.2227889
2.4561818
2.7140809
1.0000000
1.1051266
1.2213047
1.3496963
1.4915852
1.6483904
1.8216801
2.0131873
2.2248265
2.4587148
2.7171911
1.0000000
1.1051691
1.2213987
1.3498521
1.4918148
1.6487076
1.8221007
2.0137294
2.2255115
2.4595665
2.7182369
68
Solução exata
Y(x) = ex
1.0000
1.1052
1.2214
1.3499
1.4918
1.6487
1.8221
2.0138
2.2255
2.4596
2.7183
Solução exata
Y(x) = ex
1.0000000
1.1051709
1.2214028
1.3498588
1.4918247
1.6487213
1.8221188
2.0137527
2.2255409
2.4596031
2.7182818
7 -Resolução Numérica de Equações Diferenciais Ordinárias
Método de Euler
Problema: y' = 4x3; y(0) = 1
xn
0.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Solução exata
Y(x) = x4
yn
h = 0.1
0.0000
0.0000
0.0004
0.0036
0.0144
0.0400
0.0900
0.1764
0.3136
0.5184
0.8100
h = 0.05
0.0000
0.0000
0.0009
0.0056
0.0196
0.0506
0.1089
0.2070
0.3600
0.5852
0.9025
h = 0.01
0.0000
0.0001
0.0014
0.0076
0.0243
0.0600
0.1253
0.2333
0.3994
0.6416
0.9801
0.0000
0.0001
0.0016
0.0081
0.0256
0.0625
0.1296
0.2401
0.4096
0.6561
1.0000
Método de Runge-Kutta de segunda ordem
Problema: y' = 4x3; y(0) = 0
xn
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
yn
h = 0.1
0.0000000
0.0002000
0.0020000
0.0090000
0.0272000
0.0650000
0.1332000
0.2450000
0.4160000
0.6642000
1.0100000
Solução exata
Y(x) = x4
h = 0.05
0.0000000
0.0001250
0.0017000
0.0083250
0.0260000
0.0631250
0.1305000
0.2413250
0.4112000
0.6581250
1.0025000
h = 0.01
0.0000000
0.0001010
0.0016040
0.0081090
0.0256160
0.0625250
0.1296360
0.2401490
0.4096640
0.6561810
1.0001000
69
0.0000000
0.0001000
0.0016000
0.0081000
0.0256000
0.0625000
0.1296000
0.2401000
0.4096000
0.6561000
1.0000000
E.B.Hauser – Cálculo Numérico
Exercícios:
1. Utilizando o método de Euler, determinar y(X i ), se :
⎧⎪ y` = − xy − 1
⎨
, h=0,25
,
X i =1,75
a) ⎪⎩ y( 1 ) = 2
⎧⎪ y` = x + y , y > 0
b) ⎨
⎪⎩ y( 1 ) = 0
⎧⎪ y` = x 2 + 2 y
⎨
c) ⎪ y( 1 ) = 0 ,2
⎩
⎧⎪ y` = 100 y − 101e x − 100
⎨
d) ⎪⎩ y( 0 ) = 2
, h=0,1
,
X i =1,4
, h=0,1
,
X i =1,4
, h=0,1
,
⎧⎪ xy`` − y` −8 x 3 y 3 = 0
, h=0,1
e) ⎨
⎪⎩ y( 1 ) = 0.5 , y`( 1 ) = −0 ,5
⎧⎪ y`` −( 1 − y 2 ) y` + y = 0
,
f) ⎨
⎪⎩ y( 0 ) = 0 ,5 , y`( 0 ) = 0
h=0,1
X i =0,3
,
X i =1,3
,
X i =0,3
2.Utilizando o Método de Heun (Runge-Kutta de 2ª ordem), determinar y( X i ), se:
⎧ y´ = − xy
⎧⎪ y´ = 3 x 2 + y
⎨
b) ⎩ y( 0 ) = 1 , h=0,1, X i =0,2
a) ⎨
, h=0,2, X i =0,4
⎪⎩ y( 0 ) = 1
⎧⎪ y´ = 4 x 3
c) ⎨
, h=0,1 , X i =0,3
⎪⎩ y( 0 ) = 0
3.Utilizando o Método das diferenças finitas, e o valor indicado de n, resolver o PVC.
⎧ y' ' +9 y = 0
a) ⎨
,n = 4
⎩ y( 0 ) = 4 , y( 2 ) = 1
⎧⎪ x 2 y' ' +3 xy' +3 y = 0
c) ⎨
⎪⎩ y( 1 ) = 5 , y( 2 ) = 0
,n = 8
⎧ y' ' +2 y' + y = 5 x
b) ⎨
⎩ y( 0 ) = 4 , y( 1 ) = 0
,n = 5
⎧ y' ' +( 1 − x ) y' + xy = x
d) ⎨
⎩ y( 0 ) = 0 , y( 1 ) = 2
,n = 10
70
7 -Resolução Numérica de Equações Diferenciais Ordinárias
4. PVI -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. PVC - Considerar o problema de deflexão de uma viga de seção transversal retangular sujeita a
uma carga uniforme, tendo seus extremos apoiados de modo a não sofrer deflexão alguma. O
problema de valor de contorno que rege essa situação física é
d 2w
q
S
=
x( x − 1 ), 0 < x < L ,
w+
2 EI
dx 2 EI
Como não ocorre deflexão nas extremidades da viga, as condições de contorno são w(0) =0,
w(L) =0. Considerando:
• Comprimento L=120 pol;
• Intensidade de carga uniforme q=100 lb/pé;
• Módulo de elasticidade E=3.0 x 107 lb/pol2;
• Esforço nas extremidades S=1000 lb;
4
• Momento central de Inércia I=625 pol ;
aproximar a deflexão w(x) da viga a cada 20pol, utilizando diferenças finitas.
Respostas :
1.a) y(1,75)=1,2018
b) y(1,4)=0,4558
c) y(1,4)=0,7778
d) y(0,3)=-25,2206
e) y(1,3)= 0, 3647
f) y(0,3)= 0,3991
2.a) y(0,4)=1,1189
b) y(0,2)=0,9801
c) y(0,3)=0,009
y1 = 0 ,2660
y3 = 6 ,3226
y 2 = 0 ,5097
y 2 = 2 ,9640
y3 = 0 ,7357
y3 = 2 ,2064
y4 = 0 ,9471
c) y4 = 1,5826
y5 = 1,0681
d) y5 = 1,1465
y6 = 1,3353
y6 = 0 ,6430
y7 = 1,5149
y7 = 0 ,2913
y8 = 1,6855
y1 = −0 ,2259
y1 = −5 ,6774
3.a) y 2 = −2 ,5807
y1 = 3,8842
b)
y 2 = −0 ,3356
y3 = −0 ,3308
y4 = −0 ,2167
y9 = 1,8474
4) y(0.5)=0.875
71
8 – Resolução Numérica de Equações Diferenciais Parciais
8.1 -Introdução
Equação diferencial parcial (EDP) é a uma equação que envolve duas ou mais
variáveis independentes ( x , y , z ,t ,K ) e derivadas parciais de uma função incógnita(variável
dependente que queremos determinar) u ≡ u( x, y, z,t ,K) .
Um corpo é isotrópico se a condutividade térmica em cada um de seus pontos é
independente da direção do fluxo de calor através do ponto.
Em um corpo isotrópico, a temperatura, u ≡ u(x, y, z, t), é obtida resolvendo-se a
equação diferencial parcial (EDP)
∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞
∂u
⎜ k ⎟ + ⎜⎜ k ⎟⎟ + ⎜ k ⎟ = cp
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠
∂t
onde k, c e p são funções de (x,y,z), e representam respectivamente, a condutividade térmica, o
calor específico e a densidade do corpo no ponto (x,y,z).
Quando k, c e p são constantes, essa equação é denominada equação simples
tridimensional do calor, e é expressa como
∂ 2u
∂x
2
+
∂ 2u
∂y
+
2
∂ 2u
∂z
2
=
cp ∂ u
.
k ∂t
Se o domínio do problema é relativamente simples, a solução dessa equação é obtida
utilizando a série de Fourier. Na maioria das situações onde k, c e p não são constantes ou
quando o domínio é irregular, a solução da equação diferencial parcial deve ser obtida por meio
de métodos de aproximação.
Para introduzir métodos numéricos de resolução de EDP, utilizaremos as equações de
Poisson, do Calor e da Onda, as quais representam protótipos das EDP´s elípticas, parabólicas e
hiperbólicas. Será adotado um procedimento geral, seguindo os passos:
1) Construir uma malha a partir do domínio do problema;
2) Para os pontos interiores da malha, escolher a discretização das derivadas parciais;
3) Construir o sistema de equações lineares usando a discretização dos pontos interiores,
f xi , y j e as condicções de contorno.
(
)
4) Resolver o sistema de equações lineares(escolher o método masi eficiente), cuja
solução forne as aproximações da solução nos pontos interiores da malha.
8.1.1-Equação Do Potencial ou de Poisson(EDP Elíptica)
Consideremos a equação de Poisson:
∂ 2u
∂x 2
( x, y ) +
∂ 2u
∂y 2
72
( x , y ) = f ( x , y ).
E.B.Hauser – Cálculo Numérico
Nessa equação supomos que a função f descreve os dados do problema em uma região
plana R com fronteira S. Equações desse tipo aparecem durante o estudo de diversos problemas
físicos dependentes do tempo; por exemplo, a distribuição de calor para um estado estável em
uma região plana, a energia potencial de um ponto em um plano sobre o qual atuam forças
gravitacionais e os problemas bidimensionais do estado de equilíbrio que incluem fluidos não
comprimíveis.
Para se obter uma solução única para equação de Poisson é necessário impor outras
restrições. Por exemplo, o estudo da distribuição de calor no estado de equilíbrio em uma região
plana requer que f ( x , y ) ≡ 0 que é a equação de Laplace
∂ 2u
∂x 2
( x, y ) +
∂ 2u
∂y 2
( x, y ) = 0,
Se a temperatura na região é determinada por sua distribuição no limite da região, as
restrições são denominadas Condições de limite de Dirichlet, dadas por u ( x, y ) = g ( x, y ),
para todo(x,y) em S, a fronteira da região R ( ver figura 1).
Figura 1
8.1.2- Equação de Calor ou da Difusão (EDP Parabólica)
A equação do calor ou de difusão (que é uma equação diferencial parcial
parabólica)
∂u
∂ 2u
( x ,t ) +
( x ,t ) = 0 ,
∂t
∂x 2
modela matematicamente o problema físico referente ao fluxo de calor ao longo de uma barra de
comprimento l (figura 2), a qual tem uma temperatura uniforme dentro de cada elemento
transversal. Essa condição requer que a superfície lateral da barra esteja perfeitamente isolada. A
constante α é determinada pelas propriedades de condução de calor do material de que a barra é
feita e é independente da posição da barra.
Figura2
73
8 - Resolução Numérica de Equações Diferenciais Parciais
Um dos conjuntos típicos de restrições para um problema de fluxo de calor desse tipo
consiste em especificar a distribuição inicial de calor na barra: u(x,0)=f(x) e em descrever o
comportamento nas extremidades da barra. Por exemplo, se as extremidades são mantidas em
temperaturas constantes U i e U 2 , as condições de contorno têm a forma:
u( 0 ,t ) = U 1 e u( l ,t ) = U 2 ,
e a distribuição de calor se aproxima da distribuição limite de temperatura
lim u ( x , t ) = U 1 +
t→ ∞
U 2 − U1
x.
l
Se, a barra estiver isolada de modo que não flua calor por suas extremidades, as condições
de contorno serão:
∂u
∂u
(0, t ) = 0 e
(l , t ) = 0,
∂x
∂x
o que resulta em uma temperatura constante na barra como caso limite.
A equação diferencial parcial parabólica também é importante para o estudo da
difusão dos gases.
8.1.3- Equação da Onda (EDP Hiperbólica)
Consideremos a equação da Onda unidimensional , um exemplo de uma equação
diferencial parcial hiperbólica. Supomos que uma corda elástica ,de comprimento l , seja
esticada entre dois suportes no mesmo nível horizontal(figura 3)
Figura 3
Se pusermos a corda em movimento de modo que ela vibre em um plano vertical, o
deslocamento vertical u ( x , t ) de um ponto x no tempo t satisfará a equação diferencial parcial
α
2
∂ 2u
∂x 2
( x ,t ) =
∂ 2u
∂t 2
( x ,t ), para 0 < x < l , 0 < t ,
se os efeitos de amortização forem desconsiderados e a amplitude não for muito grande.
Para impor restrições a esse problema, vamos supor que a posição e a velocidade
iniciais da corda sejam dadas por
74
E.B.Hauser – Cálculo Numérico
∂u
( x ,0 ) = g( x ),
∂t
Se os pontos extremos forem fixos, teremos:
u( 0 ,t ) = 0 e u( l ,t ) = 0 .
Os outros problemas físicos envolvendo a equação diferencial parcial hiperbólica
ocorrem no estudo de vigas vibrantes com uma ou ambas as extremidades clamped e na
transmissão de eletricidade em uma linha de transmissão longa onde exista alguma perda de
corrente para o solo.
u ( x,0) = f ( x) e
8.2. -Método das Diferenças Finitas para Equação Diferencial Parcial Elíptica
Consideremos o problema de valor de contorno envolvendo a equação de Poisson (que é
uma equação diferencial parcial elíptica)
⎧ 2
∂ 2u
∂ 2u
⎪∇ u( x , y ) ≡ 2 ( x , y ) + 2 ( x , y ) = f ( x , y )
⎨
∂x
∂y
⎪
⎩ u( x , y ) = g( x , y ) para ( x , y ) ∈ S ,
onde S é a fronteira(contorno) do retângulo R = {( x , y ) / a < x < b , c < y < d } .
Se f e g são contínuas em seus domínios, então existe uma única solução para esse problema
de valor de contorno..
Utilizaremos uma adaptação do método de Diferenças Finitas .
Em R contruímos uma grade(figura4) , traçando linhas verticais e horizontais ( x = xi e ,
y = y j grid lines) pelos pontos (x i , y j ), e suas intersecções são os pontos de rede (mesh points),
onde
x i = a + ih, h=(b – a)/n , para todo i = 0,1,...,n
e
y j = c + jk , k=(d – c)/m , para todo j = 0,1,...,m
Figura 4
75
8 - Resolução Numérica de Equações Diferenciais Parciais
Para cada ponto de rede no interior da quadricula (x i , y j ), com i = 1,2,...,n-1 e com j =
1,2,...,m-1, utilizamos a série Taylor na variável x ao redor de xi para gerar a fórmula das
diferenças centrais
(
(
xi , y j ) =
2
)
(
) (
∂2u
u xi +1 , y j − 2u xi , y j + xi −1 , y j
∂x
h2
) − h2 ∂4u (ξ , y ),
12 ∂x4
i
j
(2.1)
onde ξ ∈ ( xi −1 , xi +1 ). Também utilizamos a série de Taylor na variável y ao redor de y j para gerar a
fórmula das diferenças centrais
u(xi +1 , y j ) − 2u(xi , y j ) + (xi −1 , y j ) k ∂ u
(
−
xi , y j ) =
(xi ,η j ),
2
2
12 4
∂2u
∂y
2 4
(2.3)
∂y
k
onde η ∈ ( y j −1 , y j +1 ).
A utilização dessas formulas na Equação (2.1) nos permite expressar a equação de Poisson
nos pontos (x i , y j ), como
(
)
(
) (
u xi +1 , y j − 2u xi , y j + xi −1 , y j
h2
(
)
= f xi , y j +
) + u(xi+1 , y j ) − 2u(xi , y j ) + (xi−1 , y j )
k2
(
)
(
)
h 2 ∂ 4u
k 2 ∂ 4u
+
ξ
,
y
xi ,η j ,
i j
12 ∂x4
12 ∂y4
para todo i = 1,2,...,n-1 e j = 1,2,...,m-1, e as condições de limite como
u ( x0 , y j ) = g ( x0 , y j ) e u ( x n , y j ) = g ( x n , y j ), para todo j = 0,1,...,m;
u ( xi , y 0 ) = g ( xi , y 0 ) e u ( xi , y m ) = g ( xi , y m ), para todo i = 1,2,...,n-1;
Na forma da equação de diferenças, isso resulta no método das Diferenças Finitas para a
equação de Poisson, com um erro local de truncamnto da ordem de O(h 2 + k 2 ) :
2
⎤
⎡⎛ h ⎞ 2
⎛h⎞
2 ⎢⎜ ⎟ + 1⎥ wij − wi +1, j + wi −1, j − ⎜ ⎟ wi , j + 1 + wi , j −1 = −h 2 f xi , y j ,
⎝k⎠
⎥⎦
⎢⎣⎝ k ⎠
i = 1,2,..., n - 1 e j = 1,2,..., m - 1.
(
)
(
)
(
)
(2.4)
w0 j = g ( x0 , y j ) e wnj = g ( x n , y j ) , para todo j = 0,1,...,m;
wi 0 = g ( xi , y 0 ) e wim = g ( xi , y m ), para todo i = 1,2,...,n-1;
onde wij é uma aproximação de u ( xi , y j ) .
76
(2.5)
E.B.Hauser – Cálculo Numérico
A equação em (2.4) envolve aproximações a u ( x , y ) nos pontos
(x
i −1
, y j ),
(x , y ), (x
i
j
i +1
, y j ),
(x , y ) e (x , y ) .
i
j −1
i
j +1
Reproduzindo a parte da malha na qual esses pontos estão situados (veja figura.5),
observamos que cada equação contém aproximações em uma região em forma de estrela ao redor
de (xi , y j ) .
Figura 5
Se utilizarmos a informação das condições de limite (2.5) sempre que for conveniente no
sistema dado por (12.4), poderemos dizer que, em todos os pontos (xi , y j ) adjacentes ao ponto de
rede do limite, teremos um sistema linear ( n – 1)(m – 1) x ( n – 1)(m – 1), cujas incógnitas são as
aproximações wij a u ( xi , y j ) no interior dos pontos de rede.
O sistema linear que contém essas incógnitas será expresso mais eficientemente em cálculos
matriciais se for introduzida uma remarcação dos pontos interiores da malha. Um sistema de
marcação desses pontos consiste em utilizar
Pl = xi , y j e wl = wij
(
)
onde l = i + (m − 1 − j )(n − 1), para todo i = 1,2,...,n-1 e j = 1,2,...,m-1. Isso marca
consecutivamente os ponto de rede da esquerda para direita e de cima para baixo. Por exemplo, com
n = 4 e m = 5, com a remarcação se obtém uma quadrícula cujos pontos são mostrados na Figura
12.6.Ao marcar os pontos desse modo, se garante que o sistema necessário para determinar wij seja
uma matiz de banda com uma largura de banda de, no máximo, 2n – 1.
Figura 6
77
8 - Resolução Numérica de Equações Diferenciais Parciais
Exemplo1
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
∂ 2u
∂ 2u
(
x
,
y
)
+
( x, y ) = 0,
∂x 2
∂y 2
para (x , y) no conjunto R = {( x , y ) / 0 < x < 0 ,5 , 0 < y < 0 ,5 }, 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
4wi , j − wi+1, j − wi−1, j − wi , j−i − wi , j+i = 0,
para todo i= 1,2,3 e j=1,2,3 .
Para expressar isso em função dos interiores da grade, usamos
Pl = xi , y j e wl = wij
(
e
)
wi = u ( Pi ) , e l = i+(m-1-j)(n-1).
78
E.B.Hauser – Cálculo Numérico
Também, a partir das condições de contorno e usando (2.5)
w 1 ,0 = w 2 ,0 = w 3 ,0 = w 0 ,1 = w 0 ,2 = w 0 ,3 = 0 ,
w 1 ,4 = w 4 ,1 = 25 , w 2 ,4 = w 4 ,2 = 50 , e , w 3 ,4 = w 4 ,3 = 75
Assim, para cada ponto interior da grade, Pi geramos uma equação linear:
P 1
P 2
: 4 w 1 − w 2 − w 4
: 4 w 2 − w 3 − w 1
P 3
P 4
: 4 w 3
: 4 w 4
−
P 5
P 6
: 4 w 5
: 4 w 6
−
P7
P 8
: 4 w 7
: 4 w 8
−
P 9
: 4 w 9
−
−
w 2
w 5
−
w 6
w 5
−
−
−
−
−
w 8
w 9
−
w 8
=
−
=
w 6
w 1
−
=
w 4
w 3
−
w 4
w 7
−
w 6
−
w 0 ,3 + w 1 ,4
w 5 = w 2 ,4
w 4 ,3
w 7 =
w 2 −
w 9 =
+
w 3 ,4
w 0 ,2
w 8 =
0
−
w 0 ,1
w 5 =
w 4 ,2 ,
+ w 1 ,0 ,
w 2 ,0 ,
=
w 3 ,0
+
=
w 4 ,1 ,
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.
79
8 - Resolução Numérica de Equações Diferenciais Parciais
> 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
wi
1
18.75
2
37.50
3
56.25
4
12.5
5
25.00
6
37.50
7
6.25
8
12.50
9
18.75
Assim, para cada i, wi = u ( Pi ) , isto é, wi é uma estimativa da solução em Pi = ( xi , yi ) .
Por exemplo, w6 = u( P6 ) = u( 0.375 ,0.25 ) ≅ 37.5
Exemplo2:
Consideremos o modelo para determinar a temperatura de estado estacionário para uma
placa quadrada com condições de contorno dadas:
⎧ ∂ 2u
∂ 2u
( x, y ) = 0,0 < x < 2,0 < y < 2
⎪ 2 ( x, y ) +
∂y 2
⎪ ∂x
⎪
⎨u ( 0 , y ) = 0 , u ( 2 , y ) = y( 2 − y )
⎪
x,0 < x < 1
⎪ u ( x ,0 ) = 0 , u ( x ,2 ) = ⎧⎨
⎪
⎩2 − x , 1 ≤ x < 2
⎩
80
E.B.Hauser – Cálculo Numérico
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:
ui + 1, j + ui , j = 1 + ui − 1, j + ui , j − 1 − 4ui , j = 0
Para determinar o valor de
u( 2 / 3 ,2 / 3 ) ≅ u11, u( 2 / 3,4 / 3 ) ≅ u12 , u( 4 / 3,2 / 3 ) ≅ u 21 e u( 4 / 3,4 / 3 ) ≅ u 22
Utilizando o sistema Maple para resolver o sistema linear tridiagonal:
⎧ − 4u11 + u 21 + u12 = 0
⎪u11 − 4u 21 + u 22 = −8 / 9
⎪
⎨
⎪u11 − 4u12 + u 22 = −2 / 3
⎪⎩ u 21 + u12 − 4u 22 = −14 / 9
> 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});
13
7
7
5
{ u12 = , u11 = , u22 = , u21 = }
36
36
12
12
> evalf(%);
{ u12 = .3611111111, u11 = .1944444444, u22 = .5833333333, u21 = .4166666667}
Exercícios:
1) Com h=k=10, resolver numericamente equação de Laplace
∂ 2u
∂ 2u
( x ,t ) +
( x ,t ) = 0 , 0 < x < 80 , 0 < y < 60 ,
∂x 2
∂y 2
sujeita às condições
u ( x ,0 ) = u ( x ,60 ) = u ( 80 , y ) = 0 , e
ui , j =
j↓
6
5
4
3
2
1
0
x→
u ( 0 , y ) = 100 .
ui + 1, j + ui −1, j + ui , j −i + ui , j + i
4
0
100
100
100
100
100
100
100
0
1
0
46,993
63,138
67,192
63,138
46,993
0
10
2
0
24,835
38,368
42,490
38,368
24,835
0
20
3
0
13,978
23,010
26,032
23,010
13,978
0
30
4
0
8,068
13,660
15,620
13,660
8,068
0
40
5
0
4,633
7,942
9,128
7,942
4,633
0
50
2) Resolver numericamente equação de calor
81
6
0
2,523
4,348
5,009
4,348
2,523
0
60
7
0
1,110
1,917
2,211
1,917
1,110
0
70
8
0
0
0
0
0
0
0
80
←i
60
50
40
30
20
10
0
y↓
8 - Resolução Numérica de Equações Diferenciais Parciais
∂u
∂ 2u
( x ,t ) −
( x ,t ) = 0 , 0 < x < 1, 0 ≤ t ,
∂t
∂x 2
com as condições de contorno
u ( 0 , t ) = u ( 1, t ) = 0 , 0 < t ,
e as condições iniciais
u( x ,0 ) = sen( πx ), 0 ≤ x ≤ 1.
Ajuda:
Diferenças finitas para o problema com h=0,2, k=0.01
⎞
⎛
+u
ui , j + 1 = 0 ,5 ui , j + 0 ,25⎜ u
⎟
+
−
i
1
,
j
i
1
,
j
⎠
⎝
j↓
←i
0
1
2
3
4
5
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
↑ tempo
0,2
0,4
0,6
0,8
1
Diferenças finitas para o problema com h=0,1, k=0.005
⎞
⎛
ui , j +1 = 0 ,9 ui , j + 0 ,1⎜ u
+u
⎟
i −1, j ⎠
⎝ i +1, j
j↓
6
5
4
3
2
1
0
x→
0
1
2
3
4
5
6
x
0
0
0
0
0
0
0
0
0
1
0,5189
0,4759
0,4365
0,4004
0,3673
0,3369
2
0,9869
0,9053
0,8304
0,7616
0,6986
0,6408
3
1,3584
1,2460
1,1429
1,0483
0,9616
0,8820
4
1,5969
1,4647
1,3435
1,2324
1,1304
1,0369
0,3090
0,5878
0,8090
0,9511 1,0000 0,9511
0,1
0,2
0,3
0,4
8.4 – Bibliografia para EDP`s
82
5
1,6791
1,5401
1,4127
1,2958
1,1886
1,0902
0,5
6
1,5969
1,4647
1,3435
1,2324
1,1304
1,0369
0,6
9 10 ← i
0,5189 0 0,030
0,4759 0 0,025
0,4365 0 0,020
0,4004 0 0,015
0,3673 0 0,010
0,3369 0 0,005
0,8090 0,5878 0,3090 0
0,000
↑
0,7
0,8
0,9 1 tempo
7
1,3584
1,2460
1,1429
1,0483
0,9616
0,8820
8
0,9869
0,9053
0,8304
0,7616
0,6986
0,6408
E.B.Hauser – Cálculo Numérico
--Burden, Richard L., Faires, J. Douglas.”Analise Numérica”, São Paulo, SP, 2003,
THOMSON.
--Cláudio, Dalcidio M, Marins, Jussara M, Cálculo Numérico Computacional, São Paulo,
SP,
Ed. Atlas, 1994.
--Cunha, M.Cristina, Métodos Numéricos, Campinas, SP, 2003, Ed. Unicamp.
--Zill, Denis G., Cullen, Michael R., Equações Diferenciais, Vol.2, Sâo Paulo, SP, 2002,
Makron Books.
--Schleider, Maria Amélia N, Cunha, Maria Cristina, 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, 2003,
Palgrave Macmillan.
83
Equação
Diferenças Finitas
h, k=tamanho do passo na direção x e na direção y (ou t)
Centradas, com h=k
Laplace
∂ 2u
∂ 2u
( x, y ) = 0
( x, y ) +
2
2
∂x
∂y
ui , j =
Célula Computacional
•
u i, j + 1
ui + 1, j + ui −1, j + ui , j − i + ui , j + i
4
•
u i − 1, j
o
u i, j
•
u i + 1, j
•
u i, j − 1
∂u
( x ,t ) =
α
( x ,t )
2
∂
t
∂x
Material
Prata
Cobre
Alumínio
Ferro
Concreto
o
progressivas
Calor
∂ 2u
α >0
1,7
1,5
0,85
0,15
0,005
⎛
α 2k ⎛
2α 2 k ⎞⎟
⎞
ui , j + 1 = ⎜ 1 −
ui , j +
+u
⎟
⎜u
⎜
⎟
i
1
,
j
i
1
,
j
+
−
2
2
⎠
⎝
h
h
⎝
⎠
1
k
Converge se
<
h2 2
u i, j+1
•
u
i−1, j
•
u
i, j
•
u
i+ 1, j
centradas
Onda
∂ 2u
∂ 2u
2
β
( x ,t ) =
( x ,t )
∂t 2
∂x 2
⎛
β 2 k 2 ⎞⎟
β 2k 2 ⎛
⎞
+u
ui , j + 1 = 2⎜ 1 −
ui , j +
u
⎜
⎟ − u i , j −1
⎜
i − 1, j ⎠
2 ⎟
2 ⎝ i + 1, j
h
h
⎝
⎠
u
•
u i − 1, j
Se β = 1 e
k2
= 1 , ui , j + 1 = ui −1, j + ui + 1, j − u
i , j −1
h2
84
o
i, j + 1
•
u i, j
•
u i, j − 1
•
u i + 1, j
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. São Paulo :
Thomson, 2003.
CHAPRA, Steven C., CANALE, Raymond P. Numerical Methods for Engineers with
Programming and Software Applications. 4.ed. Boston : McGraw-Hill, 2002.
CLÁUDIO,
Dalcidio
Moraes,
MARINS,
Jussara
Numérico Computacional. 2.ed. São Paulo : Atlas, 1994.
Maria.
Cálculo
CUNHA, M.Cristina, Métodos Numéricos, Campinas, SP, 2003, Ed. Unicamp.
GANDER, Walter.; HREBICEK, Jiri. Solving Problems in Scientific Computing Using
Maple and Matlab: Berlin, Springer-Verlag, 1995.
HUMES, Ana Flora P. de Castro et al. Noções de Cálculo Numérico. São Paulo: McGrawHill, 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.
numérico: aspectos teóricos e computacionais. São Paulo : McGraw-Hill, 1997.
Cálculo
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, 2003.
85
Cálculo Numérico - Formulário
⎡
⎢⎣
1a) # F = 2( β − 1 )β
⎡
⎛
xk + 1 − xk
t −1( M − m + 1 )⎤ + 1 , F ( β , t , m, M ) ; 1b) DIGSE( x , x
k k + 1 ) = − ⎢0 ,3 + log ⎜⎜ µ +
⎥⎦
xk + 1
⎢⎣
⎝
2) Seja p( x ) = a n x n + a n −1 x n −1 + ... + a3 x 3 + a 2 x 2 + a1 x + a0
p( x ) = (((...( a n x + a n −1 ) x + a n − 2 )x + ...a 2 )x + a1 ) x + a0
123
n −1
2
b) Huat: Se p(0) ≠ 0 e para algum k = 1, 2, ...,n-1, ak ≤ ak −1ak + 1 , então p tem raízes complexas.
f ( xk )
3) Método de Newton-Raphson:
xk + 1 = xk −
, k = 0 ,1,2 ,..., máx f ' ( xk ) ≠ 0
f ' ( xk )
a) Horner:
4) Sistema de n equações lineares: AX=B, A = ⎡a ij ⎤ , X = ⎡ x ⎤ e B = ⎡ b i ⎤ com i,j = 1, 2, ..., n.
⎥⎦
⎥⎦
⎢⎣
⎢⎣ i ⎥⎦
⎣⎢
4a) NORM A =
det A
onde α k = a 2 + a 2 + L a 2 , para k = 1, 2, ..., n.
kn
k2
k1
α 1α 2 Lα n
n
4b) A matriz A é Diagonal Dominante se a
>
ii
∑
∀ i , j = 1,2 ,.., n.
aij
j =1
j ≠i
4c)Gauss -Jacobi e Gauss -Seidel: ∀ i= 1, 2, ..., n e k = 1, 2, ..., máx,
⎞
⎛
⎞
⎛
⎟
⎜
⎟
⎜
n
−
i
1
n
⎟
⎟
1 ⎜
1 ⎜
x
bi −
a ij x j ⎟
x
−
=
=
bi −
a ij x j
a ij x j ⎟
⎜
⎜
i k + 1 aii ⎜
i k + 1 aii ⎜
k ⎟
k +1
k ⎟
j =1
j =1
j =i +1
⎟
⎜
⎟
⎜
j ≠i
⎠
⎝
⎠
⎝
0
5) Diferenças Finitas Ascendentes: ∀ i= 0, 1, 2, ..., n ,seja ∆ yi = yi . Para k = 1, 2, ..., n
∑
∑
a diferença finita de ordem k é
∆k yi = ∆k −1 yi + 1 − ∆k −1 yi , com i= 0, 1, 2, ..., n-k .
6) Diferenças Divididas: ∀ i= 0, 1, 2, ..., n ,seja
∆ k yi =
∆ k −1 yi +1 − ∆ k −1 yi
xi +k − xi
∑
∆0 yi = yi .
Para
k = 1, 2, ..., n, a diferença dividida de ordem k é
, com i= 0, 1, 2, ..., n-k .
7) Polinômio Interpolador de Newton para Diferenças Finitas Ascendentes:
p( x ) = y o +
Se z =
( x − xo )
h
∆y0 +
( x − x o )( x − x )L ( x − x
)
( x − x o )( x − x ) 2
1
n − 1 ∆n y
1 ∆ y +L+
o
o
2
n
2! h
n! h
∆2 yo
∆3 yo
∆n yo
( x − xo )
, p( z ) = yo + z ∆y + z( z − 1 )
.
+
−
−
+
L
+
−
L
−
−
z
(
z
1
)
(
z
(
n
1
))
z
(
z
1
)(
z
2
)
0
h
n!
3!
2!
8) Polinômio Interpolador de Newton para Diferenças Divididas
2
p( x ) = yo + ( x − xo )∆y + ( x − xo )( x − x ) ∆ yo + ...+ ( x − xo )( x − x )...(x − x
)∆ n yo
0
1
1
n−1
86
⎞⎤
⎟⎥
⎟⎥
⎠⎦
Cálculo Numérico - Formulário
p
9) Ajustamento a um Polinômio de grau p: Se Y ( x ) = ao + a x + a x 2 + L + a p x é a função que ajusta os pontos
1
2
( xi , yi ) , i = 0 ,1,...n
e
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
∑ xi2
∑ xi3
n+1
∑ xi
M
p
xi
∑
∑
n
∑ ∑
=
i =0
∑ xi
∑ xi2
M
p +1
xi
p < n , então, para
∑
M
p+2
xi
∑ xi3
∑ xi4
∑
L
L
L
9a)Ajuste Linear e Quadrático Y ( x ) = a o + a x , e
1
⎛ n+1
⎜
⎛ n+1
⎛
xi ⎞⎟ ⎛ a ⎞
yi ⎞⎟
⎜
⎜
⎜ 0⎟ = ⎜
, e ⎜ xi
⎜⎜
2
⎟
⎟
⎜
⎜
⎟
a
xi yi
xi
xi ⎟ ⎝ 1 ⎠
⎜
⎝
⎠
⎝
⎠
xi2
⎜
⎝
∑
∑
∑
∑
∑
∑
∑
⎞
⎟
⎟
⎟
⎟
M
p ⎟⎟
xi y i
⎠
∑ xip ⎞⎟⎟⎛⎜⎜ ao ⎞⎟⎟ ⎛⎜⎜ ∑ yi
∑ xip +1 ⎟⎟⎜⎜ a1 ⎟⎟ = ⎜⎜ ∑ xi yi
L
M
p+3
xi
, os parâmetros a 0 , a 1 ,..., a p constituem a solução de:
∑
M
M
2 p ⎟⎟⎜⎜ a ⎟⎟
xi
⎠⎝ p ⎠
⎜
⎜
⎝
∑
Y ( x ) = ao + a x + a 2 x 2
1
2
⎞
xi
xi ⎟ ⎛⎜ a o ⎞⎟
⎟ ⎜
⎟
xi3 ⎟ ⎜ a ⎟ =
xi2
1
⎟ ⎜
⎟
xi4 ⎟ ⎜ a ⎟
xi3
⎠ ⎝ 2⎠
∑
∑
∑
∑
∑
∑
⎛
⎜
⎜
⎜
⎜
⎝
∑ yi ⎞⎟⎟
∑ xi yi ⎟
∑ xi2 yi ⎟⎠
10) Integração: Para h = ( x n − xo ) / n , xi = xo + ih e y i = f ( xi ) , i = 0,1,...n,
∫ f ( x )dx ≅ 2 [ yo + 2( y1 + y2 + L + yn −1 ) + yn ]
xn
10a) Trapézios:
h
xo
ET ≤
h2
( xn − x0 ) max
f'' ( x )
12
x∈[ xo , xn ]
xn
10b) Simpson ( n par) :
∫
f ( x )dx ≅
ou
ET ≈
( xn − x0 )
max ∆2 yi .
12
h⎡
y o + 4( y1 + y 3 + y 5 + L + y n −1 ) + 2( y 2 + y 4 + y6 + L + y n − 2 ) + y n ⎤⎥
3 ⎢⎣
⎦
xo
ES ≤
h4
( x n − x0 )
max
f '''' ( x )
180
x∈[ xo ,x n ]
⎧⎪ y' = f ( x , y )
11) PVI: Sejam ⎨
⎪⎩ y( xo ) = yo
ou
ES ≈
( xn − x0 )
max ∆4 yi .
180
y ( xi ) ≅ y i com xi = xo + ih , i= 0,1,...,n e h = ( xn − xo ) / n .
11 a)Euler: yi + 1 = yi + hf ( xi , yi ) , i= 0,1,...,n-1.
a
h
11b) Runge-Kutta de 2 Ordem : y
= yi + (k1 + k 2 ) , k1 = f ( xi , yi ) e k 2 = f ( xi + h , yi + hk1 ) , i= 0,1,...,n-1.
i +1
2
⎧⎪ y' ' + P( x ) y' +Q( x ) y = f ( x )
, y ( xi ) ≅ y i com xi = xo + ih , i= 0,1,...,n e h = ( xn − xo ) / n .
⎪⎩ y( xo ) = yo , y( xn ) = y n .
12) PVC: Sejam ⎨
⎛
⎝
Diferenças Finitas: ⎜ 1 +
(
)
h
⎛ h
⎞
⎞
P( xi ) ⎟ yi + 1 + − 2 + h 2 Q( xi ) yi + ⎜ 1 − P( xi ) ⎟ yi −1 = h 2 f ( xi ) , i= 1,...,n-1
2
2
⎝
⎠
⎠
87
Cálculo Numérico - Formulário
Derivada
Diferença Finita
h, k=tamanho do passo na direção x e na direção y (ou t)
y´( xi )
yi + 1 − yi
avançada ,
h
y´´( xi )
yi + 1 − 2 y i + yi − 1
centrada
h2
y' ' ' ( xi )
yi + 2 − 2 yi + 1 + 2 yi − 1 − yi − 2
y IV ( xi )
yi + 2 − 4 yi + 1 + 6 y i −4 yi − 1 + yi − 2
centrada
h4
(
∂u
xi ,t j
∂t
(
2h 3
)
∂ 2u
xi ,t j
∂ x2
(
yi + 1 − yi − 1
centrada ,
2h
centrada
u ,
−u ,
i j +1
i j −1
centrada
2k
u
)
∂ 2u
xi , y j
∂ y2
)
i +1 , j
− 2u , + u
i j
i −1 , j
centrada
h2
u ,
− 2u , + u ,
i j +1
i j
i j −1
centrada
2
k
88
yi − yi − 1
atrasada
h
Cálculo Numérico - Formulário
Equação
Diferenças Finitas
h, k=tamanho do passo na direção x e na direção y (ou t)
Centradas, com h=k
Laplace
∂ 2u
∂x 2
( x, y ) +
∂ 2u
∂y 2
( x, y ) = 0
ui , j =
Célula Computacional
•
u i, j + 1
ui + 1, j + ui − 1, j + ui , j − i + ui , j + i
4
•
u
i − 1, j
o
u
i, j
•
u
i + 1, j
•
u i, j − 1
∂ 2u
∂u
α
( x ,t )
( x ,t ) =
∂t
∂x 2
Material
Prata
Cobre
Alumínio
Ferro
Concreto
o
progressivas
Calor
α >0
1,7
1,5
0,85
0,15
0,005
⎛
α 2k ⎛
2α 2 k ⎞⎟
⎞
+u
u
ui , j +
ui , j + 1 = ⎜ 1 −
⎜
⎟
⎜
i − 1, j ⎠
2 ⎝ i + 1, j
2 ⎟
h
h
⎠
⎝
k
1
<
Converge se
2
2
h
u i, j+1
•
u i−1, j
•
ui, j
•
u i+1, j
centradas
o
Onda
∂ 2u
∂ 2u
β2
( x ,t )
( x ,t ) =
∂x 2
∂t 2
⎛
β 2 k 2 ⎞⎟
β 2k 2 ⎛
⎞
+u
ui , j +
u i , j + 1 = 2⎜ 1 −
⎜ u
⎟ − u i , j −1
⎜
⎟
+
−
i
1
,
j
i
1
,
j
2
2
⎝
⎠
h
h
⎝
⎠
Se
β =1 e
k2
= 1 , ui , j + 1 = ui − 1, j + ui + 1, j − u i , j −1
h2
89
u i, j + 1
•
u
i − 1, j
•
u
i, j
•
u
i, j −1
•
u
i + 1, j
Aula de Laboratório 1 utilizando Maple
Operações, Funções
Básicas, Constantes
Adição
Subtração
Multiplicação
Divisão
Potenciação
Valor Absoluto de x
Raiz Quadrada de x
Exemplos
Notação
+
*
/
^
abs ( x )
sqrt (x)
Raiz n-ésima de x
Fatorial de n
π
Infinito
Unidade Imaginária
Número de Euler: e
Função Exponencial
x^(1/n)
n!
Pi
infinity
sqrt(-1) ou I
exp(1) ou E
exp(x) b^x
Logaritmo Natural
Logaritmo de base b
ln(x)
log[b](x) = ln(x)/ln(b).
sin(x) cos(x) tan(x)
sec(x) csc(x) cot(x)
sinh(x) cosh(x) tanh(x)
sech(x) csch(x) coth(x)
arcsin(x) arccos(x) arctan(x)
arcsec(x) arccsc(x) arccot(x)
arcsinh(x) arccosh(x) arctanh(x)
arcsech(x) arccsch(x) arccoth(x)
Funções Trigonométricas ,
Hiperbólicas e suas inversas
(argumento x em Radianos)
> 3+5;
> 758-195;
> 7.5*8;
> 39/13;
> 2^12;
> abs(-7);
>sqrt(7835);
> sqrt(7835.);
> 27^(1/3);
> 28!;
> Pi;
> infinity;
> sqrt (-1);
> exp(1);
>exp(3.5);
>7^x;
> ln (7);
> log[3](10);
> sin(Pi/2);
> cos(75.3);
>arccosh(0);
Alguns Comandos
1) > evalf(expr);
avalia expr utilizando aritmética de ponto flutuante com precisão determinada pela variável
global Digits.
> evalf(27^(1/3));
> evalf( Pi);
> evalf(arccosh(0));
2) > evalf(expr, n);
calcula expr com n dígitos de precisão.
> evalf( Pi, 21);
> evalf(ln (7),8);
3) > Digits := n;
ajusta para n o número de dígitos utilizados em ponto flutuante.
O valor "default" de Digits é 10.
> Digits:=15;
> cos(75.3);
90
4) Se usarmos : no final do comando, o mesmo é executado, porém o resultado não é
mostrado. %; mostra o conteúdo do último output e %%; do penúltimo.
> log[3](10):
> %;
> evalf(%%);
5) O Maple pode trabalhar com inteiros muito grandes
> 253!;
O comando length (%) mostra o número de dígitos de 253!
> length(%);
6) Para fazer uma atribuição à variável z e posteriormente apagar o conteúdo de z:
> z := 7;
> z;
> z := 'z';
>z;
7) O modo texto pode ser obtido clicando em
na barra de menu ou # pode ser usado
para fazer comentários
> # o conteúdo de z também pode ser deletado usando unassign ( 'z' );
> z :=7;
> z;
> unassign ( 'z' );
>z;
8) restart; reinicializa o MAPLE
Exercícios:
1) Calcular o valor de :
⎛ 30π ⎞
a) seno(30), seno(300) = seno ⎜
⎟
⎝ 180 ⎠
a + 5b 3
b)
para a = -3,02 , b = 10, c = 10,3 e d = 7,1. Armazenar o resultado na
ln( 50c / d )
variável v.
c) secante hiperbólica de v calculado no item b.
2) O valor de a = e π 163 é um número inteiro? Estimar a com 18, 29, 30, 31 e 57
dígitos. Comparar os resultados.
3) Calcular 325! pelo Maple e pela sua calculadora. Analisar os resultados. Qual o número
de dígitos de 325! ?
27
24
4) O que é maior 22 ou 25 ?
5) É possível calcular 3333333333 3333333333 no Maple?
91
Aula de Laboratório2 utilizando Maple
Objetivos : Localizar as raízes da equação f(x)=0(algébrica ou transcendente) e calculá-las
utilizando comandos do Maple e o método da Bissecção.
Aplicação: O fator de atrito λ para um duto retangular , segundo Maubach, é dado
por:
1
= 2,035 log 10 Re λ − 0,989 , onde Re é o número de Reynolds. Determinar λ para os
λ
(
)
seguintes valores de Re : 10000, 20000, 50000, 100000, 1000000.
1
2
Comandos
f:=x->f(x);
Comentários
plot(f(x),x= a..b,y=c..d);
3
plot({f(x),g(x)},x= a..b,y=c..d);
4
solve(f(x)=0, x);
fsolve(f(x)=0, x);
fsolve(f(x)=0, x=a ..b);
5
For i from io by k to in do
comandos
od;
Define a função de variável x , f(x)
Plota a função f(x) sendo:
x - variável
y - parâmetro na vertical(opcional)
a,b,c,d :parâmetros a serem especificados
Plota as funções f(x) e g(x) num mesmo
sistema de eixos.
Para opções utilizar o Help: < ?plot
Calcula as raizes da equação f(x)=0.
Para opções utilizar o Help:
< ?solve
< ?fsolve
Comando de repetição onde:
i - variável
io - valor inicial
in - valor final
k - passo
6
if condição then
comando1
else
comando2
fi;
Exemplos
Comando de teste (condicional)
1) >f:=x->sin(x/2);
> f (Pi); f (4);
2) >plot(f(x), x=-3*Pi..3*Pi);
> plot(f(x), x=-15..15, y=-2..2);
3) > plot({x^2-5*x, sqrt(x+1)}, x= -2..7, y= -7..7);
4) vide exercício 1(página2)
5) > for i from 0 by 2 to 10 do
>
print(i^2)
> od;
6) > a:=-23;
> if a<0 then
> print(-a);
> else
> print(a);
> fi;
92
Exercícios:
i) Localizar graficamente as raizes de f(x)=0
ii) Calculá-las utilizando os comandos dados em 4.
ii) Utilizando 5 e 6, implementar o Método da Bissecção para calcular a menor raiz real
positiva de f(x) = 0 em [a , b]. Para xm=(a+b)/2, a cada bissecção imprimir: a, xm, b,f(xm).
Critério de Parada: mínimo 10 bissecções.
1) xcos(x)-1=0
Solução:
> f:=x->x*cos(x)-1;
> plot(f(x), x=-20*Pi..20*Pi);
> plot({cos(x),1/x}, x=-7*Pi..5*Pi, y=-5..5, title=` cos(x) , 1/x`);
Neste caso vemos que existem infinitas raizes.
>solve(f(x)=0,x);
> evalf(%);
> fsolve(f(x)=0,x);
> fsolve(f(x)=0,x=-9..-7);
> fsolve(f(x)=0,x=4..5);
> fsolve(f(x)=0,x=13..15);
Vamos implementar o método da bisseção para calcular a raiz que está no intervalo [4 , 5]:
[> a:=4; b:=5;
⎡ > for i from 0 to 10 do
⎢
a; xm := evalf((a + b)/2); b; evalf(f(xm));
⎢
⎢>
if evalf((f(a ) * f(xm))) > 0 then a := xm;
⎢
else b := xm;
⎢>
⎢>
fi;
⎢
⎢⎣ > od;
Obs: Não foi considerado o caso f(xm)=0.
2) ln(x+11) -2x = 0
3) x 3 − 5 ,875 x 2 − 40 ,135 x − 25 ,62875 = 0
93
Aula de Laboratório3 utilizando Maple
Objetivo: Enumerar, localizar e calcular as raízes de polinômios. Analisar o comportamento
gráfico( raízes reais de multiplicidade par e ímpar, por ex). Utilizar método de NewtonRaphson para cálculo das raízes da equação f(x)=0(algébrica ou transcendente)
Aplicação: A lei para um gás ideal, PV = nRT é um conceito conhecido. Geralmente,
utilizam-se estimativas para as relações P-V-T. A equação de Beattie-Bridgeman é um
exemplo:
a
a
RT a1
P=
+
+ 2 + 3
(1)
V V 2 V3 V 4
e pode ser reescrita como um polinômio de quarto grau:
PV 4 + RTV 3 − a1V 2 − a 2 V − a = 0
(2)
3
Para um gás particular, a1=-1,06, a2 = 0,057 e a3=-0,00011. Considerando a pressão
P=25 atm , a temperatura é T = 293oK e R = 0,0,082l-atm/Kg.mol e substituindo esses
dados em (2) , o volume V é obtido de :
25V 4 + 24,03V 3 + 1,06V 2 − 0,057V + 0,00011 = 0
(3)
Determinar V utlizando o método de Newton-Raphson com V0=
Comandos
1 f:=x->f(x);
2
plot(f(x),x= a..b, y=c..d);
Comentários
Define a função de variável x , f(x)
Plota a função f(x) sendo:
x - variável
y - parâmetro na vertical(opcional)
a,b,c,d
:parâmetros
a
serem
especificados
Plota as funções f(x) e g(x) num
mesmo sistema de eixos.
Para opções utilizar o Help: < ?plot
Encontra todas as
raízes
do
polinômio p(x)
Determina a derivada de f(x)
Determina a derivada de quarta ordem
de f(x)
3 plot({f(x),g(x)},x= a..b, y=c..d);
4 fsolve(p(x)=0, x,complex);
5 diff(f(x),x);
diff(f(x),x$4);
diff(f(x),x,x,x,x)
6 subs(x=a, expr);
RT 0,082x 293
=
= 0,961 .
P
25
ou
substitui x por a na expressão expr
Exemplos
1) Raiz de multiplicidade par
> p1:=x->5*x^4-12.50*x^3-18.0375*x^2+32.31250*x+33.411125;
> plot(p1(x),x=-3..4);
> fsolve(p1(x)=0,x);
94
2) Raiz de multiplicidade ímpar
> p2:=x-> 2*x^4+6.8*x^3-16.56*x^2-86.756*x-85.1690;
> plot(p2(x), x=-5..5);
> fsolve(p2(x)=0,x);
3) Raizes complexas
> p3:=x-> x^3-2.8*x^2+8.2*x+15.6;
> plot(p3(x),x=-2..5);
> fsolve(p3(x)=0,x);
> fsolve(p3(x),x, complex);
> p4:=x-> x^4-2*x^3+11*x^2-18*x+18;
> plot(p4(x),x=-3..3,y=-5..30);
> fsolve(p4(x)=0,x);
> fsolve(p4(x),x, complex);
4) > subs( x=2, x^2+x+1 );
5) > diff(x^7,x);
> diff(x^7,x,x);
> diff(x^7,x$3);
6) Método de Newton-Raphson
>f := x ->cosh(x)*cos(x)-1;
> plot(f(x), x=-3*Pi..3*Pi);
> plot({cos(x), 1/cosh(x)}, x=-10..10);
> x[0]:=4.5;
> for i from 0 to 5 do
x[i+1]:=evalf(x[i]-(f(x[i])/subs(x=x[i],diff(f(x),x))))
od;
> abs(x[5]-x[6]);
> Digits:=20;
> abs(x[5]-x[6]);
95
Aula de Laboratório4 utilizando Maple
Objetivos : Resolver sistemas de equações lineares utilizando
diversas opções de comandos do Maple. Analisar graficamente
a solução , quando possível, no plano e no espaço.
Aplicação:
Equacionado o circuito resistivo mostrado na figura acima pelo
método dos nós, obtemos:
⎧ VA VA − VB VA −10
+
=0
⎪ 8 +
1
5
⎪
⎪ VB − VA VB −10 VB − 8 VB − VC VB − 2 VB − 0
+
+
+
+
+
=0
⎨
1
5
8
1
5
8
⎪
⎪ VC − 8 VC − VB VC − 2
+
=0
⎪ 8 +
1
5
⎩
⎧I1 = VA − VB
Aplicando a Lei de Ohm: ⎨
⎩I 2 = VC − VB
Comandos
1
array (1..m, 1..n)
2
augment( M1, M2);
3
4
5
6
7
evalm(expressão matricial);
gaussjord(M);
inverse(M);
linsolve(A,v);
matrix(m , n, [x11,x12,...xmn]);
8
solve(eqn, var);
9
10
11
12
12
transpose(M);
vector(n, [x1, x2, ..., xn]);
with(linalg);
with(plots);
A&*S;
. Determinar I1 e I2.
Comentários
Cria uma matriz cujos os elementos estão definidos.Vide
help para opções
Cria uma nova matriz colocando M1 à esquerda da M2 .
As matrizes devem ter o mesmo número de linhas.
Avalia uma expressão contendo matrizes.
Aplica o método de Gauss-Jordan
Encontra a matriz inversa de M.
Calcula um vetor x que satisfaça a equação Ax=v
Cria uma matriz com m linhas e n colunas com
elementos x11, ...,xmn.
Resolve simbolicamente equações eqn para variável var.
Para opções vide help.
Determina a matriz transposta de M.
Cria um vetor de n elementos x1, ..., xn..
Carrega a biblioteca de álgebra linear do Maple V.
Carrega o package gráfico
Expressa multiplicação de matrizes (não comutativa)
Exemplos :
> with(plots);
> with(linalg);
1) Sistema Linear Compatível Determinado (possui soluções)
> solve({2*x+3*y=18,3*x+4*y=25},{y,x});
> plot({(18-2*x)/3, (25-3*x)/4}, x=0..6);
> solve({3*x+2*y-5*z=8,-x+2*y+z=2,-x+2*y+3*z=4},{z,y,x});
> plot3d({(8-3*x-2*y)/(-5),2+2*x-2*y,(4+x-2*y)/3}, x=0..5,y=0..5, axes=box);
96
2) Sistema Linear Compatível Indeterminado (possui infinitas soluções)
> solve({2*x+y=50,4*x+2*y=100},{y,x});
> plot({50-2*x, (100-4*x)/2}, x=0..6);
> plot3d({50-2*x, (100-4*x)/2}, x=0..6,y=0..6, axes=box);
> solve({3*x+y+2*z=0,-9*x-3*y-6*z=0},{z,y,x});
> plot3d({(-3*x-y)/2,2+2*x-2*y,(9*x+3*y)/(-6)}, x=-2..2,y=-2..2, axes=box);
outra maneira:
> solvefor[t]( x+y=1, x-y+z*t=3 );
> solvefor[x]( x+y=1, x-y+z*t=3 );
3) Sistema Linear Incompatível (não admite solução)
> solve({x+y=3, x+y=-5},{y,x});
> plot({3-x, -x-5}, x=-10..10);
> plot3d({3-x, -x-5}, x=-10..10,y=-6..6, axes=box);
4) Resolver o sistema e calcular o resíduo produzido pela solução encontrada:
⎧2 x 1 + x 2 − 0.1x 3 + x 4 = 2.7
⎪
⎪0.4 x 1 + 0.5x 2 + 4 x 3 − 8.5x 4 = 21.9
⎨
⎪0.3x1 − x 2 + x 3 + 5.2 x 4 = −3.9
⎪⎩x 1 + 0.2 x 2 + 2.5x 3 − x 4 = 9.9
> A := array([ [2,1,-.1,1], [.4,.5,4,-8.5], [.3,-1,1,5.2], [1,.2,2.5,-1] ]);
> F := vector([2.7,21.9,-3.9,9.9]);
> S := linsolve(A,F);
> F1 :=evalm(A&*S);
> Resíduo := evalm(F-F1);
5) Resolver o sistema utilizando o Método de Gauss-Jordan (Matrix Inversa).
⎧2 x 1 + x 2 + 7 x 3 = b1
⎪
⎨ x 1 + 3x 2 + 2 x 3 = b 2
⎪5x + 3x + 4 x = b
2
3
3
⎩ 2
onde: a) b1 =16 b2 = -5 b3=11 b) b1 =25 b2 = -11 b3 = -5 c) b1 =3 b2 = 5 b3 = -5
> B := matrix(3,3, [2,1,7,1,3,2,5,3,4] );
> v1 :=vector(3,[16,-5,11] );
> v2 :=vector(3, [25,-11,-5] );
> v3 :=vector(3, [3,5,-5] );
> augment( B,v1,v2,v3 );
> gaussjord(%);
97
Aula de Laboratório5 utilizando Maple
Objetivos: Determinar o polinômio interpolador utilizando a resolução de sistemas de
equações lineares. Detectar se um sistema linear é mal condicionado. Gerar matrizes e vetores
utilizando comandos de teste e repetição.
Aplicação: Na modelagem de um processo de combustão é necessário expressar a entalpia
(E) como uma função da temperatura (T). Considerando os dados tabelados, estimar a
entalpia para uma temperatura de 150 o F .
E(Btu / lb) 60 80 100 120 140 160 180
T(o F) 0.0 17.2 45.2 92.9 178.8 349.4 764.3
Exemplo 1
Construir uma matriz M = [mij] , de ordem 7x7, tal que mij = i + 3j se i é igual a j e mij = 5ij
se i for diferente de j. M é simétrica? Calcular o determinante de M . M é inversível?
Solução:
> with(linalg);
>M := array(1..7,1..7);
⎡> for i to 7 do
⎢> for j to 7 do
⎢
⎢>
if i <> j then M[i, j] := 5 * i * j
⎢
else M[i, j] := i + 3 * j
⎢>
⎢>
fi
⎢
⎢> od
⎢> od,
⎣
> print(M);
> M2 :=transpose(M);
> det(M);
> M1 :=inverse(M);
> evalf(%);
> evalm(M1 &* M);
> evalm(M1 + M2);
> evalf(%);
Exemplo2:
O alongamento de uma mola foi medido em função da carga aplicada. Obteve-se:
c arg a (kg )
2 4 6 8
alongamento(cm) 1,0 2,5 5,0 6,3
Estimar o alongamento para o caso de ser aplicada uma carga de 7Kg. Analisar
graficamente.
Solução:
98
Procuramos p( x ) = a 0 + a 1x + a 2 x 2 + a 3 x 3 .
a) uma maneira de resolver no maple::
>solve ({a0+a1*2+a2*2^2+a3*2^3=1, a0+a1*4+a2*4^2+a3*4^3=2.5,
a0+a1*6+a2*6^2+a3*6^3=5, a0+a1*8+a2*8^2+a3*8^3=6.3}, {a0,a1,a2,a3});
> p:=x-> -.04583333333*x^3+.675*x^2-2.016666667*x+2.7;
> validade:= [p(2), p(4), p(6), p(8)];
> p(7);
Análise Gráfica:
> with(plots);
> plots[display]({plots[pointplot]([2,1,4,2.5,6,5,8,6.3]),plot(p(x),x=-5..15, y=0..10)});
b) outra forma de resolver o sistema:
> A := array([ [1,2,2^2,2^3], [1,4,4^2,4^3], [1,6,6^2,6^3], [1,8,8^2,8^3] ]);
> F := vector([1,2.5,5,6.3]);
> S := linsolve(A,F);
Exemplo3: Considerar o sistema linear:
⎧x1 + 1 / 2 x 2 + 1 / 3x 3 + 1 / 4 x 4 + 1 / 5x 5 = 137 / 60
⎪
⎪⎪1 / 2 x1 + 1 / 3x 2 + 1 / 4 x 3 + 1 / 5x 4 + 1 / 6 x 5 = 87 / 60
(*) ⎨1 / 3x1 + 1 / 4 x 2 + 1 / 5x 3 + 1 / 6 x 4 + 1 / 7 x 5 = 459 / 420
⎪1 / 4 x + 1 / 5x + 1 / 6 x + 1 / 7 x + 1 / 8x = 743 / 840
1
2
3
4
5
⎪
⎪⎩1 / 5x1 + 1 / 6 x 2 + 1 / 7 x 3 + 1 / 8x 4 + 1 / 9 x 5 = 1879 / 2520
a) (*) é bem condicionado ou mal condicionado? Porquê? O que isso significa?
b) calcular a solução de (*).
Solução:
> with(linalg);
> C := hilbert(5);
> cond(C);
> ?cond;
> F := vector([137/60, 87/60, 459/420, 743/840, 1879/2520]);
> L := linsolve(C,F);
> evalm(C &* L);
99
Exercício:
Consideremos a matriz de Hilbert de ordem n, Hn = [ hij ], com seus elementos genéricos definidos
por
1
hij = =
, 1≤ i , j ≤ n
i + j −1
Hn é um exemplo clássico de matriz mal condicionada.
a) Indicar se a afirmação abaixo é verdadeira ou falsa e justificar.
“ Quanto maior for n, mais mal condicionada é Hn . “
b) Resolver o sistema H5.X = B , onde B = [ bi ], i= 1, 2, 3, 4 ,5 é o vetor definido por:
n
1
bi = ∑
i + j −1
j =1
c) Calcular o determinante de Hn , a matriz inversa de Hn e o produto HnHn-1.
c) Calcular todas a s raizes do determinante da matriz sI- Hn, s numero real qualquer e
I = matriz identidade
100
Aula de Laboratório 6 utilizando Maple - Interpolação
EXEMPLO:
O alongamento de uma mola foi medido em função da carga aplicada. Obteve-se:
c arg a (kg )
2 4 6 8
alongamento(cm) 1,0 2,5 5,0 6,3
Estimar o alongamento para o caso de ser aplicada uma carga de 7Kg.
Procuramos p ( x) = a0 + a1 x + a2 x 2 + a3 x 3 Estudamos alguns algortmos para determinação dos
coeficientes de p.
1) via resolução de um sistema linear:
1.a
> with(linalg):
> A := array([ [1,2,2^2,2^3], [1,4,4^2,4^3], [1,6,6^2,6^3],
[1,8,8^2,8^3] ]);
> F := vector([1,2.5,5,6.3]);
> S := linsolve(A,F);
1.b
> solve({a0+a1*2+a2*2^2+a3*2^3=1,a0+a1*4+a2*4^2+a3*4^3=2.5,a0+a1*6
+a2*6^2+a3*6^3=5, a0+a1*8+a2*8^2+a3*8^3=6.3},{a0,a1,a2,a3});
2) Polinomio interpolador de Newtons para Diferenças Finitas Ascendentes
>
>
>
>
>
>
>
vx:=vector(4,[2,4,6,8]);
vy:=vector(4,[1,2.5,5,6.3]);
df:=array(1..3);
for i from 1 to 3 do df1[i]:=vy[i+1]-vy[i] od;
for i to 2 do df2[i]:=df1[i+1]-df1[i] od;
for i to 1 do df3[i]:=df2[i+1]-df2[i] od;
p:=x->1+((x-2)*df1[1])/(2)+((x-2)*(x-4)*df2[1])/(2!*(2^2))+((x-2
)*(x-4)*(x-6)*df3[1])/(3!*(2^3));
> simplify(p(x));
> p(7);
3) Via comando´interp´do maple
> vx:=vector(4,[2,4,6,8]):
> vy:=vector(4,[1,2.5,5,6.3]):
> pa:=interp(vx,vy,x):
> pa(7);
> paa:=unapply(pa,x):
> paa(7);
validade:
> [paa(2),paa(4),paa(6),paa(8)];
análise gráfica
> plots[display]({plots[pointplot]([2,1,4,2.5,6,5,8,6.3]),plot(paa
(x),x=-5..15, y=0..10)});
101
EXERCICIOS:
1) A tabela abaixo fornece a demanda diária máxima de energia
elétrica na Cidade A no mês de março.
x(dia )
1 11 21 31
y (demanda − MW ) 10 15 20 13
a) Determinar o polinômio interpolador p e verificar sua
validade.
b) Representar graficamente p e os dados tabelados num mesmo
sistema de eixos.
c) Estimar a demanda máxima e a data em que ocorreu.
>
2) A que temperatura a água entra em ebulição no Pico da Bandeira
(altitude de 2890m), sabendo-se que o ponto de ebuliçao da água
varia com a altitude, conforme tabela abaixo:
altitude(m)
2600 2700 2800 2900 3000
o
ponto de ebulição( C ) 91.34 91.01 90.67 90.34 90.00
>
3) A velocidade do som na água varia com a temperatura conforme
tabela:
temperatura ( oC ) 86.0 93.3 98.9 104.4 110.0
velocidade(m / s ) 1.552 1.548 1.544 1.538 1.532
a)Determinar o polinômio interpolador de Newton para diferenças
divididas,p, e verificar sua validade.
b) Estimar a velocidade do som se a temperatura da água for de
100 graus centígrados.
102
c)Representar graficamente p e os dados tabelados num mesmo
sistema de eixos.
Comandos Utilizados
diff(expr, var);deriva uma expressão " expr" na variavel "var"
diff(expr, var, var$n);deriva uma expressão expr na varivel var, var$n
siginifica derivada de ordem n.
interp([exprx1,...,exprxn+1], [expy1,...,expryn+1], var);computa um polinomio na
variável var de grau at n, no qual representa o polinomio
interpolador dos valores expx e expy .
plots[display]({plots[pointplot]([p1,p2,..,pn]),plot(f(x), x=-m..m)});plota graficos;
display, para plotar dois graficos ao mesmo tempo; pointplot, para
plotar os pontos; plot, para plotar um grafico de uma f(x)
simplify(expr);simplifica uma expressão
solve(eqn, var);resolve simbolicamente equaçoes eqn para varaiivel
var.
subs(expr velha=expr nova, expr);substitui uma expressão antiga por uma
nova expressão
unapply(P,x);converte um polinomio P em uma função na varivel x.
103
# PUCRS - Instituto de Matemática
# Cálculo Numérico -Prof. Eliete
Ajuste de Funçoes
# Exemplo1
Considerando:
i
0
1
2
3
4
5
6
x 1
2
3
4
5
6
7
y 34 45 63 88 120 159 205
Ajustar os dados a reta e a uma parábola. Comparar os resultados graficamente.
> with(stats);
[ anova, describe, fit, importdata, random, statevalf, statplots, transform ]
> xv:=[1,2,3,4,5,6,7];
xv := [ 1, 2, 3, 4, 5, 6, 7 ]
> yv:=[34,45,63,88,120,159,205];
yv := [ 34, 45, 63, 88, 120, 159, 205 ]
> g:=fit[leastsquare[[x,y],y=a*x+b,{a,b}]]([xv,yv]);
57
x − 12
2
> gl:=fit[leastsquare[[x,y]]]([[1,2,3,4,5,6,7],
[34,45,63,88,120,159,205]]);
g := y =
gl := y =
57
x − 12
2
> gll :=unapply(rhs(g),x);
57
x − 12
2
> gp:= fit[leastsquare[[x,y], y=a*x^2+b*x+c]]([xv, yv]):
> gpp:=unapply(rhs(gp), x);
gll := x →
7
1
gpp := x → x2 + x + 30
2
2
> with(plots):
> plots[display]({plots[pointplot]([1,34,2,45,3,63,4,88,5,120,6,15
9,7,205]),plot(gpp(x), x=-10..10, y=0..50),plot(gll(x), x=0..10,
y=0..210)});
104
Exemplo 2 - Ajuste por função Potência
Os dados abaixo dão a duração D de uma broca de Carborundum em função da velocidade de corte
V.
V [m/s] 100 120 150 180
D [min] 79 2 8 7.9 2.8
pede-se fazer uma tabela D = D(V) para V = 100(10) 180.
> v:=evalf([log(100),log(120),log(150),log(180)]);
v := [ 4.605170186, 4.787491743, 5.010635294, 5.192956851 ]
> d:=evalf([log(79),log(28),log(7.9),log(2.8)]);
d := [ 4.369447852, 3.332204510, 2.066862759, 1.029619417 ]
> fit[leastsquare[[x,y], y=a*x+b, {a,b}]]([v,d]);
> evalf(exp(30.52911140));
> dur :=x -> (.1813947103e14)*(x^(-5.680591334));
#validade:
> [dur(100),dur(120),dur(150),dur(180)];
#projeção:
> [dur(110),dur(130),dur(140),dur(160),dur(170)];
> plots[display]({plots[pointplot]([100,79,120,28,150,7.9,180,2.8]
),plot(dur(x), x=90..190, y=0..80)});
d 30 35 40 45 50 55 60 65 70 75
I 0.85 0.67 0.52 0.42 0.34 0.28 0.24 0.21 0.18 0.15
Os dados tabelados descrevem a intensidade da luz como uma função
Exemplo 3
105
da distância da fonte,I(d), medida num experimento.Utilizando o
sistema Maple, ajustar a uma parábola e a uma função do tipo:
1
Y(d)= 2
Ad + Bd + C
Analisar os resultados plotando num mesmo sistema de eixos os
pontos tabelados e as funções de ajuste determinadas.
> vd:=[30,35,40,45,50,55,60,65,70,75];
vd := [ 30, 35, 40, 45, 50, 55, 60, 65, 70, 75 ]
> vi:=[.85,.67,.52,.42,.34,.28,.24,.21,.18,.15];
vi := [ .85, .67, .52, .42, .34, .28, .24, .21, .18, .15 ]
> vi1:=evalf([1/.85,1/.67,1/.52,1/.42,1/.34,1/.28,1/.24,1/.21,1/.1
8,1/.15]);
vi1 := [ 1.176470588, 1.492537313, 1.923076923, 2.380952381, 2.941176471, 3.571428571,
4.166666667, 4.761904762, 5.555555556, 6.666666667 ]
> f2:= fit[leastsquare[[x,i], i=a*x^2+x*b+c]]([vd, vi1]);
f2 := i = .001329810498 x2 − .02080049421 x + .6161059331
> f2a:=unapply(rhs(f2),x);
f2a := x → .001329810498 x2 − .02080049421 x + .6161059331
> f3:=x->1/(.1329810498e-2*x^2-.2080049421e-1*x+.6161059331);
f3 := x →
1
.001329810498 x − .02080049421 x + .6161059331
> plots[display]({plots[pointplot]([30,.85,35,.67,40,.52,45,.42,50
,.34,55,.28,60,.24,65,.21,70,.18,75,.15]),plot(f2a(x),
x=10..100, i=0..3),plot(f3(x), x=10..100, i=0..3)});
2
106
> plots[display]({plots[pointplot]([30,.85,35,.67,40,.52,45,.42,50
,.34,55,.28,60,.24,65,.21,70,.18,75,.15]),plot(f2a(x),
x=-100..100, i=0..1),plot(f3(x), x=-100..100, i=0..2)});
>
Trabalho
x 0.00 0.50 1.00 1.50 2.00
y 1.00 0.50 0.30 0.20 0.20
107
Ajustar os dados tabelados a uma reta e a uma função do tipo:
1
Y(x)=
ax + b
Analisar os resultados plotando num mesmo sistema de eixos os
pontos tabelados e as funções de ajuste determinadas.
Comandos utilizados
evalf([expr1,..., exprn]);avalia numéricamente expressões, polinomios,
funções,...
fit[leastsquare[[var1,..., varn]]]; aplica o critério dos mínimos quadrados
transform[applay[f(x)]]([expr1,..., exprn]);aplica a função f(x) parra cada
expressão (expr)
transform[multiapplay[f(x)]]([list1,..., listn]);aplica os parâmetros da função
f(x) para cada elemento listado
plots[display]({plots[pointplot]([p1,p2,..,pn]),plot(f(x), x=-m..m)});plota gráficos;
display, para plotar dois gráficos ao mesmo tempo; pointplot, para
plotar os pontos; plot, para plotar um gráfico de uma f(x)
unapply(P,x);converte a expressão P em uma função na variável x
108
PUCRS-FAMAT-Cálculo Numérico-Exemplos de Integração Numérica Utilizando Sistema
Maple9
Prof. Eliete Biasotto Hauser
> restart;
> with(Student[Calculus1]):
Ex.:Utilizando tarapézios e Simpson com 10 subintervalos, estimar o valor de I e comparar
2
⌠ ( ( −x )2 )
com o valor exato. I=⎮
dx
⎮e
⌡
0
Obs: Por default, o número de subintervalos utilizados é 10. Incluindo a opção partition = n
o cálculo é obtido utilizando n subintervalos .
> ApproximateInt(exp(-x^2), x=0..2, method = trapezoid);evalf(%);
⎛ -1 ⎞
⎜ ⎟
⎛ -4 ⎞
⎜ ⎟
⎛ -9 ⎞
⎜ ⎟
⎛ -16 ⎞
⎜
⎟
⎛ -36 ⎞
⎜
⎟
⎛ -49 ⎞
⎜
⎟
⎛ -64 ⎞
⎜
⎟
⎛ -81 ⎞
⎜
⎟
1 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ( -1 ) 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠
+ e
+ e
+ e
+ e
+ e + e
+ e
+ e
+ e
10 5
5
5
5
5
5
5
5
5
+
1 ( -4 )
e
10
0.8818388107
> ApproximateInt(exp(-x^2), x=0..2, method = simpson);evalf(%);
⎛ -9 ⎞
⎜ ⎟
⎛ -9 ⎞
⎜
⎟
⎛ -16 ⎞
⎜
⎟
⎛ -49 ⎞
⎜
⎟
⎛ -81 ⎞
⎜
⎟
⎛ -81 ⎞
⎜
⎟
⎛ -4 ⎞
⎜ ⎟
⎛ -361 ⎞
⎜
⎟
1
1 ⎜⎝ 25 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠
e
e
e
e
e
e
e
e
+
+
+
+
+
+
+
+
30 15
15
15
15
15
15
15
15
⎛ -121 ⎞
⎜
⎟
⎛ -36 ⎞
⎜
⎟
⎛ -169 ⎞
⎜
⎟
⎛ -49 ⎞
⎜
⎟
2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 1 ( -4 ) 2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠ 2 ( -9 / 4 ) 1 ( -1 ) 2 ( -1 / 4 )
+
+
+
+
+
+
e
e
e +
e
e
e
e +
e
15
15
30
15
15
15
15
15
⎛ -64 ⎞
⎜
⎟
⎛ -1 ⎞
⎜
⎟
⎛ -289 ⎞
⎜
⎟
⎛ -1 ⎞
⎜ ⎟
1 ⎜⎝ 25 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠ 2 ⎜⎝ 100 ⎟⎠ 1 ⎜⎝ 25 ⎟⎠
e
e
e
e
+
+
+
+
15
15
15
15
0.8820809834
> Valor_exato:=int(exp(-x^2), x=0..2)=evalf(int(exp(-x^2),
x=0..2));
Valor_exato :=
1
erf( 2 ) π = 0.8820813910
2
109
> ApproximateInt(exp(-x^2), x=0..2, method = trapezoid, output =
plot,partition = 4);
> ApproximateInt(exp(-x^2), x=0..2, method = trapezoid, output =
plot, partition = 20):
> ApproximateInt(exp(-x^2), x=0..2, method = trapezoid, output =
animation):
110
> ApproximateInt(exp(-x^2), x=0..2, method = simpson, output =
plot, partition =2);
> ApproximateInt(ln(x), 1..100, method = simpson, output =
animation):
111
Atenção: Sempre analisar os resultados obtidos
> ApproximateInt(x*(x - 2)*(x - 3), x=0..5, method = simpson,
output = plot);
No último exemplo , área 22.916667 esta correto?
> Área:=evalf(int(x*(x - 2)*(x - 3),x=0..2))+abs(evalf(int(x*(x 2)*(x - 3),x=2..3)))+evalf(int(x*(x - 2)*(x - 3),x=3..5));
Área := 23.75000000
> evalf(int(x*(x - 2)*(x - 3),x=2..3));
-0.4166666667
112
No próximo exemplo é óbvio que a área não é nula.
> ApproximateInt(tan(x) - 2*x, x=-1..1, method = simpson, output =
plot, partition = 50);
> evalf(int(tan(x) - 2*x, x=-1..0));
0.3843735297
> evalf(int(tan(x) - 2*x, x=0..1));
-0.3843735297
> area:=evalf(int(tan(x) - 2*x, x=-1..0))+abs(evalf(int(tan(x) 2*x, x=0..1)));
area := 0.7687470594
113
PUCRS - Faculdade de Matemática
Prof. Eliete Biasotto Hauser-Cálculo Numérico ALaboratório EDO_EDP
⎧⎪ y` = − xy − 1
1) Resolver no intervalo [1 , 1.75], utilizando o método de Euler com h=0.25 , o PVI ⎨
⎪⎩ y( 1 ) = 2
Resposta: y(1.75) ≅ 1.2018
2) Com h=k=10, resolver numericamente equação de Laplace
∂ 2u
∂ 2u
( x ,t ) +
( x ,t ) = 0 , 0 < x < 80 , 0 < y < 60 ,
∂y 2
∂x 2
sujeita às condições u ( x ,0 ) = u ( x ,60 ) = u ( 80 , y ) = 0 , e
j↓
6
5
4
3
2
1
0
x→
0
100
100
100
100
100
100
100
0
1
0
46,993
63,138
67,192
63,138
46,993
0
10
2
0
24,835
38,368
42,490
38,368
24,835
0
20
3
0
13,978
23,010
26,032
23,010
13,978
0
30
4
0
8,068
13,660
15,620
13,660
8,068
0
40
5
0
4,633
7,942
9,128
7,942
4,633
0
50
u ( 0 , y ) = 100
6
0
2,523
4,348
5,009
4,348
2,523
0
60
7
0
1,110
1,917
2,211
1,917
1,110
0
70
8
0
0
0
0
0
0
0
80
←i
60
50
40
30
20
10
0
y↓
∂u
∂ 2u
( x ,t ) −
( x ,t ) = 0 , 0 < x < 1, 0 ≤ t , com as condições de contorno
∂t
∂x 2
u ( 0 , t ) = u ( 1 , t ) = 0 , 0 < t , e as condições iniciais u( x ,0 ) = sen( πx ), 0 ≤ x ≤ 1.
3 - Considerar a equação de calor
Para h=0,2 e k=0,01 obtemos as diferenças finitas para a equação do calor ui , j + 1 = 0 ,5 ui , j + 0 ,25⎛⎜ u
⎝
i + 1, j
Utilizando j=0 e i=1,2,3 e 4 obter as aproximações da temperatura na primeira linha do tempo u(x , 0.01).
←i
j↓
0
1
2
3
4
5
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
↑ tempo
0,2
0,4
0,6
0,8
1
114
+u
⎞
⎟
i −1, j ⎠