Use o botão abaixo para reportar erros ou dar sugestões.

Cálculo Numérico - Versão Python

5.1 Método de Newton para sistemas


Nesta seção, construiremos o método de Newton ou Newton-Raphson generalizado para sistemas. Assumimos, portanto, que a função F(x) é diferenciável e que existe um ponto x* tal que F(x*) = 0. Seja x(k) uma aproximação para x*, queremos construir uma nova aproximação x(k+1) através da linearização de F(x) no ponto x(k).

  • Linearização da função F(x) no ponto x(k):
    F(x) = F(x(k)) + J F x(k) x - x(k) + O x - x(k)2 (5.17)
  • A aproximação x(k) é definida como o ponto x em que a linearização F(x(k)) + J F x(k) x - x(k) é nula, ou seja:
    F(x(k)) + J F x(k) x(k+1) - x(k) = 0 (5.18)

Supondo que a matriz jacobina seja inversível no ponto x(k), temos: JF x(k) x(k+1) - x(k) = -F(x(k)) (5.19) x(k+1) - x(k) = -J F -1 x(k) F(x(k)) (5.20) x(k+1) = x(k) - J F -1 x(k) F(x(k)) (5.21)

Desta forma, o método iterativo de Newton-Raphson para encontrar as raízes de F(x) = 0 é dado por:

x(k+1) = x(k) - J F -1 x(k) F(x(k)),k > 0 x(0) = dado inicial (5.22)

Observação 5.1.1. Usamos subíndices para indicar o elemento de um vetor e superíndices para indicar o passo da iteração. Assim, x(k) se refere à iteração k e xi(k) se refere à componente i no vetor x(k).

Observação 5.1.2. A notação JF -1 x(k) enfatiza que a jacobiana deve ser calculada a cada passo.

Observação 5.1.3. Podemos definir o passo Δ(k) como

Δ(k) = x(k+1) - x(k) (5.23)

Assim, Δ(k) = -J F -1 x(k) F(x(k)), ou seja, Δ(k) resolve o problema linear:

JF x(k) Δ(k) = -F(x(k)) (5.24)

Em geral, é menos custoso resolver o sistema acima do que calcular o inverso da jacobiana e multiplicar pelo vetor F(x(k)).

Exemplo 5.1.1. Retornamos ao nosso exemplo inicial, isto é, resolver numericamente o seguinte sistema não linear: x12 3 + x22 = 1 (5.25) x12 + x22 4 = 1 (5.26)

Para tal, definimos a função F(x):

F(x) = x12 3 + x22 - 1 x12 + x22 4 - 1 (5.27)

cuja jacobiana é:

JF = 2x1 3 2x2 2x1 x2 2 (5.28)

Faremos a implementação numérica em Python. Para tal definimos as funções que implementarão F(x) e a JF (x)

>>> def F(x):  
...     y = np.zeros(2)  
...     y[0] = x[0]**2/3 + x[1]**2 - 1  
...     y[1] = x[0]**2 + x[1]**2/4 - 1  
...     return y  
...  
>>> def JF(x):  
...     y = np.zeros((2,2))  
...     y[0,0] = 2*x[0]/3  
...     y[0,1] = 2*x[1]  
...     y[1,0] = 2*x[0]  
...     y[1,1] = x[1]/2  
...     return y  
...

Desta forma, se x é uma aproximação para a raiz, pode-se calcular a próxima aproximação através dos comandos:

>>> delta = -np.linalg.inv(JF(x)).dot(F(x))  
>>> x = x + delta

Ou simplesmente

>>> x = x - np.linalg.inv(JF(x)).dot(F(x))

Observe que as soluções exatas desse sistema são ± 9 11, ± 8 11.

Exemplo 5.1.2. Encontre uma aproximação para a solução do sistema x12 = cos(x 1x2) + 1 (5.29) sen (x2) = 2 cos(x1) (5.30)

que fica próxima ao ponto x1 = 1,5 e x2 = 0,5.

Solução. Vamos, aqui, dar as principais ideias para obter a solução usando o método de Newton. Começamos definindo nossa aproximação inicial por x(1) = (1,5, 0,5). Então iteramos:

x(n+1) = x(n) - J F -1(x)F(x),n 1. (5.31)

onde

F(x) = x12 - cos(x 1x2) - 1 sen (x2) - 2 cos(x1) (5.32)

e sua jacobiana é

JF (x) = 2x1 + x2 sen (x1x2) x1 sen (x1x2) 2 sen (x1) cos(x2) (5.33)

As iterações convergem para x = (1,3468109,0,4603195).

Em Python, podemos implementá-las com o seguinte código:

def F(x):  
    y = np.zeros(2)  
 
    y[0] = x[0]**2 - np.cos(x[0]*x[1]) - 1  
    y[1] = np.sin(x[1]) - 2*np.cos(x[0])  
 
    return y  
 
def JF(x):  
    y = np.zeros((2,2))  
 
    y[0,0] = 2*x[0] + x[1]*np.sin(x[0]*x[1])  
    y[0,1] = x[0]*np.sin(x[0]*x[1])  
 
    y[1,0] =  2*np.sin(x[0])  
    y[1,1] = np.cos(x[1])  
 
    return y

E agora, basta iterar:

>>> x = np.array([1.5,0.5])  
>>> x=x-np.linalg.inv(JF(x)).dot(F(x))

5.1.1 Código Python: Newton para Sistemas


from __future__ import division  
import numpy as np  
from numpy import linalg  
 
def newton(F,JF,x0,TOL,N):  
    #preliminares  
    x = np.copy(x0).astype(double)  
    k=0  
    #iteracoes  
    while (k < N):  
       k += 1  
       #iteracao Newton  
       delta = -np.linalg.inv(JF(x)).dot(F(x))  
       x = x + delta  
       #criterio de parada  
       if (np.linalg.norm(delta,np.inf) < TOL):  
           return x  
 
    raise NameError(num. max. iter. excedido.)

Exercícios resolvidos


Esta seção carece de exercícios resolvidos. Clique em e inicie a editá-la agora mesmo. Veja outras formas de participar clicando aqui.

Exercícios


E 5.1.1. Faça o que se pede:

  • Encontre o gradiente da função
    f(x,y) = x2y + cos(xy) - 4 (5.34)
  • Encontre a matriz jacobiana associada à função
    F(x,y) = x cos(x) + y e-2x+y . (5.35)
  • Encontre a matriz jacobiana associada à função
    L(x) = a11x1 + a12x2 + a13x3 - y1 a21x1 + a22x2 + a23x3 - y2 a31x1 + a32x2 + a33x3 - y3 . (5.36)

Resposta. f = [2xy - y sen (xy),x2 - x sen (xy)]T
JF = cos(x) - x sen (x) 1 - 2e-2x+y e-2x+y (5.37)
JL ij = aij (5.38)

E 5.1.2. Encontre uma aproximação numérica para o seguinte problema não linear de três equações e três incógnitas: 2x1 - x2 = cos(x1) (5.39) -x1 + 2x2 - x3 = cos(x2) (5.40) -x2 + x3 = cos(x3) (5.41)

Partindo das seguintes aproximações iniciais:

  • x(0) = [1,1,1]T
  • x(0) = [-0,5, - 2, - 3]T
  • x(0) = [-2, - 3, - 4]T
  • x(0) = [0,0,0]T

Resposta.

Este exercício está sem resposta sugerida. Clique em e inicie a editá-la agora mesmo. Veja outras formas de participar clicando aqui.

E 5.1.3. Encontre os pontos de intersecção entre a parábola y = x2 + 1 e a elipse x2 + y24 = 1 seguindo os seguintes passos:

  • Faça um esboço das duas curvas e entenda o problema. Verifique que existem dois pontos de intersecção, um no primeiro quadrante e outro no segundo quadrante do plano xy.
  • A partir de seu esboço, encontre aproximações para x e y em cada ponto.
  • Escreva o problema na forma F x y = 0 0
  • Encontre a jacobiana JF .
  • Construa a iteração do método de Newton.
  • Implemente no computador.
  • Resolva o sistema analiticamente e compare as respostas.

Resposta. As curvas possuem dois pontos de intersecção. A posição exata destes pontos de intersecção é dada por 23 - 3,23 - 2 e -23 - 3,23 - 2. Use a solução exata para comparar com a solução aproximada obtida.

E 5.1.4. Encontre os pontos de intersecção entre a parábola y = x2 e a curva y = cos(x) seguindo os seguintes passos:

  • ) Faça um esboço das duas curvas, entenda o problema. Verifique que existem dois pontos de intersecção, um no primeiro quadrante e outro no segundo quadrante do plano xy.
  • ) A partir de seu esboço, encontre aproximações para x e y em cada ponto.
  • ) Escreva o problema na forma F x y = 0 0
  • ) Encontre a jacobiana JF .
  • ) Construa a iteração do método de Newton.
  • ) Implemente no Python.
  • ) Transforme o sistema em um problema de uma única variável e compare com a resposta do Problema 3.4.1.

Resposta. ±0.8241323, 0.6791941

E 5.1.5. Encontre uma aproximação com erro inferior a 10-5 em cada incógnita para a solução próxima da origem do sistema 6x - 2y + ez = 2 (5.42) sen (x) - y + z = 0 (5.43) sen (x) + 2y + 3z = 1 (5.44)

Resposta. x 0,259751,y 0,302736,z 0,045896

E 5.1.6. (Entenda casos particulares)

  • Considere a função L(x) = Ax - b, onde A é uma matriz n × n inversível e b um vetor coluna em n. O que acontece quando aplicamos o método de Newton para encontrar as raízes de L(x)?
  • Mostre que o método de Newton-Raphson aplicado a uma função diferenciável do tipo f : se reduz ao método de Newton estudado na primeira área.

Resposta.

Este exercício está sem resposta sugerida. Clique em e inicie a editá-la agora mesmo. Veja outras formas de participar clicando aqui.

E 5.1.7. Considere a função f(x) = sen (x) x+1 , encontre a equação da reta que tangencia dois pontos da curva y = f(x) próximos ao primeiro e segundo ponto de máximo no primeiro quadrante, respectivamente. Veja a Figura 5.1.


PIC

Figura 5.1: Reta bitangente a uma curva.


Resposta. y = mx + b com m -0.0459710 e b 0.479237 Uma metodologia possível para resolver este problema é dada a seguir:

Sejam x1 e x2 as abscissas dos dois pontos em que a reta tangencia a curva. A equação da reta bitangente assume a seguinte forma:

y = f(x1) + m(x - x1) (5.45)

onde o coeficiente angular m é dado por

m = f(x2) - f(x1) x2 - x1 (5.46)

Da condição de tangência, temos que o coeficiente angular da reta, m, deve igual à derivada da função f(x) nos dois pontos de tangência.

m = f(x 1) = f(x 2) (5.47)

E sabemos que:

f(x) = cos(x) 1 + x - sen (x) (1 + x)2. (5.48)

Assim, podemos reescrever o problema como cos(x1) 1 + x1 - sen (x1) (1 + x1)2 - cos(x2) 1 + x2 + sen (x2) (1 + x2)2 = 0 (5.49) cos(x1) 1 + x1 - sen (x1) (1 + x1)2 - f(x2) - f(x1) x2 - x1 = 0 (5.50)

Este é um sistema não linear de duas incógnitas.

Os valores iniciais para o método podem ser obtidos do gráfico buscando valores próximos aos dois primeiros pontos de máximos. Por exemplo: x1(0) = 1 e x2(0) = 8. Obtemos x1 1,2464783 e x2 8,1782997 e m pode ser obtido através desses valores.


PIC

Figura 5.2: Sistema mecânico com dois segmentos.


E 5.1.8. (Estática) Considere o sistema mecânico constituído de dois segmentos de mesmo comprimento L presos entre si e a uma parede por articulações conforme a Figura 5.2.

O momento em cada articulação é proporcional à deflexão com constante de proporcionalidade k. Os segmentos são feitos de material homogêneo de peso P. A condição de equilíbrio pode ser expressa em termos dos ângulos θ1 e θ2 conforme: kθ1 = 3PL 2 cos θ1 + k θ2 - θ1 (5.51) k θ2 - θ1 = PL 2 cos θ2 (5.52)

Considere P = 100N, L = 1m e calcule os ângulos θ1 e θ2 quando:

  • k = 1000 Nm/rad
  • k = 500 Nm/rad
  • k = 100 Nm/rad
  • k = 10 Nm/rad

Obs:Você deve escolher valores para iniciar o método. Como você interpretaria fisicamente a solução para produzir palpites iniciais satisfatórios? O que se altera entre o caso a e o caso d?

Resposta. 0.1956550; 0.2441719, 0.3694093; 0.4590564, 0.9990712; 1.1865168 e 1.4773606; 1.5552232

E 5.1.9. (estática - problemas de três variáveis) Considere, agora, o sistema mecânico semelhante ao do Problema 5.1.8, porém constituído de três segmentos de mesmo comprimento L presos entre si e a uma parede por articulações.

O momento em cada articulação é proporcional à deflexão com constante de proporcionalidade k. Os segmentos são feitos de material homogêneo de peso P. A condição de equilíbrio pode ser expressa em termos dos ângulos θ1, θ2 e θ3 conforme: kθ1 = 5PL 2 cos θ1 + k θ2 - θ1 (5.53) k θ2 - θ1 = 3PL 2 cos θ2 + k θ3 - θ2 (5.54) k θ3 - θ2 = PL 2 cos θ3 (5.55)

Considere P = 10N, L = 1m e calcule os ângulos θ1, θ2 e θ3 quando:

  • k = 1000Nm/rad
  • k = 100Nm/rad
  • k = 10Nm/rad

Resposta. 0.0449310; 0.0648872; 0.0698750, 0.3981385; 0.5658310; 0.6069019,
1.1862966; 1.4348545; 1.480127


PIC

Figura 5.3: intersecção entre duas curvas.


E 5.1.10. Considere o problema de encontrar os pontos de intersecção das curvas descritas por (ver Figura 5.3): x2 8 + (y - 1)2 5 = 1 (5.56) (5.57) tan -1(x) + x = y + y3 (5.58)

Com base no gráfico, encontre soluções aproximadas para o problema e use-as para iniciar o método de Newton-Raphson. Encontre as raízes com erro inferior a 10-5.

Resposta. -1,2085435,-1,0216674 e 2,7871115, 1,3807962

E 5.1.11. Considere o sistema de equações dado por (x - 3)2 16 + (y - 1)2 36 = 1 (5.59) tanh(x) + x = 2 sen y - 0.01y3 (5.60)

Usando procedimentos analíticos, determine uma região limitada do plano onde se encontram necessariamente todas as raízes do problema. Encontre as raízes desse sistema com pelo menos quatro dígitos significativos corretos usando o método de Newton. Você deve construir o método de Newton indicando as funções envolvidas e calculando a matriz jacobiana analiticamente. Use que d du tanh u = 1 - tanh 2u, se precisar.

Resposta. A primeira curva trata-se de uma elipse de centro (3,1) e semi-eixos 4 e 6, portanto seus pontos estão contidos no retângulo - 1 x 7 e - 5 y 7.

As soluções são -0,5384844,-1,7978634 e 2,8441544, 6,9954443.

E 5.1.12. (Otimização) Uma indústria consome energia elétrica de três usinas fornecedoras. O custo de fornecimento em reais por hora como função da potência consumida em kW é dada pelas seguintes funções C1(x) = 10 + .3x + 10-4x2 + 3.4 10-9x4 (5.61) C2(x) = 50 + .25x + 2 10-4x2 + 4.3 10-7x3 (5.62) C3(x) = 500 + .19x + 5 10-4x2 + 1.1 10-7x4 (5.63)

Calcule a distribuição de consumo que produz custo mínimo quando a potência total consumida é 1500kW. Dica: Denote por x1, x2 e x3 as potências consumidas das usinas 1, 2 e 3, respectivamente. O custo total será dado por C(x1,x2,x3) = C1(x1) + C2(x2) + C3(x3) enquanto o consumo total é x1 + x2 + x3 = 1500. Isto é, queremos minimizar a função custo total dada por:

C(x1,x2,x3) = C1(x1) + C2(x2) + C3(x3) (5.64)

restrita à condição

G(x1,x2,x3) = x1 + x2 + x3 - 1500 = 0. (5.65)

Pelos multiplicadores de Lagrange, temos que resolver o sistema dado por: C(x1,x2,x3) = λG(x1,x2,x3) (5.66) G(x1,x2,x3) = 0 (5.67)

Resposta. (x1,x2,x3) (453,62,901,94,144,43)

E 5.1.13. Encontre a função do tipo f(x) = Abx que melhor aproxima os pontos (0,3,1), (1,4,4) e (2,6,7) pelo critério dos mínimos quadrados. Dica: Você deve encontrar os valores de A e b que minimizam o resíduo dado por

R = 3,1 - f(0) 2 + 4,4 - f(1) 2 + 6,7 - f(2) 2. (5.68)

Dica: Para construir aproximações para resposta e iniciar o método, considere a função f(x) = Abx que passa pelo primeiro e terceiro ponto.

Resposta. Inicialização do método: A(0) = 3,1 e b(0) = 6,7 3,1 A 3.0297384 e b 1.4835346.

E 5.1.14. Encontre o valor máximo da função

f(x,y) = -x4 - y6 + 3xy3 - x (5.69)

na região (x,y) [-2,0] × [-2,0] seguindo os seguintes passos:

  • Defina a função z = f(x,y) = -x4 - y6 + 3xy3 - x e trace o gráfico de contorno na região.
  • Com base no gráfico, encontre valores aproximados para as coordenadas xy do ponto de máximo.
  • Sabendo que o ponto de máximo acontece quando o gradiente é nulo, escreva o problema como um sistema de duas equações não lineares e duas incógnitas.
  • Implemente o método de Newton.

Resposta. f(-1,1579702,-1,2020694) 2.376985

E 5.1.15. A função f(x,y,z) = sen (x) + sen (2y) + sen (3z) possui um máximo quando x = π2, y = π4 e z = π6. Calcule numericamente este ponto.

Resposta.

Este exercício está sem resposta sugerida. Clique em e inicie a editá-la agora mesmo. Veja outras formas de participar clicando aqui.

E 5.1.16. Encontre as raizes do problema 3x - cos(yz + z) - 12 = 0 (5.70) 4x2 - 25y2 + 0.4y + 2 = 0 (5.71) e-xy + 2x - 5z = 10 (5.72)

no cubo |x| < 2,|y| < 2,|z| < 2. Dica: Reduza a um problema de duas incógnitas e use recursos gráficos para aproximar as raízes na região.

Resposta. x 0,2982646,y -0,2990796,z -1,6620333 e x -0,0691328,y 0,2923039,z -0,8235705.

E 5.1.17. Considere o seguinte sistema de equações não lineares: x1 - x2 = 0 -xj-1 + 5(xj + xj3) - x j+1 = 10 exp(-j3),2 j 10 x11 = 1 (5.73)

  • Escreva este sistema na forma F(x) = 0 onde x = x1 x2 x 11 e calcule analiticamente a matriz jacobiana (F1,,F11) (x1,x11) . Dica: Use a regularidade nas expressões para abreviar a notação.
  • Construa a iteração para encontrar a única solução deste problema pelo método de Newton e, usando esse método, encontre uma solução aproximada com erro absoluto inferior a 10-4.

Resposta.
F x = x1 - x2
 
- x1 + 5(x2 + x23) - x 3 - 10 exp(-23)
 
- x2 + 5(x3 + x33) - x 4 - 10 exp(-33)
 
- x3 + 5(x4 + x43) - x 5 - 10 exp(-43)
 
- x9 + 5(x10 + x103) - x 11 - 10 exp(-103)
 
x 11 - 1
(5.74)
JF (x) = 1 - 1 0 0 0 0
       
- 1 5(1 + 3x22) - 1 0 0 0
       
0 - 1 5(1 + 3x32) - 1 0 0
       
0 0 - 1 5(1 + 3x42) - 1 0
       
       
0 0 0 0 0 1
(5.75)

Resposta final: 0,80447, 0,80447, 0,68686, 0,57124, 0,46535, 0,37061, 0,28883, 0,22433, 0,19443, 0,28667, 1

E 5.1.18. Considere a função

f(x,y) = e-(x-1)2-(y-2)2 1 + x2 + y2 (5.76)
  • Encontre o valor máximo desta função.
  • Usando multiplicadores de Lagrange, encontre o valor máximo desta função restrito à condição
    (x - 1)2 + (y - 2)2 = 1. (5.77)
  • Parametrize a circunferência para transformar o problema de máximo com restrição em um problema de uma única variável. Resolva usando as técnicas de equações lineares de uma variável.

Resposta. f(0,8108792, 1,6217584) 0,1950369 e f(0,5527864, 1,1055728) 0,1455298

Creative Commons License Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 15/5/2019 às 15:24:50.

Informe erros ou edite você mesmo!