Matlab - Metodo de Newton Sistemas no lineales

 
Vista:
sin imagen de perfil
Val: 43
Ha disminuido su posición en 13 puestos en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por wil (6 intervenciones) el 02/04/2020 17:35:00
Hola, copie el siguiente codigo en matlab de un articulo y no compila, donde esta el error ????

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% Newton's Method
format long;
n=5; % set some number of iterations, may need adjusting
f = @(x) [ 3*x(1)-cos(x(2)*x(3))-0.5
x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06
exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3)]; % the vector
function,3x1
% the matrix of partial derivatives
df = @(x) [3 x(3)*sin(x(2)*x(3)) x(2)*sin(x(2)*x(3))
2*x(1) -162*x(2)-16.2 cos(x(3))
-x(2)*exp(x(1)*x(2)) -x(1)*exp(x(1)*x(2)) 20];% 3x3
x = [0.1;0.1;-0.1]; % starting guess
for i = 1:n
y = -df(x)\f(x) % solve for increment, similar A\b
x = x + y % add on to get new guess
f(x) % see if f(x) is really zero
end


1
2
3
4
5
6
7
8
9
function y = F(x)
x1 = x(1);
x2 = x(2);
x3 = x(3);
y = zeros(3,1);
y(1) = 3*x(1)-cos(x(2)*x(3))-0.5; % f1(x1,x2)
y(2) = x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06;
y(3) = exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3);
end


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function x = NewtonMethod(funcF, JacobianF, n)
F = funcF;
J = JacobianF;
x = [0.1 0.1 -0.1]';
Iter = 1;
MaxIter = 100;
TOL = 1e-5;
while Iter < MaxIter
disp(['Iter = ' num2str(Iter)]);
y = J(x)\(-F(x));
x = x+y;
if norm(y,2)<TOL
break;
end
disp(['x = ' num2str(x')]);
end
if Iter >= MaxIter
disp('Maximum number of iteration exceeded!');
end
end


1
2
3
function newtonSol
x = NewtonMethod(@F,@J,2);
end
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
Imágen de perfil de JOSE JEREMIAS CABALLERO
Val: 6.975
Oro
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por JOSE JEREMIAS CABALLERO (5917 intervenciones) el 03/04/2020 01:47:45
guardar el archivo con el nombre de newtonSol.m es un solo archivo todo unido para que funcione mejor. Hice unos arreglos para mejorar la salida de la solución


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function newtonSol
x = NewtonMethod(@F,@J);
end
 
function y=J(x)
% Newton's Method
format long;
n=5; % set some number of iterations, may need adjusting
f = @(x)[ 3*x(1)-cos(x(2)*x(3))-0.5
           x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06
           exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3)    ]; % the vector function,3x1
% the matrix of partial derivatives
df = @(x)[3                      x(3)*sin(x(2)*x(3))     x(2)*sin(x(2)*x(3))
         2*x(1)                 -162*x(2)-16.2          cos(x(3))
        -x(2)*exp(x(1)*x(2))   -x(1)*exp(x(1)*x(2))    20];% 3x3
x = [0.1;0.1;-0.1]; % starting guess
y=df(x);
% for i = 1:n
% y = -df(x)\f(x) ;% solve for increment, similar A\b
% x = x + y ;% add on to get new guess
% f(x) % see if f(x) is really zero
% end
end
function y = F(x)
if nargin==0
    x=rand(1,3);
end
x1 = x(1);
x2 = x(2);
x3=x(3);
y = zeros(3,1);
y(1) = 3*x(1)-cos(x(2)*x(3))-0.5; % f1(x1,x2)
y(2) = x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06;
y(3) = exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3);
end
function x = NewtonMethod(funcF, JacobianF)
format long
F = funcF;
J = JacobianF;
x = [0.1 0.1 -0.1]';
Iter = 1;
MaxIter = 100;
TOL = 1e-5;
disp([ '      Iter                        x1              x2                     x3             f(x1)                 f(x2)                 f(x3)' ]);
while Iter < MaxIter
y = J(x)\(-F(x));
x = x+y;
M(Iter,:)=[Iter x' F(x)'];
if norm(y,2)<TOL
break;
end
Iter=Iter+1;
end
disp(M)
if Iter >= MaxIter
disp('Maximum number of iteration exceeded!');
end
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
>> newtonSol
      Iter                        x1              x2                     x3             f(x1)                 f(x2)                 f(x3)
   1.000000000000000   0.499869683245681   0.019467829453816  -0.521488532745465  -0.000339416617403  -0.344379208155604   0.032520675950323
   2.000000000000000   0.499985822222047   0.008787984589569  -0.523167915907467  -0.000031964428295  -0.148261868616623   0.004232965030093
   3.000000000000000   0.499997926468827   0.004204867055769  -0.523402648922643  -0.000003798749322  -0.069383212995024   0.001822317251152
   4.000000000000000   0.499999873586911   0.002060291408031  -0.523504585596723   0.000000202419457  -0.033639102926783   0.000854185005622
   5.000000000000000   0.500000136693704   0.001020576301393  -0.523552544340334   0.000000552832881  -0.016577529065628   0.000414467044012
   6.000000000000000   0.500000115432392   0.000508208156116  -0.523575855387564   0.000000381697929  -0.008233927398511   0.000204332359703
   7.000000000000000   0.500000069193400   0.000253721080224  -0.523587357462635   0.000000216404116  -0.004105538203387   0.000101510202068
   8.000000000000000   0.500000037449153   0.000126831188564  -0.523593073958836   0.000000114552466  -0.002050993011022   0.000050619200946
   9.000000000000000   0.500000019444241   0.000063441264060  -0.523595925144832   0.000000058884427  -0.001025586474692   0.000025288939167
  10.000000000000000   0.500000009906975   0.000031743552896  -0.523597349721348   0.000000029859052  -0.000513082423669   0.000012645888205
  11.000000000000000   0.500000005002334   0.000015885767386  -0.523598062126323   0.000000015041594  -0.000256746985301   0.000006326587295
  12.000000000000000   0.500000002514705   0.000007950516795  -0.523598418543257   0.000000007552780  -0.000128491758666   0.000003165850329

Saludos
JOSE JEREMIAS CABALLERO
Asesor de Proyectos con Matlab
Servicios de programación matlab
Servicio de Asesoría Online en Matlab
[email protected]


http://matlabcaballero.blogspot.com
https://www.facebook.com/matlabcaballero
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
sin imagen de perfil
Val: 43
Ha disminuido su posición en 13 puestos en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por wil (6 intervenciones) el 03/04/2020 03:12:50
Muy bien gracias, ahora modifique el problema para un sistema 4x4 colocando nuevas ecuaciones , pero me da error, donde estoy errando?



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
function NewtonSol
x = NewtonMethod(@F,@J);
end
 
function y=J(x)
% Newton's Method
format long;
n=30; % set some number of iterations, may need adjusting
h=1/5;
f=@(x)[(-1/h)+(2*x(1)/h)-(x(2)/h)+(1/6)+(x(1)/6)-(x(1)*x(2)/6)-((x(2)^2)/6)+(h/6)+(2*x(1)*h/3)+(x(2)*h/6)+0.3023
     -(x(1)/h)+(2*x(2)/h)-(x(3)/h)+((x(1)^2)/6)-(x(2)*x(3)/6)-((x(3)^2)/6)+(x(1)*h/6)+(2*x(2)*h/3)+(x(3)*h/6)+(x(1)*x(2)/6)+0.451
    -(x(2)/h)+(2*x(3)/h)+((x(2)^2)/6)+(x(2)*x(3)/6)-((x(3)*x(4))/6)-((x(4)^2)/6)+(x(2)*h/6)+(2*x(3)*h/3)+(x(4)*h/6)-(x(4)/h)+0.6729
     -(x(3)/h)+(2*x(4)/h)-(exp(1)/h)+((x(3)^2)/6)+(x(3)*x(4)/6)-(x(4)*exp(1)/6)-((exp(1)*exp(1))/6)+(x(3)*h/6)+((2*x(4)*h)/3)+(h*exp(1)/6)+1.004];
%f = @(x)[ 3*x(1)-cos(x(2)*x(3))-0.5
 %          x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06
  %         exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3)    ]; % the vector function,3x1
  % the matrix of partial derivatives
  df=@(x)[  (2/h)+1/6-(x(2)/6)+(2*h/3)                 -1/h -(x(1)/6)-(x(2)/3)+h/6                              0                                     0
        (-1/h)+(x(1)/3)+h/6+(x(2)/6)                  2/h-(x(3)/6)+(2*h/3) +(x(1)/6)             -1/h-(x(2)/h)-(x(3)/3) +h/6                            0
                      0                                -1/h +(x(2)/3)+(x(3)/6)  +h/6              2/h+(x(2)/6)-(x(4)/6)+ (2*h/3)             (-x(3)/6)-(x(4)/3)+(h/6)-(1/h)
                      0                                         0                                    (-1/h)+(x(3)/3)+(x(4)/6)+(h/6)              (2/h) +(x(3)/6)-(exp(1)/6)+(2*h/3) ];
%df = @(x)[3                      x(3)*sin(x(2)*x(3))     x(2)*sin(x(2)*x(3))
 %        2*x(1)                 -162*x(2)-16.2          cos(x(3))
  %      -x(2)*exp(x(1)*x(2))   -x(1)*exp(x(1)*x(2))    20];% 3x3
 x=[0;0;0;0];
%x = [0.1;0.1;-0.1]; % starting guess
y=df(x);
% for i = 1:n
% y = -df(x)\f(x) ;% solve for increment, similar A\b
% x = x + y ;% add on to get new guess
% f(x) % see if f(x) is really zero
% end
end
function y = F(x)
if nargin==0
    x=rand(1,4);
end
x1 = x(1);
x2 = x(2);
x3=x(3);
x4=x(4);
y = zeros(4,1);
y(1)=(-1/h)+(2*x(1)/h)-(x(2)/h)+(1/6)+(x(1)/6)-(x(1)*x(2)/6)-((x(2)^2)/6)+(h/6)+(2*x(1)*h/3)+(x(2)*h/6)+0.3023;
y(2)=-(x(1)/h)+(2*x(2)/h)-(x(3)/h)+((x(1)^2)/6)-(x(2)*x(3)/6)-((x(3)^2)/6)+(x(1)*h/6)+(2*x(2)*h/3)+(x(3)*h/6)+(x(1)*x(2)/6)+0.451;
y(3)=-(x(2)/h)+(2*x(3)/h)+((x(2)^2)/6)+(x(2)*x(3)/6)-((x(3)*x(4))/6)-((x(4)^2)/6)+(x(2)*h/6)+(2*x(3)*h/3)+(x(4)*h/6)-(x(4)/h)+0.6729;
 y(4)=-(x(3)/h)+(2*x(4)/h)-(exp(1)/h)+((x(3)^2)/6)+(x(3)*x(4)/6)-(x(4)*exp(1)/6)-((exp(1)*exp(1))/6)+(x(3)*h/6)+((2*x(4)*h)/3)+(h*exp(1)/6)+1.004;
%y(1) = 3*x(1)-cos(x(2)*x(3))-0.5; % f1(x1,x2)
%y(2) = x(1).^2-81*(x(2)+0.1)^2 + sin(x(3))+1.06;
%y(3) = exp(-x(1)*x(2))+20*x(3)+((10*pi-3)/3);
end
function x = NewtonMethod(funcF, JacobianF)
format long
F = funcF;
J = JacobianF;
x = [0 0 0 0]';
Iter = 1;
MaxIter = 100;
TOL = 1e-5;
disp([ '      Iter                        x1              x2                     x3             x4          f(x1)                 f(x2)                 f(x3)                f(x4)' ]);
while Iter < MaxIter
y = J(x)\(-F(x));
x = x+y;
M(Iter,:)=[Iter x' F(x)'];
if norm(y,2)<TOL
break;
end
Iter=Iter+1;
end
disp(M)
if Iter >= MaxIter
disp('Maximum number of iteration exceeded!');
end
end
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
Imágen de perfil de JOSE JEREMIAS CABALLERO
Val: 6.975
Oro
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por JOSE JEREMIAS CABALLERO (5917 intervenciones) el 03/04/2020 03:26:10
1
2
3
4
df=@(x)[ (2/h)+1/6-(x(2)/6)+(2*h/3)      -1/h-(x(1)/6)-(x(2)/3)+h/6           0                                            0
         (-1/h)+(x(1)/3)+h/6+(x(2)/6)    2/h-(x(3)/6)+(2*h/3)+(x(1)/6)    -1/h-(x(2)/h)-(x(3)/3)+h/6                   0
               0                           -1/h+(x(2)/3)+(x(3)/6)+h/6    2/h+(x(2)/6)-(x(4)/6)+ (2*h/3) (-x(3)/6)-(x(4)/3)+(h/6)-(1/h)
               0                                   0                      (-1/h)+(x(3)/3)+(x(4)/6)+(h/6)  (2/h)+(x(3)/6)-(exp(1)/6)+(2*h/3) ];
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
sin imagen de perfil
Val: 43
Ha disminuido su posición en 13 puestos en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por wil (6 intervenciones) el 03/04/2020 04:33:26
Muy bien gracias, :)

La verdad habia visto otro codigo de calcular la solucion del sistema no lineal que es el siguiente :


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
function y=newtonsistemaNL
clc;
clear all;
 
h=1/5;
x0=[0;0;0;0];
syms x y z w
fname=[(-1/h)+(2*x/h)-(y/h)+(1/6)+(x/6)-(x*y/6)-((y^2)/6)+(h/6)+(2*x*h/3)+(y*h/6)+0.3023;
         -(x/h)+(2*y/h)-(z/h)+((x^2)/6)-(y*z/6)-((z^2)/6)+(x*h/6)+(2*y*h/3)+(z*h/6)+(x*y/6)+0.451  ;
           -(y/h)+(2*z/h)+((y^2)/6)+(y*z/6)-((z*w)/6)-((w^2)/6)+(y*h/6)+(2*z*h/3)+(w*h/6)-(w/h)+0.6729 ;
        -(z/h)+(2*w/h)-(exp(1)/h)+((z^2)/6)+(z*w/6)-(w*exp(1)/6)-((exp(1)*exp(1))/6)+(z*h/6)+((2*w*h)/3)+(h*exp(1)/6)+1.004];
 fprima=jacobian(fname);
 tolerancia=1.e-10;
 maxiter = 30;
 iter = 1;
 f=inline(fname);
 jf=inline(fprima);
 error=norm(f(x0(1),x0(2),x0(3),x0(4)),2);
 fprintf(' error=%12.8f \n', error);
 while error >= tolerancia
    fx0=f(x0(1),x0(2),x0(3),x0(4));
    fpx0=jf(x0(1),x0(2),x0(3),x0(4));
    x1=x0-inv(fpx0)*fx0;
    fx1=f(x1(1),x1(2),x1(3),x1(4));
    error =norm((fx1),2);
    fprintf(' Iter  %2d raiz x=(%14.9f,%14.9f,%14.9f,%14.9f) f(x)=(%14.9f,%14.9f,%14.9f,%14.9f) \n',iter,x1(1),x1(2),x1(3),x1(4),fx1(1),fx1(2),fx1(3),fx1(4));
    if iter > maxiter
      fprintf(' Numero maximo de iteraciones excedido \n');
      return;
    end
    x0=x1;
    iter=iter+1;
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Pero lo extraño es que aqui, el valor de x(1) = 2.2267   deberia ser : 1.2224 que  se cambio por en de x(4) ,  osea solo falla en el
  %orden.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
Imágen de perfil de JOSE JEREMIAS CABALLERO
Val: 6.975
Oro
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Metodo de Newton Sistemas no lineales

Publicado por JOSE JEREMIAS CABALLERO (5917 intervenciones) el 03/04/2020 05:06:51
Estimado podría resolver en un hoja un problema y así podría entender con mayor magnitud el método de newton para sistemas de ecuaciones no lineales. Y luego de ello, hacer el código en matlab.
Saludos.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar