1.1.

SELEÇÃO PARA CONSULTA NO LABORATÓRIO 1

Tabela 1.1: Principais funções nativas no Scilab (use help nome para saber mais) Tabela 1.2: Identidades auxiliares para remover subtração catastrófica (LAB)
abs acos acosh asin asinh atan
atan2 ceil conj contour cos cosh u2 − v 2 u3 − v 3
u−v = u−v =
cotg exp expm feval floor imag u+v u2 + uv + v 2
interp intg log log10 log2 max
maxi mesh min mini modulo plot
plot2d poly power round sec sign un − v n
u−v = cos(u)2 + sin(u)2 = 1
sin sinh sqrt splin tan tanh un + un−1 v + . . . + u1 v n−1 + v n
n = 2, 3, 4, . . .

INSTITUTO DE MATEMÁTICA - UFRGS sec(u)2 = 1 + tan(u)2 cosec(u)2 = 1 + cotan(u)2

MAT01169 - CÁLCULO NUMÉRICO
MATERIAL PARA CONSULTA NAS PROVAS
1 1 v−u √
POR FAVOR NÃO ESCREVA NESTE FORMULÁRIO. − = up/q = q
up
u v uv

1.1 Seleção para consulta no laboratório log(uv) = log(u) + log(v) , u > 0, v > 0 log(uv ) = v log(u), u > 0

1.1.1 Introdução a Computação Científica

tan(u) − tan(v)
Definição: Uma medida prática e contínua da exatidão. tan(u − v) = eu−v = eu /ev
1 + tan(u) tan(v)
O digse de uma aproximação x∗ de um número x é dado por
 ∗  log(u) − log(v) = log(u/v) , u > 0, v > 0 ln(u−1 ) = − ln(u), u > 0
2|x − x|
digse(x∗ , x) = − log10 . (1.1)
|x|

1.1.2 Solução Numérica de Equações Não-lineares • M é o valor absoluto do menor dos coeficientes negativos;
Cota Básica para localização de raízes polinomiais: • ap é o último coeficiente positivo antes do primeiro coeficiente negativo.
M Observação: Considerando g(x) = −p(−x), se o grau de p é impar, ou g(x) =
|z| ≤ 1 + (1.2)
|an | p(−x), caso contrário, encontramos estimativa para as raízes negativas de p(x), comple-
tando a localização.
onde M = max{|an−1 |, |an−2 |, . . . , |a1 |, |a0 |}.
Cota de Vene para localização de raízes polinomiais1 : Algoritmo: Secante
Resultado: Toda raiz positiva α de p(x) = 0 verifica Entrada: função f (x); aproximações x0 e x−1 para a raiz procurada; parâmetro de
tolerância TOL.
M
0<α≤1+ , onde (1.3)
an + an−1 + . . . + ap Saída: aproximação xn para a raiz
1 Convenção: polinômios na forma p(x) = an xn + an−1 xn−1 + . . . a1 x + a0 , onde an > 0. Passo 1: Inicialização n = 0; segue=1,
2

Passo 2: Faça enquanto segue = 1 1.1.3 Solução Numérica de Sistemas de Equações Algébricas
f (xn ) − f (xn−1 )
Passo 3: dn =
xn − xn−1 Algoritmo: Substituição Direta.
f (xn ) Entrada: Matriz triangular inferior não-singular A = (aij ); vetor b = (bj );
Passo 4: xn+1 = xn − Saída: Vetor solução x = (xi ) de Ax = b;
dn
Passo 4: n ← n + 1 Passo 1: Para i = 1, 2, 3, . . . , n faça
Passo 5: Se TOL já foi alcançado, segue ← 0. i−1
X
Fim-faça bi − aij xj
Fim xi =
j=1
aii
Algoritmo: Newton-Raphson Modificado Retorne x = (xi );
Entrada: função f (x) com derivada f ′ (x); aproximação x0 para a raiz procurada; mul- Algoritmo: Substituição Reversa.
tiplicidade µ da raiz procurada; parâmetro de tolerância TOL. Entrada: Matriz triangular superior não-singular A = (aij ); vetor b = (bj );
Saída: Vetor solução x = (xi ) de Ax = b;
Saída: aproximação xn para a raiz
Passo 1: Para i = n, n − 1, n − 2, . . . , 2, 1 faça
Passo 1: Inicialização n = 0; segue=1, n
X
bi − aij xj
Passo 2: Faça enquanto segue = 1 j=i+1
f (xn ) xi =
Passo 3: xn+1 = xn − µ ′ aii
f (xn )
Passo 4: n ← n + 1 Retorne x = (xi );
Passo 5: Se TOL já foi alcançado, segue ← 0. Algoritmo: Eliminação Gaussiana com pivotamento parcial
Fim-faça Entrada: Matriz não-singular A = (aij ); vetor b = (bj );
Fim Saída: Vetor solução x = (xi ) de Ax = b;

Cotas para raízes polinomiais p(z) = 0 Passo 1: Para j = 1, 2, 3, . . . , n − 1 faça

Cota de Fujiwara: Sendo α qualquer raiz de p(z) = 0 então Passo 2: Encontre s, j ≤ s ≤ n tal que |a(s, j)| = max |a(i, j)|;
j≤i≤n
Passo 3: Troque linhas s e j da matriz A de posição, uma pela outra.
Passo 4: Trocar bj e bs de posição, um pelo outro.
( )
1/2 1/3 1/(n−1) 1/n
an−1 an−2 an−3 a1 a0
|α| ≤ 2 max , , ,..., , Passo 5: Para i = j + 1, . . . , n faça
an an an an an
µ = aij /ajj , bi ← bi − µbj
Passo 6: Para k = j, . . . , n faça aik ← aik − µajk
Cota de Kojima: Sendo α qualquer raiz de p(z) = 0 então |α| ≤ q1 + q2 , onde q1 e q2 Fim Para
são os dois maiores valores da lista Fim Para
Passo 7: Encontre a solução x do sistema triangular Ax = b via substituição reversa.
Retorne x.
( )
1/2 1/3 1/(n−1) 1/n
an−1 an−2 an−3 a1 a0
S= , , ,..., ,
an an an an an Algoritmo: Método iterativo de Jacobi para Sistemas Lineares
Entrada: matrizes D e N , onde A = D + N , D é diagonal; b vetor do lado direito;
aproximação x0 para a solução de Ax = b; parâmetro de tolerância TOL.
1.1. SELEÇÃO PARA CONSULTA NO LABORATÓRIO 3

Saída: aproximação xk para a solução de Ax = b. 1.1.4 Interpolação Polinomial e Ajuste de Dados

Passo 1: Inicialização: k = 0; segue = 1; Interpolação Não-Segmentada de Lagrange por {(x0 , y0 ), . . . , (xn , yn )} :
Passo 2: Enquanto segue = 1 faça
Passo 3: Resolva o sistema linear diagonal Dxk+1 = b − N xk n
X Y x − xj
Passo 5: k ← k + 1; φ(x) = yi Li (x), onde Li (x) = .
xi − xj
Passo 5: Se TOL já foi alcançado, segue ← 0; i=0 0≤j≤n,j6=i

Fim - Enquanto
Retorne xk Interpolação Não-Segmentada de Newton por (xi , f (xi )), i = 0, 1, . . . , n :
Algoritmo: Método iterativo de Gauss-Seidel para Sistemas Lineares
Entrada: matrizes L e N , onde A = L + N , L é triangular inferior; b vetor do lado φ(x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) +
direito; aproximação x0 para a solução de Ax = b; parâmetro de tolerância TOL. f [x0 , x1 , x2 , x3 ](x − x0 )(x − x1 )(x − x2 )+
Saída: aproximação xk para a solução de Ax = b. + . . . + f [x0 , . . . , xn ](x − x0 ) . . . (x − xn−1 ).

Passo 1: Inicialização: k = 0; segue = 1; onde as quantidades f [. . .] são as diferenças divididas da tabela.

Passo 2: Enquanto segue = 1 faça
Passo 3: Resolva o sistema linear triangular inferior Lxk+1 = b − N xk , Interpolação Segmentada Linear: Em cada intervalo [xi−1 , xi ], i = 1, 2, . . . , n, a
Passo 5: k ← k + 1; interpoladora φ(x) é um polinômio de primeiro grau φi (x), definido por
Passo 5: Se TOL já foi alcançado, segue ← 0;
(x − xi )yi−1 (x − xi−1 )yi
Fim - Enquanto φi (x) = + . (1.5)
xi−1 − xi xi − xi−1
Retorne xk
Algoritmo: Método de Newton para Sistemas Não-Lineares
Interpolação Cúbica Segmentada de Hermite Em cada intervalo [xi−1 , xi ], i =
Entrada: função F (x); função matriz Jacobiana JF (x); aproximação x0 para a solu-
1, 2, . . . , n, a interpoladora φ(x) é um polinômio de grau menor ou igual a 3,
ção; parâmetro de tolerância TOL.
→
Saída: aproximação xk para a solução de F (x) = 0 . φi (x) = f (xi−1 )si (x) + f (xi )ri (x) + di−1 ŝi (x) + di r̂i (x)
Passo 1: Inicialização: k = 0; segue = 1; onde hi = xi − xi−1 e
Passo 2: Enquanto segue = 1 faça
Passo 3: Resolva o sistema linear JF (xk )z k = F (xk )
 
2 2 hi 1
Passo 4: xk+1 ← xk − z k ; si (x) = 3 (x − xi ) x − xi−1 + , ŝi (x) = (x − xi−1 )(x − xi )2 (1.6)
hi 2 h2i
Passo 5: k ← k + 1;  
−2 hi 1
Passo 5: Se TOL já foi alcançado, segue ← 0; 2
ri (x) = 3 (x − xi−1 ) x − xi − , r̂i (x) = (x − xi−1 )2 (x − xi ) (1.7)
Fim - Enquanto hi 2 h2i
Retorne xk e onde os parâmetros di , i = 0, 1, 2, . . . , n devem ser prescritos.
Medida da exatidão vetorial: dado um vetor x e uma aproximação x∗ (vetor com a
Interpolação Segmentada usando Splines Cúbicos
mesma dimensão de x), definimos como digse vetorial a seguinte medida contínua da
A determinação dos di s implica na solução de um sistema linear de n − 1 equações
exatidão:
2|x∗ − x|2 φ′′i (xi ) = φ′′i+1 (xi ), i = 1, 2, . . . , n − 1 nas n + 1 incógnitas d0 , d1 , . . . , dn . Portanto 2
digsev(x∗ , x) = − log10 (1.4) equações devem ser ADICIONADAS para definir cada tipo de Spline.
|x|2
onde |u|2 denota a norma Euclidiana de um vetor, que em Scilab pode ser calculada com • Splines Cúbicos Naturais: satisfazem φ′′ (x0 ) = 0 e φ′′ (xn ) = 0. O vetor das
o comando norm(u) ou norm(u,2). incógnitas d0 , . . . , dn é calculado em Scilab via
4

vd = splin([x0 , x1 , . . . , xn ], [y0 , y1 , . . . , yn ],′ natural′ );

Tabela 1.3: Parâmetros da Quadratura de Gauss-Legendre.
n xk p wk
• Splines Cúbicos Armados: satisfazem φ′ (x0 ) = α e φ′ (xn ) = β O vetor das − 1/3 1
incógnitas d0 , . . . , dn é calculado em Scilab via 1 p
1/3 1
p
− 3/5 5/9
vd = splin([x0 , x1 , . . . , xn ], [y0 , y1 , . . . , yn ],′ clamped′ , [α, β]);
2 0p 8/9
3/5 5/9
• Splines Cúbicos Periódicos: satisfazem φ′ (x0 ) = φ′ (xn ) e φ′′ (x0 ) = φ′′ (xn ) e q p
requerem y0 = yn . O vetor das incógnitas d0 , . . . , dn é calculado em Scilab via − 3/7 + (6/7) 2/15 p
q p 1/2 − (1/18)p15/2
− 3/7 − (6/7) 2/15 1/2 + (1/18)p15/2
vd = splin([x0 , x1 , . . . , xn ], [y0 , y1 , . . . , yn ],′ periodic′ ); 3 q p
3/7 − (6/7) 2/15 1/2 + (1/18)p15/2
q 1/2 − (1/18) 15/2
Melhor solução de Ax = b no sentido dos Mínimos Quadrados, via fatoração QR
p
3/7 + (6/7) 2/15
Sendo A = QR a decomposição QR magra2 (econômica) da matriz de coeficientes A
(A deve ter posto máximo), então
onde x0 , x1 , x2 , . . . , xn são definidos conforme tabela 1.3.
Ax = b ⇔ QRx = b ⇒ QT QRx = QT b ⇔ Rx = QT b (1.8)
Quadratura de Gauss-TChebyshev para n + 1 pontos:
permite o cálculo computacional da solução x , no sentido dos Mínimos Quadrados, de
∗ Z 1
f (x) π X
n
Ax = b. √ dx ≈ f (xk ) (1.12)
−1 1 − x2 n+1
k=0

1.1.5 Integração Numérica onde x0 , x1 , x2 , . . . , xn são definidos por

Regra do Trapézio:
 
π(2k + 1)
xk = cos . (1.13)
b 3
2(n + 1)
f (a) + f (b) (b − a) (2)
Z
f (x)dx ≈ (b − a) + |f (ξ)|, (1.9)
a 2 12 Quadratura de Gauss-Laguerre para n + 1 pontos:
onde h = b − a e para algum ξ ∈ (a, b). Z ∞ Xn
f (x)e−x dx ≈ wk f (xk ) (1.14)
Regra de Simpson 0 k=0

b
f (a) + 4f (c) + f (b) (b − a)5 (4) onde x0 , x1 , x2 , . . . , xn são definidos pela tabela 1.4.
Z
f (x)dx ≈ (b − a) + |f (ξ)| (1.10)
a 6 2880 Quadratura de Gauss-Hermite para n + 1 pontos:
onde h = b − a e para algum ξ ∈ (a, b). Z ∞
2 Xn
f (x)e−x dx ≈ wk f (xk ) (1.15)
Quadratura de Gauss-Legendre para n + 1 pontos: −∞ k=0
Z 1 n
X onde x0 , x1 , x2 , . . . , xn são definidos pela tabela 1.5.
f (x)dx ≈ wk f (xk ) (1.11)
2Q tem a mesma dimensão de A, QT Q = I, enquanto R é quadrada e triangular superior.
−1 k=0
1.1. SELEÇÃO PARA CONSULTA NO LABORATÓRIO 5

Algoritmo: (Regra Composta Recursiva Trapezoidal )

Entrada: função f , números reais a e b, parâmetro de tolerância TOL.
Z b
Tabela 1.4: Parâmetros da Quadratura de Gauss-Laguerre.
n xk wk Saída: Aproximação sn para f (x)dx
a
1
0.585786437627 0.853553390593 Passo 1: n = 0; h0 = b − a; s0 = (f (a) + f (b))h0 /2; segue=1
3.414213562373 0.14644669407 Passo 2: Enquanto segue=1
0.415774556783 0.711093009929 Passo 3: n = n + 1;
2 2.294280360279 0.278517733569 hn−1
2X
n
−1
6.289945082937 0.103892565016 Passo 4: hn = , gn = f (a + i · hn )
2
0.322547689619 0.603154104342 i≥1,ímpar
1.745761101158 0.357418692438 sn−1
3 Passo 5: sn = + g n hn
4.536620296921 0.388879085150 · 10−1 2
9.395070912301 0.539294705561 · 10−3 Passo 6: Se TOL já foi alcançado, segue ← 0.
0.2635603197 0.5217556106 Fim Enquanto
1.413403059 0.39866681108 Retorne sn
4 3.59642577 0.75942244868 · 10−1 Aceleração de Richarson para Regra Composta Recursiva Trapezoidal.
u0 = s0
(
7.0858100059 0.36117586799 · 10−2
4sn − sn−1
12.64080084423 0.23369972386 · 10−4 un+1 = , n = 1, 2, 3, . . .
3

1.1.6 Solução Numérica de Equações Diferenciais

Para a equação diferencial ordinária y ′ = f (t, y), um esquema de m passos é
Tabela 1.5: Parâmetros da Quadratura de Gauss-Hermite. yi+1 + α1 yi + α2 yi−1 + . . . + αm yi+1−m
n xk wk = β0 f (ti+1 , yi+1 ) + . . . + βm f (ti+1−m , yi+1−m )
h
−0.7071067811 0.8862269255
1
0.7071067811 0.8862269255 1.1.7 Métodos Implícitos ou Corretores de Adams-Moulton
−1.22447448714 0.2954089752
2 0.00000000000 1.1816359006 São uma família de esquemas de m passos onde α1 = −1 , α2 = α3 = . . . = αm = 0
1.22447448714 0.2954089752 e os coeficientes β0 , . . . , βm são determinados de maneira a maximizar a ordem do erro
−1.6506801239 0.0813128354 de truncamento do esquema numérico. 2
−0.5246476233 0.8049140900
3
0.5246476233 0.8049140900
1.6506801239 0.0813128354 Tabela 1.6: Parâmetros para os esquemas Corretores de Adams
m β0 β1 β2 β3 T
−2.021828705 0.0199532421
−0.9585724646 0.3936193232 1 1/2 1/2 O(h2 )
2 5/12 8/12 -1/12 O(h3 )
4 0 0.9453087205
0.9585724646 0.3936193232 3 9/24 19/24 -5/24 1/24 O(h4 )
2.021828705 0.0199532421
Os esquemas explícitos de Adams-Bashforth (Preditores) são uma família de esque-
mas de m passos onde α1 = −1 , α2 = α3 = . . . = αm = 0 e os coeficientes
6

β0 , . . . , βm são determinados de maneira a maximizar a ordem do erro de truncamento k1 = f 

(ti , yi ) 
sob a restrição β0 = 0. h hk1
k2 = f ti + , y i +
 4 4 
3h 3hk1 9hk2
Tabela 1.7: Parâmetros para os esquemas Previsores de Adams k3 = f ti + , yi + +
8 32 32
m β0 β1 β2 β3 β4 T
 
12h 1932hk1 7200hk2 7296hk3
1 0 1 O(h) k4 = f ti + , yi + − +
13 2197 2197 2197
2 0 3/2 -1/2 O(h2 )
 
439hk1 3680hk3 845hk4
3 0 23/12 -16/12 5/12 O(h3 ) k5 = f ti + h, yi + − 8hk2 + −
216 513 4104
4 0 55/24 -59/24 37/24 -9/24 O(h4 )
 
h 8hk1 3544hk3 1859hk4 11hk5
k6 = f ti + , y i − + 2hk2 − + −
2 27 2565 4104 40
25hk1 1408hk3 2197hk4 hk5
Estabilidade-zero: um esquema numérico da forma yi+1 = yi + + + −
216 2565 4104 5
yn+1 + α1 yn + . . . + αm yn+1−m
= β0 f (tn+1 , yn+1 ) + . . . + βm f (tn+1−m , yn+1−m )
h k1 128k3 2197k4 k5 2k6
é zero-estável se todas as raízes λ do seu polinômio característico associado ei = − − + +
360 4275 75240 50 55
λm + α1 λm−1 + α2 λm−2 + . . . + αm−1 λ + αm = 0 Método de Euler-Cromer para y ′′ = f (t, y, v), onde v = y ′ :

yi+1 = yi + △t vi + (△t)2 f (ti , yi , vi )/2


|λ| ≤ 1 , λ raiz simples

satisfizerem (1.16)
|λ| < 1 , λ raiz múltipla. vi+1 = vi + △tf (ti , yi , vi )
Esquema de Runge-Kutta de segunda ordem, ou método de Heun:

k = f (ti , yi ) para i = 0, 1, 2, . . ..
1
Método de Verlet para y ′′ = f (t, y) :

k2 = f (ti + h, yi + hk1 )

h(k1 + k2 )
yi+1 = yi + △t vi + (△t)2 ai /2
 y (
i+1 = yi + , i = 0, 1, 2, 3, . . .

2 ai + ai+1 (1.17)
ou então, de
( outra forma (formato Previsor-Corretor), vi+1 = vi + △t
yp = yi + hf (ti , yi ) 2
h(f (ti , yi ) + f (ti+1 , yp )) para i = 0, 1, 2, . . .; onde ai = f (ti , yi ).
yi+1 = yi + , i = 0, 1, 2, 3, . . .
2 Método de Verlet Corrigido para y ′′ = f (t, y, v), onde v = y ′ :
Esquema clássico de Runge-Kutta de quarta ordem:
y = yi + △t vi + (△t)2 ai /2

k1 = f (ti , yi )

   i+1

h hk1 vp = vi + △tai

(1.18)

k2 = f ti + , y i +


ai + f (ti+1 , yi+1 , vp )

 2 2   v
i+1 = vi + △t

 

h hk2

2
k3 = f ti + , y i +

 2 2 para i = 0, 1, 2, . . .; onde ai = f (ti , yi , vi ).
k4 = f (ti + h, yi + hk3 )

Método Leapfrog para y ′′ = f (t, y):



h(k1 + 2k2 + 2k3 + k4 )


 y
i+1 = yi + , i = 0, 1, 2, 3, . . .


6 yi = yi−1 + △t vi−1/2
O método de Runge-Kutta-Fehlberg para solução numérica de y ′ = f (t, y) possui (1.19)
vi+1/2 = vi−1/2 + △tf (ti , yi )
um erro de truncamento de quarta ordem, além de fornecer uma estimativa de quinta
ordem desse erro de truncamento. FIM DO FORMULÁRIO