[Livermorio]

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

[Piockñec]

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 saber detectar exactamente por qué te falla tu programa, y eventualmente saber corregirlo, se basa en la experiencia de programar!!! ¿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

[Livermorio]

¿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.

[Livermorio]

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.

 

[Piockñec]

jajaja vaya!!! Los errores más tontos suelen ser los que mejor se esconden enhorabuena