By Newton method find the first-step approximation for the root, using the initial approximation $x_0 = y_0 = z_0 = 0.5$:
\begin{equation} \left\{\begin{matrix} x^2+y^2+z^2=0\\ 2x^2+y^2-4z=0\\ 3x^2+-4y+z^2=0 \end{matrix}\right. \end{equation}Make ten* steps of iteration method and find approximate value of the root.
Fonte: https://www.ime.usp.br/~yambar/MAE5704/Aula6Optimization-I/
Sobre o método: https://en.wikipedia.org/wiki/Newton%27s_method
O método de Newton-Raphson é um método utilizado tanto para a resolução de equações (sistemas) não lineares e otimização de funções suaves. Para resolvermos um sistema não linear de equações vamos aplicar o método Newton-Raphson. Considere a transformação de interesse $\Phi:\mathbb{R}^3 \rightarrow \mathbb{R}^3$ e o vetor $\mathbf{u}=(x, y,z)^T$ como sendo a variável do sistema abaixo:
\begin{equation} \Phi(\mathbf{u})=\begin{pmatrix} \Phi_1(\mathbf{u})\\ \Phi_2(\mathbf{u})\\ \Phi_3(\mathbf{u}) \end{pmatrix}=\begin{pmatrix} x^2+y^2+z^2\\ 2x^2+y^2-4z\\ 3x^2+-4y+z^2 \end{pmatrix}=\begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix} \end{equation}Primeiramente, para implementarmos o método de Newton-Raphson, calcularemos a matriz Jacobiana de $\Phi$, que é denominada por $J$. Temos:
\begin{equation} J(\mathbf{u})=\frac{\partial \Phi(\mathbf{u})}{\partial \mathbf{u}}=\begin{pmatrix} \frac{\partial \Phi_1}{\partial x} & \frac{\partial \Phi_1}{\partial y} & \frac{\partial \Phi_1}{\partial z}\\ \frac{\partial \Phi_2}{\partial x} & \frac{\partial \Phi_2}{\partial y} & \frac{\partial \Phi_2}{\partial z}\\ \frac{\partial \Phi_3}{\partial x} & \frac{\partial \Phi_3}{\partial y} &\frac{\partial \Phi_3}{\partial z} \end{pmatrix}=\begin{pmatrix} 2x & 2y &2z \\ 4x & 2y& -4\\ 6x&-4 & 2z \end{pmatrix} \end{equation}Para uma sequência de vetores $\big (\mathbf{u}_n \big )_{n=0:N}$, definiremos a seguinte equação de recursão dada pelo algoritmo de Newton-Raphson:
\begin{equation} \mathbf{u}_{n+1}=\mathbf{u}_{n}-J^{-1}(\mathbf{u}_{n})\Phi(\mathbf{u}_{n}) \end{equation}Vamos definir algumas funções abaixo para a implementação do método:
import numpy as np
x,y,z=.5,.5,.5
#Pega os valores de x,y,z e deixa em forma de vetor
def vec(x,y,z):
return np.array([x,y,z])
#Monta o sistema a ser resolvido
def Phi(u):
x,y,z=u[0],u[1],u[2]
return np.array([x**2+y**2+z**2, 2*x**2+y**2-4*z, 3*x**2-4*y+z**2])
#Matriz Jacobiana do Sistema (forma analítica)
def J(u):
x,y,z=u[0],u[1],u[2]
return np.array([[2*x,2*y,2*z],[4*x,2*y,-4],[6*x,-4,2*z]])
#Matriz Jacobiana do Sistema (forma numérica)
def dif2(Phi, u, e=0.001):
out=[]
for i in range(np.shape(u)[0]):
delta=[0,0,0]
delta[i]=e
delta=np.array(delta)
out.append(list((Phi(u+delta)-Phi(u-delta))/(2*e)))
return np.array(out).T
#Funções para inverter matrizes
def inv(A):
return np.linalg.inv(A)
#Algoritmo de Newton-Raphson (é possível escolher tanto a
#forma analítica ou numérica para o cálculo da matriz Jacobiana)
def n_r(u,Phi,J,it=3,numerical=False):
v=u
for i in range(it):
if numerical:
v=v-inv(J(v))@Phi(v)
else:
v=v-inv(dif2(Phi,v))@Phi(v)
return v
Se $\mathbf{u}_{0}=(0.5,0.5,0.5)^T$, então:
\begin{equation} \mathbf{u}_{1}=\mathbf{u}_{0}-J^{-1}(\mathbf{u}_{0})\Phi(\mathbf{u}_{0})= \begin{pmatrix} 0.5\\ 0.5\\ 0.5 \end{pmatrix}-\begin{pmatrix} 1 & 1& 1\\ 2& 1& -4\\ 3 & -4& 1 \end{pmatrix}^{-1}\begin{pmatrix} 0.75\\ -1.25\\ -1 \end{pmatrix}=\begin{pmatrix} 0.5\\ 0.15\\ 0.1 \end{pmatrix} \end{equation}De fato:
u=np.array([.5,.5,.5])
print('u1 é dado por:', n_r(u,Phi,J,it=1,numerical=False))
Agora aproximando a raíz de $\Phi$ com dez iterações:
u=np.array([.5,.5,.5])
solucao=n_r(u,Phi,J,it=10,numerical=False)
print('u3 é dado por:', solucao)
Para checar se a solução está correta, basta testarmos em nossa função. Se $\Phi(.)$ zerar, temos uma solução:
Phi(solucao)
Chegamos a uma solução válida.