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

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

A solução de sistemas de equações lineares é relativamente fácil e podemos tentar diversas abordagens para encontrar uma solução. Sistemas de equações não lineares são bem mais complicados de resolver. Por exemplo, o sistema $$\begin{array}{rc} x^2 + y^3 - 2 & = 0 \\ -x^2 + y & = 0 \end{array}$$ é exemplo de um sistema não linear. Para este caso específico, como temos apenas duas funções polinomiais, podemos até calcular uma solução analítica, mas, de modo geral, isso não é possível. Então, devemos recorrer a algum método numérico iterativo. Para estes casos, talvez o método mais popular (talvez "popular" não seja o termo mais adequado aqui) seja o método de Newton-Raphson. Nesta postagem iremos fazer uma pequena variação para evitar o cálculo do Jacobiano (matriz de derivadas parciais do sistema) de forma analítica, lembrando do fato: $$ \frac{df}{dx} \cong \frac{f(x+h) - f(x-h)}{2h} $$ Formalizando um pouco. Seja um sistema de equações na forma \[ \begin{array}{cc} f_1(x_1, x_2, \ldots x_n) & = 0 \\ f_2(x_1, x_2, \ldots x_n) & = 0 \\ f_3(x_1, x_2, \ldots x_n) & = 0 \\ \vdots & \vdots \\ f_n(x_1, x_2, \ldots x_n) & = 0 \\ \end{array} \] A sua solução é da forma: \begin{equation} \left[\begin{array}{c} x^{k+1}_1\\ x^{k+1}_2\\ \vdots\\ x^{k+1}_n \end{array}\right] = \left[\begin{array}{c} x^{k}_1\\ x^{k}_2\\ \vdots\\ x^{k}_n \end{array}\right] - \left[\begin{array}{cccc} \frac{\partial f_1}{\partial x^k_1} & \frac{\partial f_1}{\partial x^k_2} & \cdots & \frac{\partial f_1}{\partial x^k_n}\\ \frac{\partial f_2}{\partial x^k_1} & \frac{\partial f_2}{\partial x^k_2} & \cdots & \frac{\partial f_2}{\partial x^k_n}\\ \vdots & \vdots & \vdots & \vdots \\ \frac{\partial f_n}{\partial x^k_1} & \frac{\partial f_n}{\partial x^k_2} & \cdots & \frac{\partial f_n}{\partial x^k_n}\\ \end{array}\right]^{-1} \left[\begin{array}{c} f_1(x^k_1, \ldots, x^k_n) \\ f_2(x^k_1, \ldots, x^k_n) \\ \vdots \\ f_n(x^k_1, \ldots, x^k_n) \\ \end{array}\right] \end{equation} com x$_0$ sendo um chute inicial. Se esse chute inicial for próximo o suficiente, a solução deve convergir. Por exemplo, vamos considerar o sistema de equações não lineares $$ \begin{array}{rr} \cos(\pi x_1/2) + \ln(x_1 x_2 x_3) + x^2_3-10.8&=0\\ x_1 x_2 x_3 + \sqrt{x_1 x_3 + x_2} + \cos(\pi x_3/2) - 8.2 &=0\\ x_1 x_3 + \exp(-x_1 x_2) - \cos(\pi x_2 x_3/6) - 4.2 &=0\\ \end{array} $$ com um chute inicial igual a x$_0$ = [0,5; 1,0; 1,5]$^T$.

O código Scilab abaixo resolve esse sistema:

// função 1: entrada: "x1", "x2" e "x3":
function f1 = fx1(x1,x2,x3)
    f1 = cos(%pi*x1/2) + log(x1*x2*x3) +x3*x3-10.8; //round(5*(cos(%pi*1/2) + log(1*2*3) +3*3))/5; //(cos(%pi*1/2) + log(2*3) +3*3);
endfunction

// função 1: entrada: "x1", "x2" e "x3":
function f2 = fx2(x1,x2,x3)
  f2 = x1*x2*x3 + sqrt(x1*x3+x2) +cos(%pi*x3/2)-8.2; //round(5*(1*2*3 + sqrt(1*3+2) +cos(%pi*3/2)))/5; //(1*2*3 + sqrt(1*3) +cos(%pi*3/2));
endfunction

// função 1: entrada: "x1", "x2" e "x3":
function f3 = fx3(x1,x2,x3)
  f3 = x1*x3 + exp(-x1*x2) -cos(%pi*x2*x3/6)-4.2; //round(5*(1*3 + exp(-1*2) -cos(%pi*2*3/6)))/5; //(1*3 + exp(-1*2) -cos(%pi*2*3/6));
endfunction

clc;

// valores iniciais e passo h:
x1 = 0.5; x2 = 1; x3 = 1.5; h = 0.0001; h2 = 2*h;

// vetor inicial, erro inicial e contador do laço
x = [x1; x2; x3]; er = 1; k=0;
while abs(er)>0.0001
     // Jacobiano:
      J = [(fx1(x1+h,x2,x3)-fx1(x1-h,x2,x3))/h2, (fx1(x1,x2+h,x3)-fx1(x1,x2-h,x3))/h2, (fx1(x1,x2,x3+h)-fx1(x1,x2,x3-h))/h2
            (fx2(x1+h,x2,x3)-fx2(x1-h,x2,x3))/h2, (fx2(x1,x2+h,x3)-fx2(x1,x2-h,x3))/h2, (fx2(x1,x2,x3+h)-fx2(x1,x2,x3-h))/h2
            (fx3(x1+h,x2,x3)-fx3(x1-h,x2,x3))/h2, (fx3(x1,x2+h,x3)-fx3(x1,x2-h,x3))/h2, (fx3(x1,x2,x3+h)-fx3(x1,x2,x3-h))/h2];
      Ji = inv(J);  // inversa do Jacobiano
      Fx = [fx1(x1,x2,x3); fx2(x1,x2,x3); fx3(x1,x2,x3)];
      e = - Ji*Fx;  // erro
      x = x + e;   // atualizando x
      x1 = x(1);
      x2 = x(2);
      x3 = x(3);
      disp(x);  // mostrando o resultado
      er = max(abs(e));
      k=k+1;
      if k>10 then er = 0; end; // teste para sair do laço
end;
disp([k, er, x1, x2, x3]);  // resultado final

A resposta é:
    2.8279168
    3.187974 
    2.5561181

    0.3258810
    2.0366685
    3.7243397

    0.8161449
    1.9556474
    2.9552854

    1.0297346
    1.9228066
    3.0044189

    1.0189457
    1.9414167
    3.0077076

    1.0191206
    1.9412045
    3.0077546

    1.0191206
    1.9412045
    3.0077546

Solução:
    7.    1.163D-08    1.0191206    1.9412045    3.0077546