Projeto Computacional de Métodos Numéricos

Prévia do material em texto

UNIVERSIDADE FEDERAL DO MARANHÃO 
CENTRO DE CIÊNCIAS EXATAS E TECNOLOGIAS 
MÉTODOS NUMÉRICOS E OTIMIZAÇÃO – ENGENHARIA ELÉTRICA 
 
 
 
AMANDA JADE AQUINO JANSEN 
 
 
 
 
 
 
 
 
SOLUÇÃO DE SISTEMAS NÃO-LINEARES VIA MÉTODOS DE NEWTON-
RAPHSON E QUASI-NEWTON USANDO MATLAB 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
São Luís - MA 
2021 
1 
 
 
AMANDA JADE AQUINO JANSEN 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
SOLUÇÃO DE SISTEMAS NÃO-LINEARES VIA MÉTODOS DE NEWTON-
RAPHSON E QUASI-NEWTON USANDO MATLAB 
 
 
 
 
 
 
Projeto apresentado como 
parte dos requisitos necessários 
para a obtenção de nota na 
disciplina Métodos Numéricos e 
Otimização, ministrada pelo 
professor Anselmo Barbosa 
Rodrigues. 
 
 
 
 
 
São Luís - MA 
2021 
 
 
2 
 
 
SUMÁRIO 
 
 
1. INTRODUÇÃO ................................................................................................................ 3 
2. OBJETIVOS ..................................................................................................................... 3 
2.1 Objetivo Geral .............................................................................................................. 3 
2.2 Objetivos Específicos .................................................................................................. 4 
3. FUNDAMENTAÇÃO TEÓRICA ................................................................................... 4 
3.1 Sistemas não-lineares ................................................................................................... 4 
3.2 Método de Eliminação de Gauss ................................................................................. 5 
3.3 Método de Gauss - LU ................................................................................................. 8 
3.4 Método de Newton-Raphson ..................................................................................... 12 
3.5 Método de Quasi-Newton .......................................................................................... 13 
4. METODOLOGIA ........................................................................................................... 14 
4.1 Especificações ............................................................................................................ 14 
4.2 Descrição do Sistema Teste ....................................................................................... 14 
4.3 Problemas ................................................................................................................... 15 
5. RESULTADOS E DISCUSSÃO ................................................................................... 15 
5.1 Problema I .................................................................................................................. 15 
5.2 Problema II ................................................................................................................ 17 
5.3 Problema III ............................................................................................................... 19 
5.4 Problema IV ............................................................................................................... 20 
5.5 Problema V ................................................................................................................ 22 
6. CONCLUSÃO ................................................................................................................. 23 
REFERÊNCIAS BIBLIOGRÁFICAS ................................................................................. 23 
ANEXOS ................................................................................................................................. 24 
 
 
 
 
 
 
 
 
 
 
 
 
3 
 
 
1. INTRODUÇÃO 
 
No decorrer da disciplina, foram apresentados os métodos de Eliminação de Gauss e 
Método de Gauss-LU para a solução de sistemas lineares, assim como os métodos de Newton-
Raphson e Método de Quasi-Newton para a solução de sistemas de equações não lineares. A 
partir dos conceitos iniciais ensinados e desenvolvidos em sala de aula, é notável o quanto esses 
métodos são eficazes porém trabalhosos de se resolver manualmente logo, para otimizar o 
trabalho, é utilizado em muitos casos o software MATLAB, um programa em que o elemento 
básico de informação é uma matriz que não requer dimensionamento. Apesar da facilidade que 
um recurso computacional oferece, é preciso um domínio sobre métodos de solução e ter uma 
noção do que é custo computacional. Só assim é possível determinar qual é a maneira mais 
viável de elaborar a solução de um sistema. 
Neste projeto será desenvolvido um algoritmo para resolução de sistemas de equações 
não-lineares através dos métodos de Newton-Raphson e o método de Quasi-Newton. Foram 
utilizados como ferramentas para execução do projeto, a produção de sub-rotinas para solução 
de sistemas lineares. Para o método de Newton-Raphson foi utilizado a eliminação de Gauss e 
para o método Quasi-Newton o método de Gauss-LU. Nos dois métodos de solução de sistemas 
lineares foram implementados com a estratégia de pivoteamento parcial. 
 
2. OBJETIVOS 
 
2.1 Objetivo Geral 
 
Este projeto tem como objetivo o desenvolvimento de um programa computacional para 
solução de sistemas de equações não-lineares através dos Métodos de Newton-Raphson e 
Quasi-Newton. O principal pré-requisito para a implementação do Método de Newton-Raphson 
e suas variantes é a existência de funções (sub-rotinas) para a solução de sistemas de equações 
lineares. Desta forma, também serão implementados dois métodos para a solução de sistemas 
lineares: eliminação de Gauss (usada no Método de Newton-Raphson) e Método de Gauss-LU 
(usada no método Quasi-Newton). Ambos os métodos de solução de sistemas lineares devem 
ser implementados com a estratégia de pivoteamento parcial. O método Quasi-Newton que será 
utilizado é o Método de Newton Modificado. Neste método não há atualização da matriz 
jacobiana durante o processo iterativo, isto é, a matriz jacobiana associada com a estimativa 
𝑥(0) é usada durante todo o processo iterativo. 
4 
 
 
2.2 Objetivos Específicos 
 
• Implementar os Método de Eliminação de Gauss e Gauss-LU; 
• Desenvolver um programa computacional para solução de sistemas não-lineares através 
do Método de Newton-Raphson; 
• Descrever conceitos e códigos utilizados; 
• Apresentar e analisar resultados encontrados; 
 
3. FUNDAMENTAÇÃO TEÓRICA 
 
3.1 Sistemas não-lineares 
 
Uma equação que contenha uma expressão do tipo 𝑥2, 𝑦−2 , 𝑥𝑦, 
√𝑦
𝑧
 , 𝑠𝑒𝑛(𝑥), 𝑒𝑥+𝑧 , etc., 
é chamada não-linear em 𝑥, 𝑦, 𝑧, ..., porque ela não pode ser escrita como 𝑎𝑥 + 𝑏𝑦 +
 𝑐𝑧 + . . . = 𝑐𝑡𝑒 que é uma equação linear 𝑒𝑚 𝑥, 𝑦, 𝑧. ... 
Um sistema de 𝑛 equações e 𝑛 incógnitas 𝑥1, 𝑥2, . . . , 𝑥𝑛 é chamado de não-linear se uma 
ou mais equações é não-linear. Trazendo todos os termos diferentes de zero à esquerda de todas 
as equações, tem-se uma forma geral que pode ser usada para qualquer sistema não-linear. 
 
{
𝑓1(𝑥1, 𝑥2, … , 𝑥𝑛) = 0
𝑓2(𝑥1, 𝑥2, … 𝑥𝑛) = 0
…
𝑓𝑛(𝑥1, 𝑥2, … 𝑥𝑛) = 0
 
 
Ou simplesmente: 
 
{
𝑓1(𝑥) = 0
𝑓2(𝑥) = 0
…
𝑓𝑛(𝑥) = 0
 
 
Em que 𝑥𝑇 = (𝑥1, 𝑥2, … , 𝑥𝑛). 
Em notação vetorial, o sistema linear acima pode ser escrito como: 𝐹(𝑥) = 0, em que: 
 
5 
 
 
𝑥 = (
𝑥1
𝑥2
…
𝑥𝑛
) e 𝐹(𝑥) = (
𝑓1(𝑥)
𝑓2(𝑥)
…
𝑓𝑛(𝑥)
) 
 
Um vetor que �̅� = (�̅�1, �̅�2, … , �̅�𝑛) que satisfaz 𝐹(�̅�) = 0 é denominado raiz do sistema 
não-linear. 
Uma solução para esse sistema de equações consiste de valores para as 𝑛 variáveis, tais 
que quando esses valores são substituídos nas equações, todas elas são satisfeitas 
simultaneamente, explanaremos sobre dois métodos que apresentam soluções para essa 
configuração de sistema. 
 
3.2 Método de Eliminação de Gauss 
 
A Eliminação Gaussiana,mais conhecida como método de eliminação de Gauss, é um 
método para resolver sistemas lineares que consiste em transformar convenientemente o 
sistema linear original para se obter um sistema linear equivalente com matriz dos coeficientes 
triangular superior. 
Para modificar convenientemente um sistema linear dado de forma a obter um sistema 
equivalente, faz-se o uso do seguinte teorema: 
• Sendo 𝐴𝑥 = 𝑏, um sistema linear, aplica-se uma sequência de operações 
elementares sobre as equações deste sistema, sendo elas: 
 
1. Trocar duas equações; 
2. Multiplicar uma equação por uma constante não-nula; 
3. Adicionar um múltiplo de uma equação e uma outra equação. 
 
Onde obtém-se um novo sistema �̅�𝑥 = �̅�, e os sistemas 𝐴𝑥 = 𝑏 e �̅�𝑥 = �̅� são 
equivalentes. Deste modo, o método da Eliminação de Gauss usará este teorema para 
triangularizar a matriz 𝐴. 
Considerando que det (𝐴) ≠ 0, é possível reescrever o sistema linear de forma que o 
elemento da posição 𝑎11 seja diferente de zero, usando a operação elementar (𝑖): 
 
6 
 
 
𝐴(0)|𝑏(0) = 𝐴|𝑏 =
(
 
 
𝑎11
(0) 𝑎12
(0) … 𝑎1𝑛
(0) |𝑏1
(0)
𝑎21
(0) 𝑎22
(0) … 𝑎2𝑛
(0) ⋮
⋮ ⋮ ⋮ ⋮ ⋮
𝑎𝑛1
(0) 𝑎𝑛2
(0) ⋯ 𝑎𝑛𝑛
(0) |𝑏𝑛
(0)
)
 
 
 
 
Onde 𝑎𝑖𝑗
(𝑘)
 é a denotação para o coeficiente da linha 𝑖 e coluna 𝑗 ao final da k-ésima 
etapa, e 𝑏𝑖
(𝑘)
 é o i-ésimo elemento do vetor constante no final da etapa k. Sendo: 
 
𝑎𝑖𝑗
(0)
= 𝑎𝑖𝑗, 𝑏𝑖
(0)
= 𝑏𝑖 e 𝑎11
(0)
≠ 0. 
 
• Etapa 1: 
 
Primeiramente, é feita a eliminação da variável 𝑥1 das equações 𝑖 = 2,… , 𝑛 da seguinte 
maneira: subtrai-se da equação 𝑖, a multiplicação da primeira equação por 𝑚𝑖1, sendo este 
definido como: 
𝑚𝑖1 =
𝑎𝑖1
(0)
𝑎11
(0), 𝑖 = 2,… , 𝑛 
 
Onde os elementos 𝑚𝑖1 =
𝑎𝑖1
(0)
𝑎11
(0), 𝑖 = 2,… , 𝑛 são os multiplicadores e o elemento 𝑎11
(0)
 é 
denominado pivô da primeira etapa. Logo, ao final desta etapa, têm-se a seguinte matriz: 
 
𝐴(1)|𝑏(1) =
(
 
 
𝑎11
(1) 𝑎12
(1) … 𝑎1𝑛
(1) |𝑏1
(1)
0 𝑎22
(1) … 𝑎2𝑛
(1) |𝑏2
(1)
⋮ ⋮ ⋮ ⋮ ⋮
0 𝑎𝑛2
(1) ⋯ 𝑎𝑛𝑛
(1) |𝑏𝑛
(1)
)
 
 
 
 
Onde: 
{
𝑎1𝑗
(1) = 𝑎1𝑗
(0) 𝑝𝑎𝑟𝑎 𝑗 = 1,… , 𝑛
𝑏1
(1)
= 𝑏1
(0)
 
 
E 
 
7 
 
 
{
𝑎𝑖𝑗
(1) = 𝑎𝑖𝑗
(0) −𝑚𝑖1𝑎1𝑗
(0) 𝑝𝑎𝑟𝑎 𝑖 = 1,… , 𝑛 𝑒 𝑗 = 1,… , 𝑛
𝑏𝑖
(1)
= 𝑏𝑖
(0)
−𝑚𝑖1𝑏1
(0)
 
 
• Etapa 2: 
 
O próximo passo é reescrever a matriz 𝐴(1) sem alterar a posição da linha 1, de forma 
que o pivô 𝑎22
(1)
 seja não nulo. Para isso, deve-se ter pelo menos um elemento 𝑎12
(1) ≠ 0, para 
𝑖 = 2, … , 𝑛, caso contrário, det(𝐴(1)) = 0, o que implicará que det(𝐴) = 0, sendo que 
det(𝐴) ≠ 0, por hipótese. 
Logo, os multiplicadores dessa etapa serão os elementos: 
 
𝑚𝑖2 =
𝑎𝑖2
(1)
𝑎22
(1) para 𝑖 = 3,… , 𝑛. 
 
A variável 𝑥2 é eliminada das equações 𝑖 = 3,… , 𝑛 da seguinte forma: da equação 𝑖, 
subtrai-se a multiplicação da segunda equação por 𝑚𝑖2. Por fim, tem-se a matriz 𝐴
(2)|𝑏(2): 
 
𝐴(2)|𝑏(2) =
(
 
 
 
 
𝑎11
(2) 𝑎12
(2) 𝑎13
(2) ⋯ 𝑎1𝑛
(2) |𝑏1
(2)
0 𝑎22
(2) 𝑎23
(2) ⋯ 𝑎2𝑛
(2) |𝑏2
(2)
0 0 𝑎33
(2) ⋯ 𝑎3𝑛
(2) |𝑏3
(2)
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
0 0 𝑎𝑛3
(2) ⋯ 𝑎𝑛𝑛
(2) 𝑏𝑛
(2)
)
 
 
 
 
 
 
Onde: 
{
𝑎𝑖𝑗
(2) = 𝑎𝑖𝑗
(1) 𝑝𝑎𝑟𝑎 𝑖 = 1, 2 𝑒 𝑗 = 𝑖, 𝑖 + 1,…𝑛
𝑏𝑖
(2)
= 𝑏𝑖
(1)
 𝑝𝑎𝑟𝑎 𝑖 = 1, 2
 
 
E 
 
{
𝑎𝑖𝑗
(2) = 𝑎𝑖𝑗
(1) −𝑚𝑖2𝑎2𝑗
(1) 𝑝𝑎𝑟𝑎 𝑖 = 3,… , 𝑛 𝑒 𝑗 = 2,… , 𝑛
𝑏𝑖
(2)
= 𝑏𝑖
(1)
−𝑚𝑖2𝑏2
(1) 𝑝𝑎𝑟𝑎 𝑖 = 3,… , 𝑛
 
 
8 
 
 
Portanto, seguindo raciocínio análogo, procede-se até a etapa (𝑛 − 1) e a matriz, ao 
final desta etapa, será: 
 
𝐴(𝑛−1)|𝑏(𝑛−1) =
(
 
 
 
𝑎11
(𝑛−1) 𝑎12
(𝑛−1) 𝑎13
(𝑛−1) ⋯ 𝑎1𝑛
(𝑛−1) |𝑏1
(𝑛−1)
0 𝑎22
(𝑛−1) 𝑎23
(𝑛−1) ⋯ 𝑎2𝑛
(𝑛−1) |𝑏2
(𝑛−1)
0 0 𝑎33
(𝑛−1) ⋯ 𝑎3𝑛
(𝑛−1) |𝑏3
(𝑛−1)
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
0 0 0 ⋯ 𝑎𝑛𝑛
(𝑛−1) 𝑏𝑛
(𝑛−1)
)
 
 
 
 
 
E o sistema linear 𝐴(𝑛−1)𝑥 = 𝑏(𝑛−1) é triangular superior e equivalente ao sistema linear 
original. 
 
3.3 Método de Gauss - LU 
 
O processo de fatoração LU consiste em decompor a matriz A dos coeficientes em um 
produto de dois ou mais fatores e, em seguida, resolver uma sequência de sistemas lineares 
que conduzirá à solução do sistema linear original. Por exemplo, se 𝐴 = 𝐿𝑈, o sistema linear 
𝐴𝑥 = 𝑏 pode ser escrito como: 
(𝐿𝑈)𝑥 = 𝑏 
 
A vantagem desse processo se encontra no fato de que se pode resolver qualquer sistema 
linear que tenha uma matriz A de coeficientes. Caso o vetor b seja alterado, a resolução do novo 
sistema linear será quase imediata. 
Os fatores L e U podem ser obtidos através de fórmulas para os elementos 𝑙𝑖𝑗 e 𝑢𝑖𝑗 ou 
podem ser construídos através do método de Eliminação de Gauss, a qual utilizaremos neste 
projeto. 
Para efeito prático de aprendizagem, utilizaremos um exemplo teórico de dimensão 3 
para facilitar o entendimento deste método: 
 
{
𝑎11𝑥1 + 𝑎12𝑥2 + 𝑎13𝑥3 = 𝑏1
𝑎21𝑥1 + 𝑎22𝑥2 + 𝑎23𝑥3 = 𝑏2
𝑎31𝑥1 + 𝑎32𝑥2 + 𝑎33𝑥3 = 𝑏3
 
 
Para a matriz de coeficientes, tem-se que: 
9 
 
 
𝐴(0) = (
𝑎11
(0)
𝑎12
(0)
𝑎13
(0)
𝑎21
(0)
𝑎22
(0)
𝑎23
(0)
𝑎31
(0)
𝑎32
(0)
𝑎33
(0)
) = 𝐴 
 
Os multiplicadores dessa primeira etapa serão: 
 
𝑚21 =
𝑎21
(0)
𝑎11
(0) e 𝑚31 =
𝑎31
(0)
𝑎11
(0) (supondo que 𝑎11
(0)
≠ 0) 
 
Para eliminar 𝑥1 da linha 𝑖, 𝑖 = 2, 3, multiplica-se a linha 1 por 𝑚𝑖1 e subtrai-se o 
resultado da linha 𝑖. Logo, os coeficientes 𝑎𝑖𝑗
(0)
 serão alterados para 𝑎𝑖𝑗
(1)
, onde: 
 
{
𝑎1𝑗
(1) = 𝑎1𝑗
(0) 𝑝𝑎𝑟𝑎 𝑗 = 1, 2, 3
𝑎𝑖𝑗
(1)
= 𝑎𝑖𝑗
(0)
−𝑚𝑖1𝑎1𝑗
(0) 𝑝𝑎𝑟𝑎 𝑖 = 2, 3 𝑒 𝑗 = 1, 2, 3
 
 
Estas operações correspondem a se pré-multiplicar a matriz 𝐴(0) pela matriz 𝑀(0), dada por: 
 
𝑀(0) = (
1 0 0
−𝑚21 1 0
−𝑚31 0 1
), pois: 
 
 
𝑀(0)𝐴(0) = (
1 0 0
−𝑚21 1 0
−𝑚31 0 1
)(
𝑎11
(0) 𝑎12
(0) 𝑎13
(0)
𝑎21
(0) 𝑎22
(0) 𝑎23
(0)
𝑎31
(0) 𝑎32
(0) 𝑎33
(0)
) = 
 
= (
𝑎11
(0)
𝑎12
(0)
𝑎13
(0)
𝑎21
(0)
−𝑚21𝑎11
(0)
𝑎22
(0)
−𝑚21𝑎12
(0)
𝑎23
(0)
−𝑚21𝑎13
(0)
𝑎31
(0)
−𝑚31𝑎11
(0)
𝑎32
(0)
−𝑚31𝑎12
(0)
𝑎33
(0)
−𝑚31𝑎13
(0)
) = 
 
= (
𝑎11
(1) 𝑎12
(1) 𝑎13
(1)
0 𝑎22
(1) 𝑎23
(1)
0 𝑎32
(1) 𝑎33
(1)
) = 𝐴(1) 
 
10 
 
 
Portanto, 𝑀(0)𝐴(0) = 𝐴(1) onde 𝐴(1) é a mesma matriz obtida no final da primeira etapa 
do processo de Eliminação de Gauss. 
Supondo agora que 𝑎22
(1)
≠ 0, o multiplicador da segunda etapa será: 
 
𝑚32 =
𝑎32
(1)
𝑎22
(1)
 
 
Para eliminar 𝑥2 da linha 3, multiplica-se a linha 2 por 𝑚32 e subtrai-se o resultado da 
linha 3. Logo, os coeficientes 𝑎𝑖𝑗
(1)
 serão alterados para: 
 
{
 
 
 
 𝑎1𝑗
(2) = 𝑎1𝑗
(1) 𝑝𝑎𝑟𝑎 𝑗 = 1, 2, 3
𝑎2𝑗
(2) = 𝑎2𝑗
(1) 𝑝𝑎𝑟𝑎 𝑗 = 2, 3
𝑎3𝑗
(2) = 𝑎3𝑗
(1) −𝑚32𝑎2𝑗
(1) 𝑝𝑎𝑟𝑎 𝑗 = 2, 3
 
 
As operações efetuadas em 𝐴(1) são equivalentes a pré-multiplicar 𝐴(1) por 𝑀(1), onde: 
 
𝑀(1) = (
1 0 0
0 1 0
0 −𝑚32 1
), pois: 
 
 
𝑀(1)𝐴(1) = (
1 0 0
0 1 0
0 −𝑚32 1
)(
𝑎11
(1) 𝑎12
(1) 𝑎13
(1)
0 𝑎22
(1) 𝑎23
(1)
0 𝑎32
(1)
𝑎33
(1)
) = 
 
= (
𝑎11
(1)
𝑎12
(1)
𝑎13
(1)
0 𝑎22
(1)
𝑎23
(1)
0 𝑎32
(1)
−𝑚32𝑎22
(1)
𝑎33
(1)
−𝑚32𝑎23
(1)
) = 
 
= (
𝑎11
(2) 𝑎12
(2) 𝑎13
(2)
0 𝑎22
(2) 𝑎23
(2)
0 0 𝑎33
(2)
) 
 
11 
 
 
Portanto, 𝑀(1)𝐴(1) = 𝐴(2) onde 𝐴(2) é a mesma matriz obtida no final da segunda etapa 
do processo de Eliminação de Gauss. 
Tem-se então que: 
𝐴 = 𝐴(0) 
𝐴(1) = 𝑀(0)𝐴(0) = 𝑀(0)𝐴 
𝐴(2) = 𝑀(1)𝐴(1) = 𝑀(1)𝑀(0)𝐴(0) = 𝑀(1)𝑀(0)𝐴 
 
Onde 𝐴(2) é triangular superior. É fácil verificar que: 
 
(𝑀(0))
−1
= (
1 0 0
𝑚21 1 0
𝑚31 0 1
) e (𝑀(1))
−1
= (
1 0 0
0 1 0
0 𝑚32 1
) 
 
Assim, 
 
(𝑀(0))
−1
(𝑀(1))
−1
= (
1 0 0
𝑚21 1 0
𝑚31 𝑚32 1
) 
 
Então, 𝐴 = (𝑀(1)𝑀(0))
−1
𝐴(2) = (𝑀(0))
−1
(𝑀(1))
−1
𝐴(2): 
 
𝐴 = (
1 0 0
𝑚21 1 0
𝑚31 𝑚32 1
)(
𝑎11
(2)
𝑎12
(2)
𝑎13
(2)
0 𝑎22
(2) 𝑎23
(2)
0 0 𝑎33
(2)
) = 𝐿𝑈 
 
Ou seja, 𝐿 = (𝑀(0))
−1
(𝑀(1))
−1
 e 𝑈 = 𝐴(2). 
Paraisto, foi necessário fatorar a matriz A em duas matrizes triangulares L e U, sendo 
que o fator L é triangular inferior com diagonal unitária e seus elementos 𝑙𝑖𝑗0 para 𝑖 < 𝑗 são os 
multiplicadores 𝑚𝑖𝑗 obtidos no método de Eliminação de Gauss; o fator U é triangular superior 
e é a matriz triangular superior obtida no final da fase da triangularização do método de 
Eliminação de Gauss. 
 
 
12 
 
 
3.4 Método de Newton-Raphson 
 
O método de Newton-Raphson, mais conhecido apenas como método de Newton, é o 
método mais amplamente estudado e utilizado para resolver sistemas de equações não-lineares. 
No caso de uma equação não linear a uma variável, o método de Newton consiste em se 
tomar um modelo local linear da função 𝑓(𝑥) em torno 𝑥𝑘, e este modelo é a reta tangente à 
função em 𝑥𝑘. 
Ampliando a motivação de se construir um modelo local linear para o caso de um 
sistema de equações não lineares, tem-se conhecida a aproximação 𝑥(𝑘) ∈ 𝐷, para qualquer 
𝑥 ∈ 𝐷, existe 𝑐𝑖 ∈ 𝐷, tal que: 
 
𝑓𝑖(𝑥) = 𝑓𝑖(𝑥
(𝑘)) + ∇𝑓𝑖(𝑐𝑖)
𝑇(𝑥 − 𝑥(𝑘)) 𝑖 = 1,… , 𝑛 
 
Aproximando ∇𝑓𝑖(𝑐𝑖) por ∇𝑓𝑖(𝑥
(𝑘)), 𝑖 = 1,… , 𝑛 tem-se um modelo local linear para 
𝑓𝑖(𝑥) em torno de 𝑥
(𝑘): 
 
𝑓𝑖(𝑥) ≈ 𝑓𝑖(𝑥
(𝑘)) + ∇𝑓𝑖(𝑥
(𝑘))
𝑇
(𝑥 − 𝑥(𝑘)) 𝑖 = 1,… , 𝑛 
 
E, portanto o modelo local linear para 𝐹(𝑥) em torno de 𝑥(𝑘) é: 
 
𝐹(𝑥) ≈ 𝐿𝑘(𝑥) = 𝐹(𝑥
(𝑘)) + 𝐽(𝑥(𝑘))(𝑥 − 𝑥(𝑘)) 
A nova aproximação 𝑥(𝑘+1) será o zero do modelo local linear 𝐿𝑘(𝑥), logo: 
 
𝐿𝑘(𝑥) = 0 ↔ 𝐽(𝑥
(𝑘))(𝑥 − 𝑥(𝑘)) = −𝐹(𝑥(𝑘)) 
 
Chamando (𝑥 − 𝑥(𝑘)) de 𝑠(𝑘), tem-se que 𝑥(𝑘+1) = 𝑥(𝑘) + 𝑠(𝑘), onde 𝑠(𝑘) é a solução 
do sistema linear: 
𝐽(𝑥(𝑘))𝑠 = −𝐹(𝑥(𝑘)) 
 
Logo, dado o ponto 𝑥(𝑘), a matriz 𝐽(𝑥(𝑘)) é obtida avaliando-se 𝐽(𝑥) em 𝑥(𝑘) e, em 
seguida, o passo de Newton, 𝑠(𝑘), é obtido a partir da resolução do sistema linear. 
Portanto, uma iteração de Newton requer basicamente: 
13 
 
 
1. A avaliação da matriz Jacobiana em 𝑥(𝑘); 
2. A resolução do sistema linear 𝐽(𝑥(𝑘))𝑠 = −𝐹(𝑥(𝑘)) e, por este motivo, cada iteração é 
considerada computacionalmente cara. 
 
Em relação ao item 2, pode-se empregar métodos baseados em fatoração da matriz 
Jacobiana. Também podem ser aplicados métodos iterativos quem por serem iterativos, obtêm 
uma aproximação para a solução exata do sistema linear e, consequentemente, o passo de 
Newton, 𝑠(𝑘), não é calculado exatamente. Por esta razão, o método de Newton com a resolução 
do sistema linear através de um método iterativo é denominado método de Newton inexato. 
A característica mais atraente do método de Newton é que sob condições adequadas 
envolvendo o ponto inicial 𝑥(0), a função 𝐹(𝑥) e a matriz Jacobiana 𝐽(𝑥), a sequencia gerada 
(𝑥(𝑘)) converge a (𝑥∗) com taxa quadrática. É importante observar que os resultados de 
convergência obtidos são os locais no sentido de que a aproximação inicial 𝑥(0) deve estar 
suficientemente próxima de (𝑥∗). 
 
3.5 Método de Quasi-Newton 
 
O método de Quase-Newton a ser utilizado no decorrer deste projeto é o método de 
Newton Modificado que consiste em se tomar a cada iteração k a matriz 𝐽(𝑥(0)), em vez de 
𝐽(𝑥(𝑘)): a partir de uma aproximação inicial 𝑥(0), a sequencia (𝑥(𝑘)) é gerada através de 
𝑥(𝑘+1) = 𝑥(𝑘) + 𝑠(𝑘), onde 𝑠(𝑘) é solução do sistema linear: 
 
𝐽(𝑥(0))𝑠 = −𝐹(𝑥(𝑘)) 
 
Dessa forma, a matriz Jacobiana é avaliada apenas uma vez e, para todo k, o sistema 
linear a ser resolvido a cada iteração terá a mesma matriz de coeficientes 𝐽(𝑥(0)). 
Caso se utilize a fatoração LU para resolve-lo, os fatores L e U serão calculados apenas 
uma vez e, a partir da segunda iteração, será necessário resolver apenas dois sistemas 
triangulares para obter o vetor 𝑠(𝑘). 
Existem vários outros métodos Quase-Newton, e é importante ressaltar que em vários 
deles, além de se evitar o cálculo da matriz Jacobiana, o esforço computacional necessário para 
a resolução do sistema linear é reduzido em relação ao esforço realizado no método de Newton. 
 
14 
 
 
4. METODOLOGIA 
 
4.1 Especificações 
 
• Linguagem de programação: 
 
O programa para a solução de sistemas de equações não-lineares será implementado no 
MATLAB. 
 
• Comandos permitidos e proibidos 
 
i) Não será permitida a sobrecarga dos operadores “+”,” *”, “/”, etc. para realizar 
operações diretas com vetores e matrizes na implementação de nenhum método 
numérico. 
ii) É proibido o uso das funções “chol”, “lu”, e o operador “\” do MATLAB para 
resolver sistemas lineares. 
iii) Não é permitido o uso da função “norm” do MATLAB para calcular as normas de 
vetores e matrizes; 
 
• Sugestões 
 
i) Divida o programa em funções. Por exemplo, crie funções para realizar as 
seguintes tarefas: construção da matriz jacobiana, cálculo do vetor associado com 
as equações não lineares (F(x)); cálculo de ||F(x)||∞; etc. 
ii) Construa funções separadas para a solução dos sistemas lineares. 
iii) Defina uma função principal (main) a partir da qual as funções de solução dos 
sistemas de equações não-lineares serão chamadas 
 
4.2 Descrição do Sistema Teste 
 
Os métodos numéricos para solução de sistemas de equações não-lineares serão testados 
no sistema tridiagonal de Broyden. Este sistema é definido como se segue: 
𝐹1(𝑥) = 𝑥1(3 − 0.5𝑥1) − 2𝑥2 + 1 
𝐹𝑖(𝑥) = 𝑥𝑖(3 − 0.5𝑥𝑖) − 𝑥𝑖−1 − 2𝑥𝑖+1 + 1, 𝑖 = 2,3, … , 𝑝 − 1 
15 
 
 
𝐹𝑝(𝑥) = 𝑥𝑝(3 − 0.5𝑥𝑝) − 𝑥𝑝−1 + 1, 𝑜𝑛𝑑𝑒 𝑝 é 𝑜 𝑛ú𝑚𝑒𝑟𝑜 𝑑𝑒 𝑒𝑞𝑢𝑎çõ𝑒𝑠 𝑑𝑜 𝑠𝑖𝑠𝑡𝑒𝑚𝑎. 
O ponto de partida para a solução do sistema linear definido acima é dado por: 𝑥𝑖
(0)
=
−1 𝑝𝑎𝑟𝑎 𝑖 = 1,… , 𝑝. 
A solução do sistema tridiagonal de Broyden deve ser obtida considerando uma 
tolerância 𝜀 = 1 × 10−6 para ||𝐹(𝑥)||∞. 
 
4.3 Problemas 
 
I. O grau e o padrão de esparsidade da matriz jacobiana para 𝑝 = 1000. Com base nestes 
resultados identifique se a matriz jacobiana é densa ou esparsa. 
II. A solução do sistema de Broyden para 𝑝 = 5, 10 𝑒 20 usando o Método de Newton-
Raphson. Apresente o vetor solução 𝑥, os valores finais de ||𝐹(𝑥)||∞ e também o 
número de iterações realizadas. 
III. A solução do sistema de Broyden para 𝑝 = 5, 10 𝑒 20 usando o Método de Newton 
Modificado. Apresente o vetor solução 𝑥, os valores finais de ||𝐹(𝑥)||∞ e também o 
número de iterações realizadas. 
IV. Compare os resultados obtidos nos itens (II) e (III) com aqueles obtidos pela função 
nativa do MATLAB para a solução de sistemas não-lineares. Apresente uma descrição 
resumida do método usado pela função nativa do MATLAB para solução de sistemas 
não-lineares. 
V. Avalie o custo computacional (tempo de processamento) para o obter a solução do 
sistema de Broyden com 𝑝 = 1000 para os seguintes métodos: Método de Newton-
Raphson, Método de Newton Modificado e o método usado pela função nativa do 
MATLAB. 
 
5. RESULTADOS E DISCUSSÃO 
 
5.1 Problema I 
 
O grau de esparsidade é uma métrica que informa a percentagem dos elementos nulos 
da matriz. Daí, seja 𝐺𝑆 o grau de esparsidade. Pode-se escrever como sendo: 
 
16 
 
 
𝐺𝑆 =
𝑛2 − 𝑁
𝑛2
 
 
Onde 𝑁 é o número de elementos não-nulos e 𝑛 é a ordem da matriz quadrada cujo grau 
de esparsidade será extraído. 
Para o jacobiano do sistema de Broyden, tem-se 2 elementos não nulos na primeira e na 
última linha. Nas linhas intermediárias 𝑛 − 2, tem-se 3 elementos não nulos. Daí: 
 
𝑁 = 3 × (𝑛 − 2) + 4 = 3𝑛 − 2 
Portanto, 
𝐺𝑆 =
𝑛2 − (3𝑛 − 2)
𝑛2
=
𝑛2 − 3𝑛 + 2
𝑛2
 
 
Para o caso particular em que 𝑛 = 𝑝 = 1000, tem-se: 
 
𝐺𝑆 =
10002 − 3 × 1000 + 2
10002
× 100% = 99.7002% 
 
Desta forma, como a maior parte dos elementos do jacobiano é nula, esta matriz é 
esparsa. 
Sobre o padrão de esparsidade, pode-se afirmar que é uma representação gráfica 
bidimensional, para o caso particular de matrizes bidimensionais,dos elementos não-nulos e 
nulos da matriz. Desta forma, a região sinalizada do padrão representa elementos não-nulos e a 
região vazia é constituída de zeros. Daí, através da função spy() do MATLAB, obtém-se, para 
𝑝 = 1000, a figura a seguir. O número 𝑛𝑧 = 2998 representa o número de elementos não 
nulos, calculados através de 3 × 𝑝 − 2, quando 𝑝 = 1000. 
 
17 
 
 
 
Figura 1 - Padrão de Esparsidade do Jacobiano para 𝑝 = 1000. 
 
5.2 Problema II 
 
As soluções encontradas usando o Método de Newton-Raphson para o sistema de 
Broyden são: 
 
• Para 𝒑 = 𝟓: 
 
a) Vetor Solução: 
𝑥 =
[
 
 
 
 
 −0.9684
−1.1870
−1.1485
−0.9590
−0.5942]
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 6.6411 × 10−9 
c) Iterações: 3 
18 
 
 
• Para 𝒑 = 𝟏𝟎: 
 
a) Vetor Solução: 
𝑥 =
[
 
 
 
 
 
 
 
 
 
 −1.0301
−1.3104
−1.3799
−1.3907
−1.3796
−1.3499
−1.2907
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 6.3098 × 10−7 
c) Iterações: 3 
 
• Para 𝒑 = 𝟐𝟎: 
 
a) Vetor Solução: 
𝑥 =
[
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
−1.0324
−1.3150
−1.3887
−1.4076
−1.4125
−1.4137
−1.4139
−1.4139
−1.4136
−1.4130
−1.4119
−1.4098
−1.4055
−1.3973
−1.3813
−1.3504
−1.2908
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 1.8145 × 10−12 
c) Iterações: 4 
19 
 
 
5.3 Problema III 
 
As soluções encontradas usando o Método de Newton Modificado para o sistema de 
Broyden são: 
 
• Para 𝒑 = 𝟓: 
 
a) Vetor Solução: 
𝑥 =
[
 
 
 
 
 −0.9684
−1.1870
−1.1485
−0.9590
−0.5942]
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 3.3131 × 10−7 
c) Iterações: 9 
 
• Para 𝒑 = 𝟏𝟎: 
 
a) Vetor Solução: 
𝑥 =
[
 
 
 
 
 
 
 
 
 
 −1.0301
−1.3104
−1.3799
−1.3907
−1.3796
−1.3499
−1.2907
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 1.9219 × 10−7 
c) Iterações: 10 
 
• Para 𝒑 = 𝟐𝟎: 
 
a) Vetor Solução: 
20 
 
 
𝑥 =
[
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
−1.0324
−1.3150
−1.3887
−1.4076
−1.4125
−1.4137
−1.4139
−1.4139
−1.4136
−1.4130
−1.4119
−1.4098
−1.4055
−1.3973
−1.3813
−1.3504
−1.2908
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
b) ||𝐹(𝑥)||
∞
= 7.9266 × 10−7 
c) Iterações: 10 
 
5.4 Problema IV 
 
Para facilitar a comparação entre os resultados obtidos nos problemas II e III, foram 
levantadas as Tabelas 1, 2 e 3 que contêm a as soluções fornecidas pelos métodos de Newton-
Raphson, Quasi-Newton e pela função fsolve, nativa do Matlab, para os sistemas de Broyden, 
de acordo com as ordens de interesse. Para quatro casas de precisão, as soluções foram 
idênticas, tanto para a ordem 5, quanto 10 e 20. Porém, no sistema de ordem 5, o número de 
iterações do Método Quasi-Newton foi o triplo do Newton-Raphson e de fsolve. Mas devemos 
ter em mente que a função fsolve, apesar de utilizar a mesma tolerância dos métodos, foi 
implementada de forma mais robusta e otimizada, através do método trust region dogleg. Para 
os três sistemas a tolerância e a solução inicial foram as mesmas, contudo a norma infinita do 
método de Newton-Raphson forneceu um valor 6,6411 × 10−9, muito mais preciso que o de 
Quasi-Newton, 3,3131 × 10−7. 
 
Tabela 1 - Comparação entre os valores encontrados para 𝑝 = 5. 
𝒑 = 𝟓 Newton_Raphson Quasi_Newton fsolve 
21 
 
 
Solução 𝑥 
[
 
 
 
 
 −0.9684
−1.1870
−1.1485
−0.9590
−0.5942]
 
 
 
 
 
[
 
 
 
 
 −0.9684
−1.1870
−1.1485
−0.9590
−0.5942]
 
 
 
 
 
[
 
 
 
 
 −0.9684
−1.1870
−1.1485
−0.9590
−0.5942]
 
 
 
 
 
Número de iterações 3 9 3 
Solução inicial −𝑜𝑛𝑒𝑠(5,1) −𝑜𝑛𝑒𝑠(5,1) −𝑜𝑛𝑒𝑠(5,1) 
Tolerância 1 × 10−6 1 × 10−6 1 × 10−6 
||𝐹(𝑥)||
∞
 6.6411 × 10−9 3.3131 × 10−7 6.6392 × 10−9 
 
Na Tabela 2, foram levantados os resultados para o sistema de ordem 10, onde além das 
soluções serem idênticas, já é possível notar uma grande desvantagem no método de Quasi-
Newton em relação ao método de Newton-Raphson e da função fsolve no número de iterações 
para alcançar a solução nas mesmas condições iniciais e com a mesma tolerância, pois 
novamente o método de Quasi-Newton gastou por volta de três vezes o número de iterações das 
concorrentes. Entretanto houve um balanceamento com relação à norma infinita. 
 
Tabela 2 - Comparação entre os valores encontrados para 𝑝 = 10. 
𝒑 = 𝟏𝟎 Newton_Raphson Quasi_Newton fsolve 
Solução 𝑥 
[
 
 
 
 
 
 
 
 
 
 −1.0301
−1.3104
−1.3799
−1.3907
−1.3796
−1.3499
−1.2907
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 −1.0301
−1.3104
−1.3799
−1.3907
−1.3796
−1.3499
−1.2907
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 −1.0301
−1.3104
−1.3799
−1.3907
−1.3796
−1.3499
−1.2907
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
Número de iterações 3 10 3 
Solução inicial −𝑜𝑛𝑒𝑠(10,1) −𝑜𝑛𝑒𝑠(10,1) −𝑜𝑛𝑒𝑠(10,1) 
Tolerância 1 × 10−6 1 × 10−6 1 × 10−6 
||𝐹(𝑥)||
∞
 6.3098 × 10−7 1.9219 × 10−7 3.0453 × 10−8 
 
Na Tabela 3, foram organizados os resultados para o terceiro sistema, de ordem 20. 
Como já comprovado nos sistemas anteriores, o número de iterações do método de Quasi-
22 
 
 
Newton é bem maior e, além disso, a precisão por parte da norma infinita do método de Newton-
Raphson, 1.8145 × 10−12, é muito maior que a do método de Quasi-Newton 7.9266 × 10−7. 
 
Tabela 3 - Comparação entre os valores encontrados para 𝑝 = 20. 
𝒑 = 𝟐𝟎 Newton_Raphson Quasi_Newton fsolve 
Solução 𝑥 
[
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
−1.0324
−1.3150
−1.3887
−1.4076
−1.4125
−1.4137
−1.4139
−1.4139
−1.4136
−1.4130
−1.4119
−1.4098
−1.4055
−1.3973
−1.3813
−1.3504
−1.2908
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
−1.0324
−1.3150
−1.3887
−1.4076
−1.4125
−1.4137
−1.4139
−1.4139
−1.4136
−1.4130
−1.4119
−1.4098
−1.4055
−1.3973
−1.3813
−1.3504
−1.2908
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
−1.0324
−1.3150
−1.3887
−1.4076
−1.4125
−1.4137
−1.4139
−1.4139
−1.4136
−1.4130
−1.4119
−1.4098
−1.4055
−1.3973
−1.3813
−1.3504
−1.2908
−1.1775
−0.9675
−0.5965]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Número de iterações 3 10 3 
Solução inicial −𝑜𝑛𝑒𝑠(10,1) −𝑜𝑛𝑒𝑠(10,1) −𝑜𝑛𝑒𝑠(10,1) 
Tolerância 1 × 10−6 1 × 10−6 1 × 10−6 
||𝐹(𝑥)||
∞
 6.3098 × 10−7 1.9219 × 10−7 3.0453 × 10−8 
 
5.5 Problema V 
 
Com o intuito de melhor avaliar o custo computacional para se obter a solução do 
sistema de Broyden com 𝑝 = 1000 para os métodos de Newton-Raphson, Newton Modificado 
e a função fsolve, foi levantada mais uma tabela para efeito de organização e comparação. 
Na tabela 4, são computados os tempos de processamento de cada um dos métodos e da 
função 𝑓𝑠𝑜𝑙𝑣𝑒, onde é observado um contrabalanceamento por parte do método de Quasi-
Newton, que obteve um tempo de processamento muito menor que do método de Newton-
Raphson, devido a maior complexidade deste. Apesar da função 𝑓𝑠𝑜𝑙𝑣𝑒 também ser 
23 
 
 
implementada de forma complexa, seu tempo de processamento foi bem próximo ao do método 
de Quasi-Newton, pois fsolve resolveu o sistema através do método trust region dogleg. 
 
𝒑 = 𝟏𝟎𝟎𝟎 Tempo de processamento (s) 
Newton-Raphson 130.4324 
Quasi-Newton 42.2919 
fsolve 4.2276 
 
 
6. CONCLUSÃO 
 
Os estudos dos métodos de Newton-Raphson e Quasi-Newton podem ser uma 
ferramenta muito útil para determinar e otimizar a resolução de problemas. O método de 
Newton-Raphson apresenta a vantagem de convergência do processo iterativo mais rápida. Isso 
o torna mais vantajoso em relação aos métodos anteriormente estudados e, para encontrar as 
raízes, não é necessária a condição 𝑓 (𝑎) × 𝑓 (𝑏) < 0, como já foi mencionado anteriormente. 
O método de Newton modificado tem a vantagem de calcular uma única vez a matriz 
Jacobiana. No caso de resolver por fatoração LU, os fatores L e U serão calculados uma única 
vez. Infelizmente esse método se mostrou ineficaz para sistemas com 𝑝 ≥ 10. Por issofoi 
implementado o método de Quasi-Newton proposto por Broyden, o método de Quasi-Newton 
tem uma aplicação mais viável na grande maioria das situações. 
 
REFERÊNCIAS BIBLIOGRÁFICAS 
 
RUGGIERO, Márcia A. Gomes; LOPES, Vera Lúcia da Rocha. CÁLCULO NUMÉRICO: 
ASPECTOS TEÓRICOS E COMPUTACIONAIS. Edição 2. Editora Pearson 
Universidades, 2000. 
 
SANCHES, I. J.; FURLAN, D. C. Métodos Numéricos. Universidade Federal do Paraná. 
Departamento de Informática. Curitiba, 2004. 
 
 
24 
 
 
ANEXOS 
 
• main_NewtonRaphson.m 
 
%% Script base para execução dos testes e comparações do sistema de Broyden 
 % resolvido por Newton-Raphson e pela função fsolve do MATLAB 
 
clc 
clear all 
 
p=20; %Ordem do sistema de Broyden; 
tol=1e-6; %Tolerância (norma infinita do vetor de resíduos); 
iter_max=800; %Número máximo de iterações realizadas pelo método de Newton-
Raphson; 
time_fsolve=0; %Tempo de processamento para solução do sistema com fsolve; 
time_newtonraphson=0; %Tempo de processamento para solução do sistema com 
Newton-Raphson; 
aux=0; %Variável auxiliar; 
 
 
%% Cálculo do grau e padrão de esparsidade da matriz para p 
x0=-ones(p,1); 
spy(Jacobian_Broyden(x0),'+'); 
title('Padrão de Esparsidade do Jacobiano') 
sparsity_degree=((p^2-3*p+2)/p^2)*100; 
 
clc 
aux=cputime; 
[x_nr,norm_res,iter]=Newton_Raphson(@Broyden,@Jacobian_Broyden,x0,tol,iter_
max) 
time_newtonraphson=cputime-aux; 
 
aux = cputime; 
[x_fs] = fsolve(@Broyden,x0); 
time_fsolve = cputime-aux; 
 
[x_nr x_fs] 
 
• main_quasiNewton.m 
 
clc 
clear all 
 
p=input('Digite a ordem do sistema: '); 
x0=-ones(1, p); %Solução inicial; 
tol = 1e-6; %Tolerência (norma infinita do vetor de resíduos); 
inter_max = 800; %Número máximo de iterações realizadas pelo método de 
Quasi-Newton; 
time_fsolve = 0; %Tempo de processamento para solução do sistema com 
fsolve; 
time_quasinewton = 0; %Tempo de processamento para solução do sistema com 
Quasi-Newton; 
aux = 0; %Variável auxiliar; 
 
25 
 
 
clc 
 
aux = cputime; 
[x_qn,iter, n] = Quasi2(x0, tol, inter_max) 
time_quasinewton = cputime-aux; 
 
aux = cputime; 
[x_fs] = fsolve(@Broyden,x0); 
time_fsolve = cputime-aux; 
 
• Newton_Raphson.m 
 
function [x,norm_res,iter] = 
Newton_Raphson(system,Jacobian,x0,tol,iter_max) 
 
%% Descrição da função: retorna a solução de um sistema não-linear F(x)=0 
 %utilizando o método de Newton-Raphson 
 
%% Entrada: 
% system: handle para o sistema base a ser resolvido F(x)=0 
% jacobian: handle para o jacobiano do sistema base 
% x0: estimativa inicial de solução 
% tol: tolerância esperada, medida através da norma infinita do vetor F de 
resíduos 
% iter_max: número máximo de iterações que podem ser realizadas pela função 
 
%% Saída: 
% x: solução do sistema 
% norm_res: norma do vetor de resíduos 
% iter: número de iterações realizadas 
 
n=length(x0); 
x=x0; 
iter=0; 
res=feval(system,x); 
norm_res=norm_inf(res); 
 
while(iter<iter_max && norm_res>tol) 
 J=feval(Jacobian,x); 
 [A,P]=Lugauss_fac(J);%Fatora o jacobiano em LU com pivotamento parcial 
em armazenamento compacto 
 dx=Lugauss_sol(A,P,-res); 
 x=sum_mat(x,dx); 
 res=feval(system,x); 
 norm_res=norm_inf(res); 
 iter=iter+1; 
end 
 
 
• Quasi2.m 
 
function [x,k, n]=Quasi2(x, tol, ninter_max) 
 % O método da secante tem por finalidade encontrar a solucao do sistema 
 % não-linear de equacões da forma F(x)=0 
 
26 
 
 
k=0; %O número de iterações inicia com valor zero; 
control=0; %Atribuição de valor zero a variável de controle; 
B=Jacobian(x); %Bk recebe o valor da jacobiana de F(x) em x; 
[B,P] = Lugauss_fac(B); 
 
 while (control==0)&&(k<=ninter_max) %Especifica quando o loop finaliza; 
 
f=Broyden(x); %Determina o valor de F(x) no ponto x em questão; 
 %Bks=-f; 
 s=Lugauss_sol(B,P,-f); 
 %Armazenamento do valor de x_k 
 a=transp(s); 
 xp=sum_mat(x,a); 
%Determinando o valor de F(x) no ponto xp: 
 f2=Broyden(xp); 
 %A variável y recebe a diferença de F(xp)-F(x: 
 y=sub(f2, f); 
 %Atualizacão da matriz Bk: 
 %B=B+((y-B*s)*a)/(a*s) 
 dum1=mult(B, s); 
 dum2=sub(y, dum1); 
 dum3=mult(dum2, a); 
 dum4=mult(a,s); 
 dum=dum3/dum4; 
 B=sum_mat(B, dum); 
 
 %Teste da variável de controle: 
 if norm_inf(xp-x)<=tol 
 control=1; 
 end 
 n = norm_inf(xp-x); 
 %Atualização do x: 
 x=xp; 
 %Atualização do número de iterações: 
 k=k+1; 
 end 
 
• Broyden.m 
 
function [F] = Broyden(x) 
 
n=length(x); %Definição da ordem do sistema tridiagonal de Broyden 
if(n == 1) 
 disp('O sistema tridiagonal de Broyden não está definido para ordem 
unitária!'); 
return 
end 
F=zeros(n,1); %Alocação de memória; 
 
%Definição do sistema: 
F(1)=x(1)*(3-.5*x(1))-2*x(2)+1; 
for i=2:n-1 
 F(i)=x(i)*(3-.5*x(i))-x(i-1)-2*x(i+1)+1; 
end 
F(n)=x(n)*(3-.5*x(n))-x(n-1)+1; 
 
 
27 
 
 
• Jacobian.m 
 
function [J] = Jacobian(x) 
%Retorna o Jacobino do Sistema Tridiagonal de Broyden 
%Escolhido como Sistema Teste 
 
n=length(x); %Definição da ordem do Jacobiano; 
if(n==1) 
 disp('O Jacobiano do sistema teste não está definido para ordem 
unitária!'); 
 return 
end 
J=zeros(n); %Alocação de memória para o Jacobiano; 
 
J(1,1)=3-x(1); J(1,2) = -2; %Definição da primeira linha do Jacobiano; 
 
%Definição da i-ésima linha do Jacobiano para i $\in$ \{2,..n-1\}; 
for i=2:n-1 
 J(i,i-1)=-1; 
 J(i,i)=3-x(i); 
 J(i,i+1)=-2; 
end 
 
J(n,n-1)=-1; J(n,n)=3-x(n); %Definição da última linha do Jacobiano; 
 
return 
 
• Jacobian_Broyden.m 
 
function [J] = Jacobian_Broyden(x) 
%Retorna o Jacobino do Sistema Tridiagonal de Broyden 
%Escolhido como Sistema Teste 
 
n=length(x); %Definição da ordem do Jacobiano; 
if(n==1) 
 disp('O Jacobiano do sistema de Broyden não está definido para ordem 
unitária!'); 
 return 
end 
J=zeros(n); %Alocação de memória para o Jacobiano; 
 
J(1,1)=3-x(1); J(1,2)=-2; %Definição da primeira linha do Jacobiano; 
 
%Definição da i-ésima linha do Jacobiano para i $\in$ \{2,..n-1\}; 
for i=2:n-1 
 J(i,i-1)=-1; 
 J(i,i)=3-x(i); 
 J(i,i+1)=-2; 
end 
 
J(n,n-1)=-1; J(n,n)=3-x(n); %Definição da última linha do Jacobiano; 
 
end 
 
 
28 
 
 
• Lugauss_fac.m 
 
%% Decomposição LU através de pivotamento parcial com retorno de 
coeficientes 
 %em esquema de armazenamento compacto e esquema de permutação 
 
function [A,P] = Lugauss_fac(A) 
 
n=length(A); 
P=zeros(n,1); %Guarda as permutações do pivotamento parcial; 
 
for i=1:n 
 P(i)=i; 
end 
 
tol_piv=sqrt(eps); %Pivô deve ser maior que a raiz quadrada do menor 
%número interpretado como não nulo pelo PC; 
 
for k=1:(n-1) 
 pv=abs(A(k,k)); 
 r=k; %Guarda a linha do pivô; 
 for i=(k+1):n 
 if(abs(A(i,k))>pv) 
 pv=abs(A(i,k)); %Bloco que define o pivô; 
 r=i; 
 end 
 end 
 
 if(pv <= tol_piv) 
 fprintf('\n A matriz é singular!'); 
 fprintf('\n Parando...'); 
 return 
 else 
 if (r~=k) 
 aux = P(k); 
 P(k) = P(r); 
 P(r) = aux; 
 for j = 1:n 
 aux = A(k,j); 
 A(k,j) = A(r,j); %Permuta as linhas r e j; 
 A(r,j) = aux; 
 end 
 end 
 end 
 for i = (k+1):n 
 mik = (A(i,k))/A(k,k); 
 A(i,k) = mik; %Coeficientes da matriz L; 
 for j = (k+1):n 
 A(i,j) = A(i,j)-mik*A(k,j); 
 end 
 end 
end 
 
• Lugauss_sol.m 
 
%% Função que resolve um sistema linear, dada a matriz de decomposição 
29 
 
 
 % LU compacta A, o vetor de permutações e o vetor de coeficientes b 
 
function [x] = Lugauss_sol(A,P,b) 
 
n = length(b); 
y = zeros(n,1); 
x = zeros(n,1); 
 
for i = 1:n 
 r = P(i); 
 c(i) = b(r); 
end 
 
for i = 1:n 
 acum = 0; 
 for j = 1:(i-1) 
 acum = acum + A(i,j)*y(j); 
 end 
 y(i) = c(i)-acum; 
end 
 
for i = n:-1:1 
 acum = 0; 
 for j = (i+1):n 
 acum = acum + A(i,j)*x(j);end 
 x(i) = (y(i)-acum)/A(i,i); 
end 
 
 
• norm_inf.m 
 
function [B]=norm_inf(x) 
 
B=max(abs(x)); 
 
 
• mult.m 
 
function [ A ] = mult( X, Y ) 
[m,n]=size(X); 
[j,k]=size(Y); 
 
if(n~=j) 
 error('Dimensões de Matrizes Incompatíveis'); 
end 
A=zeros(m,k); 
for a=1:m 
 for b=1:k 
 for v=1:j 
 A(a,b)=A(a,b)+(X(a,v)*Y(v,b)); 
 end 
 
 end 
end 
30 
 
 
end 
 
• sub.m 
 
function [ A ] = sub( X, Y ) 
[m,n]=size(X); 
[j,k]=size(Y); 
 
if(m~=j)||(n~=k) 
 error('Dimensões de Matrizes Incompatíveis'); 
end 
 
A=zeros(m,n); 
 
for a=1:m 
 for b=1:n 
 A(a,b)=X(a,b)-Y(a,b); 
 end 
end 
end 
 
• sum_mat.m 
 
% Função que soma duas matrizes 
 
function C = sum_mat(A,B); 
[M,N] = size(A); 
[O,P] = size(B); 
 
if(N~=P || M ~=O) 
 error('As matrizes não podem ser somadas! Dimensões incompatíveis'); 
else 
 C = zeros(M,N); 
 for i=1:M 
 for j=1:N 
 C(i,j)=A(i,j)+B(i,j); 
 end 
 end 
end 
 
• transp.m 
 
function [ b ] = transp(v) 
[m,n]=size(v); 
b=zeros(n,m); 
 
for i=1:n 
 for j=1:m 
 b(i,j)=v(j,i); 
 end 
end 
end