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

$\{\begin{array}{c}{x}^2-y^2+2y=0\\ 2x+y^2-6=0\end{array}$

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 x₁ y la variable y por x₂. El sistema de dos ecuaciones se escribe de una forma más general

$\{\begin{array}{c}{f}_1(x_1,x_2)=0\\ f_2(x_1,x_2)=0\end{array}$

Supongamos que en la etapa k de proceso de cálculo partimos de un punto (x₁, x₂) cualesquiera y nos movemos a otro muy próximo (x₁+Δx₁, x₂+Δx₂). Los valores de las funciones son f₁ y f₂ en dicho punto son aproximadamente

$\begin{array}{c}{f}_1(x_1+\Delta x_1,x_2+\Delta x_2)\approx f_1(x_1,x_2)+\frac{\partial f_1}{\partial x_1}\Delta x_1+\frac{\partial f_1}{\partial x_2}\Delta x_2\\ f_2(x_1+\Delta x_1,x_2+\Delta x_2)\approx f_2(x_1,x_2)+\frac{\partial f_2}{\partial x_1}\Delta x_1+\frac{\partial f_2}{\partial x_2}\Delta x_2\end{array}$

Si el punto (x₁+Δx₁, x₂+Δx₂) es una solución del sistema de ecuaciones, entonces

$\begin{array}{c}{f}_1(x_1,x_2)+\frac{\partial f_1}{\partial x_1}\Delta x_1+\frac{\partial f_1}{\partial x_2}\Delta x_2=0\\ f_2(x_1,x_2)+\frac{\partial f_2}{\partial x_1}\Delta x_1+\frac{\partial f_2}{\partial x_2}\Delta x_2=0\end{array}$

Escribimos el sistema de ecuaciones en forma matricial para despejar Δx₁ y Δx₂

$\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right)+\left(\begin{array}{cc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}\end{array}\right)\left(\begin{array}{c}\Delta x_1\\ \Delta x_2\end{array}\right)=0$

Denominamos vector x al vector (x₁,x₂), el vector función F está formado por dos elementos que son las funciones (f₁,f₂) y la matriz cuadrada de dimensión dos es el Jacobiano J. Despejamos Δx₁ y Δx₂ del sistema de ecuaciones o el vector Δx.

$\begin{array}{c}F\left(x\right)+J\Delta x=0\\ \Delta x=-J^{-1}F\end{array}$

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 x_k+1, conocidas las del punto previo x_k

$x_{k+1}=x_k-J^{-1}F$

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

$\{\begin{array}{c}{f}_1(x_1,x_2\mathrm{...}{x}_n)=0\\ f_2(x_1,x_2\mathrm{...}{x}_n)=0\\ \mathrm{...}\\ f_n(x_1,x_2\mathrm{...}{x}_n)=0\end{array}$

El procedimiento se escribe

$\left(\begin{array}{c}{x}_1^{\prime }\\ x_2^{\prime }\\ \mathrm{...}\\ x_n^{\prime }\end{array}\right)=\left(\begin{array}{c}{x}_1\\ x_2\\ \mathrm{...}\\ x_n\end{array}\right)-{\left(\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \mathrm{...}& \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \mathrm{...}& \frac{\partial f_2}{\partial x_n}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& \mathrm{...}& \frac{\partial f_n}{\partial x_n}\end{array}\right)}^{-1}\left(\begin{array}{c}{f}_1(x_1,x_2\mathrm{...}{x}_n)\\ f_2(x_1,x_2\mathrm{...}{x}_n)\\ \mathrm{...}\\ f_n(x_1,x_2\mathrm{...}{x}_n)\end{array}\right)$

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

$\{\begin{array}{c}{x}_1^2-x_2^2+2x_2=0\\ 2x_1+x_2^2-6=0\end{array}$

El Jacobiano es

$J=\left(\begin{array}{cc}2x_1& -2x_2+2\\ 2& 2x_2\end{array}\right)$

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 x_my x_m+1, están muy próximos, es decir, cuando el error relativo |x_m+1-x_m|/|x_m| 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:

f₁(x,y)=2x²-xy-5x-1=0
f₂(x,y)=x+3log₁₀x-y²=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)

A fsolve tenemos que pasarle la función que hemos definido y la aproximaxión inicial (x₀, y₀) 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 f₁(x,y) y f₂(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

$\begin{array}{c}\sin \left(xy\right)+\exp (-xz)-0.9=0\\ z\sqrt{x^2+y^2}-6.7=0\\ \tan \left(\frac{y}{x}\right)+\cos z+3.2=0\end{array}$

Tomando x₀=1, y₀=2 y z₀=2 como aproximación inicial