26/01/2020, 00:38:46 am *
Bienvenido(a), Visitante. Por favor, ingresa o regístrate.
¿Perdiste tu email de activación?

Ingresar con nombre de usuario, contraseña y duración de la sesión
Noticias: Puedes practicar LATEX con el cómodo editor de Latex online
 
 
Páginas: [1]   Ir Abajo
  Imprimir  
Autor Tema: Mensaje de error en un programa para resolver sistema de ecuaciones no lineales  (Leído 1306 veces)
0 Usuarios y 1 Visitante están viendo este tema.
Livermorio
Semi pleno
***

Karma: +0/-0
Desconectado Desconectado

España España

Mensajes: 51


Ver Perfil
« : 05/05/2015, 16:10:23 pm »

Me gustaría saber por qué recibo un mensaje diciendo que la matriz que se utiliza para resolver el problema es singular para la precisión de la máquina y cómo podría arreglarlo.

El cuadradito que sale no sé cómo arreglarlo, decídmelo por favor.

% Function to solve a non-linear system of equations
function
  • =SENL(h)
n=3;
x=(rand(n,1)); % Random column vector
tol=10^-4;     % Relative tolerance of the solution
v=zeros(n,1);  % Zeros vector intended to create a canonical vector
v(1)=1;        %with this second line.
e=(2*tol*norm(x))*v; %This is the variable we use to solve the problem.
                     % I start it with double the value of tol*norm(v) so
                     % that the 'while' loop can start
i=0;                 % Variable 'i' for the iterations
while norm(e)>tol*norm(x) && i<50 % 2 different ways of exiting the loop
    i=i+1;
    f=fun(x); %----> It evaluates the function and returns a column vector
    J=funjac(x,h,n); %----> It calculates the jacobian using incremental
                     %      quotients. By numerical methods.
    e=J\f;           % This is used to calculate e=J^(-1)*f
    x=x-e;           % It updates the solution 'x'
    norm(e)          % norm(e) gives a smaller number with every iteration.
                     % With each iteration 'x' is closer to being the
                     % solution.
end
end

function [f]=fun(x) % It evaluates a f, which is a vector-valued function
                    % with vector-valued variables
f1=x(1)^2 +4*x(2) +3*x(3);
f2=x(1)   -x(2)^2 +x(3);
f3=9*x(1) +3*x(2) -x(3)^2;
f=[f1; f2; f3];
end

function [J]=funjac(x,h,n) % funjac returns the jacobian
                           % (differential matrix of the function)
I=ones(n);
for k=1:n
    alpha(k)=max(1,abs(x(k)))*sign(x(k));
    finc=fun(x+(alpha(k)*h*I(:,k)));
    f=fun(x);
    J(:,k)=(finc-f)/(alpha(k)*h);
end
end

* SENL.txt (1.63 KB - descargado 57 veces.)
En línea
Piockñec
Pleno*
*****

Karma: +1/-0
Desconectado Desconectado

Sexo: Masculino
France France

Mensajes: 1.259


Ver Perfil
« Respuesta #1 : 05/05/2015, 18:50:24 pm »

No tengo tiempo de mirarte el funcionamiento del programa, pero parece un newton raphson.

A veces existen puntos en los que el jacobiano tiene determinante 0, y el método falla (hay que buscarse otro punto inicial, o a veces ni por esas), pero me huelo que es más probable que sea por la implementación si no sabes corregirlo :cara_de_queso: saber detectar exactamente por qué te falla tu programa, y eventualmente saber corregirlo, se basa en la experiencia de programar!!! :guiño: ¿Has programado ese código tú solo?

Comprueba con el debugger (poniendo circulitos rojos donde quieras investigar al margen del código) que la función funjac te hace el jacobiano bien.
Una vez hecho eso, comprueba que las x vayan aproximándose cada vez más a una solución, para ver si es por lo que te he dicho (que en alguna x de alguna iteración el jacobiano se haga 0, y el método no se pueda usar con esa semilla).
¡Suerte!

Lo del cuadradito no lo entiendo. En matlab las funciones se definen:

function x=SENL(h)

...

end

donde esa h no veo que la uses en tu código (?!) y donde la x sería el resultado final, la solución del sistema.
Ah ya veo, h modifica la aproximación del jacobiano. Usa un número pequeño, pero no demasiado pequeño: 1e-5 ó 1e-6
En línea
Livermorio
Semi pleno
***

Karma: +0/-0
Desconectado Desconectado

España España

Mensajes: 51


Ver Perfil
« Respuesta #2 : 05/05/2015, 21:27:54 pm »

¿Has programado ese código tú solo?

Con el algoritmo al lado y teniendo la expresión para el cálculo del jacobiano por cocientes incrementales, pero sí, tampoco es tan complicado.

Gracias por contestar, aprecio mucho que te hayas tomado la molestia.

En línea
Livermorio
Semi pleno
***

Karma: +0/-0
Desconectado Desconectado

España España

Mensajes: 51


Ver Perfil
« Respuesta #3 : 09/05/2015, 14:09:49 pm »

Piockñec, ya sé dónde tengo mal el script. En funjac, puse I=ones(n), cuando en realidad quería hacer I=eye(n) para sacar los vectores canónicos.

Después de descubrir eso y haciendo dos cambios más sin importancia, la función SENL funciona bien.

 Aplauso
En línea
Piockñec
Pleno*
*****

Karma: +1/-0
Desconectado Desconectado

Sexo: Masculino
France France

Mensajes: 1.259


Ver Perfil
« Respuesta #4 : 09/05/2015, 22:32:49 pm »

jajaja vaya!!! Los errores más tontos suelen ser los que mejor se esconden :lengua_afuera: enhorabuena :guiño:
En línea
Páginas: [1]   Ir Arriba
  Imprimir  
 
Ir a:  

Impulsado por MySQL Impulsado por PHP Powered by SMF 1.1.4 | SMF © 2006, Simple Machines LLC XHTML 1.0 válido! CSS válido!