Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como calcular raízes de um polinômio de grau 3 com coeficientes reais, utilizando o método de newton

0 votos
254 visitas
perguntada Nov 5, 2020 em Economia por Philippe Azevedo (16 pontos)  
editado Dez 18, 2020 por Philippe Azevedo

Para um polinômio para o grau n, as raízes de um polinômio são os valores para que f(x) =0

\[f(x) = a_nx^n+a_{n-1}bx^{n-1}+ ...+ a_2x^2 + a_1x + a_0\]

As relações de Girard (ou teorema de Vieta) nada mais são que as equações que apresentam as relações entre as raízes de um polinômio com os seus coeficientes.

Para um polinômio de grau 2:
\[ x_1 + x_2 = -a_1/a_2\]\[x_1x_2 = a_0/a_2\]

Para um polinômio de grau 3:
\[a_3x^3 + a_2x^2+a_1x+a_0 = 0\]
Temos as seguintes equaçãoes pelas relações de girard:
\[x_1 + x_2 + x_3 = -a_2/a_3\]\[x_1x_2 + x_2x_3 + x_1x_3 = a_1/a_3\]\[x_1x_2x_3 = -a_0/a_3\]

Para a primeira equação temos as somas da raízes.
Para a segunda equação temos a soma entre todos os possíveis produtos entre duas raízes e assim sucessivamente, até que na e-nésima equação teremos o produto das n raízes.

Por fim, devemos relembrar que pelo teorema fundamental da álgebra, temos que se os coeficientes são reais e caso haja uma raiz complexa (a+bi), o conjugado dessa raiz (a-bi) também será raiz.

A ideia é utilizar uma função vetorial F(x) e verificar os valores para que F(x) =0
\[ F(x) = \begin{gather} \begin{bmatrix} x_1 + x_2 + x_3 + a_2/a_3\\ x_1x_2 + x_1x_3+ x_2x_3 - a_1/a_3 \\ x_1x_2x_3 + a_0/a_3 \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ \end{bmatrix} \end{gather} \]

Para essa questão, em virtude da existência do produto de variáveis, temos um sistema de equações não lineares.

Para resolver essa questão utilizaremos o método de newton que realiza a solução de sistemas não lineares. Ele é baseado em uma série de Taylor truncada:

\[F(x + s) ≈ F(x) + J(x)s\]

Onde s é um vetor infinitesimal, J(x) é a matriz jacobiana de F(x), formada pelas derivadas de cada
função de F(x) por cada variável, ou seja:

\[J(x)_{ij} = \{∂F_i(x)/∂x_j\}\]

Para um polinômio de grau 3, o jacobiano será:
\[ \begin{bmatrix} 1 & 1& 1\\ x_2 + x_3 & x_1 + x_3 & x_1 + x_2 \\ x_2x_3 & x_1x_3 & x_1x_2 \end{bmatrix} \]

Para que F(x+s) tenda a zero, temos que encontrar o vetor "J(x)s" que se aproxime de - F(x), cada vez mais, ou seja devemos resolver o sistema linear abaixo iterativamente:

\begin{equation}
J(x)s = −F(x)
\end{equation}
Na forma matricial:
\[ \begin{gather} \begin{bmatrix} 1 & 1& 1\\ x_2 + x_3 & x_1 + x_3 & x_1 + x_2 \\ x_2x_3 & x_1x_3 & x_1x_2 \end{bmatrix} * \begin{bmatrix} s_1\\ s_2\\ s_3\\ \end{bmatrix} = - \begin{bmatrix} x_1 + x_2 + x_3 + a_2/a_3\\ x_1x_2 + x_1x_3+ x_2x_3 - a_1/a_3 \\ x_1x_2x_3 + a_0/a_3 \end{bmatrix} \end{gather} \]

onde o vetor s, solução desse sistema linear, será o incremento para x.
Assim x + s passará a ser o novo vetor de raízes, sendo que esse processo ocorrerá até que F(x) se aproxime de 0, considerando determinada margem de erro. A verificação dessa aproximação se dá pelo meio do módulo do vetor.
Na primeira iteração devemos utilizar um vetor resposta candidato.

Compartilhe
comentou Dez 17, 2020 por danielcajueiro (5,581 pontos)  
Philipe, deve ter algum typo em uma das suas matrizes acima, pois o latex nao esta gerando.
comentou Dez 17, 2020 por Luiz Filippe (21 pontos)  
Philipe, achei sua resposta ótima! Tanto a teoria quanto os códigos, os quais consegui replicar quase que sem erro. O único erro que apareceu para mim foi quando tentei rodar a função "getCoeficientes". O erro foi que a sua função pede 2 entradas ("self" e "x") , quando na realidade só colocamos 1 entrada (o vetor inicial de 3 raízes escolhido aleatoriamente).  Essa função está fora da classe Polinomio3, certo? Então, não se deveria colocar o "self". Assim, ao invés de fazer

    def getCoeficientes(self,x):

Poderia reescrever assim:

    def getCoeficientes(x):

Eu reescrevi da forma que sugeri e consegui encontrar seus resultados. Quanto à matriz Jacobiana do polinômio de grau 3, que apareceu com um erro na questão, achei bastante curiosa. Copiei e colei seu código num outro espaço e a matriz apareceu lá sem problemas.
comentou Dez 18, 2020 por Philippe Azevedo (16 pontos)  
Verdade, do jeito que ficou modelado ficou melhor fora da classe mesmo.
Alterei conforme sua sugestão.

1 Resposta

0 votos
respondida Nov 5, 2020 por Philippe Azevedo (16 pontos)  
editado Dez 18, 2020 por Philippe Azevedo
import numpy as np
from itertools import combinations
import cmath

class Polinomio3:
    def __init__(self,coeficientes):
        if int(coeficientes[0]) != 0:
            if (len(coeficientes) == 4):
                self.__coef = list(coeficientes)
                self.__iteracoes = -1
            else:
                raise Exception("Não é polinômio de grau 3")    
        else:
            raise Exception("Coeficiente " + str(coeficientes) + " invalido!")
    #retorna o valor de f(x)
    def valor(self,x):
        valor =0
        for coef in self.__coef:
            valor = x*valor + coef
        return valor

    def __str__(self):
        i = 0 
        eq = ""
        grau = len(self.__coef) - 1
        for a in self.__coef:
            if grau - i != 0:
                if int(a) > 0:
                    if i != 0:
                        eq = eq + ' + ' + str(a) + "x^" + str(grau - i)
                    else:
                        eq = eq + str(a) + "x^" + str(grau - i)
                else:
                    if i ==0:
                        eq = eq + "-" + str(abs(a)) + "x^" + str(grau - i)
                    else:
                        eq = eq + " - " + str(abs(a)) + "x^" + str(grau - i)
            else:
                if a > 0:
                    eq = eq + ' + ' +str(a)
                elif a < 0:
                    eq = eq + ' - ' +str(abs(a))
            i= i + 1
        eq += ' = 0'
        return eq

    def grau(self):
        return len(self.__coef)-1

    # calcula o valor da funcao vetorial F
    def __F1(self,x):
        return np.array(
            [x[0] + x[1] + x[2] + self.__coef[1]/self.__coef[0],
            x[0]*x[1] + x[0]*x[2] + x[1]*x[2] - self.__coef[2]/self.__coef[0],
            x[0]*x[1]*x[2] + self.__coef[3]/self.__coef[0]])

    # calcula o jacobiano de F
    def __J1(self,x):
        return np.array(
        [[1.,1.,1.],
        [x[1] + x[2], x[0] + x[2], x[0] + x[1]],
        [x[1]*x[2],x[0]*x[2], x[0]*x[1]]])

    #calcula as raizes aproximadas para F(x) = 0
    def getRaizes(self,xx,tol): #raizes iniciais, erro
        self.__x = xx
        F = self.__F1(xx)
        modF = np.linalg.norm(F, ord=2)
        self.__iteracoes = 0
        while abs(modF) > tol and self.__iteracoes < 10000:
            delta = np.linalg.solve(self.__J1(xx), -F)
            xx += delta
            F = self.__F1(xx)
            modF = np.linalg.norm(F, ord=2)
            self.__iteracoes += 1    
        # Nao encontrou solucao ou ultrapassou o limite de iteracoes
        if abs(modF) > tol:
            self.__iteracoes = -1
        return self.__x

    # retorna a quantidade de iteracoes para obter as raizes, -1 se houve erro
    def getIteracoes(self):
        return self.__iteracoes

Definimos um procedimento para o cálculo dos coeficientes, conforme relações de girard:

#calcula coeficientes de a_2 em diante
def getCoeficientes(x):
    return np.array(
        [-(x[0] + x[1] + x[2]),
        (x[0]*x[1] + x[0]*x[2] + x[1]*x[2]),
        -x[0]*x[1]*x[2] ])

Método principal

 if __name__ == "__main__":

Devemos escolher um polinômio inicialmente, escolhemos 3 raízes aleatórias, sendo uma delas real.

print(getCoeficientes([2+4j, 2-4j,6]))

o primeiro coeficiente sempre será 1, por simplificação
[ -10.-0.j 44.+0.j -120.+0.j], ou seja os coeficientes são 10, 44 e -120

coef= np.array([1,-10,44,-120])
poli = Polinomio3(coef)
print(str(poli))

1x^3 - 10x^2 + 44x^1 - 120 = 0

print(str(poli.grau()))

3

print(poli.valor(8))

104

#tolerancia de erro desejada para o cálculo de f(x)
tol = 1e-6
#estimativa inicial para raízes
raizes_esperadas = np.array([2+3j,2-3j,4.])
x = poli.getRaizes(raizes_esperadas, tol)
print(x)

[2.+4.00000000e+00j 2.-4.00000000e+00j 6.-2.30577903e-16j]
a última raiz com parte imaginária na ordem de -16, é real :)

print(poli.getIteracoes())

5

comentou Nov 5, 2020 por Philippe Azevedo (16 pontos)  
Estava fazendo para o grau n, por meio de combinação de conjuntos, utilizando strings.

        def __getRelacoesGirard(self,n):
            input = []
            funcoes =[]
            if n > 0:
                i = 1
                while i <= n :
                    input.append('x['+ str(i-1)+ ']')
                    i = i+1
                funcoes =  sum([list(map(list, combinations(input, i))) for i in range(len(input) + 1)], [])
            else:
                print('nao é um polinomio')
            return funcoes

Estava utilizando a função eval para avaliar essas expressões em strings. Não ficou bom.
comentou Dez 17, 2020 por danielcajueiro (5,581 pontos)  
Essas coisas dao trabalho.
...