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

Cálculo Numérico - Versão Scilab

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 no Scilab. Para tal definimos as funções que implementarão F(x) e a JF (x)

function y=F(x)  
    y(1)=x(1)^2/3+x(2)^2-1  
    y(2)=x(1)^2+x(2)^2/4-1  
endfunction  
 
function y=JF(x)  
    y(1,1)=2*x(1)/3  
    y(1,2)=2*x(2)  
    y(2,1)=2*x(1)  
    y(2,2)=x(2)/2  
endfunction

Alternativamente, estas funções poderiam ser escritas como

function y=F(x)  
    y=[x(1)^2/3+x(2)^2-1; x(1)^2+x(2)^2/4-1]  
endfunction  
 
function y=JF(x)  
    y=[2*x(1)/3  2*x(2); 2*x(1) x(2)/2]  
endfunction

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

delta=-JF(x)\F(x)  
x=x+delta

Ou simplesmente

x=x-JF(x)\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).

No Scilab, podemos implementá-las com o seguinte código:

function y=F(x)  
    y(1) = x(1)^2-cos(x(1)*x(2))-1  
    y(2) = sin(x(2))-2*cos(x(1))  
endfunction  
 
function y=JF(x)  
    y(1,1) = 2*x(1)+x(2)*sin(x(1)*x(2))  
    y(1,2) = x(1)*sin(x(1)*x(2))  
 
    y(2,1) = 2*sin(x(1))  
    y(2,2) = cos(x(2))  
endfunction

E agora, basta iterar:

x=[1.5; .5]  
x=x-JF(x)\F(x) //(5 vezes)

5.1.1 Código Scilab: Newton para sistemas


function [x] = newton(F,JF,x0,TOL,N)  
  x = x0  
  k = 1  
  //iteracoes  
  while (k <= N)  
    //iteracao de Newton  
    delta = -inv(JF(x))*F(x)  
    x = x + delta  
    //criterio de parada  
    if (norm(delta,inf)<TOL) then  
      return x  
    end  
    k = k+1  
  end  
  error(Num. de iter. max. atingido!)  
endfunction

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 Scilab.
  • ) 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

Exemplo de implementação:

function z=f(x,y)  
    z=x^2/8+(y-1)^2/5-1  
endfunction  
function z=g(x,y)  
    z=atan(x)+x-y-y^3  
endfunction  
 
contour([-3:.1:3],[-2:.1:4],f,[0 0])  
contour([-3:.1:3],[-2:.1:4],g,[0 0])  
 
function y=F(x)  
    y(1)=f(x(1),x(2))  
    y(2)=g(x(1),x(2))  
endfunction  
function y=JF(x)  
    y(1,1)=x(1)/4  
    y(1,2)=2*(x(2)-1)/5  
    y(2,1)=1/(1+x(1)^2)+1  
    y(2,2)=-1-3*x(2)^2  
endfunction  
 
//primeiro ponto  
//x=[-1.2;-1.0]  
 
//segundo ponto  
//x=[2.8;1.4]  
 
x=x-JF(x)\F(x)   // 4 vezes

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.

Uma possível implementação é

 
function z=f(x,y)  
    z=(x-3)^2/16+(y-1)^2/36-1  
endfunction  
 
function z=g(x,y)  
    z=atan(x)+x-sin(y)-0.01*y^3  
endfunction  
 
contour([-1:.1:7],[-5:.1:7],f,[0 0])  
contour([-1:.1:7],[-5:.1:7],g,[0 0])  
 
function y=F(x)  
    y(1)=f(x(1),x(2))  
    y(2)=g(x(1),x(2))  
endfunction  
function y=JF(x)  
    y(1,1)=(x(1)-3)/8  
    y(1,2)=(x(2)-1)/18  
    y(2,1)=1/(1+x(1)^2)+1  
    y(2,2)=-cos(x(2))-0.03*x(2)^2  
endfunction  
 \end{resp}  
 
//primeiro ponto  
//x=[-.5;-2.0]  
 
//segundo ponto  
//x=[3;7]  
 
x=x-JF(x)\F(x)   // 4 vezes

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

Um exemplo de implementação no Scilab é:

deff(’z=f(x,y)’,’z=-x^4-y^6+3*x*y^3-x’)  
contour([-2:.01:0],[-2:.01:0],f,[ 0:.2: 3])  
deff(’z=F(x)’,’z=[-4*x(1)^3+3*x(2)^3-1;-6*x(2)^5+9*x(1)*x(2)^2]’)  
deff(’z=JF(x)’,’z=[-12*x(1)^2,9*x(2)^2;9*x(2)^2,-30*x(2)^4+18*x(1)*x(2)]’)  
x=[-1.2;-1.2]  
x=x-JF(x)\F(x)  
x=x-JF(x)\F(x)  
x=x-JF(x)\F(x)  
x=x-JF(x)\F(x)  
mprintf(’f(%f,%f)=%f’,x(1),x(2),f(x(1),x(2)))

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)

Exemplo de implementação no Scilab:

function y=F(x)  
    y(1)=x(1)-x(2)  
    for j=2:10  
        y(j)=-x(j-1)+5*(x(j)+x(j)^3)-x(j+1)-10*exp(-j/3)  
    end  
    y(11)=x(11)-1  
endfunction  
 
function y=JF(x)  
    y=zeros(11,11)  
 
    y(1,1)=1  
    y(1,2)=-1  
    for j=2:10  
 
        y(j,j-1)=-1  
        y(j,j)=5*(1+3*x(j)^2)  
        y(j,j+1)=-1  
    end  
    y(11,11)=1  
endfunction

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 30/7/2018 às 13:16:34.

Informe erros ou edite você mesmo!