Sistema não lineares – Newton-Raphson

Postado por Raony em 2 - agosto - 2011

Olá pessoal,

Primeiramente quero ~~fazer a média~~ agradecer ao pessoal do condição inicial pela confiança creditada a mim e pelo convite para fazer parte desse grupo. Espero não decepcioná-los.

Então! Vamos deixar de conversinha e trabalhar, né?

Como 1º post vou mostrar para vocês o método de Newton-Raphson, já mostrado aqui, estendido para aplicação à resolução de sistemas não lineares.

Antes vamos relembrar um pouco. Sabemos que se trata de um método iterativo onde o objetivo é encontrar a raiz de uma função partindo de uma estimativa inicial $x^0$ , e que a próxima estimativa é calculada através da tangente ao ponto $(x^0,f(x^0))$ , ou seja, a derivada em $x^0$ , certo?

No nosso problema agora teremos $n$ equações com $n$ variáveis cada:

$\left\{\begin{array}{cc} f_1(x_1,x_2,\cdots,x_n) = 0 \\ f_2(x_1,x_2,\cdots,x_n) = 0 \\ \vdots \\ f_n(x_1,x_2,\cdots,x_n) = 0\end{array}\right.$

Sendo assim o objetivo é encontrar o conjunto de $n$ variáveis $(x_i)$ que zere as $n$ equações. Como antes, tínhamos que derivar a função em torno de $x^0$ para encontrarmos $x^1$ e assim por diante, ou seja, a função era linearizada em torno de um ponto através da Série de Taylor, para sistema faremos o mesmo para todas as $n$ equações, assim vamos encontrar:

$\begin{array}{cc} f_1(x_{1}^{k+1},x_{2}^{k+1},\cdots, x_{n}^{k+1}) = f_1(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) + \dfrac{\partial f_1}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) +\\ \dfrac{\partial f_1}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) + \cdots + \dfrac{\partial f_1}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \\ f_2(x_{1}^{k+1},x_{2}^{k+1},\cdots, x_{n}^{k+1}) = f_2(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) + \dfrac{\partial f_2}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) +\\ \dfrac{\partial f_2}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) + \cdots + \dfrac{\partial f_2}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \\ f_n(x_{1}^{k+1},x_{2}^{k+1},\cdots, x_{n}^{k+1}) = f_n(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) + \dfrac{\partial f_n}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) +\\ \dfrac{\partial f_n}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) + \cdots + \dfrac{\partial f_n}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \end{array}$

Como a intenção é encontrar os valores de $x^{k+1}_{i}$ que zere as funções, fazemos $f_1(x_1^{k},x_2^{k},\cdots, x_n^{k}), f_2(x_1^{k},x_2^{k},\cdots, x_n^{k}), \cdots , f_n(x_1^{k},x_2^{k},\cdots, x_n^{k}) = 0$ e arrumando, temos:

$\begin{array}{cc} - f_1(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) = \dfrac{\partial f_1}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) + \dfrac{\partial f_1}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) +\\ \cdots + \dfrac{\partial f_1}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \\ - f_2(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) = \dfrac{\partial f_2}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) + \dfrac{\partial f_2}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) +\\ \cdots + \dfrac{\partial f_2}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \\ - f_n(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) = \dfrac{\partial f_n}{\partial x_1}\vert_{x_1 = x_{1}^{k}}\cdot (x_{1}^{k+1} - x_{1}^{k}) + \dfrac{\partial f_n}{\partial x_2}\vert_{x_2 = x_{2}^{k}}\cdot (x_{2}^{k+1} - x_{2}^{k}) +\\ \cdots + \dfrac{\partial f_n}{\partial x_n}\vert_{x_n = x_{n}^{k}}\cdot (x_{n}^{k+1} - x_{n}^{k})\\ \end{array}$

Agora o que era um sistema não linear passar a ser linear que pode ser representado na forma matricial:

$\begin{matrix} \underbrace{\begin{bmatrix} - f_1(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k})\\ \\ - f_2(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k})\\ \vdots \\ - f_n(x_{1}^{k},x_{2}^{k},\cdots, x_{n}^{k}) \end{bmatrix}}\\ \bf -F_k \end{matrix} = \begin{matrix} \underbrace{\begin{bmatrix} \dfrac{\partial f_1}{\partial x_1}\vert_{x_1 = x_{1}^{k}} & \dfrac{\partial f_1}{\partial x_2}\vert_{x_1 = x_{1}^{k}} & \cdots & \dfrac{\partial f_1}{\partial x_n}\vert_{x_n = x_{n}^{k}} \\ \\ \dfrac{\partial f_2}{\partial x_1}\vert_{x_1 = x_{1}^{k}} & \dfrac{\partial f_2}{\partial x_2}\vert_{x_1 = x_{1}^{k}} & \cdots & \dfrac{\partial f_2}{\partial x_n}\vert_{x_n = x_{n}^{k}} \\ \vdots & \vdots & \vdots & \vdots\\ \dfrac{\partial f_n}{\partial x_1}\vert_{x_1 = x_{1}^{k}} & \dfrac{\partial f_n}{\partial x_2}\vert_{x_1 = x_{1}^{k}} & \cdots & \dfrac{\partial f_n}{\partial x_n}\vert_{x_n = x_{n}^{k}} \end{bmatrix}}\\ \bf J_k \end{matrix} \cdot \begin{matrix} \underbrace{\begin{bmatrix} (x_{1}^{k+1} - x_{1}^{k})\\ \\ (x_{2}^{k+1} - x_{2}^{k})\\ \vdots \\ (x_{n}^{k+1} - x_{n}^{k}) \end{bmatrix}}\\ \bf X_{k+1} - X_k \end{matrix}$

Isolando o vetor $\bf X_{k+1}$ , vamos encontrar: $\bf X_{k+1} = X_{k} - J_{k}^{-1}\cdot F_k$

Onde $\bf J_k$ é a matriz jacobiana calculada para a iteração $k$ . Resumindo, o que precisamos é de um vetor de estimativas iniciais e uma matriz com as derivadas parciais analíticas de cada função.

Então! Embora não tenha convergência garantida, este método é sensível a estimativa inicial, porém é o mais usado devido a velocidade de convergência, ~~quando converge~~ . As mesma considerações mostradas para o caso de uma variável devem ser respeitadas aqui, lembrando que devem ser aplicadas à todas as funções e todas variáveis.

Separei uns exemplos pra ser resolvido usando como ferramenta o MATLAB, onde elaborei dois códigos, um com variáveis do tipo simbólico, onde é possível obtermos as derivadas analíticas, e outro onde usamos o fsolve, uma função interna do MATLAB que usa derivação numérica para a resolução.

CASO1: Programação usando o simbólico

Para este primeiro exemplo usaremo o toolbox simbólico do MATLAB

Dado o sistema:

$\left\{\begin{matrix} x + y = 2\\ x\cdot y = -3 \end{matrix}\right.$

Temos um sistema de 2 equações e 2 variáveis. O importante é que a função esteja escrita de forma a igualar a zero, pois o método busca o zero da função, por tanto:

$\left\{\begin{matrix} x + y - 2 = 0\\ x \cdot y + 3 = 0 \end{matrix}\right.$

E a matriz jacobiana será:

$J = \begin{bmatrix}1 & 1\\ y & x \end{bmatrix}$

Com isso o cálculo das variáveis a cada iteração será feito da seguinte forma:

$\begin{bmatrix} x\\ y \end{bmatrix}_{k+1} = \begin{bmatrix} x\\ y \end{bmatrix}_{k} + \begin{bmatrix} 1 & 1\\ y & x \end{bmatrix}^{-1} \cdot \begin{bmatrix} x + y - 2\\ x\cdot y + 3 \end{bmatrix}$

Como disse anteriormente, essa metodologia não garante convergência por isso devemos ter cuidado com a estimativa inicial. Neste exemplo, note que se as estimativas iniciais de $x$ e $y$ forem iguais, a matriz jacobiana se torna inversível, ou seja, não conseguimos resolver o sistema.

Segue o código com o algoritmo implementado para este exemplo.

^?Download NRsim.m

% Código para a resolução de sistema não lineares utilizando variaveis do
% tipo simbólico.
clear all
close all
clc
 
syms x y % Definição das variaveis simbolicas. A partir daqui,
             % toda e qualquer variavel que depedem desta se tornará
             % simbólica.
% Criação do vetor contendo as funções
f(1) = x + y - 2;
f(2) = x*y + 3;
 
% Criação do vetor com as variaveis. Dessa forma, não preciso, toda vez que
% for necessário me referencias a x e y, assim basta se referenca a var
varia = [x;y];
% Calculo da matiz jacobiana. A função jacobian, calcula a matriz
% jacobiana dado um vetor de funcoes simbolicas e o vetor das variaveis
% simbolicas que se deseja encontrar as derivadas.
jacob = jacobian(f,varia);
 
% Aqui é onde começa o calculo iterativo, onde deve ter o vetor de
% estimativa inicial.
estimativas_iniciais = [9 -9;-9 9];
intervalo = [0.5:0.01:9; -9:0.01:-0.5]; % para construcao dos graficos
tol = 1e-5; % Tolerancia de convergencia 
for ii = 1:size(estimativas_iniciais,1)
    it = 0;
    xold = estimativas_iniciais(:,ii); % Estimativa inicial
    x = xold(1);
    y = xold(2);
    % A função eval, faz o calculo numérico para o valor estabelecido para as
    % variaveis symbolicas, neste caso x e y
    F = eval(f);
    Jac = eval(jacob);
    % Aqui eremos guardar os valore de cada estimativa e os valore das funções
    % encontrada em cada interação
    vx = xold;
    vF = F';
    % O critério de parada foi a soma dos valores absolutos de F, e para previnir
    % um loop inifito complementa-se com o numero máximo de iterações  
    while sum(abs(F)) > tol && it <= 100
        xnew = xold - inv(Jac)*(F)';
        xold = xnew;
        x = xold(1);
        y = xold(2);
        F = eval(f);
        Jac = eval(jacob);
        vx = [vx xold];
        vF = [vF F'];
        it = it + 1;
    end
    % Impreção de resultados
    figure(1)
    xx = intervalo(ii,:);
    yy1 = 2-xx;
    yy2 = -3./xx;
    plot(xx,[yy1;yy2],vx(1,:),2-vx(1,:),'or',vx(1,:),-3./vx(1,:),'or','LineWidth',2)
    grid on
    xlabel('x','fontsize',16)
    ylabel('y','fontsize',16)
    legend('y = 2 - x','y = -3/x','Iterações')
    hold on
    figure(ii+1)
    plot(0:length(vF)-1,vF','-o',[0 it],[tol tol],'--','LineWidth',2)
    grid on
    xlabel('Iterações','fontsize',16)
    ylabel('Função','fontsize',16)
    legend('Função 1','Função 2',sprintf('Tolerancia = %1.1e',tol),'Location','Best')
end

Para visualizarmos a solução, construir um gráfico contendo a reta $(y = 2 - x)$ e a hipérbole $(y = -3/x)$ para $x \neq 0$ .

Vemos que há dois pontos de interseção entre as curvas, ou seja, exitem duas soluções para este sistema. Para uma estimativa de $x = 9$ e $y = -9$ o resultado encontrado é $x = 3$ e $y = -1$ , no caso de inverter as estimativas vamos encontrar $x = -1$ e $y = 3$ .

Para este exemplo, as variáveis não são grandezas físicas, portanto, ambas as soluções são válidas, porém em se tratando de grandezas físicas, é necessário conhecermos a região onde a solução se encontrar para não cometer erro na estimativa inicial.

Para ilustrar a convergências do sistema, fiz outro gráfico contendo os valores obtidos das funções cada iteração.

Vemos que há convergência de ambas as funções a partir da estimativa inicial. Nota-se que a função 1 $(x + y - 2 = 0)$ converge primeiro. No entanto como a solução tem que satisfazer ambas a equações, são necessárias mais iterações para a completa resolução do sistema.

Aparentemente ambas as funções convergem na quarta iteração, no entanto se dermos um zoom entre a quarta e quinta iteração, vemos que na quarta iteração a função 2 $(x \cdot y + 3 = 0)$ ainda não atingiu a tolerância especificada.

Caso 2: Programação usando fsolve

Para este caso usaremos o toolbox de soluções numéricas do MATLAB.

Dado duas reações que ocorrem em paralelo para gerar um produto C

$\begin{matrix} 2A + B \overset{k_1}{\rightleftharpoons} C\\ A + D \overset{k_2}{\rightleftharpoons} C \end{matrix}$

Onde as constantes de equilíbrio são dadas por:

$\left\{\begin{matrix} k_1 = \dfrac{C_C}{C_A^2\cdot C_B} = 5\times 10^{-4}\\ \\ k_2 = \dfrac{C_C}{C_A\cdot C_D} = 4\times 10^{-2} \end{matrix}\right.$

$C_A$ , $C_B$ , $C_C$ , $C_D$ , são as concentrações do de cada substância no equilíbrio calculadas em função das conversões $x_1$ e $x_2$ pelo seguinte balanço de massa:

$\left\{\begin{matrix} C_A = C_{A,0} -2 \cdot x_1\cdot C_{B,0} - x_1 \cdot C_{D,0}\\ \\ C_B = (1 - x_1)\cdot C_{B,0}\\ \\ C_C = C_{C,0} + x_1 \cdot C_{B,0} + x_2 \cdot C_{D,0}\\ \\ C_D = (1 - x_2)\cdot C_{D,0} \end{matrix}\right.$

Substituindo os balanços de massa nas equações de equilíbrio vamos encontrar:

$\left\{\begin{matrix} \dfrac{C_{C,0} + x_1 \cdot C_{B,0} + x_2 \cdot C_{D,0}}{(C_{A,0} -2 \cdot x_1 \cdot C_{B,0} - x_1 \cdot C_{D,0})^2\cdot (1-x_1)\cdot C_{B,0}} - 5\times 10^{-4} = 0\\ \\ \dfrac{C_{C,0} + x_1 \cdot C_{B,0} + x_2 \cdot C_{D,0}}{(C_{A,0} -2\cdot x_1\cdot C_{B,0} - x_1\cdot C_{D,0})\cdot (1 - x_2)\cdot C_{D,0}} - 4\times 10^{-2} = 0 \end{matrix}\right.$

Sendo as concentrações iniciais representadas por $C_{A,0}$ , $C_{B,0}$ , $C_{C,0}$ , $C_{D,0}$ .

Com isso nos sistema será resolvido para as variáveis $x_1$ e $x_2$ .

Bem! Para o uso do função fsolve, antes é preciso criar uma subrotina para o calculo do sistema, ou seja, ela deve conter as funções, as variáveis e os parâmetros que compõe o sistema. Neste caso a subrotina deve ter o seguinte formato:

function F = fname(X,param₁,param₂, … ,param_n)

A saída F deve ser definido como um vetor de comprimento igual a numero de equações contendo os valores das funções calculada para os valores contidos no vetor X, vetor este de mesmo mesmo tamnho do numero de variaveis do sistema, ou seja do mesmo tamanho de F. Pode ser criado também entradas para os parâmetros do sistema, no nosso caso, C₀ e K.

Segue o código da nossa subrotina:

^?Download SENL.m

function F = SENL(x,C0,K)
% Sistema de equacao não linear:
% Equacao 1: F1 = Cc/(Cb*Ca^2) - K1;
% Equacao 2: F2 = Cc/(Cd*Ca) - K2;
% 2 equacao para 4 variaveis, Ca, Cb, Cc, Cd. No entanto as Concentracoes
% pobem ser calculadas em funcao de conversao x1 e x2;
% Ca = Ca0 - 2*Cb0*x1 - Cd0*x2;
% Cb = Cb0*(1 - x1);
% Cc = Cc0 + Cb0*x1 + Cd0*x2;
% Cd = Cd0*(1 - x2);
% Tendo as concentrcoes de entrada (Ca0, Cb0, Cc0, e Cd0) como parametros
% do sistema. Com isso o sistema agora possui duas variaveis e duas equacoes.
% Entradas:
% Variavel do seu sistema, tem que vim com primerira entrada
%           x - Vetor contendo o valor das variaveis, no caso, duas variaveis
%                x = [x1 x2];
%           C0 - Vetor contendo as concentracoes de entrada. C0 = [Ca0 Cb0 Cc0 Cd0];
%           K - Vetor contendo as contantes de equilibrio da reacao. K =
%                [k1 k2];
% Saida:
%           F - Vetor contendo o valor de cada equacao para as respectivas
%                entradas. F = [F1
%                                    F2]
 
% Separando o vetor C0 em quatro variaveis para melhor compreensao
Ca0 = C0(1);
Cb0 = C0(2);
Cc0 = C0(3);
Cd0 = C0(4);
% Como o sistema esta em funcao das Concentracoes finais, é necessario
% cacula-las antes, pois esta estão em funcao das variaveis x1 e x2.
Ca = Ca0 - 2*Cb0*x(1) - Cd0*x(2);
Cb = Cb0*(1 - x(1));
Cc = Cc0 + Cb0*x(1) + Cd0*x(2);
Cd = Cd0*(1 - x(2));
% Montagem do sistema. As equações devem ser escritas em um vetor coluna,
% onde cada linha representa uma equacao do sistema.
F = [Cc/(Cb*Ca^2) - K(1); % Equacao 1
     Cc/(Cd*Ca) - K(2)];  % Equacao 2

Definida a subrotina vamos “chama-lá” no arquivo principal através da fsolve da seguinte forma:

Xs = fsolve(@(X) fname(X,param₁,param₂, … ,param_n),options)

As entradas da fsolve são a subrotina, para o nosso exemplo SENL, e uma variável do tipo estrutura contendo as opções de resolução (para mais detalhes digite help optimset no MATLAB), que no nosso caso não vamos alterar o default de resolução do fsolve somente a exibição das iterações. Como a subrotina pode ter mais de uma entrada é preciso informar qual deve ser considerado com variável. Isso é feito usando o @(), colocando entre parêntesis a entrada que deve ser assumida como variável. A saída é um vetor Xs contendo a solução encontrada com a mesma ordem do vetor X, ou seja, se $\textbf{X} = [x_1, x_2]$ então $\textbf{Xs} = [xs_1, xs_2]$ .

Segue o código principal para o nosso exemplo:

^?Download mainNR

% Programa principal
% Resolucao de sistema não linear.
 
% Chute inicial. Como no metodo de Newtow, para achar a raiz da equacao é
% necessario uma estimativa inicial da variavel, no caso do sitema vamos
% precisar de um vetor de estimativa inicial contendo as estimaitivas para
% cada variavel do sistema.
x0 = [0.5 0.5];
 
% Vetor de parametros.
K = [5e-4 4e-2];
C0 = [40 15 0 10];
% Foi defido na funcao 'SENL', que os parametros estao dividos em dois vetores, K = [k1 k2]
% e C0 = [Ca0 Cb0 Cc0 Cd0];
 
% Para a resolucao do sitema, basta executar o comando fsolve da seguinte forma:
x = fsolve(@(x0) SENL(x0,C0,K),x0,optimset('Display','iter')); % A saida é dada por um vetor x do mesmo tamnho de x0

Para os valores de C₀ e K, apresentados no código, a solução encontrada foi, $x_1 = 0.1203$ e $x_2 = 0.4787$ o que é coerente, pois como a constante de equilíbrio da segunda reação é maior do que a primeira, era de se esperar que mais produto fosse gerado pela reação 2.

A vantagem do fsolve é o fato de não usar derivação análica o que torna o processamento mais rápido, porém o fato de fazer uma aproximação numérica da derivada pode gerar erros numéricos o que com certeza será propagado para o valor encontrado

O fsolve pode ser usado também para o caso $f(x) = 0$ , basta definir sua subrotina onde, ao invés de entrar um vetor de variáveis, entra o valor da variável, e ao invés de sair um vetor com os valores das funções sai o valor da função.

Fica a sugestão de usar o código simbólico para este segundo caso e comparar os resultados encontrados incluindo o tempo de processamento.

Espero que vocês tenham curtido, e até a próxima!

Categorias: Destaque, Métodos Matemáticos

Tags:fsolve, MATLAB, Newton-Raphson

Sistema não lineares – Newton-Raphson

Deixe seu comentário

Inscreva-se!

Tags

Categorias

Arquivo

Mais de nós...

Recomendamos

Curtiu?