Prof. Doherty Andrade – www.metodosnumericos.com.br
- Vamos estender o método de Newton-Raphson para sistemas de equações não lineares. Deduzimos um caso particular importante do método de
Newton para sistemas de duas equações.
2 .Seja $\Omega \subset \mathbb{R}^n$ um aberto e $F:\Omega
\rightarrow \mathbb{R}^n$ uma função de classe $C^2.$ Suponha que
exista $X^* \in \Omega$ tal que $F(X^*)=0$.
Queremos determinar $X^*$ isto é, queremos resolver numericamente a equação vetorial $F(X)=0,$ onde $X=(x_1, x_2, \ldots, x_n)$.
Como $F(X) \in \mathbb{R}^n$, então $F$ é da forma $F= (f_1,f_2,
\ldots, f_n)$, onde $f_i:\Omega \rightarrow
\mathbb{R},i=1,2,\ldots,n.$ Logo, a equação $F(X)=0$ é equivalente
ao sistema de equações não lineares
$$
\begin{cases}
f_1(x_1,x_2, \ldots, x_n)=0,\cr f_2(x_1,x_2, \ldots, x_n)=0,\cr
\vdots \qquad \qquad \qquad \qquad \vdots \qquad \cr
f_n(x_1,x_2, \ldots, x_n)=0.
\end{cases}$$
Exemplo: Considere $F(x,y)= (x^2+y^2-2,x^2- \dfrac{y^2}{9}-1)$ definida em
todo $\mathbb{R}^2$. A equação $F(x,y)=0$ pode ser reescrita como
$$
\begin{cases}
x^2+y^2-2=0,\cr x^2- \dfrac{y^2}{9}-1=0.
\end{cases}$$
Geometricamente, resolver esse sistema de equações não lineares,
significa determinar onde a circunferência e a parábola de
intersectam. Veja a figura.
Note que resolver um sistema de equações lineares $AX=b$ é equivalente a
resolver uma equação $F(X)=0$, onde $F(X)= AX-b.$ Os sistema
lineares são um caso particular importante, como sabemos, existem
vários métodos diretos e iterativos para obter numericamente a
solução. No caso geral, a não linearidade de $F$ é um fator
complicador para resolver a equação $F(X)=0$.
Lembramos que a derivada de $F: \Omega \subset
\mathbb{R}^n \rightarrow \mathbb{R}^n$ dada por
$$F(x_1,x_2,\ldots, x_n)= (f_1(x_1,x_2,\ldots, x_n),
f_2(x_1,x_2,\ldots, x_n), \ldots, f_n(x_1,x_2,\ldots, x_n))$$ é a
matriz Jacobiana de $F$ que é
$$JF(X)=
\left[
\begin{array}{cccc}
\displaystyle\frac{\partial f_1}{\partial x_1} & \displaystyle\frac{\partial f_1}{\partial x_2} & \ldots & \displaystyle\frac{\partial f_1}{\partial x_n} \\
\displaystyle\frac{\partial f_2}{\partial x_1} & \displaystyle\frac{\partial f_2}{\partial x_2} & \ldots & \displaystyle\frac{\partial f_2}{\partial x_n} \\
\ldots& \ldots & \ldots & \ldots \\
\displaystyle\frac{\partial f_n}{\partial x_1} & \displaystyle\frac{\partial f_n}{\partial x_2} & \ldots & \displaystyle\frac{\partial f_n}{\partial x_n} \\
\end{array}
\right].$$
Considerando que na equação $F(X)=0$, $F$ tem derivadas até ordem
2 contínuas, tomando a expansão em série de Taylor de $F(X)$ em
torno do ponto $X=X^{(k)}$ temos
$$F(X)= F(X^{(k)})+ F^\prime (X^{(k)}) (X-X^{(k)}) +
\mbox{resto}.$$
Como $F(X)=0$, então temos
$$F^\prime (X^{(k)}) (X-X^{(k)})=- F(X^{(k)})+
\mbox{resto}.$$
Se a derivada admite inversa, então aplicando-a a ambos os lados,
e supondo que os pontos estão suficientemente próximos, podemos
desprezar o resto, obtendo
\begin{equation}X= X^{(k)}-\left[ F^\prime (X^{(k)}) \right]^{-1} F(X^{(k)})\end{equation}
uma aproximação para a raiz da equação.
Assim, obtemos uma aproximação $X=X^{(k+1)}$ para a solução dada
por
\begin{equation}\label{isso}
\displaystyle X^{(k+1)}= X^{(k)}-\left[ F^\prime (X^{(k)})
\right]^{-1} F(X^{(k)})\end{equation}
Na prática evitamos a equação \ref{isso}, pois esta escolha tem a
necessidade de calcular a inversa de $F^\prime (X^{(k)})$.
A abordagem alternativa é considerar apenas
\begin{equation}\label{mnrn1}
\displaystyle F^\prime (X^{(k)}) (X^{(k+1)}-X^{(k)})= –
F(X^{(k)}),\,\,\,\, k \geq 0,\end{equation} onde $X^{(0)}$ é
uma aproximação inicial dada.
Observe que em cada iteração resolve-se um sistema de equações
lineares. Este é o método de Newton-Raphson para sistemas de
equações não lineares.
Ao utilizarmos o método de Newton-Raphson para
sistemas não lineares, se usarmos um método i-te-ra-ti-vo
para resolver o sistema linear, então o método de Newton-Raphson é
chamado inexato.
Tomando $F^\prime (X^{(k)})=B$ constante no método de
Newton-Raphson, o método recebe o nome de método de Newton
modificado.
As condições para convergência do método de
Newton-Raphson são as
seguintes:
(1i) $F$ e suas derivadas parciais até ordem 2 são contínuas na
região $\Omega$ contendo a raiz,
(2i) O determinante da matriz Jacobiana não se anula na região
$\Omega$,
(3i) A aproximação inicial $X^{(0)}$ está suficientemente próxima
da raiz.
Como no caso de uma variável, a ordem de convergência
desse método é 2.
Os critérios de parada são os mesmos para o caso de
uma variável:
(1i) $\Vert F(X^{(k)})\Vert < \varepsilon_1,$
(2i) $\Vert X^{(k+1)}-X^{(k)} \Vert <\varepsilon_2,$ onde
$\varepsilon_1$ e $\varepsilon_2$ são dados.
Exemplo: No sistema
$$\begin{cases}
4x_1-x_1^3+x_2=0,\cr -\dfrac{x_1^2}{9}+ \dfrac{4x_2-x_2^2}{4}+1=0
\end{cases}$$
Veja a figura com a ilustração das curvas.
\section{Caso Particular: $F=(f,g)$}
Consideremos $F: \Omega \subset \mathbb{R}^2 \rightarrow
\mathbb{R}^2$ dada por $F(x,y)=(f(x,y),g(x,y))$. Aplicando a
expressão equação \ref{mnrn1}, obtemos
\begin{eqnarray}
\left[
\begin{array}{cc}
\displaystyle \frac{\partial f(x^{(k)},y^{(k)})}{\partial x} & \displaystyle\frac{\partial f(x^{(k)},y^{(k)})}{\partial y} \\
\displaystyle\frac{\partial g(x^{(k)},y^{(k)})}{\partial x} &\displaystyle \frac{\partial g(x^{(k)},y^{(k)})}{\partial y} \\ \end{array}
\right].\left[
\begin{array}{c}
\displaystyle x^{(k+1)}-x^{(k)} \\
\displaystyle y^{(k+1)}-y^{(k)}\\
\end{array}
\right] = -\left[
\begin{array}{c}
\displaystyle f(x^{(k)},y^{(k)}) \\
\displaystyle g(x^{(k)},y^{(k)}) \\
\end{array}
\right]
\end{eqnarray}
que é um sistema de equações lineares de incógnitas $ x^{(k+1)}$ e
$y^{(k+1)}$.
Esse sistema pode ser resolvido usando a Regra de Cramer. Para
isto, deve-se supor que o determinante da matriz Jacobiana não se
anula na região de busca. Neste caso tem-se (usando a regra de
Cramer),
\begin{eqnarray}
x^{(k+1)}-x^{(k)}=
\frac{(gf_y-fg_y)(x^{(k)},y^{(k)})}{(f_xg_y-g_xf_y)(x^{(k)},y^{(k)})},\,\,k
\geq 0\\
y^{(k+1)}-y^{(k)}=
\frac{(fg_x-gf_x)(x^{(k)},y^{(k)})}{(f_xg_y-g_xf_y)(x^{(k)},y^{(k)})},\,\,k
\geq 0
\end{eqnarray}
Ou equivalentemente,
\begin{equation}\displaystyle
x^{(k+1)}=x^{(k)}-
\left[\frac{fg_y-gf_y}{f_xg_y-g_xf_y}\right]
,\,\,k
\geq 0
\end{equation}
\begin{equation}\displaystyle
y^{(k+1)}=y^{(k)}-
\left[\frac{gf_x-fg_x}{f_xg_y-g_xf_y}\right],\,\,k
\geq 0
\end{equation}
em que as funções entre colchetes são aplicadas em
$(x^{(k)},y^{(k)}) $.
Uma das aplicações mais importantes para sistemas de equações não
lineares a duas dimensões consiste em determinar raízes complexas
de uma equação $f(z)=0$. Note que denotando $z=x+iy$ podemos
escrever $f(x,y)= u(x,y)+iv(x,y)=0$ se, e somente se,
$$\begin{cases} u(x,y)=0\cr
v(x,y)=0.\end{cases}
$$
Exemplo: No sistema
$$\begin{cases}
4x_1-x_1^3+x_2=0,\cr -\dfrac{x_1^2}{9}+ \dfrac{4x_2-x_2^2}{4}+1=0
\end{cases}$$
[ 3.4375 -1.4375] [ 2.56707296 -0.50535394] [ 2.11565442 -0.47031686] [ 1.9544165 -0.51310999] [ 1.93218424 -0.51813115] [ 1.93177041 -0.51822365] [ 1.93177027 -0.51822368] [ 1.93177027 -0.51822368] [ 1.93177027 -0.51822368] [ 1.93177027 -0.51822368] Obtemos assim, a aproximação para uma solução $(1.93177027,-0.51822368).$ Veja o script in Python que usamos para resolver o problema acima.
Você vai precisar dos pacotes:
from future import division
import numpy as np
from numpy import linalg