Método de Newton para sistemas não-lineares ¶ Este método determina, a cada iteração, a solução aproximada do sistema não-linear através de uma linearização das funções-alvo com a matriz Jacobiana associada ao sistema.
Passos ¶ Para o método de Newton não-linear, basicamente criamos uma espécie de “caminho” onde somamos um vetor de deslocamento s {\bf s} s
Obs.: este processo iterativo usa critérios de parada naturais em algoritmos iterativos.
Para encontrarmos o vetor solução, devemos resolver a equação matricial linearizada
J ( x ( i ) ) s ( i ) = − F ( x ( i ) ) {\bf J}({\bf x}^{(i)}){\bf s}^{(i)} = - {\bf F}({\bf x}^{(i)}) J ( x ( i ) ) s ( i ) = − F ( x ( i ) ) Em seguida, atualizamos o novo vetor da sequencia como:
x ( i + 1 ) = x ( i ) + s ( i ) . {\bf x}^{(i+1)}={\bf x}^{(i)} + {\bf s}^{(i)}. x ( i + 1 ) = x ( i ) + s ( i ) . Acima, J ( x ( i ) ) {\bf J}({\bf x}^{(i)}) J ( x ( i ) ) F {\bf F} F
No caso de um sistema em que F = [ f 1 ( x 1 , x 2 ) f 2 ( x 1 , x 2 ) ] T {\bf F} = [f_1(x_1,x_2) \ \ f_2(x_1,x_2)]^T F = [ f 1 ( x 1 , x 2 ) f 2 ( x 1 , x 2 ) ] T
[ ∂ f 1 ( x 1 , x 2 ) ∂ x 1 ∂ f 1 ( x 1 , x 2 ) ∂ x 2 ∂ f 2 ( x 1 , x 2 ) ∂ x 1 ∂ f 2 ( x 1 , x 2 ) ∂ x 2 ] [ s 0 s 1 ] = [ − f 1 ( x 1 , x 2 ) − f 2 ( x 1 , x 2 ) ] \begin{bmatrix}
\frac{\partial f_1(x_1,x_2)}{\partial x_1} & \frac{\partial f_1(x_1,x_2)}{\partial x_2} \\
\frac{\partial f_2(x_1,x_2)}{\partial x_1} & \frac{\partial f_2(x_1,x_2)}{\partial x_2} \\
\end{bmatrix}
\begin{bmatrix}
s_0 \\
s_1 \\
\end{bmatrix}
=
\begin{bmatrix}
-f_1(x_1,x_2)\\
-f_2(x_1,x_2) \\
\end{bmatrix} [ ∂ x 1 ∂ f 1 ( x 1 , x 2 ) ∂ x 1 ∂ f 2 ( x 1 , x 2 ) ∂ x 2 ∂ f 1 ( x 1 , x 2 ) ∂ x 2 ∂ f 2 ( x 1 , x 2 ) ] [ s 0 s 1 ] = [ − f 1 ( x 1 , x 2 ) − f 2 ( x 1 , x 2 ) ] %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import sympy as sy
from scipy.optimize import root
No exemplo a seguir, mostramos como podemos resolver um sistema de equações não-lineares usando o scipy.
Procuramos as soluções para o sistema não-linear
{ f 1 ( x , y ) : x 2 + y 2 = 2 f 2 ( x , y ) : x 2 − y 2 9 = 1 \begin{cases}f_1(x,y) &: x^2 + y^2 = 2 \\
f_2(x,y) &: x^2 - \frac{y^2}{9} = 1\end{cases} { f 1 ( x , y ) f 2 ( x , y ) : x 2 + y 2 = 2 : x 2 − 9 y 2 = 1 Vamos plotar o gráficos das funções:
x = np.linspace(-2,2,50,endpoint=True)
y = x[:]
X,Y = np.meshgrid(x,y)
# define funções para plotagem
F1 = X**2 + Y**2 - 2
F2 = X**2 - Y**2/9 - 1
# curvas de nível
C = plt.contour(X,Y,F1,levels=[0],colors='k')
C = plt.contour(X,Y,F2,levels=[0],colors='g')
plt.grid(True)Pela figura, vemos que existem 4 pontos de interseção entre as curvas e, portanto, 4 soluções, as quais formam o conjunto
S = { ( x 1 ∗ , y 1 ∗ ) , ( x 2 ∗ , y 2 ∗ ) , ( x 3 ∗ , y 3 ∗ ) , ( x 4 ∗ , y 4 ∗ ) ) S = \{(x_1^{*},y_1^{*}),(x_2^{*},y_2^{*}),(x_3^{*},y_3^{*}),(x_4^{*},y_4^{*})) S = {( x 1 ∗ , y 1 ∗ ) , ( x 2 ∗ , y 2 ∗ ) , ( x 3 ∗ , y 3 ∗ ) , ( x 4 ∗ , y 4 ∗ )) Agora, vamos usar a função root do scipy para computar essas soluções com base em estimativas iniciais.
# define função para o vetor F(x)
def F(x):
return [ x[0]**2 + x[1]**2 - 2,
x[0]**2 - x[1]**2/9 - 1 ]
x,y = sy.symbols('x,y')
# usa computação simbólica para determinar a matriz Jacobiana
f1 = x**2 + y**2 - 2
f2 = x**2 - y**2/9 - 1
# gradientes
f1x,f1y = sy.diff(f1,x),sy.diff(f1,y)
f2x,f2y = sy.diff(f2,x),sy.diff(f2,y)
# imprime derivadas parciais
print(f1x)
print(f1y)
print(f2x)
print(f2y)
# monta matriz Jacobiana
def jacobian(x):
return np.array([[2*x[0], 2*x[1]], [2*x[0],-2*x[1]/9]])
# resolve o sistema não-linear por algoritmo de Levenberg-Marqardt modificado
inicial = [[2,2],[-2,2],[-2,-2],[2,-2]]
S = []
i = 1
for vetor in inicial:
aux = root(F,vetor,jac=jacobian, method='lm')
S.append(aux.x)
s = 'Solução x({0})* encontrada: {1}'
print(s.format(i,aux.x))
i +=12*x
2*y
2*x
-2*y/9
Solução x(1)* encontrada: [1.04880885 0.9486833 ]
Solução x(2)* encontrada: [-1.04880885 0.9486833 ]
Solução x(3)* encontrada: [-1.04880885 -0.9486833 ]
Solução x(4)* encontrada: [ 1.04880885 -0.9486833 ]
Em seguida, vamos plotar as soluções e as curvas
# curvas de nível
C = plt.contour(X,Y,F1,levels=[0],colors='k')
C = plt.contour(X,Y,F2,levels=[0],colors='g')
plt.grid(True)
# imprime interseções
for i in range(len(S)):
plt.plot(S[i][0],S[i][1],'or')
Exercício: ¶ Resolva os sistemas não-lineares da Lista de Exercícios 4 usando a mesma abordagem acima.
Nota: Raízes de sistemas não-lineares ¶ f ( x ) = a 1 x 1 + a 2 x 2 + … + a n x n f(x) = a_1 x_1 + a_2 x_2 + \ldots + a_n x_n f ( x ) = a 1 x 1 + a 2 x 2 + … + a n x n f 2 ( x ) = a 1 x 1 x 2 + a 2 x 2 2 + a n x n x 1 f_2(x) = a_1 x_1 x_2 + a_2 x_2^2 + a_n x_nx_1 f 2 ( x ) = a 1 x 1 x 2 + a 2 x 2 2 + a n x n x 1 f 1 ( x 1 , x 2 , … , x n ) = 0 f_1(x_1,x_2,\ldots,x_n) = 0 f 1 ( x 1 , x 2 , … , x n ) = 0 f 2 ( x 1 , x 2 , … , x n ) = 0 f_2(x_1,x_2,\ldots,x_n) = 0 f 2 ( x 1 , x 2 , … , x n ) = 0 f n ( x 1 , x 2 , … , x n ) = 0 f_n(x_1,x_2,\ldots,x_n) = 0 f n ( x 1 , x 2 , … , x n ) = 0 A solução do sistema é o vetor ( x 1 ∗ , x 2 ∗ , … , x n ∗ ) (x_1^{*},x_2^{*},\ldots,x_n^{*}) ( x 1 ∗ , x 2 ∗ , … , x n ∗ ) n n n
Iteração de Ponto Fixo para sistemas não-lineares ¶ x 1 i + 1 = f ~ 1 ( x 1 i , x 2 i , … , x n i ) x_1^{i+1} = \tilde{f}_1(x_1^{i},x_2^{i},\ldots,x_n^{i}) x 1 i + 1 = f ~ 1 ( x 1 i , x 2 i , … , x n i ) x 2 i + 1 = f ~ 2 ( x 1 i , x 2 i , … , x n i ) x_2^{i+1} = \tilde{f}_2(x_1^{i},x_2^{i},\ldots,x_n^{i}) x 2 i + 1 = f ~ 2 ( x 1 i , x 2 i , … , x n i ) x 3 i + 1 = f ~ 3 ( x 1 i , x 2 i , … , x n i ) x_3^{i+1} = \tilde{f}_3(x_1^{i},x_2^{i},\ldots,x_n^{i}) x 3 i + 1 = f ~ 3 ( x 1 i , x 2 i , … , x n i ) As formas funcionais mudam porque devemos isolar a variável x i x_i x i
Exemplo: encontrar a raiz do sistema abaixo:
f 1 ( x , y ) = x 2 + x y − 10 = 0 f_1(x,y) = x^2 + xy - 10 = 0 f 1 ( x , y ) = x 2 + x y − 10 = 0 f 2 ( x , y ) = y + 3 x y 2 − 57 = 0 f_2(x,y) = y + 3xy^2 - 57 = 0 f 2 ( x , y ) = y + 3 x y 2 − 57 = 0 Solução:
Reescrevamos as equações na forma
x = f ~ 1 ( x , y ) = 10 − x y x = \tilde{f}_1(x,y) = \sqrt{10 - xy} x = f ~ 1 ( x , y ) = 10 − x y y = f ~ 2 ( x , y ) = 57 − y 3 x y = \tilde{f}_2(x,y) = \sqrt{\frac{57-y}{3x}} y = f ~ 2 ( x , y ) = 3 x 57 − y de onde temos a iteração de ponto fixo dada por
x i + 1 = 10 − x i y i x^{i+1} = \sqrt{10 - x^iy^i} x i + 1 = 10 − x i y i y i + 1 = 57 − y i 3 x i , i = 0 , 1 , 2 , … , y^{i+1} = \sqrt{\frac{57-y^i}{3x^i}}, \quad i = 0,1,2,\ldots, y i + 1 = 3 x i 57 − y i , i = 0 , 1 , 2 , … , Usando ( x 0 , y 0 ) = ( 1.5 , 3.5 ) (x^0,y^0) = (1.5,3.5) ( x 0 , y 0 ) = ( 1.5 , 3.5 )
x 1 = 10 − x 0 y 0 = 10 − 1.5 ( 3.5 ) = 2.17945 x^{1} = \sqrt{10 - x^0y^0} = \sqrt{10 - 1.5(3.5)} = 2.17945 x 1 = 10 − x 0 y 0 = 10 − 1.5 ( 3.5 ) = 2.17945 y 1 = 57 − y 0 3 x 1 = 57 − 3.5 3 ( 2.17945 ) = 2.86051 y^{1} = \sqrt{\frac{57-y^0}{3x^1}} = \sqrt{\frac{57-3.5}{3(2.17945)}} = 2.86051 y 1 = 3 x 1 57 − y 0 = 3 ( 2.17945 ) 57 − 3.5 = 2.86051 (o valor de x 1 x^1 x 1 x 0 x^0 x 0
x 2 = 10 − x 1 y 1 = 10 − 2.17945 ( 2.86051 ) = 1.94053 x^{2} = \sqrt{10 - x^1y^1} = \sqrt{10 - 2.17945(2.86051)} = 1.94053 x 2 = 10 − x 1 y 1 = 10 − 2.17945 ( 2.86051 ) = 1.94053 y 2 = 57 − y 1 3 x 2 = 57 − 2.86051 3 ( 1.94053 ) = 3.04955 y^{2} = \sqrt{\frac{57-y^1}{3x^2}} = \sqrt{\frac{57-2.86051}{3(1.94053)}} = 3.04955 y 2 = 3 x 2 57 − y 1 = 3 ( 1.94053 ) 57 − 2.86051 = 3.04955 O processo iterativo converge para a solução ( x ∗ , y ∗ ) = ( 2 , 3 ) (x^{*},y^{*}) = (2,3) ( x ∗ , y ∗ ) = ( 2 , 3 )
Notas:
A convergência por iteração de PF depende de como as equações f ~ 1 , f ~ 2 , … , f ~ n \tilde{f}_1,\tilde{f}_2,\ldots,\tilde{f}_n f ~ 1 , f ~ 2 , … , f ~ n
A iteração de PF é bastante restritiva nas soluções de sistemas não-lineares.
Newton-Raphson para sistema não-linear ¶ Depende de série de Taylor em n − n- n −
Para 2 dimensões, por exemplo, o método de Newton-Raphson pode ser escrito como:
u i + 1 = u i + ( x i + 1 − x i ) ∂ u i ∂ x + ( y i + 1 − y i ) ∂ u i ∂ y u_{i+1} = u_i + (x_{i+1} - x_i)\frac{\partial u_i}{\partial x} +(y_{i+1} - y_i)\frac{\partial u_i}{\partial y} u i + 1 = u i + ( x i + 1 − x i ) ∂ x ∂ u i + ( y i + 1 − y i ) ∂ y ∂ u i v i + 1 = v i + ( x i + 1 − x i ) ∂ v i ∂ x + ( y i + 1 − y i ) ∂ v i ∂ y v_{i+1} = v_i + (x_{i+1} - x_i)\frac{\partial v_i}{\partial x} +(y_{i+1} - y_i)\frac{\partial v_i}{\partial y} v i + 1 = v i + ( x i + 1 − x i ) ∂ x ∂ v i + ( y i + 1 − y i ) ∂ y ∂ v i ∂ u i ∂ x x i + 1 + ∂ u i ∂ y y i + 1 = − u i + ∂ u i ∂ x x i + ∂ u i ∂ y y i \frac{\partial u_i}{\partial x} x_{i+1} +
\frac{\partial u_i}{\partial y} y_{i+1} =
- u_i + \frac{\partial u_i}{\partial x}x_i
+ \frac{\partial u_i}{\partial y}y_i ∂ x ∂ u i x i + 1 + ∂ y ∂ u i y i + 1 = − u i + ∂ x ∂ u i x i + ∂ y ∂ u i y i ∂ v i ∂ x x i + 1 + ∂ v i ∂ y y i + 1 = − v i + ∂ v i ∂ x x i + ∂ v i ∂ y y i , \frac{\partial v_i}{\partial x} x_{i+1} +
\frac{\partial v_i}{\partial y} y_{i+1} =
- v_i + \frac{\partial v_i}{\partial x}x_i
+ \frac{\partial v_i}{\partial y}y_i, ∂ x ∂ v i x i + 1 + ∂ y ∂ v i y i + 1 = − v i + ∂ x ∂ v i x i + ∂ y ∂ v i y i , que é um sistema nas incógnitas x i + 1 x_{i+1} x i + 1 y i + 1 y_{i+1} y i + 1
x i + 1 = x i − J − 1 ( u i ∂ v i ∂ y − v i ∂ u i ∂ y ) x_{i+1} = x_i - J^{-1}\left(u_i\frac{\partial v_i}{\partial y} - v_i\frac{\partial u_i}{\partial y}\right) x i + 1 = x i − J − 1 ( u i ∂ y ∂ v i − v i ∂ y ∂ u i ) y i + 1 = y i − J − 1 ( v i ∂ u i ∂ x − u i ∂ v i ∂ x ) , y_{i+1} = y_i - J^{-1}\left(v_i\frac{\partial u_i}{\partial x} - u_i\frac{\partial v_i}{\partial x}\right), y i + 1 = y i − J − 1 ( v i ∂ x ∂ u i − u i ∂ x ∂ v i ) , onde
J = ∂ u i ∂ x ∂ v i ∂ y − ∂ u i ∂ y ∂ v i ∂ x J = \frac{\partial u_i}{\partial x}\frac{\partial v_i}{\partial y} - \frac{\partial u_i}{\partial y}\frac{\partial v_i}{\partial x} J = ∂ x ∂ u i ∂ y ∂ v i − ∂ y ∂ u i ∂ x ∂ v i é o determinante da matriz Jacobiana do sistema.
Exemplo: resolver o mesmo sistema do PF
Solução:
∂ u 0 ∂ x = 2 x 0 + y 0 = 2 ( 1.5 ) + 3.5 = 6.5 , ∂ u 0 ∂ y = x 0 = 1.5 \frac{\partial u_0}{\partial x} = 2x_0 + y_0 = 2(1.5) + 3.5 = 6.5, \ \ \ \ \frac{\partial u_0}{\partial y} = x_0 = 1.5 ∂ x ∂ u 0 = 2 x 0 + y 0 = 2 ( 1.5 ) + 3.5 = 6.5 , ∂ y ∂ u 0 = x 0 = 1.5 ∂ v 0 ∂ x = 3 y 0 2 = 3 ( 3.5 ) 2 = 36.75 , ∂ v 0 ∂ y = 1 + 6 x 0 y 0 = 1 + 6 ( 1.5 ) ( 3.5 ) = 32.5 \frac{\partial v_0}{\partial x} = 3y_0^2 = 3(3.5)^2 = 36.75, \ \ \ \ \frac{\partial v_0}{\partial y} = 1 + 6x_0y_0 = 1 + 6(1.5)(3.5) = 32.5 ∂ x ∂ v 0 = 3 y 0 2 = 3 ( 3.5 ) 2 = 36.75 , ∂ y ∂ v 0 = 1 + 6 x 0 y 0 = 1 + 6 ( 1.5 ) ( 3.5 ) = 32.5 u ( x 0 , y 0 ) = u 0 = − 2.5 , v ( x 0 , y 0 ) = v 0 = 1.625 u(x_0,y_0) = u_0 = -2.5, \quad v(x_0,y_0) = v_0 = 1.625 u ( x 0 , y 0 ) = u 0 = − 2.5 , v ( x 0 , y 0 ) = v 0 = 1.625 Calculamos os valores no próximo passo, i.e., ( x 1 , y 1 ) (x^1,y^1) ( x 1 , y 1 )
x 1 = 1.5 − ⋯ 156.125 = 2.03603 x^1 = 1.5 - \frac{\cdots}{156.125} = 2.03603 x 1 = 1.5 − 156.125 ⋯ = 2.03603 y 1 = 3.5 − ⋯ 156.125 = 2.84388 y^1 = 3.5 - \frac{\cdots}{156.125} = 2.84388 y 1 = 3.5 − 156.125 ⋯ = 2.84388 O processo iterativo converge para a solução ( x ∗ , y ∗ ) = ( 2 , 3 ) (x^{*},y^{*}) = (2,3) ( x ∗ , y ∗ ) = ( 2 , 3 )
Ref.: Chapra & Canale, sec. 6.6
