Solução de sistemas de equações não lineares: método iterativo.

Em uma postagem anterior (ver aqui) mostramos como resolver um sistema de equações não lineares usando o método de Newton-Raphson. Nesta postagem mostraremos uma forma mais simples que, se aplicável, também pode levar à solução do sistema. Por exemplo, para o sistema \begin{align*} x^3 + xy & = 3 \\ \sqrt{x} - y^2 & = 5 \end{align*}
Podemos escrever as seguintes equações acopladas:
\begin{align*} x_{k+1} & = \sqrt[3]{3-x_k y_k} \\ y_{k+1} & = \sqrt{5 - \sqrt{x_{k+1}}} \end{align*}
Com o chute inicial $x_0 = 0,5$ e $y_0 = 0,5$, obtemos a seguinte solução:

 1.401    1.954
 0.641    2.049
 1.190    1.977
 0.865    2.017
 1.079    1.99
 0.948    2.007
 1.031    1.996
 0.980    2.003
 1.012    1.998

Que tende a convergir para $x_f = 1$ e $y_f = 2$ depois de mais algumas iterações. Esse método é, certamente, mais lento que o método de Newton, mas é mais simples de ser usado. Exemplo 2 e código Scilab em seguida. Seja
\begin{align*} x^4 + 3xy - z^2 & = 8,00000 \\ x + y^2 + ln(z) & = 2,69315 \\ 2x - xy + e^z & = 8,38906 \end{align*}
Convergência da solução:
1.62658    1.32654    1.987
0.88206 + 0.88206i    1.12996 - 0.39031i    2.08482 - 0.13866i
0.98112 - 0.47983i    1.02295 + 0.26699i    2.02733 + 0.09638i
0.86665 + 0.1157i    1.06046 - 0.07695i    2.02627 - 0.02312i
1.02708 - 0.01715i    0.97981 + 0.01458i    1.99355 + 0.00442i
1.00151 - 0.00302i    1.00086 + 0.00040i    1.99991 + 0.00046i
0.99831 + 0.00151i    1.00087 - 0.00087i    2.00035 - 0.00032i
1.00027 - 0.00016i    0.99978 + 0.00016i    1.99993 + 0.00004i
1.00003 - 0.00004i    1. + 0.00001i    2. + 7.4D-06i
0.99998 + 0.00002i    1.00001 - 0.00001i    2.00001 - 3.7D-06i
1. - 1.3D-06i    1. + 1.6D-06i    2. + 3.9D-07i
1.00000 - 6.0D-07i    1. + 2.0D-07i    2. + 1.1D-07i

Depois de 13 iterações, o sistema converge para solução correta: $x_f = 1$, $y_f = 1$ e $z_f = 2$ com um erro pequeno (com direito a valores complexos). Código Scilab é bem simples:

---

clc;
x = 0.5; y = 0.5; z = 0.5;
v1 = 8; v2 = 2.69315; v3 = 8.38906;
for k=1:12
    x = (v1 - z*z - 3*x*y)^(1/4);
    y = sqrt(v2 - log(z) - x);
    z = log(v3 - 2*x + x*y)
    disp([x,y,z]);
end