Sistema de ecuaciones no lineales

Sea el sistema de ecuaciones

{ x 2 y 2 +2y=0 2x+ y 2 6=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

{ f 1 ( x 1 , x 2 )=0 f 2 ( x 1 , x 2 )=0

Supongamos que en la etapa k de proceso de cálculo partimos de un punto (x1, x2) cualesquiera y nos movemos a otro muy próximo (x1x1, x2x2). Los valores de las funciones son f1 y f2 en dicho punto son aproximadamente

f 1 ( x 1 +Δ x 1 , x 2 +Δ x 2 ) f 1 ( x 1 , x 2 )+ f 1 x 1 Δ x 1 + f 1 x 2 Δ x 2 f 2 ( x 1 +Δ x 1 , x 2 +Δ x 2 ) f 2 ( x 1 , x 2 )+ f 2 x 1 Δ x 1 + f 2 x 2 Δ x 2

Si el punto (x1x1, x2x2) es una solución del sistema de ecuaciones, entonces

f 1 ( x 1 , x 2 )+ f 1 x 1 Δ x 1 + f 1 x 2 Δ x 2 =0 f 2 ( x 1 , x 2 )+ f 2 x 1 Δ x 1 + f 2 x 2 Δ x 2 =0

Escribimos el sistema de ecuaciones en forma matricial para despejar Δx1 y Δx2

( f 1 f 2 )+( f 1 x 1 f 1 x 2 f 2 x 1 f 2 x 2 )( Δ x 1 Δ x 2 )=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= J 1 F

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

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

{ f 1 ( x 1 , x 2 ... x n )=0 f 2 ( x 1 , x 2 ... x n )=0 ... f n ( x 1 , x 2 ... x n )=0

El procedimiento se escribe

( x 1 ' x 2 ' ... x n ' )=( x 1 x 2 ... x n ) ( f 1 x 1 f 1 x 2 ... f 1 x n f 2 x 1 f 2 x 2 ... f 2 x n ... ... ... ... f n x 1 f n x 2 ... f n x n ) 1 ( f 1 ( x 1 , x 2 ... x n ) f 2 ( x 1 , x 2 ... x n ) ... f n ( x 1 , x 2 ... x n ) )

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

{ x 1 2 x 2 2 +2 x 2 =0 2 x 1 + x 2 2 6=0

El Jacobiano es

J=( 2 x 1 2 x 2 +2 2 2 x 2 )

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 xm y xm+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)

A fsolve tenemos que pasarle la función que hemos definido y la aproximaxión inicial (x0, y0) 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=0 z x 2 + y 2 6.7=0 tan( y x )+cosz+3.2=0

Tomando x0=1, y0=2 y z0=2 como aproximación inicial