Sistema de ecuaciones no lineales
Sea el sistema de ecuaciones
Representamos gráficamente las funciones implícitas mediante
>> 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
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 (x1+Δx1, x2+Δx2). Los valores de las funciones son f1 y f2 en dicho punto son aproximadamente
Si el punto (x1+Δx1, x2+Δx2) es una solución del sistema de ecuaciones, entonces
Escribimos el sistema de ecuaciones en forma matricial para despejar Δx1 y Δx2
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.
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
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
El procedimiento se escribe
o bien, X'=X-J\F, utilizando el operador MATLAB, división por la izquierda \. Se obtiene el nuevo vector X' de la izquierda, a partir del X previo, a la derecha de la igualdad
Para el sistema de dos ecuaciones que hemos planteado al principio de esta sección
El Jacobiano es
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=[0.5;2.2]; %punto inicial (vector columna ;) %x=[-4.6;-3.8]; %punto inicial (vector columna ;) for i=1:100 dx=-J(x)\F(x); if sqrt(norm(dx)/norm(x))<0.001 disp('Solución') disp(x+dx) break; end x=x+dx; 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
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
>> 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
Definimos la función F(x), donde x es un 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
La solución es x=3.809, y=2.356 Valores de la función = 1.18398e-10 Valores de la función = -1.65512e-11
Ejercicios
Resolver el sistema de ecuaciones no lineales
Comprobando que su solución es (1,0)
Resolver el sistema de ecuaciones no lineales
Tomando x0=1, y0=2 y z0=2 como aproximación inicial