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

Solução de um sistema de equações não lineares

0 votos
13 visitas
perguntada Nov 22 em Matemática por Luiz Filippe (11 pontos)  
editado Nov 22 por Luiz Filippe

Um modelo de combustão do propano no ar retorna o seguinte sistema de equações não lineares:

\[x_1+x_4=3 \quad,\]
\[2x_1+x_2+x_4+x_7+x_8+x_9+2x_{10}=R+10 \quad,\]
\[2x_2+2x_5+x_6+x_7=8 \quad,\]
\[2x_3+x_5=4R \quad, \]
\[x_1x_5-0,193x_2x_4=0 \quad, \]
\[x_6\sqrt{x_2}-0,002597\sqrt{x_2x_4S}=0 \quad, \]
\[x_7\sqrt{x_4}-0,003448\sqrt{x_1x_4S}=0 \quad, \]
\[x_4x_8-0,00001799x_2S=0 \quad, \]
\[x_4x_9-0,0002155x_1\sqrt{x_3S} \quad, \]
\[x^2_4(x_{10}-0,00003846S)=0 \quad, \]

no qual \(R=4,056734\) e \(S= \sum_{i=1}^{10} x_i\). Use uma biblioteca para resolver este sistema não linear. (Dica: Caso alguma raiz quadrada retorne erro devido ao input negativo, substitua tal variável por seu valor absoluto.) (Livro, página 253, questão 5.26)

REFERÊNCIA DA QUESTÃO: HEATH, M. T. Scientific Computing: An Introductory Survey. 2 edição. Ed: McGraw-Hill. Nova York: , 2002. Disponível em: < http://www-e6.ijs.si/~roman/files/tmp/M.Heath-SComputing/scientific-computing-michael-t-heath.pdf >

Compartilhe

2 Respostas

0 votos
respondida Nov 22 por Luiz Filippe (11 pontos)  
editado 2 dias atrás por Luiz Filippe

Para resolvermos o sistema de equações não lineares proposto, podemos usar a biblioteca "scipy", a qual possui diversas ferramentas elaboradas para as área de matemática (como otimização, integração, entre outros). O método empregado será o de otimização a partir da técnica "fsolve". Esta técnica aproxima a função objetivo \(f(x)\) de 0, retornando, assim, o valor de "\(x\)" que retorna \(f(x)=0\); i.e., a raiz da função.

Primeiro, define-se um possível valor para "\(x\)" e, a partir dele, chega-se ao valor verdadeiro (dentro do erro tolerado). Por padrão, a tolerância de erro é de \(1.49012e-08\), mas há a possibilidade de alterá-lo. A implementação em Python é a seguinte:

#importa as bibliotecas necessárias
import numpy as np
import scipy
from scipy import optimize
from optimize import fsolve

A técnica de fsolve aceita apenas 1 função como input. Desta forma, para podermos trabalhar com tal método num cenário com múltiplas funções, devemos usar vetores:

# definimos a função "f" que vai entrar como input no fsolve
# definimos também todas as variáveis "x", além de "R" e "S"
def f(abcdefghij):
        x1 = abcdefghijl[0]
        x2 = abcdefghijl[1]
        x3 = abcdefghijl[2]
        x4 = abcdefghijl[3]
        x5 = abcdefghijl[4]
        x6 = abcdefghijl[5]
        x7 = abcdefghijl[6]
        x8 = abcdefghijl[7]
        x9 = abcdefghijl[8]
        x10 = abcdefghijl[9]
        R = 4.056734
        S = abcdefghijl[10]

#organizamos as funções de modo a torná-las no formato f(x)=0
#cada uma dessas funções ocupará um slot do vetor de funções "f"
        f0 = x1 + x4
        f1 = 2*x1 + x2 + x4 + x7 + x8 + x9 + 2*x10 - R - 10
        f2 = 2*x2 + 2*x5 + x6 +x7 - 8
        f3 = 2*x3 + x5 - 4*R
        f4 = x1*x5 - 0.193*x2*x4
        f5 = x6*np.sqrt(abs(x2)) - 0.002597*np.sqrt(abs(x2*x4*S))
        f6 = x7*np.sqrt(abs(x4)) - 0.003448*np.sqrt(abs(x1*x4*S))
        f7 = x4*x8 - 0.00001799*x2*S
        f8 = x4*x9 - 0.0002155*x1*np.sqrt(abs(x3*S))
        f9 = (x4**2)*(x10 - 0.00003846*S)
        f10 = S - x1 - x2 - x3 - x4 - x5 - x6 - x7 - x8 - x9 - x10

#criação do vetor f       
       return np.array([f0, f1, f2, f3, f4, f5, f6, f7, f8, f9,f10])

#damos um chute inicial para a suposta raiz (para cada função)
       chute_abcdefghijl = np.array([0,0,0,0,0,0,0,0,0,0,0])

#por fim, resolvemos o sistema
       abcdefghij = fsolve(f, chute_abcdefghij)

       x1 = abcdefghijl[0]
       x2 = abcdefghijl[1]
       x3 = abcdefghijl[2]
       x4 = abcdefghijl[3]
       x5 = abcdefghijl[4]
       x6 = abcdefghijl[5]
       x7 = abcdefghijl[6]
       x8 = abcdefghijl[7]
       x9 = abcdefghijl[8]
       x10 = abcdefghijl[9]
       S=abcdefghijl[10]

Após rodarmos o código acima, podemos exibir a solução de x e o valor de f quando empregamos tal vetor "\(x\)" de soluções.

       print(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,S)

--------------
Output
--------------
[ [3.000000010727565] ,[3.999999347453441], [8.113468001378202], 
  [-1.0723254607500187e- 08], [-2.758677545712311e-09], [5.371916497930608e-11], 
  [1.3105568356247757e-06], [0.020042985282166775], [-34.3037135501549], 
  [19.170201948066815], [3.9909073325860474e-08]]


print(f(abcdefghijl))

---------------
Output
---------------
[ 4.31032987e-12  3.04822834e-12  8.17124146e-14 -2.27373675e-12
  2.31853976e-12 -1.03479307e-14  1.21667628e-11 -2.17797891e-10
  -3.37642850e-11  2.20434681e-15  2.71640488e-11]

O método foi bem sucedido porque podemos observar que as funções tendem a zero quando aplicamos o conjunto de soluções obtido.

0 votos
respondida Nov 26 por Felipe Yudi (51 pontos)  

Olá Luiz,

Interessante usar a biblioteca do Scipy, mas tenho 3 comentários:

1-) S também não é uma variável? Não deveríamos retornar o valor de S também? Da maneira que está parece que a equação \(\sum_{i=1}^{10}x_i\) não está incluída no sistema.

2-) No código, f0 está como

f0 = x1 + x4

mas deveria ser

f0 = x1 + x4 -3

3-) No código f7 está como

x7*np.sqrt(abs(x4)) - 0.003448*np.sqrt(abs(x1*x4*S))

mas deveria ser

x4*x8 - 0.00001799*x2*S

Abaixo deixo uma outra solução com a biblioteca gekko.

from gekko import GEKKO

def fGekko():
    m = GEKKO() 

    R = 4.056734

    x1 = m.Var(value = 1)
    x2 = m.Var(value = 1)
    x3 = m.Var(value = 1)
    x4 = m.Var(value = 1)
    x5 = m.Var(value = 1)
    x6 = m.Var(value = 1)
    x7 = m.Var(value = 1)
    x8 = m.Var(value = 1)
    x9 = m.Var(value = 1)
    x10 = m.Var(value = 1)
    S = m.Var(value = 1)

    m.Equation(x1 + x4 -3 == 0)
    m.Equation(2*x1 + x2 + x4 + x7 + x8 + x9 + 2*x10 - R - 10 == 0)
    m.Equation(2*x2 + 2*x5 + x6 +x7 - 8 == 0)
    m.Equation(2*x3 + x5 - 4*R == 0)
    m.Equation(x1*x5 - 0.193*x2*x4 == 0)
    m.Equation(x6*m.sqrt(m.abs(x2)) - 0.002597*m.sqrt(m.abs(x2*x4*S)) == 0)
    m.Equation(x7*m.sqrt(m.abs(x4)) - 0.003448*m.sqrt(m.abs(x1*x4*S)) == 0)
    m.Equation(x4*x8 - 0.00001799*x2*S == 0)
    m.Equation(x4*x9 - 0.0002155*x1*m.sqrt(m.abs(x3*S)) == 0)
    m.Equation((x4**2)*(x10 - 0.00003846*S) == 0)

    m.solve(disp=False)

    return (x1.Value, x2.Value, x3.Value, x4.Value, x5.Value, x6.Value, 
            x7.Value, x8.Value, x9.Value, x10.Value, S.Value)

Cujo output é:

([2.9999932241], [3.9999524149], [8.1134671282], [6.7758918215e-06], [1.7436495136e-06], [9.6845014075e-08], [9.1586058478e-05], [0.0027277107652], [4.0493323139], [0.0023183751413], [0.00020635227263])
comentou 2 dias atrás por Luiz Filippe (11 pontos)  
Obrigado pelas contribuições, Felipe! Vou corrigir essas questões.
...