% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root =  ../Metodos_Numericos.tex

\chapter[Sist. no lineales]{Solución de Ecuaciones No lineales}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introducción}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Una ecuación no lineal es aquella que tiene una forma diferente a $f(x) = a_0 + a_1 x$ en cuyo caso calcular la solución consiste en resolver despejando  de $x=-a_0/a_1$. Pero el caso es que queremos resolver un sistema de ecuaciones no lineales de la forma

\begin{eqnarray*}
f_1(x_1, x_2, \cdots, x_N) & = & 0\\ 
f_2(x_1, x_2, \cdots, x_N) & = & 0\\
\vdots & = & \vdots \\
f_n(x_1, x_2, \cdots, x_N) & = & 0
\end{eqnarray*}

donde $N$ es el número de ecuaciones y se tiene el mismo números de funciones $f_i$  que variables $x_i$. En esta sección veremos dos de estos métodos que son la iteración de punto fijo y el método de Newton-Raphson, los cuales se pueden utilizar para resolver sistemas de ecuaciones no lineales

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Métodos de punto fijo}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Considerando un sistema no lineal de ecuaciones

\begin{eqnarray*}
f_1(x_1, x_2, \cdots, x_n) & = & 0\\ 
f_2(x_1, x_2, \cdots, x_n) & = & 0\\
\vdots & = & \vdots \\
f_n(x_1, x_2, \cdots, x_n) & = & 0
\end{eqnarray*}

El método de iteración de punto fijo intentara despejar de cada una de la las ecuaciones $f_i(x_1, x_2, \cdots, x_n)$ la $i$-esima variable tal que

\begin{eqnarray*}
x_1  -  g_1(x_1, x_2, \cdots, x_n) & = & f_1(x_1, x_2, \cdots, x_n) \\ 
x_2 -  g_2(x_1, x_2, \cdots, x_n) & = & f_2(x_1, x_2, \cdots, x_n) \\
\vdots & = & \vdots \\
x_n -  g_n(x_1, x_2, \cdots, x_n) & = &f_n(x_1, x_2, \cdots, x_n) 
\end{eqnarray*}

Para resolver de manera iterativa

\begin{eqnarray*}
x_1^{(t+1)}  & = &  g_1(x_1^{(t)}, x_2^{(t)}, \cdots, x_n^{(t)})  \\ 
x_2^{(t+1)} & = &  g_2(x_1^{(t)}, x_2^{(t)}, \cdots, x_n^{(t)})  \\
\vdots & = & \vdots \\
x_n^{(t+1)} & = &  g_n(x_1^{(t)}, x_2^{(t)}, \cdots, x_n^{(t)})
\end{eqnarray*}

La implementación en Matlab es

\begin{verbatim}
function r = Punto_Fijo(g, x1) 

while 1
    x2 = g(x1)
    error = abs(norm((x2-x1)/x2));
    if(error < 0.0001) break
    else x1 = x2;
    end;
end;
r = x2;
\end{verbatim}

\subsubsection{Ejemplo 1} 

Utilice el método de iteración de punto fijo para determinar las raíces del sistema de ecuaciones dado. Considere como valores iniciales $x_1^{(0)} = 1.5$ y $x_2^{(0)}=3.5$.

\begin{eqnarray*}
x_1^2 + x_1 x_2 -10 & = & 0\\
x_2 + 3 x_1 x_2^2 - 57 & = & 0 
\end{eqnarray*}

De acuerdo con lo descrito despejamos $x_1$ y $x_2$ 

\begin{eqnarray*}
f_1(x_1, x_2) & = & x_1 -  \sqrt{10 -  x_1 x_2 }\\
f_2(x_1, x_2) & = & x_2 - \sqrt{\frac{57 - x_2}{3 x_1} }
\end{eqnarray*}

Despejando tenemos el siguiente sistema iterativo de ecuaciones

\begin{eqnarray*}
x_1^{(t+1)} & = &  \sqrt{10 -  x_1^{(t)} x_2^{(t)}  }\\
x_2^{(t+1)} & = & \sqrt{\frac{57 - x_2^{(t)}}{3 x_1^{(t+1)}} }
\end{eqnarray*}

La implementación en Matlab es

\begin{verbatim}
function x = g1(x)

x(1) = sqrt(10 - x(1)*x(2));
x(2) = sqrt((57 - x(2))/(3 * x(1)));
\end{verbatim}

El proceso iterativo será:

\textbf{Iteración 1}

Con $x^{(0)} = [1.5, 3.5]^T$ tenemos

\begin{eqnarray*}
x_1^{(1)} & = &  \sqrt{10 -  x_1^{(0)} x_2^{(0)}  } = \sqrt{10 - (1.5)(3.5)} = 2.1794\\
x_2^{(1)} & = & \sqrt{\frac{57 - x_2^{(0)}}{3 x_1^{(1)}} } = \sqrt{\frac{57 - (3.5)}{3 (2.1794)} } = 2.8605
\end{eqnarray*}

\textbf{Iteración 2}

Con $x^{(1)} = [2.1794, 2.8605]^T$ tenemos

\begin{eqnarray*}
x_1^{(2)} & = &  \sqrt{10 -  x_1^{(1)} x_2^{(1)}  } = \sqrt{10 - (2.1794)(2.8605)} = 1.9405\\
x_2^{(2)} & = & \sqrt{\frac{57 - x_2^{(0)}}{3 x_1^{(1)}} } = \sqrt{\frac{57 - (2.8605)}{3 (1.9405)} } = 3.0496
\end{eqnarray*}

El resumen el proceso iterativo se muestra en la siguiente tabla, donde podemos ver que la solución es $x^* = [2, 3]^T$ en 9 iteraciones.

\begin{center}
\begin{tabular}{|r|r|r|}
    \hline
    $k$ & $x_1^{(k)}$ & $x_2^{(k)}$ \\
    \hline
    0  &  1.5000  &  3.5000\\
    1  &  2.1794  &  2.8605\\
    2  &  1.9405  &  3.0496\\
    3  &  2.0205  &  2.9834\\
    4  &  1.9930  &  3.0057\\
    5  &  2.0024  &  2.9981\\
    6  &  1.9992  &  3.0007\\
    7  &  2.0003  &  2.9998\\
    8  &  1.9999  &  3.0001\\
    9  &   2.0000 &   3.0000\\
    \hline
\end{tabular}
\end{center}

Para realizar la ejecución dar en Matlab 

\begin{verbatim}
Punto_Fijo(@g1, [1.5; 3.5])
\end{verbatim}

\subsubsection{Ejemplo 2} 

Utilice el método de iteración de punto fijo para determinar las raíces del sistema de ecuaciones dado. Considere como valores iniciales $x_1^{(0)} = 0.0$ y $x_2^{(0)}=1.0$.

\begin{eqnarray*}
x_1^2 - 2 x_1 - x_2 + 0.5 & = & 0\\
x_1^2 + 4 x_2^2 - 4 & = & 0 
\end{eqnarray*}

De acuerdo con lo descrito despejamos $x_1$ y $x_2$ 

\begin{eqnarray*}
f_1(x_1, x_2) & = & x_1 -  \sqrt{2 x_1 + x_2 - 0.5  }\\
f_2(x_1, x_2) & = & x_2 - \sqrt{\frac{4 - x_1^2}{4} }
\end{eqnarray*}

Despejando tenemos el siguiente sistema iterativo de ecuaciones

\begin{eqnarray*}
x_1^{(t+1)} & = &  \sqrt{2 x_1^ {(t)} +  x_2^{(t)} - 0.5  }\\
x_2^{(t+1)} & = & \sqrt{\frac{4 - \left(x_1^{(t+1)} \right)^2}{4} }
\end{eqnarray*}

La implementación en Matlab es

\begin{verbatim}
function x = g2(x)

x(1) = sqrt(2*x(1) + x(2) - 0.5);
x(2) = sqrt(4 -x(1)^2)/2;  
\end{verbatim}

El proceso iterativo será:

\textbf{Iteración 1}

Con $x^{(0)} = [0,1]^T$ tenemos

\begin{eqnarray*}
x_1^{(1)} & = &  \sqrt{2 x_1^ {(0)}+  x_2^{(0)} - 0.5  } = \sqrt{2 (0) + (1) - 0.5  } = 0.7071 \\
x_2^{(1)} & = & \sqrt{\frac{4 - \left(x_1^{(1)} \right)^2}{4} }= \sqrt{\frac{4 - (0.7071)^2}{4} } = 0.9354
\end{eqnarray*}

\textbf{Iteración 2}

Con $x^{(1)} = [0.7071, 0.9354]^T$ tenemos

\begin{eqnarray*}
x_1^{(2)} & = &  \sqrt{2 x_1^ {(1)}+  x_2^{(1)} - 0.5  } = \sqrt{2 (0.7071) + (0.9354) - 0.5  } = 1.3600 \\
x_2^{(2)} & = & \sqrt{\frac{4 - \left(x_1^{(2)} \right)^2}{4} }= \sqrt{\frac{4 - (1.3600)^2}{4} } = 0.7332
\end{eqnarray*}

El resumen el proceso iterativo se muestra en la siguiente tabla, donde podemos ver que la solución es $x^* = [1.9007, 0.3112]^T$ en 8 iteraciones.

\begin{center}
\begin{tabular}{|r|r|r|}
    \hline
    $k$ & $x_1^{(k)}$ & $x_2^{(k)}$ \\
    \hline
    0  &  0.0000  & 1.0000\\
    1  &  0.7071  &  0.9354\\
    2  & 1.3600   &  0.7332\\
    3  & 1.7185  &   0.5116\\
    4  & 1.8570  &   0.3713\\
    5  & 1.8935  &  0.3220\\
    6  & 1.8997  &  0.3127\\
    7  & 1.9006  &  0.3114 \\
    8  & 1.9007  &  0.3112\\
    \hline
\end{tabular}
\end{center}

Para realizar la ejecución dar en Matlab 

\begin{verbatim}
Punto_Fijo(@g2, [0; 1])
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\section{Método iterativo de Jacobi y Método iterativo de Gauss-Seidel}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{El método de Bisecciones}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Para este método debemos considerar una función unidimensional $f:\mathbb{R} \rightarrow \mathbb{R}$ continua dentro de un intervalo $[a,b]$ tal que $f(a)$ tenga diferente signo $f(a)*f(b) <0$.

El proceso de decisión para encontrar la raíz consiste en dividir el intervalo $[inicio, fin]$ a la mitad $mitad = (inicio+fin)/2$ y luego analizar las tres posibilidades que se pueden dar.

\begin{enumerate}
\item Si $f(inicio)$ y $f(mitad)$ tienen signos opuestos, entonces hay un cero entre $[inicio, mitad]$.
\item Si $f(mitad)$ y $f(fin)$ tienen signos opuestos, entonces, hay un cero en $[mitad, fin]$.
\item Si $f(mitad)$ es igual a cero, entonces $mitad$ es un cero
\end{enumerate}

La implementación en Matlab es:

\begin{verbatim}
function r = Biseccion(f, inicio, fin)
 mitad = 0;
 while abs((fin - inicio)/fin) > 0.0001
     mitad = (fin+inicio)/2.0;

   if(f(mitad) == 0) 
       r = mitad;
       return;
   end;

   if f(inicio)*f(mitad) < 0 
     fin = mitad;
   else
     inicio = mitad;
   end;
 end;
     
 r= mitad;
 \end{verbatim}


\subsection{Ejemplo}

Calcular los ceros de la función $f(x) = x-cos(x)$ utilizando el algoritmo de Bisección en el intervalo $[0,1]$.

\begin{center}
\begin{tabular}{|r|r|r|r|r|r|r|}
\hline
iter & a & c & b & f(a) & f( c)& f(b)\\
\hline
0 & 0.0000 & 0.5000 & 1.0000 & -1.0000 & -0.3776 & 0.4597\\
1 & 0.5000 & 0.7500 & 1.0000 & -0.3776 & 0.0183 & -0.3776 \\
2 & 0.5000 & 0.6250 & 0.7500 & -0.3776 & -0.1860 & -0.3776\\
3 & 0.6250 & 0.6875 & 0.7500 & -0.1860 & -0.0853 & -0.1860\\
4 & 0.6875 & 0.7188 & 0.7500 & -0.0853 & -0.0339 & -0.0853\\
5 & 0.7188 & 0.7344 & 0.7500 & -0.0339 & -0.0079 & -0.0339\\
6 & 0.7344 & 0.7422 & 0.7500 & -0.0079 & 0.0052 & -0.0079\\
7 & 0.7344 & 0.7383 & 0.7422 & -0.0079 & -0.0013 & -0.0079\\
8 & 0.7383 & 0.7402 & 0.7422 & -0.0013 & 0.0019 & -0.0013\\
9 & 0.7383 & 0.7393 & 0.7402 & -0.0013 & 0.0003 & -0.0013\\
10 & 0.7383 & 0.7388 & 0.7393 & -0.0013 & -0.0005 & -0.0013\\
11 & 0.7388 & 0.7390 & 0.7393 & -0.0005 & -0.0001 & -0.0005 \\
12 & 0.7390 & 0.7391 & 0.7393 & -0.0001 & 0.0001 & -0.0001\\
13 & 0.7390 & 0.7391 & 0.7391 & -0.0001 & 0.0000 & -0.0001\\
\hline
\end{tabular}
\end{center}

Para ver los resultados correr 
\begin{verbatim}
	Biseccion(@f2, 0, 1)
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Método de Regula Falsi}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Una de las razones de la introducción de este método es que la velocidad de convergencia del método de Bisecciones es bastante baja. En el método de Bisección se usa el punto medio del intervalo $[a,b]$ para llevar a cabo el siguiente paso. Suele conseguirse una mejor aproximación usando el punto $[c,0]$ en el que la recta secante $L$, que pasa por los puntos $[a, f(a)]$ y $[b, f(b)]$.

En la Figura \ref{fig:Regula_Falsi} se puede ver como funciona el método. En esta figura en azul esta la función de la cual queremos calcular el cruce por cero y en negro dos líneas rectas que aproximan la solución.

 \begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/RF.eps, width=10.0cm}}
   \caption{Método Regula Falsi}
   \label{fig:Regula_Falsi}
\end{figure}

 
Para calcula la ecuación de la línea secante hacemos
\begin{eqnarray*}
p1 & = & [a, f(a)] \\
p2 &=& [b,f(b)]
\end{eqnarray*}

y sustituimos en la ecuación de la línea recta.

\[
y – f(a) = (f(b)– f(a))*(x – a)/(b – a)
\]
 

El cruce por cero de esta ecuación esta en

\[
c = a – f(a)*(b-a)/(f(b) – f(a))
\]
 

Entonces el método de Bisecciones puede ser modificado, en lugar de calcular $c = (a+b)/2$ hacemos $c = a – f(a)*(b-a)/(f(b) – f(a))$ y aplicamos las mismas tres reglas de la Bisección. La implementación en Matlab es:

\begin{verbatim}
function r = Regula_Falsi(f, a, b)

 while abs((b - a)/a) > 0.0001
     c = a - f(a)*(b-a)/(f(b) - f(a));

   if(f(c) == 0) 
       r = c;
       return;
   end;

   if f(a)*f(c) < 0 
     b = c;
   else
     a = c;
   end;
 end;
     
 r= c;
 \end{verbatim}

Para ejecutar hacer

\begin{verbatim}
Regula_Falsi(@f2, 0, 1)
\end{verbatim}

\subsection{Ejemplo}

Calcular los ceros de la función $f(x) = x-cos(x)$ utilizando el algoritmo de regula falsi en el intervalo $[0,1]$.


\begin{center}
\begin{tabular}{|r|r|r|r|r|r|r|}
\hline
iter. & a & c & b & f(a) & f(c ) &f(b) \\
\hline
0 & 0.0000 & 0.6851 & 1.0000 & -1.0000 & -0.0893 & 0.4597 \\
1 & 0.6851 & 0.7363 & 1.0000 & -0.0893 & -0.0047 & -0.0893 \\
2 & 0.7363 & 0.7389 & 1.0000 & -0.0047 & -0.0002 & -0.0047 \\
3 & 0.7389 & 0.7391 & 1.0000 & -0.0002 & 0.0000 & -0.0002\\
4 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.0000\\
5 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.0000\\
6 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.0000\\
\hline
\end{tabular}
\end{center}

\subsection{Solución de un circuito con un diodo}

Consideremos un circuito de corriente alterna formado por una fuente, una resistencia y un diodo tal como se muestra en la figura \ref{fig:Cto_diodo}.

 \begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/diodo.eps, width=10.0cm}}
   \caption{Circuito de corriente alterna con un diodo}
   \label{fig:Cto_diodo}
\end{figure}

La ecuación que modela al diodo esta dada por \eqref{eq:diodo} y la curva se muestra en la figura \ref{fig:Comp_diodo}:

\begin{equation}
i_d = I_s(e^{v_d/V_t} -1)
\label{eq:diodo}
\end{equation}

\noindent
donde $I_s = 1^{-12}$ y  $V_t = 25.85^{-3}$

 \begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/comp_diodo.eps, width=10.0cm}}
   \caption{Comportamiento del diodo}
   \label{fig:Comp_diodo}
\end{figure}

Para este ejemplo la ecuación que hay que resolver es:

\begin{equation*}
		V(t) - R i(t) - V_d(t) = 0 \\
\end{equation*}

\begin{equation*}
		V_d(t) = \left \{
			\begin{array}{ccc}
			log(I_d(t)/I_s +1) * V_t & Si &I_d(t) \geq 0 \\
			\\
			V & si no & Id(t)< 0 \\ 
			\end{array}
		\right .
\end{equation*}

La solución implementada se muestra en el siguiente código y la solución se muestra en la Fig. \ref{fig:Sol_cto_diodo}

\begin{verbatim}
function Simula_Circuito_Diodo()
    % Paramentros del Cicuito
    
    Vmax = 10;  % Volts;
    R    =  1.8; % ohms
    
    %tiempo de evaluación
    t = [0:0.001:0.1];
    %Voltaje en función del tiempo
    v = Vmax*sin(377*t);
    
    for k = 1:length(t)
        I(k) = ReglaFalsa(@Circuito_Diodo, -Vmax/R, Vmax/R, [v(k), R]);
    end
    
    plot(t, v, t, I);  
    title('Solución de un circuito con Diodo');
    xlabel('tiempo s');
    legend('Voltaje', 'Corriente');
end


function di = Circuito_Diodo(I, p)
   V = p(1);
   R = p(2);
   di = V - R*I - Vd(V, I);
end

function r = ReglaFalsa(f, a, b, p)
%ReglaFalsa(f, a, b, p)
%   f : función a evaluar
%   a : inicio del intervalo
%   b : fin frl intervalo
%   p : parametros

 while 1 
     c = a - f(a,p)*(b-a)/(f(b,p) - f(a,p));
     if(abs(f(c,p)) <= 1e-06)
         break;
     end;
     if f(a,p)*f(c,p) < 0 
         b = c;
     else
         a = c;
   end;
 end;
 r = c;
end

function V = Vd(Vf, I)
% V = Vd(Vf, I)
%   Vf voltaje de la Fuente
%   I Corriente en el circuito

  Is = 1e-12;
  Vt = 25.85e-3;
  
  if I>=0
      V = log(I/Is+1)*Vt;
  else
      V = Vf;
  end
end
\end{verbatim}

 \begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/sol_cto_diodo.eps, width=10.0cm}}
   \caption{Solución al circuito con diodo}
   \label{fig:Sol_cto_diodo}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{El método Newton Raphson}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Consideremos el sistema de ecuaciones no lineales

\begin{eqnarray*}
f_1(x_1, x_2, \cdots, x_n) & = & 0\\ 
f_2(x_1, x_2, \cdots, x_n) & = & 0\\
\vdots & = & \vdots \\
f_n(x_1, x_2, \cdots, x_n) & = & 0
\end{eqnarray*}

Utilizando la serie de Taylor podemos hacer una aproximación lineal para una función $f_i(x_1^{(t)} +\delta x_1^{(t)}, x_2^{(t)} +\delta x_2^{(t)}, \cdots, x_n^{(t)} + \delta x_n^{(t)})$ en un incremento $\delta x_i^{(t)} = x_i^{(t+1)}- x_i^{(t)}$ como:

\begin{eqnarray*}
f_1(x_1^{(t)} +\delta x_1^{(t)}, x_2^{(t)} +\delta x_2^{(t)}, \cdots, x_n^{(t)} + \delta x_n^{(t)}) \\
= f_1(x_1^{(t)}, x_2^{(t)} + \cdots, + x_n^{(t)})  + \frac{\partial f_1}{\partial x_1} \delta x_1^{(t)} + \frac{\partial f_1}{\partial x_2} \delta x_2^{(t)} + \cdots + \frac{\partial f_1}{\partial x_n} \delta x_n^{(t)}\\
f_2(x_1^{(t)} +\delta x_1^{(t)}, x_2^{(t)} +\delta x_2^{(t)}, \cdots, x_n^{(t)} + \delta x_n^{(t)}) \\
= f_2(x_1^{(t)}, x_2^{(t)} + \cdots, + x_n^{(t)})  + \frac{\partial f_2}{\partial x_1} \delta x_1^{(t)} + \frac{\partial f_2}{\partial x_2} \delta x_2^{(t)} + \cdots + \frac{\partial f_2}{\partial x_n} \delta x_n^{(t)}\\
\vdots \\
f_n(x_1^{(t)} +\delta x_1^{(t)}, x_2^{(t)} +\delta x_2^{(t)}, \cdots, x_n^{(t)} + \delta x_n^{(t)}) \\
= f_n(x_1^{(t)}, x_2^{(t)} + \cdots, + x_n^{(t)})  + \frac{\partial f_n}{\partial x_1} \delta x_1^{(t)} + \frac{\partial f_n}{\partial x_2} \delta x_2^{(t)} + \cdots + \frac{\partial f_n}{\partial x_n} \delta x_n^{(t)}\\
\end{eqnarray*}

Si escribimos el sistema en forma matricial tenemos

\begin{eqnarray*}
\left[
\begin{array}{c}
f_1(x_1^{(t)} + \delta x_1^{(t)}, x_2^{(t)}+ \delta x_2^{(t)}, \cdots,  x_n^{(t)}+ \delta x_n^{(t)} )  \\
f_2(x_1^{(t)} + \delta x_1^{(t)}, x_2^{(t)}+ \delta x_2^{(t)}, \cdots,  x_n^{(t)}+ \delta x_n^{(t)} )  \\
\vdots\\ 
f_n(x_1^{(t)} + \delta x_1^{(t)}, x_2^{(t)}+ \delta x_2^{(t)}, \cdots,  x_n^{(t)}+ \delta x_n^{(t)} )  \\
\end{array}
\right] = \\
\left[
\begin{array}{c}
f_1(x_1^{(t)}, x_2^{(t)}, \cdots,  x_n^{(t)} )  \\
f_2(x_1^{(t)}, x_2^{(t)}, \cdots,  x_n^{(t)} )  \\
\vdots\\ 
f_n(x_1^{(t)}, x_2^{(t)}, \cdots,  x_n^{(t)} )  \\
\end{array}
\right] +
\left[
\begin{array}{cccc}
\frac{\partial f_1(x)}{\partial x_1}  &  \frac{\partial f_1(x)}{\partial x_2}  & \cdots & \frac{\partial f_1(x)}{\partial x_n}   \\
\frac{\partial f_2(x)}{\partial x_1}  &  \frac{\partial f_2(x)}{\partial x_2}  & \cdots & \frac{\partial f_2(x)}{\partial x_n}   \\
\vdots  & \vdots  & \cdots & \vdots \\
\frac{\partial f_n(x)}{\partial x_1}  &  \frac{\partial f_n(x)}{\partial x_2}  & \cdots & \frac{\partial f_n(x)}{\partial x_n}   \\  
\end{array}
\right]
\left[
\begin{array}{c}
\delta x_1^{(t)} \\
\delta x_2^{(t)} \\
\vdots \\
\delta x_n^{(t)} \\
\end{array}
\right]
\end{eqnarray*}

En forma compacta

\begin{eqnarray*}
f(x^{(t)} + \delta x^{(t)}) &=& f(x^{(t)}) + J(x^{(t)}) \delta x^{(t)} \\
f(x^{(t)} + \delta x^{(t)}) &=& f(x^{(t)}) + J(x^{(t)}) (x^{(t+1)} - x^{(t)} ) \\
\end{eqnarray*}


donde $J(x)$ es la matriz de primeras derivada o Jacobiano.

Como queremos encontrar el cero de la función,  la aproximación lineal que debemos resolver es:

\begin{equation*}
0 = f(x^{(t)}) + J(x^{(t)})  (x^{(t+1)} - x^{(t)} )\\
\end{equation*} 

donde las actualizaciones las hacemos como

\begin{equation*}
x^{(t+1)}  =  x^{(t)}  -  \left[J(x^{(t)})\right]^{-1} f(x^{(t)})
\end{equation*} 

La implementación del algoritmo en Matlab es:

\begin{verbatim}
function r = Newton_Raphson(f, J, x1)

while 1
    x2 = x1 - inv(J(x1))* f(x1);
    if norm((x1-x2)./x2) < 0.0001 break;
    else x1 = x2;
    end;
end;
r = x2;
\end{verbatim}

\subsection{Ejemplo 1}

Resolver el siguiente sistema de ecuaciones dado por y cuya gráfica se muestra en la Figura \ref{fig:SNL-01}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/SNL-01.eps,width=8.0cm}}
   \caption{Curvas con la solución del sistema del ejemplo 1}
   \label{fig:SNL-01}
\end{figure}

\begin{eqnarray*}
f_1(x_1,x_2) &=& x_1^2 + x_1 x_2 -10\\
f_2(x_1,x_2) &= & x_2 + 3 x_1 x_2^2 -57
\end{eqnarray*}

El Jacobiano es
\[
J(x) = 
\left[
\begin{array}{c|c}
2 x_1 + x_2  &  x_1   \\
\hline
3 x_2^2    &   1 + 6 x_1 x_2    
\end{array}
\right]
\]
y el arreglo de funciones

\[
f(x) = 
\left[
\begin{array}{c}
x_1^2 + x_1 x_2 -10\\
x_2 + 3 x_1 x_2^2 -57
\end{array}
\right]
\]

Considerando como valores iniciales $x^{(0)} = [1.5, 3.5]^T$ tenemos:
 
\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(1)}    \\
x_2^{(1)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    1.5000 \\
    3.5000     
\end{array}
\right] -
\left[
\begin{array}{rr}
    6.5000  &  1.5000 \\
   36.7500 &  32.5000
\end{array}
\right]^{-1}
\left[
\begin{array}{l}
   -2.5000 \\
    1.6250   
\end{array}
\right] =
\left[
\begin{array}{r}
    2.0360 \\
    2.8439
\end{array}
\right]
\end{equation*}


\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(2)}    \\
 x_2^{(2)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    2.0360 \\
    2.8439     
\end{array}
\right] -
\left[
\begin{array}{rr}
    6.9159  &  2.0360 \\
   24.2629 &  35.7413\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.0644 \\
   -4.7562 
\end{array}
\right] =
\left[
\begin{array}{r}
    1.9987\\
    3.0023
\end{array}
\right]
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(3)}    \\
x_2^{(3)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    1.9987 \\
    3.0023     
\end{array}
\right] -
\left[
\begin{array}{rr}
    6.9997  &  1.9987 \\
   27.0412  & 37.0041
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.0045 \\
    0.0496
\end{array}
\right] =
\left[
\begin{array}{r}
    2.0000 \\
    3.0000
\end{array}
\right]
\end{equation*}

\textbf{Cuarta iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(4)}    \\
 x_2^{(4)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    2.0000\\
    3.0000    
\end{array}
\right] -
\left[
\begin{array}{rr}
    7.0000  &  2.0000 \\
   27.0000 &  37.0000\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.00000129\\
   -0.00002214
\end{array}
\right] =
\left[
\begin{array}{r}
    2.0000 \\
    3.0000
\end{array}
\right]
\end{equation*}

Note que solamente 4 iteraciones son necesarias para llegar a la solución $x = [2,3]^T$

La implementación en Matlab para este ejemplo son:

Para la función tenemos

\begin{verbatim}
function f = f1(x)

n = length(x);
f = zeros(n,1);

f(1) = x(1)^2 + x(1)*x(2) - 10; 
f(2) = x(2) + 3*x(1)*x(2)^2 - 57;
end
\end{verbatim}

El Jacobiano

\begin{verbatim}
function J = J1(x)
n = length(x);
J = zeros(n, n);

J(1,1) = 2*x(1) + x(2);
J(1,2) = x(1);
J(2,1) = 3*x(2)^2;
J(2,2) = 1 + 6*x(1)*x(2);
end
\end{verbatim}

y la ejecución

\begin{verbatim}
>> Newton_Raphson(@f1, @J1, [1.5,3.5]')

ans =

    2.0000
    3.0000

\end{verbatim}

\subsection{Ejemplo 2}

Resolver el siguiente sistema de ecuaciones dado por las ecuaciones y cuya solución gráfica se muestra en la figura \ref{fig:SNL-02}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/SNL-02.eps,width=8.0cm}}
   \caption{Curvas con la solución del sistema del ejemplo 2}
   \label{fig:SNL-02}
\end{figure}


\begin{eqnarray*}
f_1(x_1,x_2) &=& x_1^2 – 2x_1 – x_2 + 0.5\\
f_2(x_1,x_2) &= & x_1^2 + 4x_2^2 – 4
\end{eqnarray*}

El Jacobiano es
\[
J(x) = 
\left[
\begin{array}{c|c}
2 x_1 -2  &  -1   \\
\hline
2 x_1    &   8x_2    
\end{array}
\right]
\]
y el arreglo de funciones

\[
f(x) = 
\left[
\begin{array}{c}
x_1^2 – 2x_1 – x_2 + 0.5\\
x_1^2 + 4x_2^2 – 4
\end{array}
\right]
\]

Considerando como valores iniciales $x^{(0)} = [0, 1]^T$ tenemos:
 
\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(1)}    \\
x_2^{(1)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
0    \\
1     
\end{array}
\right] -
\left[
\begin{array}{rr}
-2  & -1     \\
 0  &  8  
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
-0.5    \\
0.0   
\end{array}
\right] =
\left[
\begin{array}{r}
-0.25    \\
1.00  
\end{array}
\right]
\end{equation*}


\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(2)}    \\
x_2^{(2)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
-0.2500    \\
1.0000      
\end{array}
\right] -
\left[
\begin{array}{rr}
2 .5000 & -1.0000     \\
-0.5000  &  8.0000  
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
0.0625    \\
0.0625  
\end{array}
\right] =
\left[
\begin{array}{r}
-0.2226   \\
 0.9939  
\end{array}
\right]
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(3)}    \\
x_2^{(3)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.2226\\
    0.9939    
\end{array}
\right] -
\left[
\begin{array}{rr}
   -2.4451  & -1.0000 \\
   -0.4451  &  7.9512\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.0008\\
    0.0009 
\end{array}
\right] =
\left[
\begin{array}{r}
   -0.2222 \\
    0.9938
\end{array}
\right]
\end{equation*}


\textbf{Cuarta iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(4)}    \\
x_2^{(4)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.2222 \\
    0.9938   
\end{array}
\right] -
\left[
\begin{array}{rr}
   -2.4444  &  -1.0000\\
   -0.4444   &  7.9505
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.00002300\\
    0.00000038
\end{array}
\right] =
\left[
\begin{array}{r}
   -0.2222 \\
    0.9938
\end{array}
\right]
\end{equation*}

La implementación en Matlab para este ejemplo son:


Para la función tenemos

\begin{verbatim}
function f = f2(x) 

n = length(x);
f = zeros(n,1);

f(1) = x(1)*x(1) - 2*x(1) - x(2) + 0.5; 
f(2) = x(1)*x(1) + 4*x(2)*x(2) - 4;
end
\end{verbatim}

El Jacobiano

\begin{verbatim}
function y = J2(x)
n = length(x);
y = zeros(n, n);
y(1,1) = 2*x(1) -2;
y(1,2) = -2;
y(2,1) = 2*x(1);
y(2,2) = 8*x(2);
end
\end{verbatim}

y la ejecución

\begin{verbatim}
>> Newton_Raphson(@f2, @J2, x)

ans =

   -0.2222
    0.9938

\end{verbatim}

\subsection{Ejemplo}

Para el circuito que se muestra en la figura \ref{fig:sistema}, esta constituido por dos mallas y dos cargas. Se sea calcular la corriente que circula por cada uno de los elementos y la potencia que debe suministrar la fuente de voltaje.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Sistemas_No_Lineales/figuras/sistema.eps,width=8.0cm}}
   \caption{Circuito eléctrico de dos mallas}
   \label{fig:sistema}
\end{figure}

La ecuaciones de Voltaje para cada una de las mallas es:

\begin{eqnarray*}
V_1 - R_1 (i_1 + i_2) - L_1/i_1 = 0\\
V_1 - R_1 (i_1 + i_2) - R_2 i_2 - L_2/i_2 =0
\end{eqnarray*}

Reorganizando las funciones tenemos

\begin{equation*}
F = \left[
\begin{array}{c}
V i_1 - R_1 (i_1 + i_2) i_1 - L_1\\
V i_2 - R_1 (i_1 + i_2) i_2- R_2 i_2^2 - L_2
\end{array}
\right] = 0
\end{equation*}

\begin{verbatim}
function f = Fun(i)
N = length(i);
f = zeros(N,1);

V = 100;
R1 = 2;
R2 = 3;
L1 = 100;
L2 = 200;

f(1) = V*i(1) - R1*(i(1)+i(2))*i(1) - L1;
f(2) = V*i(2) - R1*(i(1)+i(2))*i(2) - R2*i(2)^2 - L2;
end
\end{verbatim}

El Jacobiano para este sistema es:

\begin{equation*}
J = \left[
\begin{array}{c|c}
V  - R_1 (2 i_1 + i_2) \  & - i_1 R_1\\
\hline
-i_2 R1 & V_1  - R_1 (i_1 + 2 i_2)  - 2 R_2 i_2   \\
\end{array}
\right]
\end{equation*}

\begin{verbatim}
function J = JacSistema(i)
V = 100;
R1 = 2;
R2 = 3;
%L1 = 100;
%L2 = 200;

J(1,1) = V - R1*(2*i(1) + i(2)) ;
J(1,2) = -i(1)*R1;
J(2,1) = -i(2)*R1;
J(2,2) = V - (i(1) + 2*i(2))*R1 - 2*i(2)*R2 ;
end
\end{verbatim}

La solución considerando un valor inicial de corrientes $I^{(0)} = [0,0]^T$:

\begin{verbatim}

>> I = NewtonRaphson(@FunSistema, @JacSistema, [0,0]')
1.-      1     2

2.-     1.0720    2.3114

3.-     1.0728    2.3185

4.-     1.0728    2.3185


I =

    1.0728
    2.3185
\end{verbatim}

Potencia de la fuente $P_V =  V(i_1+i_2) = 100(1.0728 +  2.3185) = + 339.1282$\\ 
Potencia consumida por $R_1$  es $P_{R_1} = -R_1(i_1+i_2)^2 = -2(1.0728 +  2.3185)^2  = - 23.0016$\\
Potencia consumida por $R_2$  es $P_{R_1} = -R_2i_2^2 = - 3 (2.3185) ^2 = - 16.1266$\\
La potencia consumida por las cargas es $L1 + L2 = -300$.\\
La suma de la potencia es cero.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section{El método de Secantes}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Convergencia del método de Newton Raphson}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Aunque el método de Newton-Raphson en general es muy eficiente, hay situaciones en que se comporta en forma deficiente. Un caso especial, raíces múltiples.

\subsubsection{Ejemplo}

Determinar la raíz de la función $f(x) = x^{10} –1$.

La solución utilizando el método de Newton-Raphson queda:

\begin{equation*}
x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})} {f'(x^{(k+1)})}
\end{equation*}

Sustituyendo valores tenemos

\begin{equation*}
x^{(k+1)} = x^{(k)} - \frac{\left(x^{(k)} \right)^{10}-1} {10 \left(x^{(k)} \right)^9}
\end{equation*}

Y la solución numérica es:

\begin{center}
\begin{tabular}{|r|r|r|}
\hline
$x^{(k)}$ & $f(x^{(k)})$ & $f'(x^{(k)})$ \\
\hline
0.5000 & -0.9990 & 0.0195 \\
51.6500 & 135114904483914000.0000 & 26159710451871000.0000 \\
46.4850 & 47111654129711500.0000 & 10134807815362300.0000 \\
41.8365 & 16426818072478500.0000 & 3926432199748670.0000 \\
37.6529 & 5727677301318310.0000 & 1521180282851980.0000 \\
33.8876 & 1997117586819850.0000 & 589336409039672.0000\\
30.4988 & 696351844868619.0000 & 228320999775654.0000\\
27.4489 & 242802875029547.0000 & 88456233382052.8000\\
24.7040 & 84660127717097.5000 & 34269757191973.2000\\
22.2336 & 29519161271064.1000 & 13276806089225.7000\\
20.0103 &10292695105054.7000 & 5143706707446.1600\\
18.0092 & 3588840873655.1100 & 1992777367871.5700\\
16.2083 & 1251351437592.9200 & 772042782329.1500\\
14.5875 & 436319267276.5290 & 299105192259.1190\\
13.1287 & 152135121499.2910 & 115879479847.7330\\
11.8159 & 53046236848.5329 & 44894084747.9692\\
10.6343 & 18496079117.2577 & 17392888266.5936\\
9.5708 & 6449184014.3077 & 6738361277.7304\\
8.6138 & 2248691421.7628 & 2610579221.6818\\
7.7524 & 784070216.9426 & 1011391879.0870\\
\hline
\end{tabular}
\end{center}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section{Raíces múltiples y el método de Newton modificado}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

