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
14 visitas
perguntada Nov 5 em Economia por Philippe Azevedo (16 pontos)  
editado Nov 5 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\\
x2 + x3 & x1 + x3 & x1 + x2 \\
x2x3 & x1x3 & x1x2
\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

1 Resposta

0 votos
respondida Nov 5 por Philippe Azevedo (16 pontos)  
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(self,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 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.
...