Métodos Iterativos para Soluções de Equações Algébricas
não-Lineares
Fernando Deeke Sasse
[email protected]
Departamento de Matemática, UDESC - Joinville
Introdução
Consideraremos aqui métodos para resolver equações algébricas não-lineares f x = 0. Em partitular,
determinaremos métodos interativos para determinar os valores x que satisfazem a equação,
possivelmente dentro de um erro pré-determinado. Denominamos este problema como o da
determinação da raiz ou zero de f x .
A existência e unicidade de soluções de equações não-lineares é frequentemente difícil de determinar,
havendo uma enorme variedade de comportamentos possíveis. Além disso, uma equação não-linear
possui um número de raízes que pode ir de zero a infinito.
Consideremos, por exemplo, a função não-linear
f x = x2 K 4 sin x
Certamente um zero desta função é x = 0. No entanto, outros possíveis zeros desta função não podem
ser determinados exatamente. Há vários modos de ao menos estimar o número de zeros. No entanto,
um dos mais simples é através de um gráfico:
O f d x2 K 4 sin x :
O plot f, x =K1 ..2.3
4
3
2
1
K1
0
K1
1
2
x
K2
Como
limx/GN f x =CN
vemos que, além de x = 0 há uma única outra raiz, que deve estar próxima a 1.8 ou 1.9. Para tentar
determinar esta raiz com maior acurácia, poderíamos fazer outro gráfico numa região mais próxima a
este valor. Por exemplo:
O plot f, x = 1.7 ..2
0.2
0
1.8
1.9
2
x
K0.2
K0.4
K0.6
K0.8
K1
Aqui já vemos que a solução deve ser próxima a 1.93.
O f d sin x Kcos x$$2
f := sin x K cos x2
(1.1)
O plot f, x = 4 ..6
0.5
0
4.5
K0.5
K1
K1.5
A equação a seguir possui raízes múltiplas (6):
5
x
5.5
6
O P6 d expand
1Kx
6
P6 := 1 K 6 x C 15 x2 K 20 x3 C 15 x4 K 6 x5 C x6
(1.2)
O plot P6 , x = 0.995 ..1.005
1.5 # 10 - 14
1. # 10 - 14
5. # 10 - 15
0
0.996
K5. # 10 - 15
0.998
1
x
1.002
1.004
Exercício: Use gráficos e argumentos heurísticos para mostrar que
(i) ex C 1 = 0 não tem solução real.
(ii) eKx K x = 0 tem uma solução real.
(iii) x2 K 9 x2 K x3 K 4 sin x =0 tem três soluções reais.
(iv) 8 x3 K 4$x2 K 3 x C 2 = 0 tem uma solução real.
(v) sin x K cos x2 = 0 tem 5 soluções no intervalo [0..4] e 6 soluções no intervalo [4..6].
Obviamente tal processo gráfico, embora útil para estimativas iniciais em técnicas que estudaremos
mais adiante, não é eficiente para a determinação da raiz com uma acurácia de 18 dígitos.
Uma função não-linear f x pode ter múltiplas raízes, onde f x = 0 e f ' x = 0. Isso signica que a
curva tem uma tangente horizontal no eixo x. Se f x = 0 e f ' x s 0, então diz-se que x tem uma raiz
simples.
Mal condicionamento
Suponhamos que xap é uma solução aproximada de uma equação não-linear f x = 0. Isso significa
que
f xap
z0
ou
xap K xt z 0 ,
onde xt é a verdadeira solução. A primeira relação corresponde ao fato de termos um pequeno resíduo,
enquanto que a segunda dá uma medida da proximidade da solução verdadeira (desconhecida). Tais
critérios não são simultaneamente pequenos. Isso é ilustrado na Fig. 1, onde ambas as funções têm a
mesma incerteza nos seus valores, mas muito diferentes valores para as respectivas localizações de suas
raízes. O gráfico em azul representa uma função associada a uma equação bem condicionada. O gráfico
em cor preta representa uma função associada a uma equação mal condicionada.
0.10
0.05
0
2.5
3
x
K0.05
3.5
4
K0.10
Fig. 1 Condicionamento
O condicionamento de um problema de determinação de raiz para uma dada função é o oposto do
problema de avaliação da função: se o valor da função é insensitivo ao valor do argumento então a raiz
será sensitiva. Por outro lado, se o valor da função é sensitivo ao argumento, então a raiz é insensitiva.
Esta propriedade faz sentido, pois os dois problemas são inversos um do outro: se y = f x , então
encontrar x, dado y é tem o condicionamento oposto ao do problema de encontrar y, dado x.
Convergência de métodos iterativos
Ao contrário do que acontece com equações lineares, a maior parte da equações não-lineares não pode
ser resolvida en número finito de passos. Portanto, usualmente devemos recorrer a um método iterativo
que produz soluções mais acuradas a cada passo, a ser terminado quando a acurácia desejada for
alcançada.
O custo computacional total da resolução do problema depende do:
1. Custo por iteração.
2. Número de iterações necessárias para a convergência.
Frequentemente há um "trade-off" entre esses dois fatores. Para comparar a efetividade dos métodos
iterativos necessitamos caracterizar suas taxas de convergência. Definimos o erro ek numa iteração por
ek = xk K x) ,
onde xk é a solução aproximada na iteração k e x) é a solução verdadeira. Métodos iterativos para
equações não-lineares não produzem diretamente uma solução aproximada xk, mas sim um intervalo
que contém a solução, sendo que o comprimento deste intervalo diminui a cada iteração. Para tais
métodos tomamos ek como sendo o comprimento deste intervalo na iteração k. Um método iterativo é
dito convergente com taxa r se
ek C 1
limk/N
=C,
ek r
onde C é uma constante finita, não nula. Alguns casos de interesse são:
• Se r = 1 e C ! 1, a taxa de convergência é linear.
• Se r O 1, a taxa de convergências é superlinear.
• Se r = 2, a taxa de convergência é quadrática.
Um modo de interpretar a distinção entre convergêncial linear e superlinear é que, assintoticamente,
uma sequência que converge linearmente ganha um número constante de dígitos acurados por iteração,
enquanto que uma sequência que converge superlinearmente ganha um número sempre maior de dígitos
de acurados a cada iteração. Em particular, um método quadraticamente convergente dobra o número de
dígitos de acurácia a cada iteração.
Bissecção
Na aritmética de precisão finita pode não haver um número de ponto flutuante x tal que f x seja
exatamente zero. Uma alternativa é buscar um por um pequeno intervalo a, b no qual f tem uma
mudança de sinal, de modo que garantidamente a correspondente função contínua deve ser zero em
algum ponto deste intervalo.
Um intervalo para o qual o sinal de f difere nos seus extremos é chamado um bracket. O método da
bissecção começa com bracket inicial e sucessivamente reduz seu tamanho até que uma solução tenha
sido isolada com a acurácia desejada. A cada iteração a função é avaliada no ponto médio do intervalo
corrente, e metade do intervalo é descartada, dependendo do sinal da função no ponto médio. Mais
formalmente, o pseudocódigo é o seguinte:
while ((b-a)>tol) do
c=a+(b-a)/2
if sign(f(a))=sign(f(c)) then
a=c
else
b=c
end
end
O fluxograma do algoritmo é dado abaixo: (Shaharuddin Salleh, Albert Y. Zomaya, Sakhinah Abu
Bakar. Computing for numerical methods using Visual c++, John Wiley & Sons, Inc., 2008):
A evolução do processo é ilustrada na figura abaixo ( ibid.):
Exemplo Tentemos encontrar a menor raiz positiva da equação
x2 K 4$sin x2 C 2 = 0
Solução
Uma análise gráfica preliminar é útil para a determinação do intervalo inicial:
O restart :
O f d x2 K 4$sin x2 C 2 :
O plot f, x = 0.7 ..1
0.6
0.5
0.4
0.3
0.2
0.1
0
K0.1
0.8
0.9
1.0
x
K0.2
K0.3
Vemos do gráfico que a menor raiz positiva está, por exemplo, no intervalo 0.7, 1 . Implementaremos
o pseudocódigo em Maple da seguinte forma:
O G d unapply f, x
G := x/x2 K 4 sin x2 C 2
(4.1)
O f(3);
x 3 2 K 4 sin x2 3 C 2
(4.2)
O G(3.);
9.351526059
(4.3)
O Digits := 11;
Digits := 11
(4.4)
O a:=0.7:
b:=1.:
i:=0:
tol:=10^(-10):
while abs((b-a))>tol do
m:=(a+b)/2.;
if sign(G(a))=sign(G(m)) then
a:=m;
else
b:=m;
fi:
i:=i+1
od:
O m, i
0.87308372119, 32
(4.5)
Portanto, com precisão de 11 dígitos decimais, 32 iterações foram necessárias.
Transformemos o loop anterior num procedimento que tem como entradas uma expressão, o intervalo e
a tolerância. A saída é a solução da equação e o número de iterações.
O bis1:=proc(f,a,b,tol)
local i,G,m,af,bf:
i := 0:
af:=evalf(a):
bf:=evalf(b):
G:=unapply(f,x):
while abs(bf-af) > tol do
m := (af+bf)/2.;
if sign(G(af)) = sign(G(m)) then
af := m
else
bf := m
end if;
i := i+1
end do:
return([m,i]);
end:
Para testar o procedimento resolvamos o mesmo problema novamente.
O Digits d 11 :
O bis1(f,0.7,1,10^(-10));
0.87308372120, 32
(4.6)
O pseudocódigo a seguir implementa o algoritmo de bissecção na forma de um procedimento, de tal
forma que o processo de convergência pode ser acompanhado. A entrada consiste na função, intervalo,
número máximo de iterações e tolerância de erro.
procedure Bisseccao(f,a,b,nmax,tol)
integer n, nmax; real a,b,c, fa, fb, fc, error
fa ) f a
fb ) f b
fa ) f a
if sign(fa) = sign(fb) then
output a,b, f a , f b
output "funcao tem os mesmos sinais em a e b"
return
end if
error ← b K a
for n = 0 to nmax do
error ) erro/2
c ) a + error
fc ) f c
output n, c, fc error
if error ! tol then
output "convergencia"
return
end if
if sign(fa) s sign(fc) then
b)c
fb ) fc
end if
end do
end procedure Bisseccao
Temos no algoritmo acima três critérios de parada. Uma implentação deste procedimento no Maple é a
seguinte:
O bis2 := proc (f, a, b, nmax, tol)
local af, bf, F,c, erro,n;
af := evalf(a); bf := evalf(b);
F := unapply(f, x);
if sign(F(af)) = sign(F(bf)) then
printf(" a = %11.9f, b = %11.9f, F(a) = %11.9f, F(b) =
%11.9f \n", a, b, F(af), F(bf));
error "Função tem os mesmos sinais em a e b. Tente
novamente"
end if ;
erro:=bf-af:
for n from 1 to nmax do
erro:=erro/2.;
c:=af+erro;
printf(" n = %2.0f, c = %13.19f, f(c) = %13.9f, erro =
%13.9f \n", n, c, F(c), erro); #dados em ponto flutuante, com 7
dígitos.
if abs(erro)<tol then
print(`Convergência foi alcançada`);
break;
end if;
if sign(F(af))<> sign(F(c)) then
bf:=c
else
af:=c
fi:
end do:
end :
Exemplo: Façamos um teste para este procedimento, para encontrar a primeira raiz positiva da equação
x4 K 4 sin x K 3 = 0.
Façamos um gráfico para ter uma idéia da localização da raiz:
O f d x4 K 4$sin x C 2
f := x4 K 4 sin x C 2
O plot f, x =K2 ..2
(4.7)
20
15
10
5
K2
K1
0
1
x
2
A primeira raiz positiva está entre 0 e 1, próxima a 0.5. Aplicando nosso procedimento temos:
O bis2(f, 0, 1, 40, 10^(-9));
n = 1, c = 0.5000000000000000000, f(c) =
0.144797846, erro
= 0.500000000
n = 2, c = 0.7500000000000000000, f(c) = -0.410148790, erro
= 0.250000000
n = 3, c = 0.6250000000000000000, f(c) = -0.187801201, erro
= 0.125000000
n = 4, c = 0.5625000000000000000, f(c) = -0.033097779, erro
= 0.062500000
n = 5, c = 0.5312500000000000000, f(c) =
0.053206013, erro
= 0.031250000
n = 6, c = 0.5468750000000000000, f(c) =
0.009362052, erro
= 0.015625000
n = 7, c = 0.5546875000000000000, f(c) = -0.012044834, erro
= 0.007812500
n = 8, c = 0.5507812500000000000, f(c) = -0.001385136, erro
= 0.003906250
n = 9, c = 0.5488281250000000000, f(c) =
0.003977584, erro
= 0.001953125
n = 10, c = 0.5498046875000000000, f(c) =
0.001293498, erro
= 0.000976562
n = 11, c = 0.5502929687500000000, f(c) = -0.000046502, erro
= 0.000488281
n = 12, c = 0.5500488281200000000, f(c) =
0.000623328, erro
= 0.000244141
n = 13, c = 0.5501708984300000000, f(c) =
0.000288370, erro
= 0.000122070
n = 14, c = 0.5502319335900000000, f(c) =
0.000120924, erro
= 0.000061035
n = 15, c = 0.5502624511700000000, f(c) =
0.000037208, erro
= 0.000030518
n = 16, c = 0.5502777099600000000, f(c) = -0.000004647, erro
= 0.000015259
n = 17, c = 0.5502700805600000000, f(c) =
0.000016280, erro
= 0.000007629
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
= 18, c = 0.5502738952600000000, f(c) =
0.000003815
= 19, c = 0.5502758026100000000, f(c) =
0.000001907
= 20, c = 0.5502767562800000000, f(c) =
0.000000954
= 21, c = 0.5502762794500000000, f(c) =
0.000000477
= 22, c = 0.5502760410300000000, f(c) =
0.000000238
= 23, c = 0.5502759218200000000, f(c) =
0.000000119
= 24, c = 0.5502759814200000000, f(c) =
0.000000060
= 25, c = 0.5502760112200000000, f(c) =
0.000000030
= 26, c = 0.5502760261200000000, f(c) =
0.000000015
= 27, c = 0.5502760186700000000, f(c) =
0.000000007
= 28, c = 0.5502760149500000000, f(c) =
0.000000004
= 29, c = 0.5502760168100000000, f(c) =
0.000000002
= 30, c = 0.5502760158800000000, f(c) =
0.000000001
Convergência foi alcançada
0.55027601588
A segunda raiz está em
O bis2 f, 1, 1.5, 40, 10^ K9 ;
n = 1, c = 1.2500000000000000000, f(c) =
= 0.250000000
n = 2, c = 1.1250000000000000000, f(c) =
= 0.125000000
n = 3, c = 1.1875000000000000000, f(c) =
= 0.062500000
n = 4, c = 1.1562500000000000000, f(c) =
= 0.031250000
n = 5, c = 1.1406250000000000000, f(c) =
= 0.015625000
n = 6, c = 1.1328125000000000000, f(c) =
= 0.007812500
n = 7, c = 1.1289062500000000000, f(c) =
= 0.003906250
n = 8, c = 1.1269531250000000000, f(c) =
= 0.001953125
n = 9, c = 1.1259765625000000000, f(c) =
= 0.000976562
n = 10, c = 1.1264648438000000000, f(c) =
= 0.000488281
n = 11, c = 1.1267089844000000000, f(c) =
= 0.000244141
n = 12, c = 1.1268310547000000000, f(c) =
= 0.000122070
n = 13, c = 1.1267700196000000000, f(c) =
= 0.000061035
n = 14, c = 1.1268005372000000000, f(c) =
= 0.000030518
n = 15, c = 1.1268157960000000000, f(c) =
= 0.000015259
0.000005816, erro
0.000000584, erro
-0.000002031, erro
-0.000000723, erro
-0.000000069, erro
0.000000258, erro
0.000000094, erro
0.000000012, erro
-0.000000028, erro
-0.000000008, erro
0.000000002, erro
-0.000000003, erro
-0.000000000, erro
(4.8)
0.645467773, erro
-0.007263736, erro
0.278792980, erro
0.126142639, erro
0.057089806, erro
0.024332539, erro
0.008390138, erro
0.000527242, erro
-0.003377223, erro
-0.001427236, erro
-0.000450559, erro
0.000038201, erro
-0.000206214, erro
-0.000084015, erro
-0.000022909, erro
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
= 16, c = 1.1268234254000000000, f(c) =
0.000007646,
0.000007629
= 17, c = 1.1268196107000000000, f(c) = -0.000007631,
0.000003815
= 18, c = 1.1268215180000000000, f(c) =
0.000000007,
0.000001907
= 19, c = 1.1268205644000000000, f(c) = -0.000003812,
0.000000954
= 20, c = 1.1268210412000000000, f(c) = -0.000001902,
0.000000477
= 21, c = 1.1268212796000000000, f(c) = -0.000000948,
0.000000238
= 22, c = 1.1268213988000000000, f(c) = -0.000000470,
0.000000119
= 23, c = 1.1268214584000000000, f(c) = -0.000000232,
0.000000060
= 24, c = 1.1268214882000000000, f(c) = -0.000000112,
0.000000030
= 25, c = 1.1268215031000000000, f(c) = -0.000000053,
0.000000015
= 26, c = 1.1268215106000000000, f(c) = -0.000000023,
0.000000007
= 27, c = 1.1268215143000000000, f(c) = -0.000000008,
0.000000004
= 28, c = 1.1268215162000000000, f(c) = -0.000000000,
0.000000002
= 29, c = 1.1268215171000000000, f(c) =
0.000000003,
0.000000001
Convergência foi alcançada
1.1268215171
Resolvendo o mesmo problema com o procedimento bis1, temos
O L1 d bis1 f, 0, 1, 10^ K9 ; L2 d bis1 f, 1, 1.5, 10^ K9
L1 := 0.55027601595, 30
L2 := 1.1268215171, 29
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
(4.9)
(4.10)
O R1 d L1 1 ; R2 d L2 1
R1 := 0.55027601595
R2 := 1.1268215171
(4.11)
O método da bissecção não faz uso das magnitudes dos valores das funções, mas somente dos seus
sinais. Como resultado, a bissecção converge com certeza, mas o faz de forma lenta. Mais
especificamente, a cada iteração o comprimento do intervalo contendo a solução (e consequentemente o
erro) é reduzido pela metade. Isso significa que o método da bissecção é linearmente convergente, com
r = 1 e C = 0.5.
Notemos ainda que dado um intervalo inicial a, b , o comprimento do intervalo após k iterações é
bKa
,
2k
Se um erro de tolerância for prescrito é possível determinar o número de passos necessários no método
da bissecção. Suponha que desejamos x) K ck ! e. Então devemos ter um número k de iterações tal
que
bKa
x) K ck %
%e
2k
O número de iterações necessário para que o erro esteja dentro dentro do intervalo de tolerância e é
bKa
k R log2
e
iterações, independente de f. No presente caso, de fato,
1.
O log 2
10^ K9
29.897352854
(4.12)
Exemplo: Determinemos o número de iterações necessárias para resolver a equação para resolver
f x = x3 C 4 x2 K 10 = 0
com uma tolerância de 10K3 , supondo a = 1 e b = 2.
O f := x^3+4*x^2-10;
f := x3 C 4 x2 K 10
(4.13)
O plot f, x =K3.5 ..2.5
30
20
10
K3
K2
0
K1
1
2
x
K10
Necessitamos então um número de iterações
k R log2
O evalf log2
bKa
= log2
e
3 ln 10
ln 2
2K1
10K3
(4.14)
2K1
10K3
9.9657842847
(4.15)
K3
Portanto, 10 iterações garantirão uma aproximação com erro que não ultrapassa 10 . De fato,
O bis1 f, 1, 2, 10K3
1.3642578125, 10
bis2 f, 1, 2, 13, 10K3
= 1, c = 1.5000000000000000000, f(c) =
2.375000000,
0.500000000
= 2, c = 1.2500000000000000000, f(c) = -1.796875000,
0.250000000
= 3, c = 1.3750000000000000000, f(c) =
0.162109375,
0.125000000
= 4, c = 1.3125000000000000000, f(c) = -0.848388672,
0.062500000
= 5, c = 1.3437500000000000000, f(c) = -0.350982666,
0.031250000
= 6, c = 1.3593750000000000000, f(c) = -0.096408844,
0.015625000
= 7, c = 1.3671875000000000000, f(c) =
0.032355786,
0.007812500
= 8, c = 1.3632812500000000000, f(c) = -0.032149970,
0.003906250
= 9, c = 1.3652343750000000000, f(c) =
0.000072025,
0.001953125
= 10, c = 1.3642578125000000000, f(c) = -0.016046691,
0.000976562
Convergência foi alcançada
1.3642578125
Vejamos como o problema converge com uma precisão maior e menor tolerância:
O Digits d 19 :
O
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
n
=
(4.16)
erro
erro
erro
erro
erro
erro
erro
erro
erro
erro
(4.17)
O bis1 f, 1, 2, 10K18
1.365230013414096845, 60
Podemos comparar este resultado com o exato:
O evalf solve f, x
1.365230013414096846, K2.682615006707048423 C 0.3582593599240429912 I,
K2.682615006707048423 K 0.3582593599240429912 I
O dígito incorreto no processo iterativo é marcado a seguir:
1.365230013414096845
1.365230013414096846
(4.18)
(4.19)
Notemos que limitante inferior para o número de iterações k no método da bissecção não significa que
a convergência possa ser alcançada muito antes.
Exemplo: Consideremos novamente o problema de buscar a raiz positiva de
x3 K 0.165 x2 C 0.000003993 = 0
O Digits d 19 :
O f d x3 K 0.165 x2 C 0.000003993;
f := x3 K 0.165 x2 C 0.000003993
O plot f, x =K0.02 ..0.02
(4.20)
K0.02
K0.01
0
0.01
x
K0.00001
0.02
K0.00002
K0.00003
K0.00004
K0.00005
K0.00006
K0.00007
O bis1 f, 0, 0.01, 10^ K19
;
0.004995553672659094264, 57
O evalf solve f = 0, x
0.004995553672659094289, 0.1648530717782957632, K0.004848625450954857531
Os dígitos incorretos estão marcados abaixo:
0.004995553672659094264
0.004995553672659094289
Notemos que uma escolha diferente de intervalo pode mudar um pouco o resultado:
O bis1 f, 0.0049, 0.1, 10^ K19 ;
0.004995553672659094312, 60
0.004995553672659094312
Equações com raízes múltiplas não podem ser resolvidas por este método. Por exemplo:
O f d 1Kx
(4.21)
(4.22)
(4.23)
3
f := 1 K x
O plot f, x = 0.6 ..1.4
3
(4.24)
0.06
0.04
0.02
0
0.7
0.8
0.9
1.0
x
1.1
1.2
1.3
1.4
K0.02
K0.04
K0.06
Vantagens do método de bissecção
(i) O método é sempre convergente, podendo ser usado como um iniciador para métodos de maior
acurácia.
(ii) Como os intervalos contendo a raiz são reduzidos à metade a cada iteração, podemos predizer o erro
na solução.
Desvantagens do método de bissecção
(i) A convergência é lenta (linear), pois a cada passo o intervalo contendo a raiz é somente reduzido à
metade.
(ii) O método não funciona para a determinar raízes múltiplas
Exercícios:
Resolva todos os problemas ou todos os problemas ímpares (Exercícios 2.1) de Burden-Faires, Análise
Numérica, 8 ed. páginas 51 e 52.
Métodos de ponto fixo
Vamos considerar aqui o problema de determinar os valores de x tal que f x = 0, onde f x é uma
função não-linear, bem comportada no intervalo a, b . Seja
g x = xCc x f x ,
onde c x s 0 em a, b . Então o problema de determinar os valores x = b, tais que f b = 0, é
equivalente ao de encontrar valores de x para os quais b = g b .
O método iterativo consiste em partir de uma solução aproximada e determinar as outras aproximações
por
xi C 1 = g xi
De acordo com as diferentes escolhas para g x temos diferentes métodos de ponto fixo.
d
g x ! 1, em [a,b]. Além
dx
disso, pode-se mostrar que neste caso há somente um x = b neste intervalo tal que g b = b.
A convergência do processo iterativo no intervalo [a,b] é garantida se
Se fizermos c x = 1 temos o chamado método da iteração linear. Neste caso, dada a função f x
devemos escolher um g x tal que
g x = xCf x .
Normalmente temos diversas escolhas possíveis para g x .
Exemplo. Determinemos as raízes reais de
f x = x3 K 5 x C 3.
Façamos um gráfico para obter a localização aproximada das raízes:
O restart;
O f := x^3-5*x+2;
f := x3 K 5 x C 2
O plot(f, x = -3 .. 3);
(5.1)
10
5
K3
K2
K1
0
1
x
2
3
K5
K10
Algumas possíveis escolhas de g x para formar as relações de iteração
xi C 1 = g xi
são obtidas isolando x de f x = 0 de diferentes formas:
(i)
x3 C 2
g x =
,
5
(ii)
g x = 5 xK2
(iii)
g x =
1
3
,
5 xK2
x2
Tomemos a escolha (i) de g x para tentar encontrar a raiz negativa que está no intervalo K3,K2 .
Analisemos a condição suficiente para convergência neste intevalo, g' x ! 1:
O g := (x^3+2)*(1/5);
1 3
2
g :=
x C
(5.2)
5
5
O gp := abs(diff(f, x));
gp := 3 x2 K 5
(5.3)
O plot(gp, x = -3 .. -2);
22
20
18
16
14
12
10
8
K3
K2.8
K2.6
K2.4
x
K2.2
K2
Este gráfico mostra que a convergência não é garantida. De fato:
O gf := unapply(g,x);
gf := x/
1 3
2
x C
5
5
(5.4)
O r[1] := -3;
for i to 7 do
r[i+1] := evalf(gf(r[i]));
end do;
r1 := K3
r2 := K5.
r3 := K24.60000000
r4 := K2976.987200
r5 := K5.276681702 109
r6 := K2.938411998 1028
r7 := K5.074205616 1084
r8 := K2.612968538 10253
O
Esta divergência pode ser ilustrada graficamente:
O path d NULL :
(5.5)
O
for i from 1 to 7 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x =K10 ..2, y =K10 ..6, color = red, green,
blue ;
6
4
2
K10
K8
K6
K4
0
K2
x
K2
K4
y
K6
K8
K10
No intervalo 0.5, 1 , que contém a primeira raiz positiva, temos
O plot gp, x = 0 ..0.5
2
5.0
4.9
4.8
4.7
4.6
4.5
4.4
4.3
0
0.1
0.2
x
0.3
0.4
0.5
de modo que a convergência não é garantida.
No intervalo 1.5, 2.5 , que contém a segunda raiz positiva, temos
O plot gp, x = 1.5 ..2.5
12
10
8
6
4
2
1.6
1.8
2.0
x
2.2
2.4
e a convergência também não é garantida. Notemos, no entanto, o gráfico da iterações sugere que talvez
a convergência para a primeira raiz positiva seja possível:
O r[1] := -1; for i to 12 do
r[i+1] := evalf(gf(r[i]))
end do;
r1 := K1
r2 := 0.2000000000
r3 := 0.4016000000
r4 := 0.4129542152
r5 := 0.4140843142
r6 := 0.4142002612
r7 := 0.4142121931
r8 := 0.4142134214
r9 := 0.4142135479
r10 := 0.4142135609
r11 := 0.4142135622
r12 := 0.4142135624
r13 := 0.4142135624
Ilustremos graficamente esta convergência:
O path := NULL;
for i to 7 do
path := path, [r[i], r[i+1]], [r[i+1], r[i+1]]
end do;
plot([g(x), x, [path]], x = -1.3 .. 2.2, y = -.5 .. 2.2, color
= [red, green, blue]);
path :=
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000 , 0.4016000000, 0.4129542152 ,
0.4129542152, 0.4129542152
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000 , 0.4016000000, 0.4129542152 ,
0.4129542152, 0.4129542152 , 0.4129542152, 0.4140843142 , 0.4140843142,
0.4140843142
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000 , 0.4016000000, 0.4129542152 ,
0.4129542152, 0.4129542152 , 0.4129542152, 0.4140843142 , 0.4140843142,
(5.6)
0.4140843142 , 0.4140843142, 0.4142002612 , 0.4142002612, 0.4142002612
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000 , 0.4016000000, 0.4129542152 ,
0.4129542152, 0.4129542152 , 0.4129542152, 0.4140843142 , 0.4140843142,
0.4140843142 , 0.4140843142, 0.4142002612 , 0.4142002612, 0.4142002612 ,
0.4142002612, 0.4142121931 , 0.4142121931, 0.4142121931
path := K1, 0.2000000000 , 0.2000000000, 0.2000000000 , 0.2000000000,
0.4016000000 , 0.4016000000, 0.4016000000 , 0.4016000000, 0.4129542152 ,
0.4129542152, 0.4129542152 , 0.4129542152, 0.4140843142 , 0.4140843142,
0.4140843142 , 0.4140843142, 0.4142002612 , 0.4142002612, 0.4142002612 ,
0.4142002612, 0.4142121931 , 0.4142121931, 0.4142121931 , 0.4142121931,
0.4142134214 , 0.4142134214, 0.4142134214
2
1.5
y
1
0.5
K1
0
1
2
x
K0.5
A convergência para esta mesma raiz também ocorre se partirmos de de x = 1, por exemplo:
O r 1 d 1;
for i to 12 do
r i C 1 d evalf gf r i ;
end do;
r1 := 1
r2 := 0.6000000000
r3 := 0.4432000000
r4 := 0.4174112219
r5 := 0.4145452891
r6 := 0.4142477389
r7 := 0.4142170809
r8 := 0.4142139246
r9 := 0.4142135997
r10 := 0.4142135662
r11 := 0.4142135628
r12 := 0.4142135624
r13 := 0.4142135624
(5.7)
O path d NULL :
for i from 1 to 7 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x =K1.3 ..2.2, y =K0.2 ..2.2, color = red, green, blue ;
2
1.5
y
1
0.5
K1
0
1
x
Este mesmo gráfico, no entanto, sugere divergência se escolhermos x = 2.5:
O r 1 d 2.5;
for i to 6 do
r i C 1 d evalf gf r i ;
end do;
r1 := 2.5
r2 := 3.525000000
2
r3 := 9.160065624
r4 := 154.1183630
r5 := 7.321387530 105
r6 := 7.848925036 1016
r7 := 9.670758526 1049
(5.8)
O path d NULL :
for i from 1 to 6 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x = 1.5 ..10, y =K0.2 ..12, color = red, green, blue ;
12
10
8
y
6
4
2
0
2
3
4
5
6
x
Utilizemos neste intervalo, a escolha (iii) para g x :
g x =
5 xK2
x2
7
8
9
10
O gd
5 xK2
x2
g :=
5 xK2
(5.9)
x2
O gf d unapply g, x :
O gp d abs diff g, x
gp := K
5
C
2 5 xK2
x2
Analisemos a convergência no intervalo 1.9, 4
1.9, 4
:
O plot gp, x = 1.9 ..4
(5.10)
x3
(5.11)
0.8
0.7
0.6
0.5
0.4
0.3
2
2.5
x
3
3.5
Portanto, a convergência é assegurada se partirmos de x = 3, por exemplo. De fato,
O r 1 d 3;
for i to 30 do
r i C 1 d evalf gf r i
end do;
;
r1 := 3
r2 := 1.444444444
r3 := 2.502958580
r4 := 1.678391988
r5 := 2.269066631
4
r6 := 1.815098886
r7 := 2.147613933
r8 := 1.894536937
r9 := 2.081950997
r10 := 1.940181397
r11 := 2.045771884
r12 := 1.966188863
r13 := 2.025646566
r14 := 1.980928458
r15 := 2.014395021
r16 := 1.989255351
r17 := 2.008087426
r18 := 1.993950749
r19 := 2.004546101
r20 := 1.996595584
r21 := 2.002556212
r22 := 1.998084474
r23 := 2.001437562
r24 := 1.998922345
r25 := 2.000808532
r26 := 1.999393764
r27 := 2.000454769
r28 := 1.999658974
r29 := 2.000255799
r30 := 1.999808169
r31 := 2.000143882
Notamos que a convergência é lenta. Uma explicação para este fato pode ser obtido do gráfico de
interações:
O path d NULL :
for i from 1 to 12 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x = 1.5 ..3.5, y =K0.2 ..3, color = red, green, blue ;
(5.12)
3
2
y
1
0
2
2.5
x
3
3.5
Vemos aqui a convergência oscilante.
Exemplo. Determinar as raízes reais de
f x = 5 sin x K ex
pelo método da iteração linear .
Façamos um gráfico para ter noção da localização das raízes:
O restart :
O f := 5*sin(x)-exp(x);
f := 5 sin x K ex
O plot(f, x = 0 .. 2);
(5.13)
1
0
0.5
1
x
1.5
2
K1
K2
Busquemos inicialmente a raiz contida no intervalo 0, 0.5 . Isolando x da função exponencial temos
ex = 5 sin x ,
x = ln 5 sin x .
O processo iterativo é, então,
xi C 1 = g xi = ln 5 sin xi
O g d ln 5$sin x
g := ln 5 sin x
(5.14)
O gf d unapply g, x
gf := x/ln 5 sin x
Analisemos a convergência no intervalo 0.1, 0.5 :
O gp d abs diff g, x
cos x
gp :=
sin x
O plot gp, x = 0.1 ..0.5
(5.15)
(5.16)
9
8
7
6
5
4
3
2
0.1
0.2
0.3
x
0.4
0.5
Ou seja, a convergência não é garantida neste caso. De fato,
O r 1 d 0.1;
for i to 7 do
r i C 1 d evalf gf r i ;
end do;
r1 := 0.1
r2 := K0.6948144032
r3 := 1.163530238 C 3.141592654 I
r4 := 4.059164941 C 0.4059094416 I
r5 := 1.500765599 K 2.855062860 I
r6 := 3.774628152 K 0.06956983537 I
r7 := 1.091382760 C 3.047215575 I
r8 := 3.964802847 C 0.4775698774 I
No intervalo que contém a segunda raiz positiva, 1.5, 1, 8 , por exemplo, temos
O plot gp, x = 1.5 ..1.8
(5.17)
0.20
0.15
0.10
0.05
0
1.6
1.7
1.8
x
A convergência é garantida. De fato,
O r 1 d 1.8;
for i to 8 do
r i C 1 d evalf gf r i ;
end do;
r1 := 1.8
r2 := 1.582937488
r3 := 1.609364207
r4 := 1.608693987
r5 := 1.608719624
r6 := 1.608718652
r7 := 1.608718689
r8 := 1.608718687
r9 := 1.608718687
Num gráfico de iterações,
O path d NULL :
for i from 1 to 12 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x = 1.4 ..2, y = 1.5 ..1.7, color = red, green, blue ;
(5.18)
1.70
1.65
y 1.60
1.55
1.50
1.4
1.5
1.6
1.7
x
1.8
1.9
2.0
Exemplo. Determinar os zeros de
O g d arctan
exp x
5
g := arctan
1 x
e
5
(5.19)
Transformemos a expressão numa função:
O gf d unapply g, x
gf := x/arctan
1 x
e
5
Analisemos a convergência no intervalo 0.1, 0.5 :
O gp d abs diff g, x
1
eR x
gp :=
5
1
1C
ex
25
O plot gp, x = 0.1 ..0.5
(5.20)
2
(5.21)
0.29
0.28
0.27
0.26
0.25
0.24
0.23
0.22
0.1
0.2
0.3
x
0.4
A convergência é portanto garantida no intervalo 0.1, 0.5 , pois aí gp ! 1. De fato,
O r 1 d 0.11;
for i to 16 do
r i C 1 d evalf gf r i
end do;
;
r1 := 0.11
r2 := 0.2196534917
r3 := 0.2441587360
r4 := 0.2499694594
r5 := 0.2513657636
r6 := 0.2517023538
r7 := 0.2517835532
r8 := 0.2518031454
r9 := 0.2518078729
r10 := 0.2518090136
r11 := 0.2518092888
r12 := 0.2518093552
r13 := 0.2518093713
0.5
r14 := 0.2518093753
r15 := 0.2518093762
r16 := 0.2518093764
r17 := 0.2518093764
(5.22)
O path d NULL :
for i from 1 to 10 do
path d path, r i , r i C 1 , r i C 1 , r i C 1 ;
end do:
plot g x , x, path , x = 0.1 ..0.3, y = 0.1 ..0.3, color = red, green, blue ;
0.30
0.25
y 0.20
0.15
0.10
0.10
0.15
0.20
x
0.25
Exercícios
1. Utilizando o método da iteração linear, encontre todas as raízes reais de
(i) f x = x5 K x2 K 5 x K1,
(ii) f x = x2 C tan x K x C 1
Faça a análise de convergência e ilustre graficamente a convergência.
Método de Newton-Raphson
0.30
Os métodos de ponto fixo consistem em definir um g x tal que
g x = xCc x f x ,
tal que x = g x é equivalente a f x = 0. O processo iterativo é então
xi C 1 = g xi .
Veremos agora que o método de Newton é um caso particular de método de iteração fixa, onde c x
tem uma forma que otimiza o processo de convergência.
Suponhamos que x) é uma raiz de f x . Como
g' x = 1 C c' x f x C c x f ' x ,
temos
g' x) = 1 C c x) f ' x) .
A convergência mais rápida possível ocorre quando g' x) = 0, ou seja, quando escolhemos
1
c x) =K
.
f ' x)
Como x) não é conhecido vamos definir
1
c x =K
,
f' x
de modo que
f x
g x =xK
.
f' x
Error, unable to match delimiters
f x
g x =xK
.
f' x
O processo iterativo é então
xi C 1 = xi K
f xi
f ' xi
Este é o chamado método da iteração de Newton ou Newton-Raphson.
Exemplo. Busquemos as raízes reais de
f x = x K sin x C 2
O restart;
O f := x-> 4*x^5-3*x^3-3*x^2-5*x+3 ;
f := x/4 x5 K 3 x3 K 3 x2 K 5 x C 3
O plot(f(x), x = -1.5 .. 1.5);
(6.1)
5
K1.5
K1
K0.5
0
0.5
1
1.5
x
K5
K10
K15
O Df := D(f);
Df := x/20 x4 K 9 x2 K 6 x K 5
(6.2)
Busquemos a raiz negativa. O processo iterativo pode agora ser implementado da seguinte forma:
O Digits := 40;
Digits := 40
(6.3)
O it := x-> evalf(x-f(x)/Df(x));
x1 := -1.5;
for i to 8 do x1 := it(x1)
end do;
f x
it := x/evalf x K
Df x
x1 := K1.5
x1 := K1.305882352941176470588235294117647058824
x1 := K1.216146061136080173760996927982363843430
x1 := K1.197775193827272108126242798711852317755
x1 := K1.197076944817166324720774663089980288863
x1 := K1.197075966364972296272551048218807791636
x1 := K1.197075966363053370759790832262420180433
x1 := K1.197075966363053370759783451616526918495
x1 := K1.197075966363053370759783451616526918495
(6.4)
As outras raízes são dadas por
O it := x-> evalf(x-f(x)/Df(x));
x2 := .3;
for i to 8 do
x2 := it(x2)
end do;
it := x/evalf x K
f x
Df x
x2 := 0.3
x2 := 0.4555746509129967776584317937701396348013
x2 := 0.4434827618709262763416769299902774224404
x2 := 0.4434256635022914148743827692660945472992
x2 := 0.4434256621831649768671517698982280196579
x2 := 0.4434256621831649761629765071208659677392
x2 := 0.4434256621831649761629765071208659675387
x2 := 0.4434256621831649761629765071208659675386
x2 := 0.4434256621831649761629765071208659675386
O it := x-> evalf(x-f(x)/Df(x));
x3 := 1.1;
for i to 10 do
x3 := it(x3)
end do;
f x
it := x/evalf x K
Df x
x3 := 1.1
x3 := 1.641955241460541813898704358068315665489
x3 := 1.442808304247545711669117624894270413942
x3 := 1.340743032099843726747024916568099491291
x3 := 1.312419896681541966304776443656985621841
x3 := 1.310362371603789310443717441415657363606
x3 := 1.310351946682090752194964357667833976967
x3 := 1.310351946415413215202161408801849091767
x3 := 1.310351946415413215027657244683799723504
x3 := 1.310351946415413215027657244683799723429
x3 := 1.310351946415413215027657244683799723429
Portanto, as raízes são
O x1, x2, x3;
K1.197075966363053370759783451616526918495,
0.4434256621831649761629765071208659675386,
1.310351946415413215027657244683799723429
(6.5)
(6.6)
(6.7)
Definamos agora um procedimento que tem como entrada, a função, um ponto inicial e o número de
iterações.
O restart;
O Digits := 20;
Digits := 20
(6.8)
O newton := proc (g, x0, N)
local it, i, Dg, x;
(6.9)
O
Dg := D(g);
it := x-> evalf(x-g(x)/Dg(x));
x := x0; for i to N do x := it(x) end do; x end proc;
newton := proc g, x0, N
local it, i, Dg, x;
Dg := D g ; it := x/evalf x K g x / Dg x ; x := x0; for i to N do
x := it x
end do;
x
end proc
(6.9)
Exemplo. Determinemos as raízes reais de
f x = exp x K 1.5 K arctan x .
O f := exp(x)-1.5-arctan(x);
f := ex K 1.5 K arctan x
O plot(f, x = -1 .. 1.5);
(6.10)
1.5
1
0.5
K1
K0.5
0
0.5
1
1.5
x
K0.5
Usando nosso procedimento, temos
O F := unapply(f, x);
F := x/ex K 1.5 K arctan x
O newton(F, .5, 10);
0.76765326620127889818
(6.11)
(6.12)
Procedimento Recursivo
Uma outra maneira de implementar o processo de Newton-Raphson é através de um procedimento
recursivo
O x := proc (n)
option remember;
global f, x0;
if n <= 0 then x0
else x(n-1)-f(x(n-1))/(D(f))(x(n-1))
end if
end proc:
Por exemplo, vamos encontrar a primeira raiz positiva de f x = x sin x Kln x .
O f := x-> x*sin(x)-ln(x);
f := x/x sin x K ln x
O plot(f(x), x = 1 .. 3);
(6.13)
1
0.8
0.6
0.4
0.2
0
1.5
K0.2
2
x
2.5
3
K0.4
K0.6
O ponto inicial pode ser x = 2.5.
O x0 := 2.5;
Com 6 iterações,
O x(6);
x0 := 2.5
2.7649217502367386442
Para ver a sequência de pontos gerados,
O seq(x(n), n = 0 .. 6);
2.5, 2.8213776546862923541, 2.7664994616157348773, 2.7649230797328876762,
2.7649217502376848093, 2.7649217502367386442, 2.7649217502367386442
(6.14)
(6.15)
(6.16)
Exercícios
Nos exercícios 1 a 3 encontre todas raízes reais das equações abaixo, usando o método de NewtonRaphson, com tolerância de 10K8 . Determine também os pontos críticos (extremos relativos) no
intervalo K4, 4 :
1. f x = arctan x K sin x C 1
2. f x = sin x2 C 1 C x K 1
3. f x = x5 K 4 x4 C 3 x3 K x2 C x K 1
4. Verifique se o método de Newton-Raphson funciona para o caso de raízes múltiplas.
Polinômios e raízes complexas
Embora não tenhamos estudado a teoria dos processos iterativos para a obtenção de raízes complexas, o
método de Newton é capaz de encontrar todas as raízes de um polinômio, reais e complexas.
Consideremos um polinômio de grau n, escrito na seguinte forma
Pn x := a0 C a1 x C a2 x2 + . . . + an xn . an s 0.
Então é possível mostrar que os zeros de P
x
n
estão localizados, no plano complexo, dentro de uma
região circular definida por
x % 1 C max
ak
an
, 1 % k , k % n.
Exemplo. Determinemos todas as raízes complexas do polinômio
P x = 4 x4 K 3 x3 C x2 K 5 x K 8.
O f d x/4 x^ 4 K 3 x^ 3 C x^ 2 K 5 x K 8;
f := x/4 x4 K 3 x3 C x2 K 5 x K 8
Inicialmente faremos um procedimento para coletar os coeficientes do polinômio:
O n:=degree(f(x));
n := 4
O c[0]:=f(0);
c0 := K8
O
(6.17)
(6.18)
(6.19)
for k from 1 to n do
c[k]:=coeff(f(x),x^k);
od;
c1 := K5
c2 := 1
c3 := K3
c4 := 4
As raízes devem estar dentro de círculo, no plano complexo, de raio r, dado por
O r0 := 1+max(seq(c[i]/c[n], i = 0 .. n-1));
5
r0 :=
4
Podemos resumir estes comandos em um procedimento:
O Raio := proc (f)
local n, c, k, r0;
n := degree(f(x));
c[0] := f(0);
for k to n do
c[k] := coeff(f(x), x^k)
end do; r0 := 1+max(seq(c[i]/c[n], i = 0 .. n-1));
(6.20)
(6.21)
(6.22)
eval(r0)
end proc;
Raio := proc f
local n, c, k, r0;
n := degree f x ;
c 0 := f 0 ;
for k to n do c k := coeff f x , x^ k end do;
r0 := 1 C max seq c i / c n , i = 0 ..n K 1 ;
eval r0
end proc
O Raio(f);
5
4
(6.22)
(6.23)
Geramos ângulos aleatórios
O r := rand(0 .. 1);
r := proc
proc
option builtin = RandNumberInterface; end proc 6, 2, 1 end proc (6.24)
O theta := rand(0 .. 24);
q := proc
proc
option builtin = RandNumberInterface; end proc 6, 25, 5 end proc (6.25)
O N := 0;
N := 0
(6.26)
O while N<degree(f(x)) do
L:=[]:
for m to 30 do
R[m]:=(r()*(cos(theta())+I*sin(theta()))):
x[m]:=newton(f,R[m],20):
if abs(Im(x[m]))<10^(-30) then
x[m]:=Re(x[m])
fi;
if not evalf[5](x[m]) in evalf[5](L) then
L:=[op(L),x[m]];
fi;
od:
N:=nops(L);
od:
op(L);
K0.014468405394238195539 C 1.2477276721654362863 I, K0.80893860824194532739,
K0.014468405394238195539 K 1.2477276721654362863 I, 1.5878754190304217185
O mesmo resultado pode ser obtido através do comando fsolve do Maple:
O printlevel:=1:
O fsolve(f(x),x,complex);
K0.80893860824194532739, K0.014468405394238195540 K 1.2477276721654362863 I,
K0.014468405394238195540 C 1.2477276721654362863 I, 1.5878754190304217185
Vamos agora fazer um procedimento que plota as raízes de um polinômio no plano complexo:
O raizplot:=proc(p::polynom(constant,x))
(6.27)
(6.28)
O
local R, points;
O
R:=[fsolve(p,x,complex)];
O
points:=map(z->[Re(z),Im(z)],R);
O
plot(points,style=point,symbol=circle);
O end:
Usando o comando randpoly gera um polinômio randômico temos
O y:=randpoly(x,degree=24);
y := 98 x24 C 77 x16 K 60 x10 K 28 x7 K 47 x5 C 34 x2
O yf:=unapply(y,x);
yf := x/98 x24 C 77 x16 K 60 x10 K 28 x7 K 47 x5 C 34 x2
O RR:=Raio(yf);
25
RR :=
14
O with(plots):
O g1:=raizplot(y):
O g2:=polarplot([RR],t=0..2*Pi):
O display([g1,g2]);
(6.29)
(6.30)
(6.31)
p
2
3p
4
p
p
4
0 0.20.40.60.8 1 1.2 1.41.6
5p
4
0
7p
4
3p
2
1. A partir dos comandos definidos para obter as raízes de um polinômio, construa um procedimento
que determina todas as raízes de um polinômio e as plota no plano complexo. .
2. Utilize o método de Newton-Raphson para encontrar as 5 primeiras raízes reais positivas de
f := x7 sin x K 3 x2 .
3. Determine todas as raízes de f = 47 x12 K 16 x11 C 38 x6 K 53 x4 K 42 x2 K 37 x +9 e plote-as no
plano complexo.
Método da Secante
Esta introdução é baseada em Peter Stone. Neste método devemos partir de duas aproximações x0 e
x1 para uma raiz de f x = 0. Sejam y0 = f x0 e y1 = f x1 . Podemos considerar a linha que junta
os pontos (x0 , y0 ) e (x1 , y1 ) como uma aproximação linear para a função f x próximo a x0 e x1 , como
mostra a figura abaixo:
y = f(x)
(x 1 ,y1 )
x0
x2
x1
(x 0 ,y0 )
Esta reta tem equação
y K y1 =
y0 K y1
x0 K x1
x K x1 .
Ela encontra o eixo x no ponto
x2 = x1 K
x0 K x1
y1
y0 K y1
Podemos tomar o valor x2 como uma nova aproximação para a raiz r .
O processo pode então ser repetido com x2 e x1 substituindo x1 e x0 respectivamente.
A fórmula de iteração é dada por:
xn C 1 = xn K
xn K 1 K xn
f xn ,
f xn K 1 K f xn
O procedimento secap use pode ser usado para um passo do método da secante:
O restart;
O secap := proc(a,b)
local fa,fb;
fa := f(a);
fb := f(b);
evalf(b - fb*(a - b)/(fa - fb));
end proc:
Testemos este procedimento para determinar 2 :
O f := x-> x^2-2;
f := x/x2 K 2
O plot(f(x), x = -1 .. 2);
(7.1)
2
1
K1
0
1
x
2
K1
K2
Tomemos dois pontos quaisquer próximos à raiz:
O a := 0.; b := 1.;
a := 0.
b := 1.
O c := secap(a, b);
c := 2.000000000
(7.2)
(7.3)
Renomeamos a e b e repetimos o processo:
O a := b;
b := c;
a := 1.
b := 2.000000000
(7.4)
Portanto,
O c := secap(a, b);
c := 1.333333333
Automatizaremos o processo daqui para frente:
O while abs(b-c) > 10^(-7) do
a := b;
b := c;
c := secap(a, b)
end do;
a := 2.000000000
b := 1.333333333
c := 1.400000000
a := 1.333333333
b := 1.400000000
c := 1.414634146
a := 1.400000000
b := 1.414634146
c := 1.414211438
a := 1.414634146
b := 1.414211438
c := 1.414213562
a := 1.414211438
b := 1.414213562
c := 1.414213562
De fato,
O sqrt(2.);
1.414213562
(7.5)
(7.6)
(7.7)
Construiremos agora um procedimento que toma como entrada uma expressão algébrica, dois pontos
iniciais, tolerância e número máximo de iterações.
O Secante:= proc(expr::algebraic,x,x1,x2,Delta,it)
local x3,nit,delta,f,x1t,x2t, xdg3,g3,f3,x_2,x_1,dg3;
f:=unapply(expr,x);
nit:=0: delta:=1:
dg3:=f(x2t)-f(x1t):
g3:=(x1t*f(x2t)-x2t*f(x1t))/dg3:
f3:=unapply (g3,x1t,x2t);
x_2:=x2: x_1:=x1:
while nit < it and delta > Delta do
x3:=f3(x_1,x_2);
x_1:=x_2:
x_2:=x3:
nit:=nit+1:
delta:=abs(x_2-x_1);
od;
print(x_2);
end:
Apliquemos o método da secante ao problema de encontrar as raízes reais de
f x = x5 K 3 x4 C 5 x3 C 7 x2 K x C 2 .
O f := x^5-3*x^4+5*x^3+7*x^2-x+2;
f := x5 K 3 x4 C 5 x3 C 7 x2 K x C 2
O plot(f, x = -2 .. 2);
(7.8)
40
20
K2
K1
0
K20
1
x
2
K40
K60
K80
O Digits := 19;
Digits := 19
O Secante(f, x, -1.5, -.5, 10^(-8), 10);
K1.052704929237755516
(7.9)
(7.10)
O processo iterativo da secante está implementado na biblioteca Student do Maple:
O with(Student[NumericalAnalysis]):
O f := x^2-2;
f := x2 K 2
(7.11)
O Secant(f, x = [1, 8], tolerance = 10^(-7));
1.414213562372920999
(7.12)
O Secant(f, x = [1.1, 8], tolerance = 10^(-7), output = sequence)
;
1.1, 8., 1.186813186813186813, 1.251196172248803828, 1.429418678628810507,
(7.13)
1.413288887062976678, 1.414208616456125317, 1.414213563990556179,
1.414213562373092221
O processo pode visualizado geometricamente:
O Secant(f, x = [0, 1], tolerance = 10^(-3), output = plot);
5 iteration(s) of the secant method applied to
f x = x2 K 2
with initial points a = 0. and b = 1.
b
p5
x
f x
Iniciando o processo com o intervalo K1, 2 , temos
O Secant(f, x = [-1, 2], tolerance = 10^(-3), output = plot);
5 iteration(s) of the secant method applied to
f x = x2 K 2
with initial points a = K1. and b = 2.
The stopping criterion is not met
a
b
x
f x
Suponhamos que agora queremos que o critério de parada seja f xk ! e
O Secant(f(x), x = [-1, 2], tolerance = 10^(-8),
stoppingcriterion = function_value);
1.414213562057320463
(7.14)
O processo iterativo abaixo está animado:
O Secant(f, x = [-1, 2], output = animation, stoppingcriterion =
function_value);
7 iteration(s) of the secant method applied to
f x = x2 K 2
with initial points a = K1. and b = 2.
a
p7
b
x
f x
Exercícios
Nos exercícios 1 a 3 encontre todas raízes reais das equações abaixo, usando o método da secante, com
tolerância de 10K8 . Determine também os pontos críticos (extremos relativos) no intervalo K4, 4 :
1. f x = arctan x K sin x C 1
2. f x = sin x2 C 1 C x K 1
3. f x = x5 K 4 x4 C 3 x3 K x2 C x K 1
4. Verifique se o método da secante funciona para o caso de raízes múltiplas.
5. Verifique a diferença de tempo para os métodos de Newton-Raphson e secante. Por exemplo,
3
comparamos o tempo para a determinação de 5 usnado o método da secante e o comando fsolve do
Maple:
O g := x^3-5;
g := x3 K 5
(7.15)
O Digits := 19;
Digits := 19
(7.16)
O st := time();
st := 2.218
(7.17)
O Secante(g, x, 1., 2., 10^(-9), 20);
(7.18)
O time() - st;
O st := time();
1.709975946676696989
(7.18)
0.
(7.19)
st := 2.219
(7.20)
O fsolve(g, x = 1 .. 2);
1.709975946676696989
O time()-st;
0.017
(7.21)
(7.22)