Matemática computacional
MEFT - Instituto Superior Técnico (IST)
João Cunha no 67921 Miguel Prata no 67933 Natália Morais no 67936
Trabalho de Casa com Componente Computacional No 26
19 de Novembro de 2010
Objectivos
Na primeira parte deste trabalho pretende-se:
1. Usar o método de Steffensen para calcular valores aproximados da solução z da equação f (x) = 0, onde f : I ⊂ R → R
é uma função continuamente diferenciável numa vizinhança da solução z, sendo aplicado este método, primeiramente
a uma função w que escolhemos e em segundo lugar à função f do exercı́cio [3d.2].
2. Usar o método de Newton generalizado para calcular um valor aproximado da solução z de um sistema de equações
f (x) = 0 onde F : D ⊂ Rd → Rd , d ≥ 2, é uma função continuamente diferenciável numa vizinhança de z. Usar o
método de eliminação de Gauss com pesquisa parcial de pivot para resolver o sistema linear que determina a diferença
de duas iteradas sucessivas do método de Newton. Um dos sistemas f é nos dado e outro é escolhido por nós (exercı́cio
[5a.7]).
1
1.1
Método de Steffensen
Note-se que as condições dadas por este teorema são suficientes mas não necessárias.
Introdução
Teorema 2. Seja z ponto fixo de g, suponha-se que g ∈ C 1
0
Existem vários métodos de cálculo de zeros de uma função. numa vizinhança de z e que |g (x)| < 1, então a sucessão
gerada por (2) converge para z e são válidas as conclusões do
Grande parte deles reduzem-se a considerar
Teorema 1 desde que x0 esteja suficientemente perto de z
z = g(z)
(1)
Há também mais um teorema que interessa referir para a
seguir
estudar o Método de Steffensen. Note-se que se tem
uma vez que g (a função iteradora) transforma z em si mesmo,
e
=
z
− xm e em+1 = z − xm+1 .
z é denominado ponto fixo de g e desta forma pode-se obter m
Teorema 3. Seja z um ponto fixo de g ∈ C p numa vizinhança de z e g (p) contı́nua nessa vizinhança. Seja ainda
o método que consiste em construir a sucessão x0 , x1 , x2 , ...
tal que
xm+1 = g(xm ), m ∈ N0
(2)
g 0 (z) = ... = g (p−1) (z) = 0
g (p) 6= 0
sendo este método apelidado de Método do Ponto Fixo.
Então para qualquer x0 suficientemente perto de z, a sucessão
O Método do ponto fixo possui naturalmente condições de
(2) converge para z e
convergência. Vejamos primeiro a seguinte definição:
(−1)(p+1) (p)
em+1
lim
=
g (z)
(5)
Definição 1. Uma função g diz-se Lipschitziana no intervalo
p
m→∞ em
p!
I = [a, b] se existir um constante L ≥ 0 tal que
e portanto a convergência é supralinear de ordem p e tem-se
|g(y) − g(x)| ≤ L|y − x|,
∀y, x ∈ I
(3) que o coeficiente assintótico de convergência K∞ = | g(p) (z) | .
p!
O Método de Steffensen pode tratar-se como um caso parNo caso em que L < 1 a função g diz-se contractiva. A L
ticular do método do ponto fixo já exposto. Pode-se escrever
chama-se constante de contractividade.
a forma iterativa do método de Steffensen como
Nas condições da definição anterior
[f (xm )]2
xm+1 = xm −
,m ∈ N
(6)
f (xm + f (xm )) − f (xm )
L = max |g 0 (x)|
(4)
x∈I
que se pode considerar um método de ponto fixo bastando
Sabendo isto pode-se agora analisar o teorema seguinte, escolher
[f (x)]2
que chamamos do ponto fixo.
g(x) = x −
(7)
f (x + f (x)) − f (x)
Teorema 1. Se existir um intervalo I = [a, b] no qual a
Desta forma, são válidas para este método, as condições
função g é contractiva e g(I) ⊂ I, i. e., a ≤ g(x) ≤ b, então: de convergência impostas à função g pelo Teorema 1.
Para estudar a ordem de convergência deste método va1. g possui pelos menos um ponto fixo z em I
mos reescrevê-lo como
2. z é o único ponto fixo em I
[f (x)]2
F (x) = x −
(8)
f (x + f (x)) − f (x)
3. para qualquer x0 em I a sucessão xm gerada por (2)
converge para z
Note-se que xm , f (xm ) e f (f (xm )) são 3 iteradas consecutivas
de um método iterativo cuja função iteradora é
1
|xm+1 − xm | ,
4. |z − xm | ≤ 1−L
f
.
Se
f
∈
C 1 numa vizinhança de z, pode-se verificar que
L
|z − xm+1 | ≤ 1−L |xm+1 − xm |
limx→z F (x) = z e portanto podemos definir por continuidade
F
em z escrevendo F (z) = z. Se f 0 (z) 6= 1 então F 0 (z) = 0, e
Lm
5. |z − xm | ≤ Lm |z − x0 | ≤ 1−L
|x1 − x0 |
portanto a convergência do método de Steffesen é pelo menos
1
quadrática (pelo Teorema 3). Note-se que para f 0 (x) = 0 não
se deve aplicar o método, pois conduz a um enorme número
de iteradas.
Verifica-se, em comparação com o método do ponto fixo,
uma aceleração da convergência do método do ponto fixo
quando a primeira derivada da função iteradora é diferente de
zero (o método iterativo original converge linearmente, p = 1,
ou não converge). Se a primeira derivada da função iteradora
é zero, então (pelo Teorema 3) o método tem convergência
supralinear (p > 1).
Conclui-se assim que, em termos da expressão (7), o
método de Steffensen, tendo em conta os Teoremas 1 e 2,
converge quadraticamente (p = 2) para |g 0 (z)| < 1 (tendo
em atenção que para o caso de g 0 (z) = 0 não se deve aplicar
este método) e converge linearmente (p=1) quando g 0 (z) = 1.
Tem-se também (do Teorema 3) um coeficiente assintótico
00 g (z) (9)
K∞ = 2! de modo a que o novo subintervalo contenha um zero de f ; se
pelo contrário forem f (bm ) e f (xm+1 ) a ter sinais diferentes
então faça-se
am+1 = xm+1
1.2.1
bm+1 = bm
Deste modo garante-se que o novo subintervalo Im+1 =
[am+1 , bm+1 ] contém pelo menos um zero de f .
No ciclo (ver CICLO1, no código steffensen.cpp) de
cálculo relativo à determinação do subintervalo final é utilizada a condição de paragem bm − am < 0.01, ou seja,
consegue-se assim um intervalo (que é impresso, para informação do utilizador) suficientemente pequeno em que se
vai aplicar o Método de Steffensen.
Assim, fazendo, por exemplo, x0 = bm aplica-se o método
de Steffensen já exposto. O ciclo (ver CICLO2, no código
steffensen.cpp) de aplicação do método de Steffensen imprime a iteração, o xm , o em ”exacto”e em estimado, os valores da sucessão |em |/|em−1 |2 e ainda uma estimativa de g 0 .
Utilizando as conclusões do Teorema 1 pode-se concluir
que é admissı́vel a estimativa
xm+1 − xm
(12)
g 0 (xm ) ≈
xm − xm−1
(Pode ainda convergir para |g 0 (x)| > 1, sendo que assim, g
não é contractiva)
1.2
e
, possı́vel quando m ∈ N0 : m ≥ 2 e do mesmo teorema
pode-se concluir também que para o caso de convergência
supralinear (caso do Método de Steffensen)
Aplicação do método
Programa steffensen.cpp
|em | ≈ |xm+1 − xm | e
Em resposta ao TCCC foi desenvolvido o programa steffensen.cpp (ver Anexo 1 e 2) escrito na linguagem de programação conhecida por C++.
Para compilar e executar o programa aconselha-se a
utilização do sistema operativo Windows com o programa
Dev-cpp da Bloodshed disponı́vel para download no link
http://sourceforge.net/projects/dev-cpp/files/Binaries/DevC%2B%2B%204.9.9.2/devcpp-4.9.9.2 setup.exe/download .
O programa pretende aplicar o método de Steffensen num
intervalo I = [a, b] fornecido pelo utilizador do programa.
Depois de inserido o intervalo o programa verifica se existe mais que um zero (da função f ) nesse intervalo e verifica se f (a)f (b) > 0, condição na qual não se pode aplicar
o Corolário do Teorema de Bolzano o que não garante a existência de zero nesse intervalo. Caso se verifique uma das
duas situações anteriores o programa é encerrado.
Caso contrário o programa tenta encontrar um novo intervalo I 0 que verifique as condições de convergência do Teorema
1. Para isso aplica-se o método da bissecção à função f no
intervalo I.
O método da bissecção aplica-se a uma função que verifique as condições do Corolário do Teorema de Bolzano (f
é contı́nua em I e f (a)f (b) < 0). Consiste em construir
subintervalos Im = [am , bm ] ⊂ I = [a, b] por divisões sucessivas a meio e relativamente aos quais também se verifica que
f (am )f (bm ) < 0. Deste modo, confina-se o zero da função
f a intervalos tão pequenos quanto se queira. Concretizando
tem-se
Im = [am , bm ] m = 0, 1, 2, . . .
(10)
|em+1 | |xm+1 − xm |
(13)
00
Tira-se também de (8) que K∞ = | g 2!(z) |, o que se pode
verificar observando a convergência da sucessão |em |/|em−1 |2 .
Durante o ciclo são recolhidos dados para avaliar as
condições suficientes de convergência exigidas pelo Teorema
1. Se alguma não se verificar no final do ciclo o utilizador
é informado. Como essas condições são apenas suficientes e
não necessárias apresenta-se sempre a tabela para que outras
conclusões possam ser retiradas. O critério de paragem é tal
que |xm − xm−1 | < sendo dado.
Caso o método não convirja em 1000 iteradas é apresentada essa informação.
No final do ciclo é apresentado o resultado encontrado para o zero, assim como a informação do coeficiente
assintótico de convergência.
Menciona-se a utilização de uma particular função externa
à função principal main, a função quicksort. O seu objectivo
é ordenar por ordem crescente os componentes do vector g_d
que se traduz nos valores de |g 0 (xm )|. Desda forma, conseguese identificar L = maxx∈I 0 |g 0 (x)|, como sendo o último elemento do vector g_d. Avalia-se assim se g é contractiva.
Outras funções externas à função principal main têm o
seguinte objectivo: a função f retorna o valor de f num determinado ponto dado como argumento; f_d retorna a derivada
de f num determinado ponto dado como argumento; f_d2 retorna a segunda derivada de f num determinado ponto dado
como argumento.
Informação referentes a comandos próprios da linguagem
C++ encontram-se nos comentários ao código.
em que a0 = a e b0 = b e seja xm+1 o ponto médio do intervalo
Im , i. e.,
1.2.2 Função [3d.2]
am + bm
xm+1 =
(11)
Neste TCCC pretendia-se aplicar o método exposto à
2
função
2
Considera-se então 3 situações: se f (xm+1 ) = 0 então
f (x) = ex −1 + 10 sin(2x) − 5
(14)
xm+1 é um zero e o processo iterativo termina; se f (am ) e
que
se
trata
de
uma
equação
transcendente.
f (xm+1 ) tiverem sinais diferentes faça-se
Para averiguar a posição dos zeros, fez-se uma repream+1 = am e bm+1 = xm+1
sentação gráfica da função f que se pode observar a seguir
2
método, a imagem dos extremos são de sinais diferentes e a
f é monótona crescente/decrescente nesse intervalo.
Outras análises
prendem-se com a comparação
ao output
em entre a sucessão em−1 2 e o valor de K∞ , tal não se verifica
novamente devido a erros de arredondamento.
1.3
Através da informação dos Outputs, do conhecimento
teórico e do conhecimento adquirido na produção do programa, podem-se concluir algumas vantagens/desvantagens
em relação ao Método de Steffensen.
Do lado das vantagens insere-se a rapidez de convergência,
visto que o método converge muito rapidamente sendo
necessárias poucas iteradas; em comparação com outros
métodos tais como o método de Newton ganha vantagem uma
vez que tem uma ordem de convergência similar e não necessita de qualquer informação acerca da derivada da função.
Do lado das desvantagens inserem-se o tempo de computação para funções f complicadas que é elevado, principalmente devido ao termo f (x + f (x)); e a principal desvantagem é o facto do intervalo de convergência local ser bastante
reduzido, sendo necessário inserir um x0 muito próximo do
zero z, caso contrário o método não converge para a solução.
Quanto a este último assunto verifica-se que o método converge mais rapidamente quanto mais perto estiver x0 da
solução z, i. e., é necessário um menor número de iteradas.
Gráfico 1: Representação gráfica da função f
obtém-se da representação gráfica os zeros
z1
z2
z3
z4
= −1.59019078359636 . . .
= 0.2395826576263948 . . .
= 1.489630021514447 . . .
= 1.801987287952526 . . .
Note-se que a função f é contı́nua pois é a soma de funções
contı́nuas, a primeira a composição de uma exponencial com
uma quadrática (contı́nua em R), a segunda uma função
trigonométrica (contı́nua em R) e a última uma função constante.
Desta forma, o Corolário do Teorema de Bolzano garante a
existência de pelo menos um zero de f num intervalo I = [a, b]
em que f (a)f (b) < 0. Se para além disso se verificar que f 0
existe, é contı́nua e não se anula em (a, b) conclui-se, pelo 1.4
Teorema de Rolle que a solução z é única em (a, b).
Note-se também que f é de classe C ∞ em R pois é a soma
de duas funções de classe C ∞ em R.
1.4.1
1.2.3
Conclusões
Aplicação na resolução de um problema
de Fı́sica
Problema
Para aplicação do método a um problema de fı́sica, foi
enunciado o seguinte:
Análise de output
Para analisar os resultados do programa, escolheu-se
primeiro uma função w(x) =√x2 − 2. Verifica-se analiticamente que w(x) = 0 ⇒ x = ± 2 ≈ ±1.41421.
Analisemos os seguintes outputs presentes no Anexo 1.
Repare-se que e_m1 é a estimativa do erro obtida pela estimativa (13) e e_m2 é o erro. Pode observar-se nos outputs
o que foi deduzido teoricamente.
Os resultados obtidos coincidem, tendo em conta precisão da máquina, com os resultados analı́ticos. A precisão
da máquina dita também a máxima precisão do resultado e
sendo que em C++ só é possı́vel obter resultados com no
máximo 6 casas decimais, é considerada a condição de paragem |xm+1 − xm| < 1E − 6.
em Nota-se também que a sucessão em−1
2 converge para
algo similar a K∞ , o facto dos valores diferirem substancialmente deve-se a erros de arredondamento.
Problema 1. Uma partı́cula desloca-se num espaço a 3 dimensões numa trajectória definida pela função P (t) :]0, 4π[→
R3 , que descreve num referencial cartesiano Oxyz as posições
da partı́cula em função do tempo t. A função P (t) é dada por:
P (t) = (x(t), y(t), z(t)) =
t
t
t
cos (t), 2 + cos
cos (t), sin
=
2 + cos
2
2
2
(15)
Em que x(t) :]0, 4π[→ R, y(t) :]0, 4π[→ R e z(t) :]0, 4π[→
R.
A trajectória da partı́cula pode ser ilustrada pelo seguinte
gráfico:
Analise-se agora o output do programa para o caso de f .
Podem apreciar-se os outputs presentes no Anexo 2.
Nos Outputs 3,4,5 e 6 observa-se o output para cada um
dos zeros da função f e que estão coerentes com os zeros
encontrados geometricamente. Nos Outputs 7 e 8 observamse situações em que o programa identifica mais que um zero
no intervalo ou não existe zero nesse intervalo (podem ser
observados intervalos nessa situação quando considerando o
gráfico1).
Embora não tenha sido feito até agora pode-se confirmar
que no intervalo que o programa encontra para aplicação do
Grafico 2: Gráfico das posições da partı́cula no referencial Oxyz
3
Determine os extremos de posição para cada função x(t), modo similar ao exemplos anteriores. Vai-se tentar por isso
y(t) e z(t), ou seja as coordenadas x, y e z máximas e aplicar o método de Steffensen para encontrar os zeros das
mı́nimas da trajectória da partı́cula.
funções x0 (t), y 0 (t) e z 0 (t).
1.4.2
Resolução
1.4.3
Sabe-se que os máximos e mı́nimos de uma função de
I → R em que I ⊂ R correspondem aos zeros da sua derivada,
se a função for diferenciável.
t
2 é uma função diferenciável em R, pois é o produto de
1
duas funções
diferenciáveis em R, t e 2 .
t
cos 2 é uma função diferenciável em R pois é a composta
t
de duas funções
diferenciáveis em R, cos (t) e 2 .
t
2 + cos 2 é diferenciável em R,
pois é a soma de duas
funções diferenciáveis
em R, cos 2t e uma função constante.
x(t) = 2 + cos 2t cos (t) é uma função diferenciável em
R, pois é o produto de duas funções diferenciáveis em R,
2 + cos 2t e cos (t), logo a função x(t) é diferenciável no
intervalo ]0, 4π[ em particular.
y(t) = 2 + cos 2t sin (t) é uma função diferenciável em
R, pois é o produto de duas funções diferenciáveis em R,
2 + cos 2t e sin (t), logo a função y(t) é diferenciável no intervalo ]0, 4π[ em particular.
z(t) = sin 2t é uma função diferenciável em R, pois é a
composta de duas funções diferenciáveis em R, sin t e 2t , logo
a função z(t) é diferenciável no intervalo ]0, 4π[ em particular.
Analogamente, verifica-se que as três funções são de classe
C ∞.
Pode-se então calcular a derivada de cada uma das funções
x(t), y(t) e z(t).
O programa steffensen aplicacao.cpp é semelhante ao das
alı́neas anteriores, tendo sido adaptado para a resolução do
problema em questão. As diferenças mais significativas são
as seguintes:
A função pos é a função das posições da partı́cula, sendo
a função vel a velocidade, dvel a primeira derivada da velocidade e *dvel a segunda derivada da mesma.
A função graf escreve no ficheiro ”trajectoria.txt”, as coordenadas x, y e z em cada instante t com dt =h=0.001 e desenha os gráficos da trajectória da partı́cula (Figura 1) e das
velocidades segundo os eixos x, y e z (Figura 2), já mostrados
acima.
A função steffensen aplica o método de Steffensen às
funções x0 (t), y 0 (t) ou z 0 (t), como já foi descrito nas alı́neas
anteriores.
Na função main é chamada a função graf para desenhar
os gráficos já falados. pergunta-se ao utilizador se pretende
encontrar os extremos das posições no eixos dos x’s, y ou z
e consoante a escolha, é chamada a função steffensen para
aplicar o método.
Mais informação sobre o código do programa encontra-se
nos comentários ao código.
1.4.4
t
t
cos (t) − sin (t) cos
+2
2
2
t
t
1
sin (t) + cos (t) cos
+2
y 0 (t) = − sin
2
2
2
t
1
z 0 (t) = cos
2
2
1
x0 (t) = − sin
2
Que são representadas no seguinte gráfico:
Os outputs correspondentes ao programa steffensen aplicacao.cpp encontram-se no Anexo 3 (Outputs 9
a 18).
(17)
Verifica-se que os zeros encontrados estão de acordo com
os previstos, com precisão até à sexta casa decimal. (Outputs
(18) 9 a 17)
em A sucessão em−1
2 converge para K∞ (tendo em conta os
erros de arredondamento).
Repare-se também que a derivada da função g, respeita
sempre a condição |g 0 (x)| < 1 ∀x ∈ I, sendo I um intervalo
válido fornecido pelo utilizador.
O Output 18 corresponde a uma situação em que existe
mais do que um zero no intervalo, tendo sido dado erro pelo
programa.
Grafico 3: Gráfico das velocidades da partı́cula em x, y e zem função
do tempo
= 1.4506848...
= 4.4342238...
= 8.1321468...
= 11.115686...
Método de Newton Generalizado e
Eliminação de Gauss com pesquisa
parcial de pivot
2.1
Analizando os gráficos do Gráfico 3, verifica-se que no intervalo pretendido, x0 (t) tem 3 zeros, y 0 (t) 4 zeros e z 0 (t) 2
zeros, sendo estes:
y1
y2
y3
y4
Resultados e análise
(16)
2
x1 = 2.9109431...
x2 = 6.2831853... = 2π
x3 = 9.6554275...
Programa steffensen aplicacao.cpp
2.1.1
Introdução
Método de Newton Generalizado
Consideremos um sistema de equações não lineares do
z1 = 3.1415927... = π
z2 = 9.4247780... = 3πtipo:

2

x1 + 3 sin x2 = 4
x41 − 6x21 x22 + x43 = −1


x1 x2 x3 − x21 + x32 = 2
Já vimos que a todas as funções são diferenciáveis no intervalo ]0, 4π[, sendo consequentemente constı́nuas nesse intevalo, podendo assim ser aplicado o Teorema de Bolzano de
4
(19)
(iii) sendo 0 = 2kx(1) − x(0) k e K =
sigualdade 2K0 < 1;
Podemos definir:

2

f1 (x1 , x2 , x3 ) = x1 + 3 sin x2 − 4
f2 (x1 , x2 , x3 ) = x41 − 6x21 x22 + x43 + 1


f3 (x1 , x2 , x3 ) = x1 x2 x3 − x21 + x32 − 2
M2
2M1 ,
verifica-se a de-
(0)
(20) (iv) B0 (x ) ⊂ D
Então conclui-se
Visto que as equações são não lineares, não é possı́vel escrever
o sistema em notação matricial. Contudo, é possı́vel fazê-lo
em notação vectorial:


f1 (x1 , x2 , x3 )
F (x1 , x2 , x3 ) =  f2 (x1 , x2 , x3 ) 
(21)
f3 (x1 , x2 , x3 )
1. f tem um único zero z em B0 (x(0) )
2. a sucessão de Newton com condição inicial x(0) é bem
definida, permanece em B0 (x(0) ) converge para z
3. verifica-se a estimativa de erro a priori
1
kz − x(m) k ≤ (K0 )2m
Podemos assim, exprimir o sistema não linear como a
K
equação F (x1 , x2 , x3 ) = 0 (sendo 0 o vector nulo).
Assim, o Teorema de Kantorovich garante-nos a conPara a resolver desejamos aplicar o método de Newton
vergência do Método de Newton Generalizado para algum
Generalizado.
O método iterativo de Newton é utilizado para calular x(0) . Pode-se mostrar, que, para f de classe C 2 numa vizinhança do zero z, tal que det(Jf (z)) 6= 0. Então o método
zeros de funções de uma variável e é dado por:
de Newton converge para z desde que x(0) esteja suficientef (xn )
, n∈N
(22) mente perto de z, ou seja, garante-se a convergência local do
xn+1 = xn − 0
f (xn )
método.
Pode-se também obter a proposição seguinte
Sejam F = [f1 , f2 , . . . , fn ] e X = [x1 , x2 , . . . , xn ] vectores
de dimensão n. Para aplicar o método de Newton à equação
Proposição 5. Seja f ∈ C 2 (Vz ), onde Vz é uma vizinhança
F (X) = 0 podemos começar por escrever:
de z, zero de f , tal que det(Jf (z)) 6= 0. Então o método
de Newton, quando converge para z, tem convergência pelo
(n)
F (X )
X (n+1) = X (n) − 0 (n)
(23) menos quadrática, i. e., existe K > 0 tal que
F (X )
kz − x(m+1) k ≤ Kkz − x(m) k2 ,
Com



0
F (x1 , x2 , . . . , xn ) = 


∂f1
∂x1
∂f2
∂x1
∂f1
∂x2
∂fn
∂x1
...
...
...
...
..
.
∂f1
∂xn
∂f2
∂xn
..
.
...
∂fn
∂xn






ou
kz − x(m) k ≤
(24)
m∈N
1
(Kkz − x(0) k)2m
K
A proposição anterior garante-nos que o método converge
pelo menos quadraticamente.
Dado que F 0 (X) é uma matriz, mais concretamente a ma2.1.2 Eliminação de Gauss
triz Jacobiana, J(X), do sistema, a equação 23 toma a forma:
O método de eliminação de Gauss permite resolver sisX (n+1) = X (n) − F (X (n) )[J(X (n) )]−1
(25) temas de equações lineares através da sucessiva eliminação
de incógnitas.
Contudo, inverter uma matriz é computacionalmente
Considerando um sistema Ax = b eliminam-se sucessivaproblemático e trabalhoso, pelo que rearranjando a equação
mente
incógnitas de modo a obter um sistema equivalente
25 se obtém o método de Newton generalizado:
U x = g, sendo U uma matriz triagular superior.
J(X (n) )∆X (n) = −F (X (n) )
(26)
Designemos o sistema original por A(1) x = b(1) , e o siscom
tema após o passo κ da eliminação de Gauss por A(κ) x = b(κ) .
X (n+1) = X (n) + ∆X (n)
(27)
Supunhamos que a11 6= 0. Devemos obter coeficientes
m
i1 tal que as entradas a21 a an1 sejam nulas quando lhes
Partindo de pontos iniciais x0 , x1 , ..., xn , a equação 26
for
subtraı́da a primeira equação multiplicada por mi1 . Os
é uma equação linear, podendo ser resolvida utilizando o
coeficientes
mi1 são dados por:
método de eliminação de Gauss com pesquisa parcial de pivot,
de modo a encontrar ∆X. Iterando ambas as equações obtemos aproximações cada vez melhores da solução, devendo este
processo parar quando ∆X < p, sendo p a precisão desejada.
Repare-se agora no seguinte teorema suficiente de convergência do método, apelidado de Teorema de Kantorovich.
mi1 =
(2)
(1)
(1)
aij = aij − mi1 aij
(29)
É também necessário actualizar as entradas do vector b:
(2)
bi
1
M1 ,
(28)
E as novas entradas por:
Teorema 4. Seja D ⊂ Rn um conjunto aberto convexo e
f ∈ C 1 (D). Suponha-se que para alguma norma em Rn e
algum x(0) ∈ D são verificadas as condições
(i) det(Jf (x)) 6= 0, ∀x ∈ D;
−1
∃M1 > 0 : k [Jf (x)] k ≤
a1i1
a111
∀x ∈ D;
(1)
= bi
Com i = 2, . . . , n e j = 1, . . . , n.
(ii) ∃M2 > 0 : kJf (x) − Jf (y)k ≤ M2 kx − yk, ∀x, y ∈ D;
5
(1)
− mi1 b1
(30)
Assim, num passo κ, tal que 1 6 κ 6 n − 1 e supondo que
aκκ =
6 0 (elemento pivot):
2.2.2
Função [5a.7]
O sistema de equações não lineares a resolver é:
(κ)
miκ =
aiκ
(31)
(κ)

2

xy − z + 3 = 0
xyz − x2 + y 2 − 2 = 0

 x
e − ey + z − 3 = 0
aκκ
(κ)
(κ)
aκ+1
= aij − miκ aκj
ij
(32)
bκ+1
= bκi − miκ b(κ)
κ
i
(33)
(36)
Com uma matriz J:
Com: i = j = k + 1, . . . , n

y
 yz − 2x
ex
Pesquisa Parcial de Pivot
x
xz + 2y
−ey

−2z
xy 
1
(37)
A pesquisa parcial de pivot num dado passo κ da eliminação de Gauss tem como objectivo evitar as divisões por
As soluções exactas foram determinadas utilizando
0 e também por números muito pequenos que desiquilibrem
o software Wolfram Alpha disponı́vel em http://www.
a matriz devido a erros de arredondamento, devido ao facto
wolframalpha.com/ sendo estas (ver Figura 1 no Anexo 4):
dos coeficientes m serem muito grandes (note-se que quando
o pivot → 0, m → ∞).
x = 1.29279
Assim, num dado passo κ da eliminação de Gauss devemos
encontrar o máximo absoluto numa dada coluna κ e trocar
a linha desta entrada com a linha κ, fazendo o mesmo às
y = 0.998577
mesmas linhas de x e b.
z = 2.07146
Resolução do sistema triangular
Estando o sistema transformado na forma U x = g sendo
este é facilmente resolúvel por subsituição ascendente:
gn
unn nn
Pn
(gκ − j=κ+1 ukj xj )
xn =
xk =
uκκ
Com os pontos iniciais:
(34)
x0 = 1
(35)
y0 = 1.5
z0 = 2
Com κ = n − 1, n − 2, ..., 1
2.2
2.2.1
Aplicação do método
2.2.3
Função escolhida
Programa newtongeneralizado.c
Modificámos o código do programa newtongeneralizado.c
O programa escrito para resolver sistemas de equações não de modo a que resolvesse o seguinte sistema de equações não
lineares foi escrito em linguagem C, podendo ser compilado e lineares escolhido por nós:
e executado com software apropriado (como o Dev-C++).
Primeiro são pedidos ao utilizador os pontos iniciais
(
10x2 + sin y − 20 = 0
x1 , x2 , . . . , xn .
(38)
É então inicializado um vector b e uma matriz A, ambos
x4 + 5y − 6 = 0
de dimensão n, que servirão, respectivamente, para guardar
a matriz Jacobiana e o vector F (X).
Com uma matriz J:
É agora possı́vel começar o método de Newton, utilizando
os pontos x1 , x2 , . . . , xn dados pelo utilizador.
20x cos y
Após esta etapa o programa inicia um ciclo de modo a
(39)
4x3
5
iterar sucessivamente as equações 26 e 27.
Para resolver 26 são utilizadas as funções externas
E soluções exactas (ver Figura 2):
forward_elim e solve.
A função forward_elim aplica eliminação de Gauss com
pesquisa parcial de pivot ao sistema. E a função solve resolve
x = 1.39929
o sistema por substituição ascendente.
As soluções deste sistema são ∆x1 , ∆x2 , . . . , ∆xn , que são
y = 0.433232
utilizados na equação 27.
Se a condição ∆X < 1 × 10−6 for satisfeita o ciclo é inCom os pontos iniciais:
terrompido e o programa acaba. A cada ciclo são impressos
no ecrã os valores de x1 , x2 , . . . , xn , ||X (κ) − X (κ−1) ||, ||e(κ) ||
x0 = 1
e(κ)
e || e(κ−1)
||. O erro e(κ) é dado pela maior diferença entre as
soluções exactas e as respectivas soluções determinadas nuy0 = 1
mericamente.
6
2.2.4
Análise de resultados e Conclusões
Verificámos que quanto mais próximos da solução estiverem os pontos iniciais mais rápida é a convergência.
Para analisar os resultados devemos comparar as soluções
obtidas com as consideradas exactas (ver Output 19 e 20).
O programa estima com uma precisão bastante elevada a
solução, desde que tal como foi visto no Teorema 4, os pontos iniciais estejam numa vizinhança relativamente perto da
solução. Caso contrário, surgem resultados inesperados.
O output do programa newtongeneralizado.c mostra
uma convergência após 8 iterações. Note-se que os pontos
iniciais utilizados são diferentes dos que foram utilizados no
Wolfram, contudo a solução obtida é a mesma.
3
3.1
No programa newtongeneralizado2.c são utilizados os
mesmos pontos iniciais que no Wolfram e verifica-se a convergência após 4 iteradas. Note-se que o erro e(κ) da última
iteração deste programa é menor que o do anterior. Isto pode
estar relacionado com o os erros de arredondamento introduzidos pelo método de eleminação de Gauss com pesquisa parcial
de pivot serem maiores no programa newtongeneralizado.c,
pois a matriz Jacobiana é de dimensão 3 × 3 enquanto que em
newtongeneralizado2.c a Jacobiana é uma matriz 2 × 2.
Anexos
Anexo 1
Output 5: Output do programa steffensen.cpp
Output 1: Output do programa steffensen_sqrt2.cpp
Output 6: Output do programa steffensen.cpp
Output 2: Output do programa steffensen_sqrt2.cpp
Output 7: Output do programa steffensen.cpp
3.2
Anexo 2
Output 8: Output do programa steffensen.cpp
3.3
Anexo 3
Output 3: Output do programa steffensen.cpp
Output 9: Output do programa steffensen_aplicacao.cpp
Output 4: Output do programa steffensen.cpp
7
Output 10: Output do programa steffensen_aplicacao.cpp
Output 14: Output do programa steffensen_aplicacao.cpp
Output 15: Output do programa steffensen_aplicacao.cpp
Output 11: Output do programa steffensen_aplicacao.cpp
Output 16: Output do programa steffensen_aplicacao.cpp
Output 12: Output do programa steffensen_aplicacao.cpp
Output 17: Output do programa steffensen_aplicacao.cpp
Output 13: Output do programa steffensen_aplicacao.cpp
Output 18: Output do programa steffensen_aplicacao.cpp
8
3.4
Anexo 4
Figura 1: Resultados obtidos para o sistema 36
Figura 2: Resultados obtidos para o sistema 38
Output 19: Resultados obtidos para o sistema 36 com o programa desenvolvido
Output 20: Resultados obtidos para o sistema 38 com o programa desenvolvido
9
3.5
3.5.1
Códigos dos programas
steffensen.cpp
1 /∗PROG QUE APLICA O METODO DE STEFFENSEN∗/
/∗IST , MEFT, MC 2010 ∗/
3 /∗JOAO CUNHA 67921 , NATALIA MORAIS 67936 , MIGUEL PRATA 67933 ∗/
5 #include <i o s t r e a m >
#include <cmath>
7 #include <c s t d l i b >
// h e a d e r que p o s s u i f u n ç o e s de i n p u t / o u t p u t
// h e a d e r que p o s s u i f u n ç o e s s t a n d a r t
// h e a d e r com f u n ç o e s matematicas
9 using namespace s t d ; // namespace com informaçao de f u n ç o e s
11 /∗ f u n ç a o f que r e t o r n a o v a l o r da f u n ç a o para um d et er mi n ad o ponto x dado
como argumento ∗/
13 double f ( double x )
{ return ( exp ( x∗x −1) + 10∗ s i n ( 2 ∗ x ) −5) ; }
15
/∗ f u n ç a o f d que r e t o r n a o v a l o r da p r i m e i r a d e r i v a d a de f para um
17 determi n ad o ponto x dado como argumento ∗/
double f d ( double x )
19
{ return ( 2 ∗ exp ( x∗x−1) ∗x+20 ∗ c o s ( 2 ∗ x ) ) ; }
21 /∗ f u n ç a o f d 2 que r e t o r n a o v a l o r da p r i m e i r a d e r i v a d a de f para um
determi n ad o ponto x dado como argumento ∗/
23 double f d 2 ( double x )
{ return ( 4 ∗ exp ( x∗x−1) ∗x∗x+2 ∗ exp ( x∗x−1)−40 ∗ s i n ( 2 ∗x ) ) ; }
25
/∗ f u n ç a o t r o c a v a l o r e s que t r o c a d o i s e l e m e n t o s de um v e c t o r , r e c e b e n d o
27 como argumentos p o n t e i r o s para os e l e m e n t o s do v e c t o r ∗/
void t r o c a v a l o r e s ( double ∗n1 , double ∗n2 )
29 {
double i 1 ;
31
i 1 = ( ∗ n1 ) ;
33
( ∗ n1 ) = ( ∗ n2 ) ;
( ∗ n2 ) = i 1 ;
35 }
37 /∗ f u n ç a o q u i c k s o r t que ordena por ordem c r e c e n t e os e l e m e n t o s de um
v e c t o r r1 dado como argumento e os e l e m e n t o s i n i c i a l e f i n a l do v e c t o r ∗/
39 void q u i c k s o r t ( double ∗ r 1
, int e s q , int d i r )
{
41
int i , j , m ;
double rm ;
43
i = esq ;
45
j = dir ;
m = ( esq + d i r ) / 2 ;
47
rm = r 1 [m] ;
while ( i <= j )
49
{
while ( r 1 [ i ] < rm) ++i ;
51
while (rm < r 1 [ j ] ) −−j ;
i f ( i <= j )
53
{
t r o c a v a l o r e s (& r 1 [ i ] , &r 1 [ j ] ) ;
55
++i ;
−−j ;
57
}
}
59
i f ( e s q < j ) q u i c k s o r t ( r1 , esq , j ) ;
i f ( i < d i r ) q u i c k s o r t ( r1 , i , d i r ) ;
61 }
10
63 //−−−−−−−−−−−−−−−−−−−−−−−−−−−FUNÇAO PRINCIPAL MAIN−−−−−−−−−−−−−−−−−−−−−−−−−−−
int main ( )
65
{
// i n i c i a l i z a ç a o de v a r i a v e i s u t i l i z a d a s no programa
67
/∗NOTA: ’ d o u b l e ’ i n d i c a que a v a r i a v e i s s u p o r t a p r e c i s a o d u p l a
(6 a l g a r i s m o s s i g n i f i c a t i v o s / 6 c a s a d e c i m a i s ;
69
’ i n t ’ i n d i c a que a v a r i a v e l é um i n t e i r o ∗/
double c =0, f d c =0, x m , z , k , s s , a m , b m ;
71
int m;
73
75
77
79
// v a l o r dos z e r o s da f u n ç ã o f c o n h e c i d o s g e o m e t r i c a m e n t e
double z1= −1.590190783596;
double z2= 0 . 2 3 9 5 8 2 6 5 7 6 2 6 3 9 4 ;
double z3= 1 . 4 8 9 6 3 0 0 2 1 5 1 4 4 ;
double z4= 1 . 8 0 1 9 8 7 2 8 7 9 5 2 5 2 ;
// i n i c i a l i z a ç a o de v e c t o r e s de 1000 e l e m e n t o s cada
double g d [ 1 0 0 0 ] , x [ 1 0 0 0 ] , e m [ 1 0 0 0 ] , ex m [ 1 0 0 0 ] ;
81
83
/∗NOTA: O comando ’ c o u t ’ imprime no ecran ; o comando ’ e n d l ’ l im p a
o b u f f e r e c o l o c a p a s s a a proxima l i n h a ∗/
c o u t << ”\tMETODO DE STEFFENSEN” << e n d l << e n d l ;
85
87
89
// l e r os v a l o r e s a m e b m do t e c l a d o
c o u t << ” I n s i r a o extremo ’ a ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >> a m ;
/∗NOTA: ’ c i n ’ p e r m i t e r e c e b e r a informaçao do t e c l a d o e g u a r d a r numa
v a r i a v e l j a i n i c i a l i z a d a ∗/
91
93
95
97
99
101
103
105
107
109
111
113
c o u t << ” I n s i r a o extremo ’ b ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >> b m ;
/∗NOTA: ’ i f ’ é um comando c o n d i c i o n a l que e f e c t u a a i n s t r u ç a o s e a c o n d i ç a o
s e f o r v e r d a d e i r a ∗/
// v e r i f i c a q u a n t o s z e r o s tem no i n t e r v a l o
i f ( z1<b m && z1>a m ) { z = z1 ;
f d c ++;}
i f ( z2<b m && z2>a m ) { z = z2 ; f d c ++;}
i f ( z3<b m && z3>a m ) { z = z3 ; f d c ++;}
i f ( z4<b m && z4>a m ) { z = z4 ; f d c ++;}
// s e e x i s t i r e m mais que um z er o , ou f ( a m ) ∗ f ( b m )>0 imprime e s a i
// s e nao c o n t i n u a a c o r r e r o programa
i f ( f ( a m ) ∗ f ( b m )>0 | | f d c >1)
{
c o u t << e n d l ;
c o u t << ”ATENCAO! Nao s e pode a p l i c a r o metodo n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” Causas p r o v a v e i s : ” << e n d l ;
c o u t << ” −ha mais que um z e r o n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” −nao ha z e r o s n e s t e i n t e r v a l o ” << e n d l ;
system ( ”PAUSE” ) ; // i n s t r u ç ã o para p a u s ar o s i s t e m a
return 0 ; // i n s t r u ç a o de s a i d a do programa
}
115
double x b i s 1 , x b i s=b m ; // i n i c i a l i z a ç a o de v a r i a v e i s do t i p o ’ d o u b l e ’
117
119
121
//CICLO1 : a p l i c a
// v a i de b i s =0 a
f o r ( int b i s =0;
{
xbis1 = (a
o metodo da b i s s e c ç a o ao i n t e r v a l o [ a m , b m ]
b i s <1000 incrementando a v a r i a v e l b i s
b i s < 1 0 0 0 ; b i s ++)
m+b m ) / 2 ; // f o r m u l a i t e r a t i v a do metodo
123
125
i f ( f ( a m ) ∗ f ( x b i s 1 ) <0) b m=x b i s 1 ;
e l s e i f ( f ( b m ) ∗ f ( x b i s 1 ) <0) a m=x b i s 1 ;
11
i f ( abs ( b m−a m ) <0.01) { x [ 0 ] = b m ; break ; }
// c a s o | b m−a m | <0.01 e n t a o guarda em x [ 0 ] o v a l o r de b m e s a i
/∗NOTA: a i n s t r u ç a o ’ b r e a k ’ p r o v o c a a s a i d a do c i c l o ∗/
127
129
}
131
133
135
137
139
141
143
// imprime o novo i n t e r v a l o I ’ e d i z o v a l o r e s c o l h i d o para x 0
c o u t << e n d l << ” Foi e n c o n t r a d o o i n t e r v a l o [ ” << a m << ” , ” << b m << ” ] . ” ;
c o u t << e n d l << ” Escolheu −s e x 0= ” << x [ 0 ] ;
c o u t << e n d l << e n d l ;
// c a l c u l a os e r r o s i n i c i a i s
e m [0]= z − x [ 0 ] ;
ex m [ 0 ] = z − x [ 0 ] ;
// imprime o c a b e ç a l h o da t a b e l a e a p r i m e i r a l i n h a d e s t a
c o u t << ” I t . x m \ te m1 \ t \ te m2 \ t
| e m /{ e {m−1}}ˆ2| \ t g ’ ” << e n d l ;
c o u t << ” 0
” << x [ 0 ] << ” \ t ” << e m [ 0 ] << ” \ t ” << ex m [ 0 ]
<< ”\ t−−−
\t
−−−” << e n d l ;
145
147
149
151
153
155
/∗CICLO2 : a p l i c a o metodo de s t e f f e n s e n ,
v a i de i =1 a t e i <1000 incrementando i ∗/
f o r ( int i =1; i <1000; i ++)
{
// c a l c u l a r x {m+1} p e l o metodo de s t e f f e n s e n
x [ i ] = x [ i −1] − f ( x [ i −1]) ∗ f ( x [ i −1]) / ( f ( x [ i −1]+ f ( x [ i −1]) )−f ( x [ i −1]) ) ;
s s = ( x [ i ]−x [ i −1]) / ( x [ i −1]−x [ i −2]) ; // c a l c u l a r a d e r i v a d a g ’
g d [ i ] = abs ( s s ) ;
// c a l c u l a r | g ’ |
ex m [ i ] = z − x [ i − 1 ] ;
// c a l c u l a r e r r o
e m [ i ] = x [ i ] − x [ i −1];
// c a l c u l a r e r r o e s t i m a d o
157
// c a l c u l a r a s u c e s s a o que t e n d e para o c o e f i c i e n t e a s s i n t o t i c o
k = abs ( e m [ i ] / ( e m [ i −1]∗ e m [ i −1]) ) ;
159
// s e g ( x m ) nao p e r t e n c e ao i n t e r v a l o , c=1
i f ( x [ i ]<a m && x [ i ]>b m ) c =1;
161
163
// i m p r e s s a o da t a b e l a
c o u t << i << ”
” << x [ i ] << ” \ t ” << e m [ i ] << ” \ t ”
<< ex m [ i ] << ” \ t ” << k << ”
\t” ;
i f ( i >= 2 ) c o u t << s s <<e n d l ;
e l s e c o u t << ”−−−” << e n d l ;
165
167
169
// c o n d i ç a o de paragem , s e o e r r o e s t i m a d o f o r menor que 1E−6
i f ( i >1 && abs ( e m [ i ] )< 1E−6)
{
m=i ;
// guarda a i t e r a ç a o em m
x m = x[ i ];
// guarda o v a l o r de x m
break ;
// s a i do c i c l o
}
171
173
175
}
177
179
181
// c a s o nao c o n v i r j a a p r e s e n t a e s s a informaçao
i f ( abs ( x [999] − x [ 9 9 8 ] ) > 1E−6) c o u t << ”O metodo nao c o n v e r g i u ” << e n d l ;
// c r i a um novo v e c t o r de m+1 e l e m e n t o s
double ∗ cgd= new double [m+ 1 ] ;
183
185
f o r ( int j =0; j <m+1; j ++)
cgd [ j ]= g d [ j ] ;
187
q u i c k s o r t ( cgd , 0 , m) ;
189
// c i c l o que p r e e n c h e o novo v e c t o r
// u t i l i z a q u i c k s o r t para o r d e n a r o novo v e c t o r
// imprimir o v a l o r do c o e f i c i e n t e a s s i n t o t i c o , c a l c u l a d o com r e c u r s o
12
// a segunda d e r i v a d a e o v a l o r do z e r o e n c o n t r a d o
c o u t << e n d l << ”K {\\ i n f t y } = ” << abs ( f d 2 ( x m ) / 2 ) << e n d l ;
c o u t << e n d l << ” Foi e n c o n t r a d o o z e r o ” << x m << ” . ” << e n d l << e n d l ;
191
193
i f ( c ) // s e c é d i f e r e n t e de ze r o , imprime
{
c o u t << e n d l << ”ATENCAO! A c o n d i ç a o de que g ( I ) nao p e r t e n c e a I nao s e
verifica”
<< e n d l ;
c o u t << ”O método pode não c o n v e r g i r ! ” << e n d l ;
}
195
197
199
201
i f ( cgd [m] >1) // s e o u l t i m o e l e m e n t o de cgd é maior que 1 , g nao e c o n t r a c t i v a
c o u t << ”ATENCAO! A f u n c a o g nao e c o n t r a c t i v a , o metodo pode nao c o n v e r g i r ”
<< e n d l ;
203
205
207
system ( ”PAUSE” ) ; // pausa o s i s t e m a
return 0 ;
// s a i do programa
}
3.5.2
steffensen sqrt2.cpp
1 /∗PROG QUE APLICA O METODO DE STEFFENSEN∗/
/∗IST , MEFT, MC 2010 ∗/
3 /∗JOAO CUNHA 67921 , NATALIA MORAIS 67936 , MIGUEL PRATA 67933 ∗/
5 #include <i o s t r e a m >
#include <cmath>
7 #include <c s t d l i b >
// h e a d e r que p o s s u i f u n ç o e s de i n p u t / o u t p u t
// h e a d e r que p o s s u i f u n ç o e s s t a n d a r t
// h e a d e r com f u n ç o e s matematicas
9 using namespace s t d ; // namespace com informaçao de f u n ç o e s
11 /∗ f u n ç a o f que r e t o r n a o v a l o r da f u n ç a o para um d et er mi n ad o ponto x dado
como argumento ∗/
13 double f ( double x )
{ return ( x∗x −2) ; }
15
/∗ f u n ç a o f d que r e t o r n a o v a l o r da p r i m e i r a d e r i v a d a de f para um
17 determi n ad o ponto x dado como argumento ∗/
double f d ( double x )
19
{ return 2∗ x ; }
21 /∗ f u n ç a o f d 2 que r e t o r n a o v a l o r da p r i m e i r a d e r i v a d a de f para um
determi n ad o ponto x dado como argumento ∗/
23 double f d 2 ( double x )
{ return 2 ; }
25
/∗ f u n ç a o t r o c a v a l o r e s que t r o c a d o i s e l e m e n t o s de um v e c t o r , r e c e b e n d o
27 como argumentos p o n t e i r o s para os e l e m e n t o s do v e c t o r ∗/
void t r o c a v a l o r e s ( double ∗n1 , double ∗n2 )
29 {
double i 1 ;
31
i 1 = ( ∗ n1 ) ;
33
( ∗ n1 ) = ( ∗ n2 ) ;
( ∗ n2 ) = i 1 ;
35 }
37 /∗ f u n ç a o q u i c k s o r t que ordena por ordem c r e c e n t e os e l e m e n t o s de um
v e c t o r r1 dado como argumento e os e l e m e n t o s i n i c i a l e f i n a l do v e c t o r ∗/
39 void q u i c k s o r t ( double ∗ r 1
, int e s q , int d i r )
{
41
int i , j , m ;
double rm ;
43
13
45
47
49
51
53
55
57
59
i = esq ;
j = dir ;
m = ( esq + d i r ) / 2 ;
rm = r 1 [m] ;
while ( i <= j )
{
while ( r 1 [ i ] < rm) ++i ;
while (rm < r 1 [ j ] ) −−j ;
i f ( i <= j )
{
t r o c a v a l o r e s (& r 1 [ i ] , &r 1 [ j ] ) ;
++i ;
−−j ;
}
}
i f ( e s q < j ) q u i c k s o r t ( r1 , esq , j ) ;
i f ( i < d i r ) q u i c k s o r t ( r1 , i , d i r ) ;
61 }
//−−−−−−−−−−−−−−−−−−−−−−−−−−−FUNÇAO PRINCIPAL MAIN−−−−−−−−−−−−−−−−−−−−−−−−−−−
63 int main ( )
{
65
// i n i c i a l i z a ç a o de v a r i a v e i s u t i l i z a d a s no programa
/∗NOTA: ’ d o u b l e ’ i n d i c a que a v a r i a v e i s s u p o r t a p r e c i s a o d u p l a
67
(6 a l g a r i s m o s s i g n i f i c a t i v o s / 6 c a s a d e c i m a i s ;
’ i n t ’ i n d i c a que a v a r i a v e l é um i n t e i r o ∗/
69
double c =0, f d c =0, x m , z , k , s s , a m , b m ;
int m;
71
// v a l o r dos z e r o s da f u n ç ã o f c o n h e c i d o s g e o m e t r i c a m e n t e
73
double z1= −1.41421;
double z2= 1 . 4 1 4 2 1 ;
75
// i n i c i a l i z a ç a o de v e c t o r e s de 1000 e l e m e n t o s cada
77
double g d [ 1 0 0 0 ] , x [ 1 0 0 0 ] , e m [ 1 0 0 0 ] , ex m [ 1 0 0 0 ] ;
79
81
83
85
87
89
/∗NOTA: O comando ’ c o u t ’ imprime no ecran ; o comando ’ e n d l ’ l im p a
o b u f f e r e c o l o c a p a s s a a proxima l i n h a ∗/
c o u t << ”\tMETODO DE STEFFENSEN” << e n d l << e n d l ;
// l e r os v a l o r e s a m e b m do t e c l a d o
c o u t << ” I n s i r a o extremo ’ a ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >> a m ;
/∗NOTA: ’ c i n ’ p e r m i t e r e c e b e r a informaçao do t e c l a d o e g u a r d a r numa
v a r i a v e l j a i n i c i a l i z a d a ∗/
c o u t << ” I n s i r a o extremo ’ b ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >> b m ;
91
93
95
/∗NOTA: ’ i f ’ é um comando c o n d i c i o n a l que e f e c t u a a i n s t r u ç a o s e a c o n d i ç a o
s e f o r v e r d a d e i r a ∗/
// v e r i f i c a q u a n t o s z e r o s tem no i n t e r v a l o
i f ( z1<b m && z1>a m ) { z = z1 ;
f d c ++;}
i f ( z2<b m && z2>a m ) { z = z2 ; f d c ++;}
97
99
101
103
105
107
// s e e x i s t i r e m mais que um z er o , ou f ( a m ) ∗ f ( b m )>0 imprime e s a i
// s e nao c o n t i n u a a c o r r e r o programa
i f ( f ( a m ) ∗ f ( b m )>0 | | f d c >1)
{
c o u t << e n d l ;
c o u t << ”ATENCAO! Nao s e pode a p l i c a r o metodo n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” Causas p r o v a v e i s : ” << e n d l ;
c o u t << ” −ha mais que um z e r o n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” −nao ha z e r o s n e s t e i n t e r v a l o ” << e n d l ;
system ( ”PAUSE” ) ; // i n s t r u ç ã o para p a u s ar o s i s t e m a
14
109
return 0 ; // i n s t r u ç a o de s a i d a do programa
}
111
double x b i s 1 , x b i s=b m ; // i n i c i a l i z a ç a o de v a r i a v e i s do t i p o ’ d o u b l e ’
113
//CICLO1 : a p l i c a
// v a i de b i s =0 a
f o r ( int b i s =0;
{
xbis1 = (a
115
117
o metodo da b i s s e c ç a o ao i n t e r v a l o [ a m , b m ]
b i s <1000 incrementando a v a r i a v e l b i s
b i s < 1 0 0 0 ; b i s ++)
m+b m ) / 2 ; // f o r m u l a i t e r a t i v a do metodo
i f ( f ( a m ) ∗ f ( x b i s 1 ) <0) b m=x b i s 1 ;
e l s e i f ( f ( b m ) ∗ f ( x b i s 1 ) <0) a m=x b i s 1 ;
119
121
i f ( abs ( b m−a m ) <0.01) { x [ 0 ] = b m ; break ; }
// c a s o | b m−a m | <0.01 e n t a o guarda em x [ 0 ] o v a l o r de b m e s a i
/∗NOTA: a i n s t r u ç a o ’ b r e a k ’ p r o v o c a a s a i d a do c i c l o ∗/
123
125
127
129
}
// imprime o novo i n t e r v a l o I ’ e d i z o v a l o r e s c o l h i d o para x 0
c o u t << e n d l << ” Foi e n c o n t r a d o o i n t e r v a l o [ ” << a m << ” , ” << b m << ” ] . ” ;
c o u t << e n d l << ” Escolheu −s e x 0= ” << x [ 0 ] ;
c o u t << e n d l << e n d l ;
131
133
// c a l c u l a os e r r o s i n i c i a i s
e m [0]= z − x [ 0 ] ;
ex m [ 0 ] = z − x [ 0 ] ;
135
137
139
141
143
145
147
149
// imprime o c a b e ç a l h o da t a b e l a e a p r i m e i r a l i n h a d e s t a
c o u t << ” I t . x m \ te m1 \ t \ te m2 \ t
| e m /{ e {m−1}}ˆ2| \ t g ’ ” << e n d l ;
c o u t << ” 0
” << x [ 0 ] << ” \ t ” << e m [ 0 ] << ” \ t ” << ex m [ 0 ]
<< ”\ t−−−
\t
−−−” << e n d l ;
/∗CICLO2 : a p l i c a o metodo de s t e f f e n s e n ,
v a i de i =1 a t e i <1000 incrementando i ∗/
f o r ( int i =1; i <1000; i ++)
{
// c a l c u l a r x {m+1} p e l o metodo de s t e f f e n s e n
x [ i ] = x [ i −1] − f ( x [ i −1]) ∗ f ( x [ i −1]) / ( f ( x [ i −1]+ f ( x [ i −1]) )−f ( x [ i −1]) ) ;
s s = ( x [ i ]−x [ i −1]) / ( x [ i −1]−x [ i −2]) ; // c a l c u l a r a d e r i v a d a g ’
g d [ i ] = abs ( s s ) ;
// c a l c u l a r | g ’ |
ex m [ i ] = z − x [ i − 1 ] ;
// c a l c u l a r e r r o
e m [ i ] = x [ i ] − x [ i −1];
// c a l c u l a r e r r o e s t i m a d o
151
// c a l c u l a r a s u c e s s a o que t e n d e para o c o e f i c i e n t e a s s i n t o t i c o
k = abs ( e m [ i ] / ( e m [ i −1]∗ e m [ i −1]) ) ;
153
155
// s e g ( x m ) nao p e r t e n c e ao i n t e r v a l o , c=1
i f ( x [ i ]<a m && x [ i ]>b m ) c =1;
157
// i m p r e s s a o da t a b e l a
c o u t << i << ”
” << x [ i ] << ” \ t ” << e m [ i ] << ” \ t ”
<< ex m [ i ] << ” \ t ” << k << ”
\t” ;
i f ( i >= 2 ) c o u t << s s <<e n d l ;
e l s e c o u t << ”−−−” << e n d l ;
159
161
163
// c o n d i ç a o de paragem , s e o e r r o e s t i m a d o f o r menor que 1E−6
i f ( i >1 && abs ( e m [ i ] )< 1E−6)
{
m=i ;
// guarda a i t e r a ç a o em m
x m = x[ i ];
// guarda o v a l o r de x m
break ;
// s a i do c i c l o
}
165
167
169
171
}
15
173
// c a s o nao c o n v i r j a a p r e s e n t a e s s a informaçao
i f ( abs ( x [999] − x [ 9 9 8 ] ) > 1E−6) c o u t << ”O metodo nao c o n v e r g i u ” << e n d l ;
175
// c r i a um novo v e c t o r de m+1 e l e m e n t o s
double ∗ cgd= new double [m+ 1 ] ;
177
179
f o r ( int j =0; j <m+1; j ++)
cgd [ j ]= g d [ j ] ;
// c i c l o que p r e e n c h e o novo v e c t o r
181
q u i c k s o r t ( cgd , 0 , m) ;
// u t i l i z a q u i c k s o r t para o r d e n a r o novo v e c t o r
183
// imprimir o v a l o r do c o e f i c i e n t e
// a segunda d e r i v a d a e o v a l o r do
c o u t << e n d l << ”K {\\ i n f t y } = ”
c o u t << e n d l << ” Foi e n c o n t r a d o o
185
187
189
a s s i n t o t i c o , c a l c u l a d o com r e c u r s o
zero encontrado
<< abs ( f d 2 ( x m ) / 2 ) << e n d l ;
z e r o ” << x m << ” . ” << e n d l << e n d l ;
i f ( c ) // s e c é d i f e r e n t e de ze r o , imprime
{
c o u t << e n d l << ”ATENCAO! A c o n d i ç a o de que g ( I ) nao p e r t e n c e a I nao s e
verifica”
<< e n d l ;
c o u t << ”O método pode não c o n v e r g i r ! ” << e n d l ;
}
191
193
195
i f ( cgd [m] >1) // s e o u l t i m o e l e m e n t o de cgd é maior que 1 , g nao e c o n t r a c t i v a
c o u t << ”ATENCAO! A f u n c a o g nao e c o n t r a c t i v a , o metodo pode nao c o n v e r g i r ”
<< e n d l ;
197
199
system ( ”PAUSE” ) ; // pausa o s i s t e m a
return 0 ;
// s a i do programa
201
}
3.5.3
steffensen aplicacao
/∗PROG QUE APLICA O METODO DE STEFFENSEN∗/
2 /∗IST , MEFT, MC 2010 ∗/
/∗JOAO CUNHA 67921 , NATALIA MORAIS 67936 , MIGUEL PRATA 67933 ∗/
4
#include <i o s t r e a m > // h e a d e r que p o s s u i f u n ç o e s de i n p u t / o u t p u t
6 #include <f s t r e a m >
// h e a d e r com f u n c o e s para e s c r i t a e l e i t u r a de f i c h e i r o s
#include <cmath>
// h e a d e r com f u n ç o e s matematicas
8
using namespace s t d ; // namespace com informaçao de f u n ç o e s
10
double h = 0 . 0 0 1 ; // v a r i a v e l g l o b a l , s u b s t i t u i t o d a s as o c o r r e n c i a s de h por 0 . 0 0 1
12
/∗ f u n c a o que r e t o r n a um v e c t o r com o v a l o r da p o s i c a o em R3 da p a r t i c u l a para
14 um i n s t a n t e t ∗/
double ∗ pos ( double t )
16
{
double ∗ r=new double [ 3 ] ; /∗ c r i a o e s p a ç o na memoria c o r r e s p o n d e n t e a um v e c t o r
18
t r i d i m e n s i o n a l ∗/
r [0]=(2+ c o s ( t / 2 ) ) ∗ c o s ( t ) ; // f u n c a o das p o s i c o e s no e i x o dos xx
20
r [1]=(2+ c o s ( t / 2 ) ) ∗ s i n ( t ) ; // f u n c a o das p o s i c o e s no e i x o dos yy
r [2]= s i n ( t /2) ;
// f u n c a o das p o s i c o e s no e i c o dos z z
22
return r ; // r e t o r n a o p o n t e i r o para o v e c t o r p o s i c a o
24
}
26 /∗ f u n c a o que r e t o r n a um v e c t o r com o v a l o r da v e l o c i d a d e em R3 da p a r t i c u l a para
um i n s t a n t e t ∗/
28 double ∗ v e l ( double t )
{
30
double ∗v=new double [ 3 ] ; /∗ c r i a o e s p a ç o na memoria c o r r e s p o n d e n t e a um v e c t o r
16
t r i d i m e n s i o n a l ∗/
v [0]=0 −1∗ s i n ( t / 2 ) ∗ c o s ( t ) /2− s i n ( t ) ∗ ( c o s ( t / 2 ) +2) ; // f u n c a o das v e l o c i d a d e s em x
v [ 1 ] = ( c o s ( t / 2 ) +2)∗ c o s ( t )−s i n ( t / 2 ) ∗ s i n ( t ) / 2 ;
// f u n c a o das v e l o c i d a d e s em y
v [2]= cos ( t /2) /2;
// f u n c a o das v e l o c i d a d e s em z
32
34
36
return v ; // r e t o r n a o p o n t e i r o para o v e c t o r v e l o c i d a d e
}
38
/∗ f u n c a o que r e t o r n a um v e c t o r com o v a l o r da d e r i v a d a da v e l o c i d a d e em R3 da
40 p a r t i c u l a para um i n s t a n t e t ∗/
double ∗ d v e l ( double t )
42
{
double ∗dv=new double [ 3 ] ; /∗ c r i a o e s p a ç o na memoria c o r r e s p o n d e n t e a um v e c t o r
44
t r i d i m e n s i o n a l ∗/
dv [ 0 ] = s i n ( t / 2 ) ∗ s i n ( t )−c o s ( t / 2 ) ∗ c o s ( t ) /4−( c o s ( t / 2 ) +2)∗ c o s ( t ) ;
// f u n c a o em x
46
dv[1]=0 − s i n ( t / 2 ) ∗ c o s ( t )−s i n ( t ) ∗ c o s ( t / 2 ) /4−( c o s ( t / 2 ) +2)∗ s i n ( t ) ; // f u n c a o em y
dv[2]= − s i n ( t / 2 ) / 4 ;
// f u n c a o em z
48
return dv ; // r e t o r n a o p o n t e i r o para o v e c t o r d e r i v a d a da v e l o c i d a d e
50
}
52 /∗ f u n c a o que r e t o r n a um v e c t o r com o v a l o r da segunda d e r i v a d a da v e l o c i d a d e em
R3 da p a r t i c u l a para um i n s t a n t e t ∗/
54 double ∗ d d v e l ( double t )
{
56
double ∗ddv=new double [ 3 ] ; /∗ c r i a o e s p a ç o na memoria c o r r e s p o n d e n t e a um v e c t o r
t r i d i m e n s i o n a l ∗/
58
ddv [ 0 ] = 1 3 ∗ s i n ( t / 2 ) ∗ c o s ( t ) /8+3∗ s i n ( t ) ∗ c o s ( t / 2 ) /4+( c o s ( t / 2 ) +2)∗ s i n ( t ) ; // f u n c a o em x
ddv [ 1 ] = 1 3 ∗ s i n ( t / 2 ) ∗ s i n ( t ) /8−3∗ c o s ( t / 2 ) ∗ c o s ( t ) /4−( c o s ( t / 2 ) +2)∗ c o s ( t ) ; // f u n c a o em y
60
ddv[2]= − c o s ( t / 2 ) / 8 ;
// f u n c a o em z
62
return ddv ; // r e t o r n a o p o n t e i r o para o v e c t o r segunda d e r i v a d a da v e l o c i d a d e
}
64
/∗ f u n ç a o t r o c a v a l o r e s que t r o c a d o i s e l e m e n t o s de um v e c t o r , r e c e b e n d o
66 como argumentos p o n t e i r o s para os e l e m e n t o s do v e c t o r ∗/
void t r o c a v a l o r e s ( double ∗n1 , double ∗n2 )
68
{
double i 1 ;
70
i 1 = ( ∗ n1 ) ;
72
( ∗ n1 ) = ( ∗ n2 ) ;
( ∗ n2 ) = i 1 ;
74
}
76 /∗ f u n ç a o q u i c k s o r t que ordena por ordem c r e c e n t e os e l e m e n t o s de um
v e c t o r r1 dado como argumento e os e l e m e n t o s i n i c i a l e f i n a l do v e c t o r ∗/
78 void q u i c k s o r t ( double ∗ r 1
, int e s q , int d i r )
{
80
int i , j , m ;
double rm ;
82
i = esq ;
84
j = dir ;
m = ( esq + d i r ) / 2 ;
86
rm = r 1 [m] ;
while ( i <= j )
88
{
while ( r 1 [ i ] < rm) ++i ;
90
while (rm < r 1 [ j ] ) −−j ;
i f ( i <= j )
92
{
t r o c a v a l o r e s (& r 1 [ i ] , &r 1 [ j ] ) ;
94
++i ;
17
−−j ;
}
96
}
i f ( e s q < j ) q u i c k s o r t ( r1 , esq , j ) ;
i f ( i < d i r ) q u i c k s o r t ( r1 , i , d i r ) ;
98
100
}
102 /∗ f u n c a o que e s c r e v e os v a l o r e s , para cada i n s t a n t e t , da p o s i ç a o da p a r t i c u l a
em R3 num f i c h e i r o , a v e l o c i d a d e segundo x para cada i n s t a n t e t e o i n s t a n t e t
104 n o u t r o f i c h e i r o , e i g u a l m e n t e para as v e l o c i d a d e s segundo y e z ∗/
void g r a f ( double ∗ r , double ∗v )
106
{
o f s t r e a m p ( ” t r a j e c t o r i a . t x t ” ) ; // a b r e o f i c h e i r o t r a j e c t o r i a . t x t
108
o f s t r e a m vx ( ” vxx . t x t ” ) ;
// a b r e o f i c h e i r o v x x . t x t
o f s t r e a m vy ( ” vyy . t x t ” ) ;
// a b r e o f i c h e i r o vyy . t x t
110
o f s t r e a m vz ( ” vzz . t x t ” ) ;
// a b r e o f i c h e i r o v z z . t x t
double t =0, t f =4∗M PI ; // v a r i a v e i s que guardam o i n s t a n t e f i n a l e i n i c i a l
112
while ( t<t f ) // c i c l o que s e r e a l i z a enquanto o i n t a n t e a c t u a l f o r menor que o
114
{
// i n s t a n t e f i n a l
r=pos ( t ) ; // guarda no v e c t o r r a p o s i c a o no i n s t a n t e t
116
v=v e l ( t ) ; // guarda no v e c t o r r a v e l o c i d a d e no i n s t a n t e t
p<<r [0]<< ’ ’<<r [1]<< ’ ’<<r [2]<< e n d l ; /∗ e s c r e v e as p o s i c o e s x , y e z no
118
f i c h e i r o t r a j e c t o r i a . t x t ∗/
vx<<t<< ’ ’<<v[0]<< e n d l ; /∗ e s c r e v e o i n s t a n t e t e a v e l o c i d a d e segundo o
120
e i x o dos xx no f i c h e i r o v x x . t x t ∗/
vy<<t<< ’ ’<<v[1]<< e n d l ; /∗ e s c r e v e o i n s t a n t e t e a v e l o c i d a d e segundo o
122
e i x o dos yy no f i c h e i r o vyy . t x t ∗/
vz<<t<< ’ ’<<v[2]<< e n d l ; /∗ e s c r e v e o i n s t a n t e t e a v e l o c i d a d e segundo o
124
e i x o dos z z no f i c h e i r o v z z . t x t ∗/
t+=h ; // i n c r e m e n t a o i n s t a n t e de 0 . 0 0 1
126
}
// f e c h a os f i c h e i r o s
128
p . close () ;
vx . c l o s e ( ) ;
130
vy . c l o s e ( ) ;
vz . c l o s e ( ) ;
132
/∗ as l i n h a s a b a i x o contem o e n d e r e c o do e x e c u t a v e l do g n u p l o t , que s e r v i r a
134
para d e s e n h a r os g r a f i c o s p r e t e n d i d o s . Se s e p r e t e n d e r d e s e n h a r os g r a f i c o s
t e r a de s e g a r a n t i r que os f i c h e i r o s p l o t . t x t e p l o t 2 . t x t e s t a o na mesma
136
p a s t a que o programa e que o e n d e r e c o para o e x e c u t a v e l do g n u p l o t s e
e n c o n t r a c o r r e c t o ∗/
138
system ( ” \” c : \ \ g n u p l o t \\ b i n a r y \\ g n u p l o t . exe \” p l o t . t x t ” ) ;
system ( ” \” c : \ \ g n u p l o t \\ b i n a r y \\ g n u p l o t . exe \” p l o t 2 . t x t ” ) ;
140
}
142 // f u n c a o que a p l i c a o metodo de s t e f f e n s e n
void s t e f f e n s e n ( char e )
144
{
// i n i c i a l i z a ç a o de v a r i a v e i s u t i l i z a d a s na f u n c a o
146
double a m , b m , ∗ z1 , z , f d c =0, c =0, k , s s , x m ;
double x [ 1 0 0 0 ] , e m [ 1 0 0 0 ] , g d [ 1 0 0 0 ] , ex m [ 1 0 0 0 ] ;
148
int n , nz , m;
150
cout<<endl <<” P o s i c o e s ”<<e<<” maximas ou minimas ”<<endl <<e n d l ;
152
i f ( e==’ x ’ ) // s e a e s c o l h a do u t i l i z a d o r f o r x ( v e r a b a i x o f u n c a o main )
{
z1=new double [ 3 ] ; // a l o c a o e s p a c o n e c e s s a r i o para um v e c t o r t r i d i m e n s i o n a l
// guarda no v e c t o r os v a l o r e s dos z e r o s
z1 [ 0 ] = 2 . 9 1 0 9 4 3 1 ; z1 [ 1 ] = 2 ∗ M PI ; ; z1 [ 2 ] = 9 . 6 5 5 4 2 7 5 ;
/∗ v a r i a v e i s com u t i l i z a d e mais a f r e n t e na f u n c a o . indicam que as
i n f o r m a c o e s r e f e r e n t e s ao e i x o dos xx e s t a o g u a r d a d a s na p o s i c a o 0 dos
154
156
158
18
v e c t o r e s p o s i c a o , v e l o c i d a d e , p r i m e i r a e segunda d e r i v a d a da v e l o c i d a d e
e que o numero de z e r o s da f u n c a o v e l o c i d a d e no i n t e r v a l o de 0 a 4 p i é 3 ∗/
n=0; nz =3;
160
184
}
e l s e i f ( e==’ y ’ ) // s e a e s c o l h a do u t i l i z a d o r f o r y ( v e r a b a i x o f u n c a o main )
{
z1=new double [ 4 ] ; // a l o c a o e s p a c o n e c e s s a r i o para um v e c t o r t r i d i m e n s i o n a l
// guarda no v e c t o r os v a l o r e s dos z e r o s
z1 [ 0 ] = 1 . 4 5 0 6 8 4 8 ; z1 [ 1 ] = 4 . 4 3 4 2 2 3 8 ; z1 [ 2 ] = 8 . 1 3 2 1 4 6 8 ; z1 [ 3 ] = 1 1 . 1 1 5 6 8 6 ;
/∗ v a r i a v e i s com u t i l i z a d e mais a f r e n t e na f u n c a o . indicam que as
i n f o r m a c o e s r e f e r e n t e s ao e i x o dos yy e s t a o g u a r d a d a s na p o s i c a o 0 dos
v e c t o r e s p o s i c a o , v e l o c i d a d e , p r i m e i r a e segunda d e r i v a d a da v e l o c i d a d e
e que o numero de z e r o s da f u n c a o v e l o c i d a d e no i n t e r v a l o de 1 a 4 p i é 4 ∗/
n=1; nz =4;
}
e l s e i f ( e==’ z ’ ) // s e a e s c o l h a do u t i l i z a d o r f o r z ( v e r a b a i x o f u n c a o main )
{
z1=new double [ 2 ] ; // a l o c a o e s p a c o n e c e s s a r i o para um v e c t o r t r i d i m e n s i o n a l
// guarda no v e c t o r os v a l o r e s dos z e r o s
z1 [ 0 ] = M PI ; z1 [ 1 ] = 3 ∗ M PI ;
/∗ v a r i a v e i s com u t i l i z a d e mais a f r e n t e na f u n c a o . indicam que as
i n f o r m a c o e s r e f e r e n t e s ao e i x o dos z z e s t a o g u a r d a d a s na p o s i c a o 0 dos
v e c t o r e s p o s i c a o , v e l o c i d a d e , p r i m e i r a e segunda d e r i v a d a da v e l o c i d a d e
e que o numero de z e r o s da f u n c a o v e l o c i d a d e no i n t e r v a l o de 2 a 4 p i é 2 ∗/
n=2; nz =2;
}
186
do
162
164
166
168
170
172
174
176
178
180
182
188
190
192
194
196
198
200
202
204
{ // c i c l o que p e r m i t e o b t e r os l i m i t e s a e b do i n t e r v a l o I =[a , b ]
do
{ // c i c l o que p e r m i t e o b t e r os l i m i t e s a e b do i n t e r v a l o I =[a , b ]
cout<<”A f u n c a o das p o s i c o e s de ”<<e<<” tem ”<<nz<<” extremos ”<<e n d l ;
cout<<” I n s i r a o l i m i t e ’ a ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >>a m ; // l e o v a l o r de a
cout<<” I n s i r a o l i m i t e ’ b ’ do i n t e r v a l o I =[a , b ] : ” ;
c i n >>b m ; // l e o v a l o r de b
/∗ Se os v a l o r e s de a e b nao s e encontrarem no i n t e r v a l o de 0 a 4 p i
imprime mensagem de e r r o ∗/
i f ( a m<=0 | | b m<=0 | | a m>=4∗M PI | | b m>=4∗M PI )
cout<<” Erro : a f u n c a o s o e s t a d e f i n i d a em ] 0 , 4 p i [ ”<<endl <<e n d l ;
}
// Enquanto os v a l o r e s de a e b nao s e encontrarem no i n t e r v a l o de 0 a 4 p i
while ( a m<=0 | | b m<=0 | | a m>=4∗M PI | | b m>=4∗M PI ) ;
f o r ( int i =0; i <nz ; ++i ) // v e r i f i c a q u a n t o s z e r o s tem no i n t e r v a l o
{ i f ( z1 [ i ]<b m && z1 [ i ]>a m ) { z = z1 [ i ] ;
f d c ++;}}
206
208
210
212
214
216
i f ( v e l ( a m ) [ n ] ∗ v e l ( b m ) [ n]>0 | | f d c >1)
{ /∗ s e e x i s t i r mais que um z er o , ou f ( a ) ∗ f ( b )>0 imprime mensagem de
e r r o ∗/
c o u t << e n d l ;
c o u t << ”ATENCAO! Nao s e pode a p l i c a r o metodo n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” Causas p r o v a v e i s : ” << e n d l ;
c o u t << ” −ha mais que um extremo n e s t e i n t e r v a l o ” << e n d l ;
c o u t << ” −nao ha extremos n e s t e i n t e r v a l o ” << endl <<e n d l ;
}
} // enquanto e x i s t i r mais que um z e ro , ou f ( a ) ∗ f ( b )>0
while ( v e l ( a m ) [ n ] ∗ v e l ( b m ) [ n]>0 | | f d c >1) ;
218
double x b i s 1 , x b i s=b m ; // i n i c i a l i z a ç a o de v a r i a v e i s
220
222
//CICLO1 : a p l i c a o metodo da b i s s e c ç a o ao i n t e r v a l o [ a m , b m ]
// v a i de b i s =0 a b i s <1000 incrementando a v a r i a v e l b i s
19
224
for ( int b i s =0; b i s < 1 0 0 0 ; b i s ++)
{
x b i s 1 = ( a m+b m ) / 2 ; // f o r m u l a i t e r a t i v a do metodo
226
228
230
232
234
236
i f ( v e l ( a m ) [ n ] ∗ v e l ( x b i s 1 ) [ n] <0) b m=x b i s 1 ;
e l s e i f ( v e l ( b m ) [ n ] ∗ v e l ( x b i s 1 ) [ n] <0) a m=x b i s 1 ;
i f ( abs ( b m−a m ) <0.01) { x [ 0 ] = b m ; break ; }
// c a s o | b−a | <0.01 e n t a o guarda em x [ 0 ] o v a l o r de b e s a i
}
// imprime o novo i n t e r v a l o I ’ e d i z o v a l o r e s c o l h i d o para x 0
cout<<endl <<” Foi e n c o n t r a d o o i n t e r v a l o [ ”<<a m<<” , ”<<b m<<” ] . ” ;
c o u t << e n d l << ” Escolheu −s e x 0= ” << x [ 0 ] ;
cout<<endl <<e n d l ;
238
240
// c a l c u l a os e r r o s i n i c i a i s (1 e 2)
e m [ 0 ] = z−x [ 0 ] ;
ex m [ 0 ] = z − x [ 0 ] ;
242
244
246
248
250
252
254
256
258
// imprime o c a b e ç a l h o da t a b e l a e a p r i m e i r a l i n h a
c o u t << ” I t . x m \ te m1 \ t \ te m2 \ t
| e m /{ e {m−1}}ˆ2| \ t g ’ ” << e n d l ;
c o u t << ” 0
” << x [ 0 ] << ” \ t ” << e m [ 0 ] << ” \ t ” << ex m [ 0 ]
<< ”\ t−−−
\t
−−−” << e n d l ;
/∗CICLO2 : a p l i c a o metodo de s t e f f e n s e n , v a i de i =1 a t e i <1000 incrementando i ∗/
for ( int i =1; i <1000; i ++)
{
// c a l c u l a r x {m+1} p e l o metodo de s t e f f e n s e n
x [ i ]=x [ i −1]− v e l ( x [ i −1]) [ n ] ∗ v e l ( x [ i −1]) [ n ] / ( v e l ( x [ i −1]+ v e l ( x [ i −1]) [ n ] ) [ n]− v e l ( x [ i −1])
[n]) ;
s s =(x [ i ]−x [ i −1]) / ( x [ i −1]−x [ i −2]) ; // c a l c u l a r a d e r i v a d a g ’
g d [ i ]= abs ( s s ) ;
// c a l c u l a r | g ’ |
e m [ i ]=x [ i ]−x [ i − 1 ] ;
// c a l c u l a r e r r o
ex m [ i ] = z − x [ i − 1 ] ;
// c a l c u l a r e r r o e s t i m a d o
// c a l c u l a r a s u c e s s a o que t e n d e para o c o e f i c i e n t e a s s i n t o t i c o
k=abs ( e m [ i ] / ( e m [ i −1]∗ e m [ i −1]) ) ;
260
262
264
266
// s e g ( x m ) nao p e r t e n c e ao i n t e r v a l o , c=1
i f ( x [ i ]<a m && x [ i ]>b m ) c =1;
// i m p r e s s a o da t a b e l a
cout<<i <<”
”<<x [ i ]<<” \ t ”<<e m [ i ]<<” \ t ”<<ex m [ i ]<<” \ t ”<<k<<”
i f ( i >= 2 ) c o u t << s s <<e n d l ;
e l s e c o u t << ”−−−” << e n d l ;
268
270
272
274
276
278
// c o n d i ç a o de paragem , s e o e r r o e s t i m a d o f o r menor que 1E−6
i f ( i >1 && abs ( x [ i ]−x [ i −1])< 1E−6)
{
m=i ;
// guarda a i t e r a ç a o em m
x m = x [ i ] ; // guarda o v a l o r de x m
break ;
// s a i do c i c l o
}
}
// c a s o nao c o n v i r j a a p r e s e n t a e s s a informaçao
i f ( abs ( x [999] − x [ 9 9 8 ] ) >1E−6) cout<<”O metodo nao c o n v e r g i u ”<<e n d l ;
280
282
284
// c r i a um novo v e c t o r de m+1 e l e m e n t o s
double ∗ cgd= new double [m+ 1 ] ;
for ( int j =0; j <m+1; j ++) // c i c l o que p r e e n c h e o novo v e c t o r
cgd [ j ]= g d [ j ] ;
20
\t” ;
286
q u i c k s o r t ( cgd , 0 , m) ; // u t i l i z a q u i c k s o r t para o r d e n a r o novo v e c t o r
288
/∗ imprimir o v a l o r do c o e f i c i e n t e a s s i n t o t i c o , c a l c u l a d o com r e c u r s o
a segunda d e r i v a d a e o v a l o r do z e r o e n c o n t r a d o ∗/
c o u t << e n d l << ”K {\\ i n f t y } = ” << abs ( d d v e l ( x m ) [ n ] / 2 ) << e n d l ;
c o u t << e n d l << ” Foi e n c o n t r a d o o extremo ” << x m << ” que c o r r e s p o n d e a um ” ;
290
292
294
// s e o v a l o r da p o s i c a o f o r maior que z e r o e n t a o é um maximo
i f ( pos ( x m ) [ n] >0) cout<<”maximo” ;
// s e o v a l o r da p o s i c a o f o r menor que z e r o e n t a o é um minimo
i f ( pos ( x m ) [ n] <0) cout<<”minimo” ;
296
298
cout<<” de p o s i c a o ”<<endl <<e n d l ;
300
c o u t << e n d l ;
302
if (c)
{ // s e c é d i f e r e n t e de ze r o , imprime
c o u t << e n d l << ”ATENCAO! A c o n d i ç a o de que g ( I ) nao p e r t e n c e a I nao s e v e r i f i c a ”
<< e n d l ;
c o u t << ”O método pode não c o n v e r g i r ! ” << e n d l ;
}
304
306
308
i f ( g d [m] >1) // s e o u l t i m o e l e m e n t o de cgd é maior que 1 , g nao e c o n t r a c t i v a
c o u t << ”ATENCAO! A f u n c a o g nao e c o n t r a c t i v a , o metodo pode nao c o n v e r g i r ” << e n d l ;
310
}
312
//−−−−−−−−−−−−−−−−−−−−−−−−−−−FUNÇAO PRINCIPAL MAIN−−−−−−−−−−−−−−−−−−−−−−−−−−−
314 int main ( )
{
316
// i n i c i a l i z a c a o das v a r i a v e i s
double ∗ r , ∗v ;
318
char e ;
320
322
// Imprime um t i t u l o com informacao do programa
cout<<” \ t T r a j e c t o r i a de uma p a r t i c u l a em R3”<<endl <<e n d l ;
cout<<” r =((2+ c o s ( t / 2 ) ) ∗ c o s ( t ) ,(2+ c o s ( t / 2 ) ) ∗ s i n ( t ) , s i n ( t / 2 ) ) ”<<e n d l ;
cout<<” t em ] 0 , 4 p i [ ”<<endl <<e n d l ;
324
326
//chama a f u n c a o para d e s e n h a r os g r a f i c o s
graf (r , v) ;
328
do
330
332
334
{ /∗ c i c l o que p e r m i t e ao u t i l i z a d o r e s c o l h e r s e p r e t e n d e
e x t r e m o s das p o s i c o e s no e i x o dos xx , yy ou z z ∗/
cout<<” P ret en de : ”<<e n d l ;
cout<<” Enc ontrar o s extremos das p o s i c o e s no e i x o dos
cout<<” Enc ontrar o s extremos das p o s i c o e s no e i x o dos
cout<<” Enc ontrar o s extremos das p o s i c o e s no e i x o dos
c i n >>e ; // guarda a e s c o l h a na v a r i a v e l e ( x , y ou z )
e n c o n t r a r os
( x ) ’ s ”<<e n d l ;
( y ) ’ s ”<<e n d l ;
( z ) ’ s ”<<e n d l ;
336
338
340
342
344
// s e a e s c o l h a nao f o r v a l i d a imprime mensagem de e r r o
i f ( e != ’ x ’ && e != ’ y ’ && e != ’ z ’ )
cout<<” Erro : E s c o l h a ( x ) , ( y ) ou ( z ) ”<<e n d l ;
}
// enquanto a e s c o l h a nao f o r v a l i d a
while ( e != ’ x ’ && e != ’ y ’ && e != ’ z ’ ) ;
s t e f f e n s e n ( e ) ; /∗chama a f u n c a o que a p l i c a o metodo de s t e f f e n s e n , e n v i a n d o
para e s t a a v a r i a v e i que contem a e s c o l h a do u t i l i z a d o r ∗/
346
348
system ( ”PAUSE” ) ; // pausa o s i s t e m a
return 0 ;
// s a i do programa
21
}
3.5.4
newtongeneralizado.c
1 /∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
∗
TRABALHO DE MATEMATICA COMPUTACIONAL
∗
3 ∗
∗
∗
METODO DE NEWTON GENERALIZADO
∗
5 ∗
∗
∗
JOAO CUNHA, MIGUEL PRATA, NATALIA MORAIS
∗
7 ∗
∗
∗
IST − MEFT − 2010/2011
∗
9 ∗
∗
∗
∗
11 /∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗HEADERS∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗/
#include <s t d i o . h>
13 #include < s t d l i b . h>
#include <math . h>
15 #define max( a , b ) ( a > b ) ? a : b
/∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗FUNCOES EXTERNAS∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗/
17 f l o a t maxi ( f l o a t a , f l o a t b , f l o a t c )
{
19 f l o a t t 1 ;
i f ( a<b ) t 1=b ;
21 e l s e t 1=a ;
i f ( t1<c ) return c ;
23 e l s e return t 1 ;
}
25 // Funcao F − 1 a l i n h a do s i s t e m a
float funcf ( float x , float y , float z )
27 {
return ( x∗y − z ∗ z +3) ;
29 }
// Derivada p a r c i a l de F em ordem a x
31 f l o a t d x f ( f l o a t x , f l o a t y , f l o a t z )
{
33 return y ;
}
35 // Derivada p a r c i a l de F em ordem a y
float dyf ( float x , float y , float z )
37 {
return x ;
39 }
// Derivada p a r c i a l de F em ordem a z
41 f l o a t d z f ( f l o a t x , f l o a t y , f l o a t z )
{
43 return −2∗z ;
}
45 // Funcao G − 2 a l i n h a do s i s t e m a
float funcg ( float x , float y , float z )
47 {
return ( x∗y∗ z − x∗x + y∗y −2) ;
49 }
// Derivada p a r c i a l de G em ordem a x
51 f l o a t dxg ( f l o a t x , f l o a t y , f l o a t z )
{
53 return ( y∗ z − 2∗ x ) ;
}
55 // Derivada p a r c i a l de G em ordem a y
f l o a t dyg ( f l o a t x , f l o a t y , f l o a t z )
57 {
return ( x∗ z + 2∗ y ) ;
59 }
// Derivada p a r c i a l de G em ordem a z
61 f l o a t dzg ( f l o a t x , f l o a t y , f l o a t z )
22
{
63
65
67
69
71
73
75
77
79
81
83
85
87
89
return ( x∗y ) ;
}
// Funcao H − 3 a l i n h a do s i s t e m a
f l o a t funch ( f l o a t x , f l o a t y , f l o a t z )
{
return ( exp ( x ) − exp ( y ) + z − 3 ) ;
}
// Derivada p a r c i a l de H em ordem a x
f l o a t dxh ( f l o a t x , f l o a t y , f l o a t z )
{
return exp ( x ) ;
}
// Derivada p a r c i a l de H em ordem a y
f l o a t dyh ( f l o a t x , f l o a t y , f l o a t z )
{
return −exp ( y ) ;
}
// Derivada p a r c i a l de H em ordem a z
f l o a t dzh ( f l o a t x , f l o a t y , f l o a t z )
{
return 1 ;
}
// f u n c a o que f a z e l i m i n a c a o de g a u s s com p e s q u i s a p a r c i a l de p i v o t
void f o r w a r d e l i m ( int m a t r i x s i z e , f l o a t ∗∗A, int ∗L)
{
int i , j , k , tempi , tempk ;
f l o a t ∗S ;
f l o a t xmult , smax , rmax , r a t i o ;
91
S = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
93
95
97
99
101
103
105
107
109
111
113
115
117
119
121
for ( i = 1 ; i <= m a t r i x s i z e ; i ++)
{
L[ i ] = i ;
smax = 0 . 0 ;
fo r ( j = 1 ; j <= m a t r i x s i z e ; j ++)
smax = 1 ;
S [ i ] = smax ;
}
for ( k = 1 ; k < m a t r i x s i z e ; k++)
{
rmax = 0 . 0 ;
for ( i = k ; i <= m a t r i x s i z e ; i ++)
{
tempi = L [ i ] ;
r a t i o = f a b s (A[ tempi ] [ k ] / S [ tempi ] ) ;
i f ( r a t i o > rmax )
{
rmax = r a t i o ;
j = i;
}
}
tempk = L [ j ] ;
L[ j ] = L[ k ] ;
L [ k ] = tempk ;
123
for ( i = k+1; i <= m a t r i x s i z e ; i ++)
{
tempi = L [ i ] ;
125
xmult = A[ tempi ] [ k ] /A[ tempk ] [ k ] ;
23
A[ tempi ] [ k ] = xmult ;
fo r ( j = k+1; j<= m a t r i x s i z e ; j ++)
A[ tempi ] [ j ] −= xmult ∗ A[ tempk ] [ j ] ;
127
}
129
}
131 }
// f u n c a o que r e s o l v e o s i s t e m a de e q u a c o e s l i n e a r e s
133 void s o l v e ( int m a t r i x s i z e , f l o a t ∗∗A, int ∗L , f l o a t ∗B, f l o a t ∗X)
{
135
int i , j , k , tempi , tempk , tempn ;
137
f l o a t sum ;
139
for ( k = 1 ; k < m a t r i x s i z e ; k++)
{
tempk = L [ k ] ;
for ( i = k+1; i <= m a t r i x s i z e ; i ++)
{
tempi = L [ i ] ;
141
143
145
B [ tempi ] −= A[ tempi ] [ k ] ∗ B [ tempk ] ;
147
149
}
}
tempn = L [ m a t r i x s i z e ] ;
X[ m a t r i x s i z e ] = B [ tempn ] / A[ tempn ] [ m a t r i x s i z e ] ;
151
153
155
157
159
for ( i = m a t r i x s i z e −1; i>= 1 ; i −−)
{
tempi = L [ i ] ;
sum = B [ tempi ] ;
for ( j = i +1; j <= m a t r i x s i z e ; j ++)
sum −= A[ tempi ] [ j ] ∗ X[ j ] ;
X[ i ] = sum / A[ tempi ] [ i ] ;
}
}
161 // f u n c a o MAIN
void main ( )
163 {
// i n i c i a l i z a c a o de p o n t e i r o s
165 f l o a t ∗∗A, ∗B, ∗X;
int ∗L ;
167 // i n i c i a l i z a c a o de v a r i a v e i s
int i , j , m a t r i x s i z e , n ;
169 f l o a t x0 , y0 , z0 , x n , y n , z n , d e l t a x , d e l t a y , d e l t a z , e n ;
// d e f i n i c a o dos p o n t o s i n i c i a i s
171 x0=−21;
y0 =10;
173 z0 =2;
175 // d e f i n i c a o do tamanho das m a t r i z e s
matrix size = 3;
177 // a l o c a o de memoria para as m a t r i z e s
A = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ∗ ) ) ;
179 B = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
X = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
181 L = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
// c r i a c a o de m a t r i z e s
183
for ( i = 0 ; i <= m a t r i x s i z e ; i ++)
A[ i ] = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
185 // i m p r e s s a o de t a b e l a s c a b e c a l h o e i n t e r f a c e
p r i n t f ( ” \n\n\ t \ t \ t \ t \ t \ tMetodo de Newton G e n e r a l i z a d o \n\n” ) ;
187
p r i n t f ( ” | \ n | xy − z +3 = 0\n | xyz − x ˆ2 + y ˆ2 − 2 = 0\n | e ˆx − e ˆy + z − 3 = 0\n | \ n\n\n” )
;
24
189
191
193
195
197
199
201
203
205
207
209
211
213
215
217
219
221
223
225
227
p r i n t f ( ” I n s i r a o ponto i n i c i a l x0 ” ) ;
s c a n f ( ”%f ” ,&x0 ) ;
fflush ( stdin ) ;
p r i n t f ( ” \ n I n s i r a o ponto i n i c i a l y0 ” ) ;
s c a n f ( ”%f ” ,&y0 ) ;
fflush ( stdin ) ;
p r i n t f ( ” \ n I n s i r a o ponto i n i c i a l z0 ” ) ;
s c a n f ( ”%f ” ,& z0 ) ;
fflush ( stdin ) ;
p r i n t f ( ” \ t \ t \ t \ t \ t \ t \ t R e s u l t a d o s : \ n\n\n” ) ;
printf (” Iterada k
xˆk
yˆk
z ˆk
| | Xˆk − Xˆ ( k−1) | |
| | e ˆk | |
| | e ˆk | | / | | e ˆ ( k−1) | | \ n\n” ) ;
// metodo de newton
for ( n=0; n < 1 0 0 ; n++)
{
// guarda c o p i a s de x0 , y0 , z0 para comparacao com os p o n t o s s e g u i n t e s
x n=x0 ;
y n=y0 ;
z n=z0 ;
e n = maxi ( f a b s (1.29279 − x n ) , f a b s (0.998577 − y n ) , f a b s (2.07146 − z n ) ) ;
// p r e e n c h i m e n t o da m a t r i z J a c o b i a n a com as d e r i v a d a s p a r c i a i s
A [ 1 ] [ 1 ] = d x f ( x0 , y0 , z0 ) ; A [ 1 ] [ 2 ] = d y f ( x0 , y0 , z0 ) ; A [ 1 ] [ 3 ] = d z f ( x0 , y0 , z0 ) ;
A [ 2 ] [ 1 ] = dxg ( x0 , y0 , z0 ) ; A [ 2 ] [ 2 ] = dyg ( x0 , y0 , z0 ) ; A [ 2 ] [ 3 ] = dzg ( x0 , y0 , z0 ) ;
A [ 3 ] [ 1 ] = dxh ( x0 , y0 , z0 ) ; A [ 3 ] [ 2 ] = dyh ( x0 , y0 , z0 ) ; A [ 3 ] [ 3 ] = dzh ( x0 , y0 , z0 ) ;
// p r e e n c h i m e n t o do v e c t o r com as f u n c o e s
B [ 1 ] = −f u n c f ( x0 , y0 , z0 ) ;
B [ 2 ] = −f u n c g ( x0 , y0 , z0 ) ;
B [ 3 ] = −funch ( x0 , y0 , z0 ) ;
// a p l i c a as f u n c o e s de e l e m i n c a o de g a u s s e r e s o l u c a o do s i s t e m a l i n e a r
f o r w a r d e l i m ( m a t r i x s i z e , A, L) ;
s o l v e ( m a t r i x s i z e , A, L , B,X) ;
// guarda os r e s u l t a d o s
x0 = x0 + X [ 1 ] ;
y0 = y0 + X [ 2 ] ;
z0 = z0 + X [ 3 ] ;
// v e r i f i c a c a o da c o n v e r g e n c i a
i f ( f a b s (X [ 1 ] ) < 0 . 0 0 0 0 0 0 1
| | f a b s (X [ 2 ] ) < 0 . 0 0 0 0 0 0 1 | | f a b s (X [ 3 ] ) < 0 . 0 0 0 0 0 0 1 ) break
;
// i m p r e s s a o de r e s u l t a d o s
i f ( n==0) p r i n t f ( ” %d
%f
%f
%f
−
%f
−\n” , n , x0 , y0 , z0 , maxi ( f a b s (1.29279 − x0 ) , f a b s (0.998577 − y0 ) , f a b s
(2.07146 − z0 ) ) ) ;
i f ( n>0) p r i n t f ( ” %d
%f
%f
%f
%f
%f
%f \n” , n , x0 , y0 , z0 , maxi ( f a b s ( x0−x n ) , f a b s ( y0−y n ) , f a b s ( z0−z n ) ) , maxi (
f a b s (1.29279 − x0 ) , f a b s (0.998577 − y0 ) , f a b s (2.07146 − z0 ) ) , maxi ( f a b s (1.29279 − x0 ) , f a b s
(0.998577 − y0 ) , f a b s (2.07146 − z0 ) ) / e n ) ;
229
231 }
// l i b e r t a a memoria a l o c a d a
233
f r e e (A) ;
f r e e (X) ;
235
f r e e (B) ;
f r e e (L) ;
237 p r i n t f ( ” \ n V e r i f i c o u −s e c o n v e r g e n c i a apos %d i t e r a d a s \n\n” , n ) ;
i f ( n==99) ( ”Não s e v e r i f i c o u c o n v e r g e n c i a apos 100 i t e r a d a s , t e n t e o u t r a s c o n d i c o e s i n i c a i s
\n\n” ) ;
239
system ( ”PAUSE” ) ;
return 0 ;
241 }
3.5.5
newtongeneralizado2.c
1 /∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
25
∗
TRABALHO DE MATEMATICA COMPUTACIONAL
∗
3 ∗
∗
∗
METODO DE NEWTON GENERALIZADO
∗
5 ∗
∗
∗
JOAO CUNHA, MIGUEL PRATA, NATALIA MORAIS
∗
7 ∗
∗
∗
IST − MEFT − 2010/2011
∗
9 ∗
∗
∗
∗
11 /∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗HEADERS∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗/
#include <s t d i o . h>
13 #include < s t d l i b . h>
#include <math . h>
15 #define max( a , b ) ( a > b ) ? a : b
/∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗FUNCOES EXTERNAS∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗/
17 f l o a t maxi ( f l o a t a , f l o a t b )
{
19 f l o a t t 1 ;
i f ( a<b ) return b ;
21 e l s e return a ;
}
23 // Funcao F − 1 a l i n h a do s i s t e m a
float funcf ( float x , float y)
25 {
return ( 1 0 ∗ x∗x + s i n ( y ) −20) ;
27 }
// d e r i v a d a p a r c i a l de f em ordem a x
29 f l o a t d x f ( double x )
{
31 return 20∗ x ;
}
33 // d e r i v a d a p a r c i a l de f em ordem a y
float dyf ( float y )
35 {
return c o s ( y ) ;
37 }
// Funcao G − 2 a l i n h a do s i s t e m a
39 f l o a t f u n c g ( f l o a t x , f l o a t y )
{
41 return ( pow ( x , 4 ) + 5∗ y −6) ;
}
43 // d e r i v a d a p a r c i a l de g em ordem a x
f l o a t dxg ( f l o a t x )
45 {
return ( 4 ∗ pow ( x , 3 ) ) ;
47 }
// d e r i v a d a p a r c i a l de g em ordem a y
49 f l o a t dyg ( f l o a t y )
{
51 return 5 ;
}
53 // f u n c a o que f a z e l i m i n a c a o de g a u s s com p e s q u i s a p a r c i a l de p i v o t
void f o r w a r d e l i m ( int m a t r i x s i z e , f l o a t ∗∗A, int ∗L)
55 {
int i , j , k , tempi , tempk ;
57 f l o a t ∗S ;
f l o a t xmult , smax , rmax , r a t i o ;
59
S = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
61
for ( i = 1 ; i <= m a t r i x s i z e ; i ++)
63 {
L[ i ] = i ;
65
smax = 0 . 0 ;
26
fo r ( j = 1 ; j <= m a t r i x s i z e ; j ++)
smax = 1 ;
S [ i ] = smax ;
67
}
69
71
for ( k = 1 ; k < m a t r i x s i z e ; k++)
{
rmax = 0 . 0 ;
for ( i = k ; i <= m a t r i x s i z e ; i ++)
{
tempi = L [ i ] ;
r a t i o = f a b s (A[ tempi ] [ k ] / S [ tempi ] ) ;
i f ( r a t i o > rmax )
{
rmax = r a t i o ;
j = i;
}
}
73
75
77
79
81
83
85
tempk = L [ j ] ;
L[ j ] = L[ k ] ;
L [ k ] = tempk ;
87
89
for ( i = k+1; i <= m a t r i x s i z e ; i ++)
{
tempi = L [ i ] ;
91
93
xmult = A[ tempi ] [ k ] /A[ tempk ] [ k ] ;
A[ tempi ] [ k ] = xmult ;
fo r ( j = k+1; j<= m a t r i x s i z e ; j ++)
A[ tempi ] [ j ] −= xmult ∗ A[ tempk ] [ j ] ;
95
}
97
}
99
}
101 // f u n c a o que r e s o l v e o s i s t e m a de e q u a c o e s l i n e a r e s
void s o l v e ( int m a t r i x s i z e , f l o a t ∗∗A, int ∗L , f l o a t ∗B, f l o a t ∗X)
103 {
int i , j , k , tempi , tempk , tempn ;
105
f l o a t sum ;
107
for ( k = 1 ; k < m a t r i x s i z e ; k++)
109
{
tempk = L [ k ] ;
111
for ( i = k+1; i <= m a t r i x s i z e ; i ++)
{
113
tempi = L [ i ] ;
B [ tempi ] −= A[ tempi ] [ k ] ∗ B [ tempk ] ;
115
}
117
119
}
tempn = L [ m a t r i x s i z e ] ;
X[ m a t r i x s i z e ] = B [ tempn ] / A[ tempn ] [ m a t r i x s i z e ] ;
121
123
125
127
for ( i = m a t r i x s i z e −1; i>= 1 ; i −−)
{
tempi = L [ i ] ;
sum = B [ tempi ] ;
for ( j = i +1; j <= m a t r i x s i z e ; j ++)
sum −= A[ tempi ] [ j ] ∗ X[ j ] ;
X[ i ] = sum / A[ tempi ] [ i ] ;
}
129 }
27
// f u n c a o MAIN
131 void main ( )
{
133 // i n i c i a l i z a c a o de p o n t e i r o s
f l o a t ∗∗A, ∗B, ∗X;
135 int ∗L ;
// i n i c i a l i z a c a o de v a r i a v e i s
137 int i , j , m a t r i x s i z e , n ;
f l o a t x0 , y0 , x n , y n , d e l t a x , d e l t a y , e n ;
139 // d e f i n i c a o dos p o n t o s i n i c i a i s
x0 =1;
141 y0 =1;
// d e f i n i c a o do tamanho das m a t r i z e s
143 m a t r i x s i z e = 2 ;
// a l o c a o de memoria para as m a t r i z e s
145 A = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ∗ ) ) ;
B = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
147 X = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
L = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
149 // c r i a c a o de m a t r i z e s
for ( i = 0 ; i <= m a t r i x s i z e ; i ++)
151
A[ i ] = c a l l o c ( ( m a t r i x s i z e +1) , s i z e o f ( f l o a t ) ) ;
// i m p r e s s a o de t a b e l a s e c a b e c a l h o
153 p r i n t f ( ” \n\n\ t \ t \ t \ t \ t \ tMetodo de Newton G e n e r a l i z a d o \n\n” ) ;
p r i n t f ( ” | \ n | 1 0 x ˆ2 + s i n ( y ) − 20 = 0\n | x ˆ4 + 5y − 6 = 0\n | \ n\n\n\n” ) ;
155
p r i n t f ( ” I n s i r a o ponto i n i c i a l x0 ” ) ;
157
s c a n f ( ”%f ” ,&x0 ) ;
fflush ( stdin ) ;
159
p r i n t f ( ” \ n I n s i r a o ponto i n i c i a l y0 ” ) ;
161
s c a n f ( ”%f ” ,&y0 ) ;
fflush ( stdin ) ;
163
p r i n t f ( ” \ t \ t \ t \ t \ t \ t \ t R e s u l t a d o s : \ n\n\n” ) ;
165
printf (” Iterada k
xˆk
yˆk
| | Xˆk − Xˆ ( k−1) | |
| | eˆ
k||
| | e ˆk | | / | | e ˆ ( k−1) | | \ n\n” ) ;
// metodo de newton
167 for ( n=0; n < 1 0 0 ; n++)
{
169 // guarda c o p i a s de x0 e y0 para comparacao com os p o n t o s s e g u i n t e s
x n=x0 ;
171 y n=y0 ;
e n = maxi ( f a b s (1.39929 − x n ) , f a b s (0.433232 − y n ) ) ;
173 // p r e e n c h i m e n t o da m a t r i z J a c o b i a n a com as d e r i v a d a s p a r c i a i s
A [ 1 ] [ 1 ] = d x f ( x0 ) ; A [ 1 ] [ 2 ] = d y f ( y0 ) ;
175 A [ 2 ] [ 1 ] = dxg ( x0 ) ; A [ 2 ] [ 2 ] = dyg ( y0 ) ;
// p r e e n c h i m e n t o do v e c t o r com as f u n c o e s
177 B [ 1 ] = −f u n c f ( x0 , y0 ) ;
B [ 2 ] = −f u n c g ( x0 , y0 ) ;
179 // a p l i c a as f u n c o e s de e l e m i n c a o de g a u s s e r e s o l u c a o do s i s t e m a l i n e a r
f o r w a r d e l i m ( m a t r i x s i z e , A, L) ;
181
s o l v e ( m a t r i x s i z e , A, L , B,X) ;
// guarda os r e s u l t a d o s
183
x0 = x0 + X [ 1 ] ;
y0 = y0 + X [ 2 ] ;
185
// v e r i f i c a c a o da c o n v e r g e n c i a
i f ( f a b s (X [ 1 ] ) < 0 . 0 0 0 0 0 0 1
| | f a b s (X [ 2 ] ) < 0 . 0 0 0 0 0 0 1 ) break ;
187 // i m p r e s s a o de r e s u l t a d o s
i f ( n == 0 ) p r i n t f ( ” %d
%f
%f
−
%f
−\n” , n , x0 , y0 , maxi ( f a b s (1.39929 − x0 ) , f a b s (0.433232 − y0 ) ) ) ;
189 i f ( n > 0 ) p r i n t f ( ” %d
%f
%f
%f
%f
%f \n” , n , x0 , y0 , maxi ( f a b s ( x0−x n ) , f a b s ( y0−y n ) ) , maxi ( f a b s (1.39929 − x0 ) ,
f a b s (0.433232 − y0 ) ) , maxi ( f a b s (1.39929 − x0 ) , f a b s (0.433232 − y0 ) ) / e n ) ;
28
191 }
// l i b e r t a a memoria a l o c a d a
193
f r e e (A) ;
f r e e (X) ;
195
f r e e (B) ;
f r e e (L) ;
197 // i m p r e s s a o de r e s u l t a d o s
p r i n t f ( ” \nForam n e c e s s a r i a s %d i t e r a d a s a t e a t i n g i r a c o n v e r g e n c i a . \ n\n” , n ) ;
199
system ( ”PAUSE” ) ;
return 0 ;
201 }
29