Método de Newton para sistemas não lineares

Em notação vetorial, um sistema não linear é representado por $$ F(x) = 0, $$ onde $F:\R^n\rightarrow \R^n$ e $x=(x_1,\ldots,x_n)^T\in\R^n.$ Assim $F(x) = (f_1(x),\ldots,f_n(x))^T,$ onde cada $f_j$ é uma função de $\R^n$ em $\R.$ Por exemplo, o sistema a seguir $$ \left\{ \begin{array}{l} (x_1-p)^2 + (x_2-q)^2 = 1\\ ax_1 + bx_2 = c \end{array} \right. $$ seria representado por $$ F(x) = \left(\begin{array}{c}(x_1-p)^2 + (x_2-q)^2 - 1 \\\ ax_1 + bx_2 - c \end{array}\right) = 0.$$ Note que o "$0$" na equação acima representa o vetor nulo.

A construção do método de Newton para sistemas não lineares é feita de forma muito semelhante à derivação algébrica realizada na construção do método de Newton para uma única variável. A ideia é construir uma aproximação linear para a função $F$ em torno do vetor $x^{(k)},$ e resolver o problema mais simples de encontrar o zero da aproximação linear. Esse zero é definido como nova aproximação para a solução do sistema não linear.

Para fazer isso, precisamos saber como construir o polinômio de Taylor para uma função vetorial. Vamos começar relembrando como é polinômio de Taylor de primeira ordem para uma função $f:\R^n\rightarrow \R.$ Se $s\in\R^n,$ então \begin{align*} f(x+s) &= f(x) + \frac{\partial f}{\partial x_1}(x) s_1 + \cdots + \frac{\partial f}{\partial x_n}(x) s_n + {\cal O}(\|s\|^2)\\[2mm] & = f(x) + [\nabla f(x)]^Ts + {\cal O}(\|s\|^2), \end{align*} onde ${\cal O}(\|s\|^2)$ representa o erro da aproximação da função $f$ pelo polinômio. Para $F:\R^n\rightarrow \R^n,$ temos que \begin{align*} F(x+s) &= \begin{pmatrix}f_1(x+s)\\ \vdots \\ f_n(x+s)\end{pmatrix}\\[2mm] &= \begin{pmatrix}f_1(x)\\\ \vdots\\\ f_n(x) \end{pmatrix} + \begin{pmatrix}[\nabla f_1(x)]^Ts \\\ \vdots\\\ [\nabla f_n(x)]^Ts \end{pmatrix} + {\cal O}(\|s\|^2)\\[2mm] & = F(x) + J(x)s + {\cal O}(\|s\|^2), \end{align*} onde $J(x)$ é a matriz Jacobina de $F,$ dada por $$ J(x) = \left( \begin{array}{c} \nabla^T f_1(x)\\ \vdots\\ \nabla^T f_n(x) \end{array} \right), $$ isto é, a linha $j$ de $J(x)$ é o gradiente de $f_j(x),$ como um vetor linha. Por exemplo, \begin{equation}\label{exemplo} F(x) = \left( \begin{array}{c} x_1^2 - e^{-x_1x_2}\\ x_1x_2 + \sin x_1 \end{array} \right), \qquad J(x) = \left( \begin{array}{cc} 2x_1 +x_2e^{-x_1x_2} & x_1 e^{-x_1x_2}\\ x_2 + \cos x_1 & x_1 \end{array} \right). \end{equation}

Se $x^{(k)}$ for a aproximação atual para o vetor $x^*,$ solução do sistema não linear, então $x^{(k+1)}$ é definido de modo a ser um zero da aproximação linear para $F$ em torno de $x^{(k)},$ ou seja $x^{(k+1)} = x^{(k)} + s^{(k)}$, onde $s$ satisfaz $$ F(x^{(k)}) + J(x^{(k)})s^{(k)} = 0. $$ Desta forma, $s^{(k)}$ é solução do sistema linear \begin{equation}\label{sl} J(x^{(k)})s^{(k)} = - F(x^{(k)}). \end{equation}

Assim, cada iteração do método de Newton consiste em computar o vetor $F(x^{(k)}),$ a matriz $J(x^{(k)})$, e resolver o sistema linear \eqref{sl} para determinar o vetor $s^{(k)}.$ Feito isso, $x^{(k+1)} = x^{(k)} + s^{(k)}.$ Um protótipo de algoritmo seria:

Como exemplo, podemos aplicar o método de Newton à $F$ definida em \eqref{exemplo}. Isto seria feito com

F = @(x) [x(1)^2-exp(-x(1)*x(2)); x(1)*x(2)+sin(x(1))];
J = @(x) [2*x(1)+x(2)*exp(-x(1)*x(2)), x(1)*exp(-x(1)*x(2));
            x(2)+cos(x(1)), x(1)];

x = [2;1];  # ponto inicial

# resolução do sistema linear e atualização de x
s = J(x)\(-F(x)); x = x + s;

# norma infinito do passo e de F(x)
[norm(s,inf) norm(F(x),inf)]
ans =

   1.20485   0.67601

# 2ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   0.67440   1.52392

# 3ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   0.22960   0.32557

# 4ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   0.095528   0.021026

# 5ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   0.006059981   0.000076539

# 6ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   0.0000246324086   0.0000000011395

# 7ª iteração
s = J(x)\(-F(x)); x = x + s;
[norm(s,inf) norm(F(x),inf)]
ans =

   4.0684e-10   4.4409e-16

Observe que a partir da terceira iteração podemos perceber a taxa de convergência quadrática característica do método de Newton. Note também que ao computar a norma de $s^{(k)}$ ou de $F(x^{(k)}),$ usamos a norma infinito ao invés da norma Euclidiana usual (ou norma 2), isso porque a norma infinito de um vetor, definida como $$ \|x\|_\infty = \max_k \{ |x_k| \}, $$ é mais barata de ser computacionalmente calculada.

Claro que no exemplo acima, o melhor seria ter implementado um laço para repetir as iterações e a definição de um critério de parada. Os critérios de paradas usuais são similares aos empregados no método de Newton para uma variável, baseando-se ou em $\|F(x^{(k)})\|_\infty$ ou em $\|s^{(k)}\|_\infty.$

Assim como em uma variável, quando método de Newton falha se a derivada for nula em algum iterando, aqui a iteração não pode prosseguir se a matriz Jacobina for singular. Nesse caso, não é possível resolver o sistema linear e portanto não há como progredir. A discussão sobre como superar ou mitigar essa dificuldade foge ao escopo deste curso.

Por fim, perceba que a parte mais custosa em cada iteração do método é a resolução do sistema linear. Não estamos explorando este aspecto aqui, mas diferentes implementações do método de Newton podem surgir quando diferentes métodos numéricos são empregados na resolução do sistema linear.

Por falar em sistemas lineares, esse é o próximo tópico deste curso.

Referências

J. E. Dennis e R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM, 1996.

P. Deuflhard. Newton Methods for Nonlinear Problems: Affine Invariance and Adaptative Algorithms. Springer, 2004.

C. T. Kelley. Solving Nonlinear Equations with Newton's Method. SIAM, 2003.

1. Estime a solução dos sistemas não lineares abaixo. Como critério de parada para o método de Newton utilize $\|F(x^{(k)})\| < 10^{-8}.$ Considere o ponto inicial indicado. Tente esboçar graficamente a trajetória que os iterandos computados pelo método de Newton descreveram no plano.

1. $\left\{ \begin{array}{ccc} x_1^2+x_2^2&=&2 \\ e^{x_1-1}+x_2^3&=&2 \end{array} \right.\quad$ com aproximação inicial $x^{(0)} = \left( \begin{array}{c}1.5\\\ 2.0\end{array}\right).$
2. $\left\{ \begin{array}{ccl} y - \sin x &=& 0.5 \\ y^2 - \cos 2x &=&2 \end{array} \right.\quad$ com aproximação inicial $x^{(0)} = \left( \begin{array}{c}3\\\ 1\end{array}\right).$