viernes, 29 de marzo de 2019

Sistema de ecuaciones no lineales graficas, algoritmo y fsolve con matlab u octave


http://www.sc.ehu.es/sbweb/fisica3/numerico/raices/raices_5.html

Sistema de ecuaciones no lineales

Sea el sistema de ecuaciones
{x2y2+2y=02x+y26=0
Representamos gráficamente las funciones implícitas mediante ezplot
>> hold on
>> ezplot('x^2-y^2+2*y',[-8,8,-8,8])
>> ezplot('2*x+y^2-6',[-8,8,-8,8])
>> hold off
>> title('Funciones implícitas')
>> grid on
Observamos que el sitema de ecuaciones tiene cuatro raíces, vamos a obtener estas raíces aplicando el procedimiento de Newton-Raphson

Procedimiento de Newton-Raphson

Para aplicar este procedimiento para un sistema de n ecuaciones con n incógnitas, representamos la variable x por x1 y la variable y por x2. El sistema de dos ecuaciones se escribe de una forma más general
{f1(x1,x2)=0f2(x1,x2)=0
Supongamos que en la etapa k de proceso de cálculo partimos de un punto (x1x2) cualesquiera y nos movemos a otro muy próximo (x1x1x2x2). Los valores de las funciones son f1 y f2 en dicho punto son aproximadamente
f1(x1+Δx1,x2+Δx2)f1(x1,x2)+f1x1Δx1+f1x2Δx2f2(x1+Δx1,x2+Δx2)f2(x1,x2)+f2x1Δx1+f2x2Δx2
Si el punto (x1x1x2x2) es una solución del sistema de ecuaciones, entonces
f1(x1,x2)+f1x1Δx1+f1x2Δx2=0f2(x1,x2)+f2x1Δx1+f2x2Δx2=0
Escribimos el sistema de ecuaciones en forma matricial para despejar Δx1 y Δx2
(f1f2)+(f1x1f1x2f2x1f2x2)(Δx1Δx2)=0
Denominamos vector x al vector (x1,x2), el vector función F está formado por dos elementos que son las funciones (f1,f2) y la matriz cuadrada de dimensión dos es el Jacobiano J. Despejamos Δx1 y Δx2 del sistema de ecuaciones o el vector Δx.
F(x)+JΔx=0Δx=J1F
J-1 es la matriz inversa de J y Δx es el vector diferencia entre el vector que nos da las coordenadas del nuevo punto xk+1, conocidas las del punto previo xk
xk+1=xkJ1F
Hemos obtenido una expresión similar al procedimiento de Newton-Raphson que utilizamos para calcular una raíz de la ecuación f(x)=0
Para un sistema de n ecuaciones
{f1(x1,x2...xn)=0f2(x1,x2...xn)=0...fn(x1,x2...xn)=0
El procedimiento se escribe
(x1'x2'...xn')=(x1x2...xn)(f1x1f1x2...f1xnf2x1f2x2...f2xn............fnx1fnx2...fnxn)1(f1(x1,x2...xn)f2(x1,x2...xn)...fn(x1,x2...xn))
o bien, X=X-J\F, utilizando el operador MATLAB, división por la izquierda \. Se obtiene el nuevo punto, el vector X de la izquierda, a partir del punto X previo, a la derecha de la igualdad
Para el sistema de dos ecuaciones que hemos planteado al principio de esta sección
{x12x22+2x2=02x1+x226=0
El Jacobiano es
J=(2x12x2+222x2)
El código MATLAB para calcular las raíces es similar al empleado para calcular una raíz de la ecuación f(x)=0
F=@(x) [x(1)^2-x(2)^2+2*x(2);2*x(1)+x(2)^2-6];
J=@(x) [2*x(1),-2*x(2)+2;2,2*x(2)];
x=[-4.6;-3.8]; %punto inicial
i=0;
while i<100
    y=-J(x)\F(x);
    if sqrt(norm(y)/norm(x))<0.001
        disp('Solución')
        disp(x)
        break;
    end
    x=x+y;
    i=i+1;
end
if i>=100
    disp('Se ha soprepasado el número de iteracciones');
end
En la primera línea se define el vector columna F de las funciones. En la segunda, la matriz cuadrada J que representa el Jacobiano. En la tercera, las coordenadas del punto inicial. En la representación gráfica de las dos funciones (figura más arriba), utilizamos el icono denominado Data cursor del menú de la ventana gráfica, para conocer el valor aproximado de una de las raíces (0,5, 2.2), que tomaremos como vector inicial de partida para el procedimento de cálculo.
La función norm calcula el módulo un vector. El procedimiento concluye cuando los puntos xmxm+1, están muy próximos, es decir, cuando el error relativo |xm+1-xm|/|xm| sea menor que una cantidad prefijada
Si el procedimiento no convergiera, el programa entraría en un bucle indefinido. Limitando el número de iteracciones se asegura que siempre va a terminar.
    0.6252
    2.1794
>> F(x)
ans =   1.0e-09 *
    0.3991
    0.0536
Cuando introducimos el valor de la raíz buscada x en la función F, vemos que se obtiene un valor próximo a cero
Con Data cursor buscamos un valor próximo a la segunda raíz (-4.6,-3.8). Modificamos en la tercera línea de código el vector inicial que tomamos para calcular la raíz y obtenemos al correr el script modificado
   -4.8642
   -3.9659
>> F(x)
ans =   1.0e-09 *

    0.2845
    0.1122
El problema más importante para aplicar este procedimiento es la elección del vector inicial, cuando el número de ecuaciones es tres o más.

La función fsolve

MATLAB dispone de una función denominada fsolve que realiza un cálculo similar
>> F=@(x) [x(1)^2-x(2)^2+2*x(2);2*x(1)+x(2)^2-6];
>> fsolve(F,[0.5,2.2])
ans =    0.6252    2.1794
>> fsolve(F,[-4.6,-3.8])
ans =   -4.8642   -3.9659
Ejemplo
Calcular las raíces del sistema de dos ecuaciones:
f1(x,y)=2x2-xy-5x-1=0
f2(x,y)=x+3log10x-y2=0
En la figura vemos el punto de intersección y mediante Tools/Data Cursor sus coordenadas aproximadas..
x=linspace(2.5,5.5,50);
y1=(2*x.^2-5*x-1)./x;
y2=sqrt(x+3*log10(x));
plot(x,y1,'b',x,y2,'r')
Vamos a utilizar la función fsolve de MATLAB para obtener el punto de intersección.
Definimos la función F(x), donde x es un vector, x(1) representa a la variable x y x(2) representa la variable y. La primera ecuación, es el primer elemento del vector F y la segunda ecuación, es el segundo elemento de dicho vector
F=@(x) [2*x(1)^2-x(1)*x(2)-5*x(1)-1,  x(1)*3*log10(x(1))-x(2)^2];
x0 =[3 2]; %valor inicial
[x,fval] = fsolve(F,x0); 
fprintf('La solución es x=%1.3f, y=%1.3f\n',x(1),x(2));
fprintf('Valores de la función = %g\n',fval)
fsolve tenemos que pasarle la función que hemos definido y la aproximaxión inicial (x0y0) que hemos encontrado anteriormente de forma gráfica. Esta función nos devuelve las coordenadas (x,y) del punto de intersección buscado y los valores de las funciones f1(x,y) y f2(x,y) en dicho punto.
La solución es x=3.958, y=2.664
Valores de la función = -1.33692e-010
Valores de la función = -3.07305e-009
Ejercicio
Resolver el sistema de ecuaciones no lineales
sin(xy)+exp(xz)0.9=0zx2+y26.7=0tan(yx)+cosz+3.2=0
Tomando x0=1, y0=2 y z0=2 como aproximación inicial

1 comentario:

  1. Your post is very great.I read this post. It’s very helpful. I will definitely go ahead and take advantage of this. You absolutely have wonderful stories. Cheers for sharing with us your blog. For more learning about data science visit at data science course in bangalore

    ResponderEliminar