% 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[SE No Lineales]{Solución de Ecuaciones no-lineales}

\section{Método Gráfico}

Un método simple para obtener, visualmente la raíz de la ecuación $f(x)=0$, consiste en graficar la función y observar en donde cruza el eje $x$. Este punto que representa el valor de $x$ para el cual $f(x)=0$, proporciona una aproximación inicial de la raíz.

\subsection{Ejemplo}

Utilice el método gráfico para observar algunas de las raíces de la función $f(x) = sen 10x + cos 3x$, en el intervalo $[0,5]$

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_3/figuras/ejemplo1.eps, width=10.0cm}}
   \caption{Método Gráfico}
   \label{fig:eje_1}
\end{figure}

En la Figura \ref{fig:eje_1} se muestran los cruces por cero de la función en el intervalo dado.

\subsection{Implementación}

La implementación en Matlab es

\begin{verbatim}
function Metodo_Grafico(funcion, x)
   f = funcion;
   y = f(x);
   plot(x, y);
   xlabel('x');
   ylabel('f(x)')
   grid on
\end{verbatim}

Ver Metodo\_Grafico.m. La implementación de la función 

\begin{verbatim}
function y = f1(x) 
   y  = sin (10*x) + cos(3*x);
\end{verbatim}

La corrida del método es 

\begin{verbatim}
Metodo_Grafico(@f1, 0:0.01:5)
\end{verbatim}

\section{Método de Bisección}

Para este método debemos considerar una función 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.000000 & 0.500000 & 1.000000 & -1.000000 & -0.377583 & 0.459698 \\ 
1 & 0.500000 & 0.750000 & 1.000000 & -0.377583 & 0.018311 & 0.459698 \\ 
2 & 0.500000 & 0.625000 & 0.750000 & -0.377583 & -0.185963 & 0.018311 \\ 
3 & 0.625000 & 0.687500 & 0.750000 & -0.185963 & -0.085335 & 0.018311 \\ 
4 & 0.687500 & 0.718750 & 0.750000 & -0.085335 & -0.033879 & 0.018311 \\ 
5 & 0.718750 & 0.734375 & 0.750000 & -0.033879 & -0.007875 & 0.018311 \\ 
6 & 0.734375 & 0.742188 & 0.750000 & -0.007875 & 0.005196 & 0.018311 \\ 
7 & 0.734375 & 0.738281 & 0.742188 & -0.007875 & -0.001345 & 0.005196 \\ 
8 & 0.738281 & 0.740234 & 0.742188 & -0.001345 & 0.001924 & 0.005196 \\ 
9 & 0.738281 & 0.739258 & 0.740234 & -0.001345 & 0.000289 & 0.001924 \\ 
10 & 0.738281 & 0.738770 & 0.739258 & -0.001345 & -0.000528 & 0.000289 \\ 
11 & 0.738770 & 0.739014 & 0.739258 & -0.000528 & -0.000120 & 0.000289 \\ 
12 & 0.739014 & 0.739136 & 0.739258 & -0.000120 & 0.000085 & 0.000289 \\ 
13 & 0.739014 & 0.739075 & 0.739136 & -0.000120 & -0.000017 & 0.000085 \\ 
\hline
\end{tabular}
\end{center}

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

\section{Método de la falsa posición (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=Capitulo_3/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*}
p_1 & = & [a, f(a)] \\
p_2 &=& [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.4597 \\
2 & 0.7363 & 0.7389 & 1.0000 & -0.0047 & -0.0002 & 0.4597 \\
3 & 0.7389 & 0.7391 & 1.0000 & -0.0002 & 0.0000 & 0.4597\\
4 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.4597\\
5 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.4597\\
6 & 0.7391 & 0.7391 & 1.0000 & 0.0000 & 0.0000 & 0.4597\\
\hline
\end{tabular}
\end{center}

\section{Método de iteración de punto fijo}

Los métodos de punto fijo son métodos que predicen la raíz de la ecuación a partir de una formula dada. La formula puede ser desarrollada por una sustitución sucesiva al rearreglar la ecuación $f(x)=0$ de tal modo que $x$ quede del lado izquierdo de la ecuación

\begin{equation*}
f(x) \equiv x - g(x) =0
\end{equation*} 

Lo que es igual a 

\begin{equation*}
x = g(x)
\end{equation*} 

Esta transformación se puede llevar a cabo mediante operaciones algebraicas o simplemente agregando $x$ a cada lado de la ecuación original. Por ejemplo:

\begin{eqnarray*}
x^2 - 2x + 3 =0 \\
x = \frac{x^2 +3}{2} \\
x = g(x)
\end{eqnarray*}

mientras que $sen(x) =0$ puede transformarse haciendo

\begin{equation*}
x = sen(x) +x
\end{equation*}

Las iteraciones del método dado un punto inicial $x^{(0)}$, se pueden llevar a cabo calculando $x^{(k+1)}$

\begin{equation*}
x^{(k+1)} = g(x^{(k)})
\end{equation*}

La implementación en Matlab es:

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

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

\subsection{Ejemplo}

Utilizando la iteración de punto fijo calcular las raíces $f(x) = e^{-x} - x$, para un valor inicial $x^{(0)} =0$.

De acuerdo con lo revisado la función $g(x) = e^{-x}$ y las iteraciones hasta convergencia son:

\begin{center}
\begin{tabular}{|r|r|r|}
\hline
$k$ & $x^{(k)}$ & $error(x^{(k)})$\\
\hline
1 & 1.0000 & 1.0000 \\
2 & 0.3679 & 1.7183 \\
3 & 0.6922 & 0.4685 \\
4 & 0.5005 & 0.3831 \\
5 & 0.6062 & 0.1745 \\
6 & 0.5454 & 0.1116 \\
7 & 0.5796 & 0.0590 \\
8 & 0.5601 & 0.0348 \\
9 & 0.5711 & 0.0193 \\
\vdots & \vdots & \vdots \\
19 & 0.5672 & 6.7008e-05 \\
\hline
\end{tabular}
\end{center}

Para hacer la ejecución correr:

\begin{verbatim}
Punto_Fijo(@g1, 0)
\end{verbatim}

con

\begin{verbatim}
function y = g4(x)

y = exp(-x);
\end{verbatim}

\subsection{Ejemplo}

Calcular utilizando iteración de punto fijo, un cero de la función $f(x) = sen(10x)+cos(3x)$

Para llevar a cabo la iteración de punto fijo hacemos

\begin{equation*}
x = sen(10x) + cos(3x) + x
\end{equation*}

La implementación en Matlab es

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

y = sin(10*x)+cos(3*x) +x;
\end{verbatim}

Para ejecutar hacer:

\begin{verbatim}
Punto_Fijo(@g1, 0)
\end{verbatim}

\section{Método de Newton-Raphson}

Dada una función $f(x): \mathbb{R} \rightarrow \mathbb{R}$ continua, con un valor inicial $x^{(0)}$ cerca de una raíz $x^*$ y si la derivada $f’(x)$ existe, entonces puede utilizarse esta información para desarrollar algoritmos que produzcan sucesiones $x^{(k)}$ que converjan a $x^*$ más rápidamente que los algoritmos como los vistos en la secciones anteriores.

El método considera que dada una función $f(x)$, podemos hacer una representación lineal utilizando la serie de Taylor. Así la representación lineal en Serie de Taylor es:

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

Si de la ecuación anterior despejamos el valor de $x^{(k+1)}$, suponemos que $f(x^{(k+1)}) = 0$ y despreciamos el residuo $R_1$,  tenemos que la sucesión $x^{(0)}, x^{(1)}, x^{(2)}, …, x^{(k)}, x^{(k+1)}$  puede ser calculada  utilizando la formula:

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

La implementación del Método en Matlab es:

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

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


\subsection{Ejemplo}

Calcular los ceros de la función $f(x) = x-cos(x)$ utilizando el algoritmo de Newton Raphson con  $x_0 = 0$.

\begin{center}
\begin{tabular}{|r|r|}
\hline
$k$ & $x^{(k)}$\\
\hline
0 & 0.0000 \\
1 & 1.0000 \\
2 & 0.7504 \\
3 & 0.7391 \\
4 & 0.7391 \\
5 & 0.7391 \\
6 & 0.7391 \\
7 & 0.7391 \\
\hline
\end{tabular}
\end{center}

En estas Figura \ref{fig:NR} se muestra tres iteraciones del método. En la figura de la izquierda mostramos la línea recta que es tangente a la función $f(x)$ (en negro) en $x^{(0)} = 0$, note que la línea recta cruza el eje $x$ en $x^{(1)} = 1$. En el centro la tenemos la segunda iteración tomando como valor inicial $x^{(1)}=1$ y solución $x^{(2)} = 0.7504$. Finalmente a la izquierda tenemos la linea recta y la solución para $x^{(3)} = 0.7391$. Note como poco a poco se acerca a la solución.

\begin{figure}[!htb]
\begin{center}
	\subfigure[$x^{(0)} = 0$ ] {\psfig{figure=Capitulo_3/figuras/NR_0.eps,width=4.6cm}}
	\subfigure[$x^{(1)} = 1$ ] {\psfig{figure=Capitulo_3/figuras/NR_1.eps,width=4.6cm}}	
	\subfigure[$x^{(2)} = 0.7391$] {\psfig{figure=Capitulo_3/figuras/NR_2.eps,width=4.6cm}}	
\end{center}
\caption{Iteraciones del Método de Newton-Raphson} \label{fig:NR}
\end{figure}

Para correr este ejemplo dar:

\begin{verbatim}
Newton_Raphson(@f2, @df2, 0)
\end{verbatim}

\section{Aplicaciones}

\subsection{Cálculo de la Raíz Cuadrada }

Hacer un algoritmo iterativo que permita hacer el cálculo de la raíz cuadrada de A.

Para este caso nuestra función a resolver es $f(x) = x^2-A$. La solución cuando $f(x)=0$ es $x = \sqrt{A}$. Para nuestro calculo hacemos

\begin{eqnarray*}
f(x) & = & x^2-A \\
f’(x)  & = & 2 x
\end{eqnarray*}

Los cálculos numéricos suponiendo para calcular la raíz cuadrada de 5 son:

\begin{center}
\begin{tabular}{|r|r|}
\hline
$k$ & $x^{(k)}$\\
\hline
0 & 2.0000 \\
1 & 2.2500 \\
2 & 2.2361 \\
3 & 2.2361 \\
4 & 2.2361 \\
5 & 2.2361 \\
6 & 2.2361 \\
7 & 2.2361 \\
8 & 2.2361 \\
9 & 2.2361 \\
10 & 2.2361 \\
\hline
\end{tabular}
\end{center}

Para ver la ejecución dar
\begin{verbatim}
Newton_Raphson(@f3, @df3, 2)
\end{verbatim}

Si para este ejemplo ponemos sustituimos en la formulación original los valores de la derivada y de la función tenemos:

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

Reduciendo términos tenemos

\begin{equation*} 
x^{(k+1)} =   \frac{{x^{(k)}}^2+A}{2 x^{(k)}}
\end{equation*}

La función en Matlab queda:

\begin{verbatim}
function y = Raiz_Cuadrada(A)
   r = A;
   while abs(r*r - A) > 0.00001 
      r = (r*r+A)/(2*r);
   end;
end
\end{verbatim}

\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=Capitulo_3/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=Capitulo_3/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=Capitulo_3/figuras/sol_cto_diodo.eps, width=10.0cm}}
   \caption{Solución al circuito con diodo}
   \label{fig:Sol_cto_diodo}
\end{figure}

\subsection{Solución del problema de flujos de Potencia}

Se tiene un circuito formado por tres elementos: una fuente de voltaje $V$,  una resistencia $R$ y una carga $P_L$ en una combinación en serie. Dados los valores $V = 10$ volts, $R = 2$ ohms y $P_L=12$ Wats, calcular la corriente que circula por el circuito.

La ecuación de voltaje del circuito es:

\begin{equation*}
	V - R i - V_L =0	
\end{equation*}

multiplicado por la corriente $i$ en el circuito tenemos 

\begin{equation*}
	Vi - R i^2 - P_L =0	
\end{equation*}

La solución implementada es:

\begin{verbatim}
function Simula_Flujos_Potencia()
    % Parametros del Circuito
    
    V  =  20; % Volts
    R  =  2; % ohms
    P  = 12; % Watts
    
    fprintf('La solución es %10.6f\n', NR(@fp, @der_fp, V/R, [V, R, P]));
end

function d = fp(I, p)
   V = p(1);
   R = p(2);
   PL = p(3);
   d = V*I - R*I*I - PL;
end

function der = der_fp(I, p)
   V = p(1);
   R = p(2);
   PL = p(3);
   der = V-2*R*I;
end

function r = NR(f, df, x1, p)
k=1;
while 1
    x2 = x1 - f(x1, p)/df(x1, p);
    if abs(f(x2, p)) < 1e-6 break;
    else x1 = x2;
    end;
    k = k+1;
end;
r = x2;
end
\end{verbatim}

La solución obtenida de la corrida es:

\begin{verbatim}
>> Simula_Flujos_Potencia
La solución es   3.000000

\end{verbatim}


\section{Evaluación de Polinomios}

Consideremos como ejemplo el siguiente polinomio.

\begin{equation*}
f_3(x) = a_0 + a_1x + a_2 x^2 + a_3 x^3
\end{equation*}


En este caso los coeficientes del polinomio los podemos poner en un vector de datos dado como $a = [a_0, a_1, a_2, a_3]$, y realizar la evaluación polinomial utilizando ciclos. Adicionalmente necesitamos de una función que nos permita calcular la potencias de $x$.

Manipulando el polinomio, podemos dar una forma eficiente de realizar la evaluación. Así el mismo polinomio, lo podemos representar como:

\begin{equation*}
f_3(x) = a_0 + x(a_1+ x(a_2 + x (0 x + a_3)))
\end{equation*}
 
De manera general podemos expresar como

\begin{equation*}
f_N(x) = a_0 + x(a_1 + x(a_2  +\cdots +  x(a_{N-1} + x( 0 x + a_N))))
\end{equation*}
 
Note que ahora ya no tenemos potencia de x y su implementación resulta mas eficiente. La implementación en Matlab es:

\begin{verbatim}
function y = Evalua_Polinomio(x, a) 

n = length(a);
p = 0.0;

for i=n:-1: 1
    p = p.*x + a(i);
end
y = p;
\end{verbatim}

\subsection{Ejemplo}

Consideremos el siguiente polinomio $f_2(x) = 1+2x+3x^2$, el cual deseamos evaluar en $x = 20$. Para implementarlo hacemos


\begin{center}
\begin{tabular}{|c|c|c|}
\hline
$i$ & $p = p*x + a_i$ & $p=0$\\
\hline
$2$ & $p = 0*20+3$ & $p = 3$ \\
$1$ & $p = 3*20+2$ & $p = 62$\\
$0$ & $p = 62*20+1$ & $p = 1241$\\
\hline
\end{tabular}
\end{center}

\subsection{Ejemplo}

Hacer la evaluación del polinomio $f_5(x) =  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25$ en $x=2.5$

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
$i$ & $p = p*x + a_i$ & $p=0$\\
\hline
5 & $p = 0.000000 * 2.500000  + 1.000000$ & 1.000000 \\ 
4 & $p = 1.000000 * 2.500000  - 3.500000$ & -1.000000 \\ 
3 & $p = -1.000000 * 2.500000  + 2.750000$ & 0.250000 \\ 
2 & $p = 0.250000 * 2.500000  + 2.125000$ & 2.750000 \\ 
1 & $p = 2.750000 * 2.500000  - 3.875000$ & 3.000000 \\ 
0 & $p = 3.000000 * 2.500000  + 1.250000$ & 8.750000 \\ 
\hline
\end{tabular}
\end{center}

\subsection{Cálculo de derivadas}

Cuando se buscan los ceros de una función, como es el caso del método de Newton, es necesario no solo hacer la evaluación del polinomio sino calcular también su derivada. En este caso la derivada de cualquier polinomio la podemos calcular como:

 
\begin{equation*}
f_n(x) =  a_0  +  a_1x + a_2x^2 + a_3x^3  +\cdots  + a_{n-1}x^{n-1} +  a_nx^n 
\end{equation*}
 
\begin{equation*}
f'_n(x) =  a_1 + 2a_2x + 3 a_3x^2 +  \cdots  +  (n-1)a_{n-1}x^{n-2} +  na_nx^{n-1} 
\end{equation*}
 

Note que el polinomio $f’n(x)$ es un polinomio de menor grado.

Los coeficientes $a$  del polinomio original pueden tomarse para calcular los coeficiente del nuevo polinomio, para esto tomamos en cuenta lo siguiente

\begin{equation*}
\begin{array}{c}
a = [a_0, a_1, a_2, a_3, a_4, \cdots, a_{n-1}, a_n] \\
\\
a' = [ a_1, 2a_2, 3  a_3, 4 a_4, \cdots, (n-1)  a_{n-1}, n a_n]
\end{array}
\end{equation*}

Con esto tenemos que dado $a$ los coeficientes de del polinomio de derivadas pueden ser calculados y la evaluación de la derivada se reduce a una simple evaluación polinomial.

Así por ejemplo considerando el polinomio $f_3(x) = 1+2x+3x^2+ 4x^3$,  tenemos que $a = [1, 2, 3, 4]$  y los coeficiente de $f'_3(x)$ son $a' = [2, 6,12]$. La evaluación del polinomio resultante lo podemos hacer para $x = 20$ utilizando:

\begin{center}
\begin{tabular}{|c|r|r|}
\hline
$i$ & $p = p*x + a'_i$ & $p$\\
\hline
2 &     0*20 +12 & 12\\
1 &   12*20 + 6  &  246\\
0 & 246*20 + 2 & 4922 \\
\hline
\end{tabular}
\end{center}

\subsection{Ejemplo}

Dado el polinomio de la serie de Taylor correspondiente a la  función $f(x) = seno(x)$, hacer la evaluación de la derivada en $x=0.5$.

El polinomio de la serie de Taylor para seno es:

\begin{equation*}
f(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} + \cdots
\end{equation*}

Los coeficientes de este polinomio son $a = [0, 1,  0, -1/3!, 0, 1/5!, 0, -1/7!, 0, 1/9!]$, por lo tanto los coeficientes para la derivada son $a^* = [1,  0, -3/3!, 0, 5/5!, 0, -7/7!, 0, 9/9!]$

\begin{center}
\begin{tabular}{|c|r|r|}
\hline
$i$ & $p = p*x + a_i^*$ & $p$ \\
\hline
8 & $p = 0.000000 * 0.500000  + 0.000025$ & 0.000025 \\ 
7 & $p = 0.000025 * 0.500000  + 0.000000$ & 0.000012 \\ 
6 & $p = 0.000012 * 0.500000  - 0.001389$ & -0.001383 \\ 
5 & $p = -0.001383 * 0.500000  + 0.000000$ & -0.000691 \\ 
4 & $p = -0.000691 * 0.500000  + 0.041667$ & 0.041321 \\ 
3 & $p = 0.041321 * 0.500000  + 0.000000$ & 0.020660 \\ 
2 & $p = 0.020660 * 0.500000  - 0.500000$ & -0.489670 \\ 
1 & $p = -0.489670 * 0.500000  + 0.000000$ & -0.244835 \\ 
0 & $p = -0.244835 * 0.500000  + 1.000000$ & 0.877583 \\ 
\hline
\end{tabular}
\end{center}


La implementación en Matlab es
\begin{verbatim}
function y = Evalua_Derivada(x, a)

    N = length(a);
    b = a(2:N).*[1:N-1];
    
    y = Evalua_Polinomio(x, b);

end
\end{verbatim}

\subsection{Calculo de la integral}

Dado un polinomio de grado $f_n(x)$

\begin{equation*}
f_n(x) =  a_0  +  a_1x + a_2x^2 + a_3x^3 \cdots  + a_{n-1}x^{n-1} +  a_nx_n 
\end{equation*}

La integral de $f_n(x)$ se calcula de manera analítica como

\begin{equation*}
\int f(x) dx = CTE +  a_0 x +  \frac{a_1}{2}x^2 + \frac{a_2}{3}x^3 + \frac{a_3}{4}x^4 \cdots  +  \frac{a_{n-1}}{n}x^n +  \frac{a_n}{n+1}x^{n+1} 
\end{equation*}
 
Note que los coeficiente del nuevo polinomio $a^*$ pueden calcularse a partir de los coeficientes $a$ de polinomio original :

 \begin{equation*}
\begin{array}{c}
a = [a_0, a_1, a_2, a_3, a_4, \cdots, a_{n-1}, a_n] \\
\\
a^* = [ CTE, \frac{a_0}{1}, \frac{a_1}{2}, \frac{a_2}{3}, \frac{ a_3}{4}, \frac{a_4}{5}, \cdots, \frac{  a_{n-1}}{n}, \frac{a_n}{n+1}]
\end{array}
\end{equation*}

Así la evaluación de la integral del polinomio se reduce a calcular los nuevos coeficientes del nuevo polinomio y hacer la evaluación del nuevo polinomio con coeficientes $a^*$.

En el caso de tener una integral definida debemos notar:

\begin{equation*}
\int_{x_0}^{x_1} f(x) dx = \left . \left (CTE + a_0 x +  \frac{a_1}{2}x^2 + \frac{a_2}{3}x^3 + \frac{a_3}{4}x^4 \cdots  +  \frac{a_{n-1}}{n}x^n +  \frac{a_n}{n+1}x^{n+1} \right)  \right |_{x_0}^{x_1}
\end{equation*}

Lo que significa que debemos evaluar el polinomio con coeficientes $a^*$ en $x_1$ y restar la evaluación en $x_2$. 

La implementación en Matlab para este queda

\begin{verbatim}
function y = Evalua_Integral(x1, x2, a) 
    N = length(a);
    b = [0, a./[1:N]];
    y = Evalua_Polinomio(x2, b) - Evalua_Polinomio(x1, b);
end
\end{verbatim}

\subsection{Ejemplo}

Para la serie de Taylor de la función seno hacer la evaluación de la integral:

\begin{equation*}
\int_0^{\pi/2} sin(x) dx
\end{equation*}

El polinomio de la serie de Taylor para seno es:

\begin{equation*}
f(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} + \cdots
\end{equation*}

Los coeficientes de este polinomio son $a = [0, 1,  0, -1/3!, 0, 1/5!, 0, -1/7!, 0, 1/9!]$, por lo tanto los coeficientes para la integral serán :

\begin{equation*}
\begin{array}{c}
a^* = [0, 0, \frac{1}{2},  0, \frac{-1}{4 \times 3!}, 0, \frac{1}{6\times 5!}, 0, \frac{-1}{8 \times 7!}, 0, \frac{1}{10 \times 9!}] \\
\\
a^* = [0,   0,   0.5,  0,   -0.041667,   0,   0.001389,   0,   -0.000025,   0]
\end{array}
\end{equation*}

La evaluación del polinomio es $x=\pi/2$ con los coeficientes $a^*$ se muestra a continuación y la evaluación con $x=0$ .

\begin{center}
\begin{tabular}{|c|r|r|}
\hline
$i$ & $p = p*x + a_i^*$ & $p$ \\
\hline
10 & p = 0.000000 * 1.570796  + 0.000000 & 0.000000 \\ 
9 & p = 0.000000 * 1.570796  + 0.000000 & 0.000000 \\ 
8 & p = 0.000000 * 1.570796  - 0.000025 & -0.000024 \\ 
7 & p = -0.000024 * 1.570796  + 0.000000 & -0.000038 \\ 
6 & p = -0.000038 * 1.570796  + 0.001389 & 0.001329 \\ 
5 & p = 0.001329 * 1.570796  + 0.000000 & 0.002088 \\ 
4 & p = 0.002088 * 1.570796  - 0.041667 & -0.038387 \\ 
3 & p = -0.038387 * 1.570796  + 0.000000 & -0.060297 \\ 
2 & p = -0.060297 * 1.570796  + 0.500000 & 0.405285 \\ 
1 & p = 0.405285 * 1.570796  + 0.000000 & 0.636620 \\ 
0 & p = 0.636620 * 1.570796  + 0.000000 & \bf{1.000000} \\ 
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|r|r|}
\hline
$i$ & $p = p*x + a_i^*$ & $p$ \\
\hline
10 & p = 0.000000 * 0.000000  + 0.000000 & 0.000000 \\ 
9 & p = 0.000000 * 0.000000  + 0.000000 & 0.000000 \\ 
8 & p = 0.000000 * 0.000000  - 0.000025 & -0.000025 \\ 
7 & p = -0.000025 * 0.000000  + 0.000000 & 0.000000 \\ 
6 & p = 0.000000 * 0.000000  + 0.001389 & 0.001389 \\ 
5 & p = 0.001389 * 0.000000  + 0.000000 & 0.000000 \\ 
4 & p = 0.000000 * 0.000000  - 0.041667 & -0.041667 \\ 
3 & p = -0.041667 * 0.000000  + 0.000000 & 0.000000 \\ 
2 & p = 0.000000 * 0.000000  + 0.500000 & 0.500000 \\ 
1 & p = 0.500000 * 0.000000  + 0.000000 & 0.000000 \\ 
0 & p = 0.000000 * 0.000000  + 0.000000 & \bf{0.000000} \\ 
\hline
\end{tabular}
\end{center}

Con esto tenemos que el valor de la integral es igual a $1-0 = 1$.

La solución en Matlab queda

\begin{verbatim}
>> a = [0, 1, 0, -1/factorial(3), 0, 1/factorial(5), 
0, -1/factorial(7), 0, 1/factorial(9)]

>> Evalua_Integral(0, pi/2, a)

ans =

    1.0000
\end{verbatim}

\subsection{Ejemplo cálculo de raíces de un polinomio}

 Dado un polinomio de grado $n$
 
\begin{equation*}
f_n(x) = a_0 + a_1x+ a_2x^2 + ... + a_nx^n.
\end{equation*}

es posible realizar una factorización de la siguiente forma

\begin{equation*}
f_n(x) =(x+s)*(a'_0 + a'_1x+ a'_2x^2 + ... + a'_{n-1}x^{n-1}).
\end{equation*}

y podemos comprobar que tenemos un monomio de la forma $x-s$, si evaluando la función $f_n(x)$  en $x=-s$, el resultado es $f_n(-s)=0$. Una manera de resolverlo es utilizar un método de prueba que sería la manera mas tediosa.

El problema lo podemos solucionar utilizando el método de Newton-Raphson, sin embargo hay que hacer un pequeño cambio en la formulación dado que el argumento de la función es negativo. A partir de la serie de 
Taylor tenemos:

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

haciendo  $x^{(k)} = -s^{(k)}$ tenemos

\begin{equation*}
f(-s^{(k+1)}) \approx f(-s^{(k)}) + f'(-s^{(k)})(-s^{(k+1)} +s^{(k)}) =0
\end{equation*}

Despejando para $s^{(k+1)}$ tenemos la iteración de Newton-Raphson para este problema.

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

Para hacer la solución más eficiente evaluaremos la función utilizando evaluación polinomial tanto para la función como para la derivada.

\begin{verbatim}
function sn = factor_NR(s0, a)
    % factor_NR(s0, a) Calcula el factor (x+s) de un polinomio con 
    % coeficiente a dado un valor inicial s0 mediante el 
    % metodo de Newton-Raphson


    for iter = 1: 100000
        sn = s0 + Evalua_Polinomio(-s0, a)/Evalua_Derivada(-s0, a);
    
        if abs(Evalua_Polinomio(sn, a)) < 1e-6
        break;
        else
            s0 = sn
        end
    end;
end
\end{verbatim}

Como ejemplo numérico calculamos la raíz del polinomio $f_6(x) = x^6-1$ con un valor inicial $s^{(0)} = 10$

\begin{equation*}
\begin{array}{l}
s_{1} = 10.000000 + \frac{999999.000000}{-600000.000000} = 8.333335 \\ 
\\
s_{2} = 8.333335 + \frac{334897.378558}{-241126.784337} = 6.944450 \\ 
\\
s_{3} = 6.944450 + \frac{112156.191251}{-96903.735989} = 5.787052 \\ 
\\
s_{4} = 5.787052 + \frac{37560.618299}{-38943.785362} = 4.822569 \\ 
\\
s_{5} = 4.822569 + \frac{12578.711852}{-15651.050576} = 4.018871 \\ 
\\
s_{6} = 4.018871 + \frac{4212.321940}{-6290.306217} = 3.349218 \\ 
\\
s_{7} = 3.349218 + \frac{1410.434918}{-2528.533032} = 2.791411 \\ 
\\
s_{8} = 2.791411 + \frac{472.088718}{-1016.880869} = 2.327159 \\ 
\\
s_{9} = 2.327159 + \frac{157.838757}{-409.526161} = 1.941741 \\ 
\\
s_{10} = 1.941741 + \frac{52.597923}{-165.618133} = 1.624156 \\ 
\\
s_{11} = 1.624156 + \frac{17.355481}{-67.809321} = 1.368210 \\ 
\\
s_{12} = 1.368210 + \frac{5.560199}{-28.768380} = 1.174936 \\ 
\\
s_{13} = 1.174936 + \frac{1.630779}{-13.434500} = 1.053548 \\ 
\\
s_{14} = 1.053548 + \frac{0.367497}{-7.787952} = 1.006360 \\ 
\\
s_{15} = 1.006360 + \frac{0.038774}{-6.193251} = 1.000100 \\ 
\\
s_{16} = 1.000100 + \frac{0.000598}{-6.002990} = 1.000000 \\
\end{array}
\end{equation*}

Como resultado tenemos que $x^6 -1 = (x+1) (a_5 x^5 + a_4 x^4 + a_3 x^3 + a_2 x^2 + a_1 x + a_0)$.

\section{Deflación de polinomios y División de Polinomios}

La deflación polomial se define como el cálculo de un polinomio de grado $f_{n_1}(x)$ tal que $f_n(x) = (x-s)*f_{n-1}(x)$.  Para esto suponga que conoce, una de las raíces de un polinomio y que podemos realizar la división de este polinomio para obtener un polinomio de grado menor. Así por ejemplo si tenemos

\begin{equation*}
f_n(x) = a_0 + a_1x+ a_2x^2 + ... + a_nx^n.
\end{equation*}
 
Y una de sus raíces es $s$, entonces podemos escribir

\begin{equation*}
f_n(x) =(x+s)*(a'_0 + a'_1x+ a'_2x^2 + ... + a'_{n-1}x^{n-1}).
\end{equation*}
 
En este caso el residuo de la división es cero y podemos calcular un polinomio de grado $n-1$. Para un polinomio de orden 3 tenemos que

\begin{center}
\begin{tabular}{r|rrrr}
           & $a_3x^2$ & $(a_2-a_3s) x$      & $a_1-(a_2-a_3s)s$ & \\
\hline
$x+s$ & $a_3x^3$ & $+a_2x^2$             & $+ a_1x$                 & $+a_0$\\
           &$-a_3x^3$ & $- a_3 s x^2$         &                                &        \\	
           & 0               &$(a_2-a_3s) x^2$   &  $+ a_1x$ & $+a_0$\\

           &                   & $-(a_2-a_3s) x^2$ & $-(a_2-a_3s)sx$ & \\

            &                  &  $0$          & $(a_1-(a_2-a_3s)s)x$& $+a_0$  \\
            &                  &                  & $-(a_1-(a_2-a_3s)s)x$ & $-(a_1-(a_2-a_3s)s)s$ \\
            &                  &                  & 0          &  $a_0-(a_1-(a_2-a_3s)s)s$\\
\end{tabular}
\end{center}

La división sintética, es una manera de hacer lo mismo, pero de forma compacta. Así tenemos para el ejemplo anterior que podemos hacer:


\begin{center}
\begin{tabular}{ccccc}
$x+s$ & $a_3$  & $a_2$             & $a_1$                        & $a_0$ \\
           &             &$- a_3 s$        & $-(a_2- a_3 s)s $      & $-(a_1-(a_2-s a_3)s)s$\\
\hline
            & $a_3$ & $a_2- a_3 s$  & $a_1-(a_2- a_3 s)s$ & $a_0-(a_1-(a_2- a_3 s)s)s$ \\
\end{tabular}
\end{center}

y de manera general suponiendo que el grado del polinomio es $n$:

\begin{center}
\begin{tabular}{ccccc}
$x+s$ & $a_n$  & $a_{n-1}$             & $a_{n-2}$                        & $a_{n-3}$ \\
           &             & $- a_{n} s$        & $-(a_{n-1}- a_n s)s $      & $-(a_{n-2}-(a_{n-1}-s a_{n})s)s$\\
\hline
            & $a_n$ & $a_{n-1}- a_{n} s$  & $a_{n-2}-(a_{n-1}- a_n s)s$ & $a_{n-3}-(a_{n-2}-(a_{n-1}- a_n s)s)s$ \\
\end{tabular}
\end{center}

Renombrado términos tenemos

\begin{eqnarray*}
b_n & = & a_n\\
b_{n-1} &=& a_{n-1}- a_{n} s = a_{n-1}- b_{n} s\\
b_{n-2} &=& a_{n-2}-(a_{n-1}- a_n s)s = a_{n-2}-b_{n-1}s\\
\vdots & = \vdots\\
b_2 &=&  a_2 - b_3 s\\
b_1 & = &  a_1 - b_2s\\
b_0 & = & a_0 - b_1 s
\end{eqnarray*}

En general 

\begin{eqnarray*}
b_n & = & a_n\\
b_{k} &=& a_{k}- b_{k+1} s \\
\\
\forall k &\in& [n-1, n-2, \cdots, 2, 1, 0]
\end{eqnarray*}


La implementación en Matlab es:

\begin{verbatim}
function y = Division_Sintetica (a, s)

n = length(a);
b = zeros(1,n);
 
b(n) = a(n);
 

for k= n-1: -1: 1
    b(k) = a(k) - b(k+1) * s;
end;

disp('residuo ');
disp(b(1));

y = b(2:n);
\end{verbatim}

Y para realizar la ejecución hacemos

\begin{verbatim}
Division_Sintetica([-120, -46, 79, -3, -7, 1], -4)
\end{verbatim}

\subsection{Ejemplo}

Dado el polinomio $f_5(x) = x^5 - 7x^4 – 3x^3 + 79x^2 – 46 x –120$ encontrar la división sintética con el monomio $(x-4)$.

\begin{center}
\begin{tabular}{ccccccc}
$x-4$ & 1 & -7 & -3   &  79 &-46 &-120 \\
          &    &  4 & -12 & -60 & 76 & 120 \\
\hline          
          & 1 & -3 & -15 & 19 & 30 & 0\\
\end{tabular}
\end{center}

Para este ejemplo dar en Matlab

\begin{verbatim}
Division_Sintetica([-120, -46, 79, -3, -7, 1], -4)

El resultado es

30    19   -15    -3     1
\end{verbatim}

El polinomio resultante, en este caso, es $f_4(x) = x^4 - 3x^3 – 15x^2 + 19x +30$. Note que el residuo es cero y en este caso tenemos una deflación polinomial.

En algunos casos en deseable hacer la división sintética de la función $f_n(x)$ entre un polinomio de segundo orden de la forma $f_2(x) = x^2 + rx + s$. Aplicando la técnica anterior podemos realizar la división $f_n(x)/f_2(x)$ como


\begin{center}
\begin{tabular}{c|cccccccc}
                         & $a_5 x^3$        & $+(a_4 - a_5 r) x^2$                 &     $+((a_3 - s a_5) - (a_4 -a_5 r)r)x$             &                  &              & \\
\hline                    
$x^2 + rx + s$  & $a_5 x^5$ & $+ a_4 x^4$ & $+  a_3 x^3$ & $+ a_2 x^2$ & \dots \\
                        & $-a_5 x^5$ & $-a_5 r x^4$ & $-a_5 s x^3$ \\ 
\cline{2-5}
                        &                   & $(a_4 -a_5 r)x ^4$& $+  (a_3 -s a_5)x^3$ & $+ a_2 x^2$ &\dots \\
                        &                   & $-(a_4 - a_5 r)x ^4$& $- (a_4 - a_5 r)rx^3$ & $-(a_4 - a_5 r)sx^2$ & \dots\\
\cline{3-5}
                        &                    &                                & $((a_3 -s a_5) - (a_4 - a_5 r)r)x^3$ & +$(a_2 - (a_4 - a_5r)s)x^2$ \\
           & & & \vdots & \vdots & \dots
\end{tabular}
\end{center}

Si en la ecuación anterior sustituimos 

\begin{eqnarray*}
b_5&=& a_5\\
b_4 = a_4 - a_5 r &=& a_4 - r b_5 \\
b_3 = a_3 -s a_5 - (a_4 - a_5 r)r &=& a_3 -s b_5 - rb_4 \\
b_2 &= & a_2 - s b_4 - r  b_3\\
b_1 &= & a_1 - s b_3 - r  b_2\\
b_0 &= & a_0 - s b_2 - r  b_1
\end{eqnarray*}


En general para un polinomio de grado $f_n(x)$ tenemos que podemos hacer la división sintética aplicando la siguiente recurrencia

\begin{eqnarray*}
b_n &=& a_n\\
b_{n-1} &=& a_{n-1} - r b_n  \\
b_{k} &=& a_k - rb_{k+1} - sb_{k+2} \\
\forall k &\in& \{n-2, n-3, \cdots 1\}
\end{eqnarray*}

Finalmente el residuo de la división se calcula como

\begin{equation*}
R = b_1(x+r) + b_0
\end{equation*}

La implementación en Matlab es:

\begin{verbatim}
function [c, r] = Division_Sintetica2(a, r, s)

n = length(a);

b=zeros(n,1);
 
b(n) = a(n);
b(n-1) = a(n-1) - r*b(n);

for k=n-2: -1:1
      b(k) = a(k) - r*b(k+1) - s*b(k+2);
end;

fprintf(1,'Residuo %6.4f x + %6.4f\n', b(2), b(1)+r*b(2)); 

c = b(3:n);
r = [b(2); b(1)+r*b(2)];
end
\end{verbatim}

\subsection{Ejemplo}

Hacer la división sintética del polinomio $f_3(x) =  x^3 - 1$, entre el polinomio $f_2(x) = x^2+3x+4$

De acuerdo con la sucesión 
\begin{eqnarray*}
b_3 &=& a_3 = 1.0\\
b_2 &=& a_2 - r b_3 = 0 - (1)(3) = -3\\
b_1 &=& a_1 - r b_2 - s b_3 = 0 - (3)  (-3) -( 4)(1) = 5\\
b_0 &=& a_0 - r b_1 - s b_2 = -1 - (3)(5) - (4) (-3) = -4
\end{eqnarray*}

El residuo de acuerdo con lo expuesto es
\[
R = b_1(x+r) + b_0 = 5(x+3)-4 = 5x + 11 
\]

La corrida correspondiente en Matlab queda

\begin{verbatim}
>> [b, r] = Division_Sintetica2([-1, 0, 0, 1], 3, 4)
Residuo 5.0000 x + 11.0000

b =

    -3
     1


r =

    11
     5
\end{verbatim}

\subsection{Ejemplo}

Hacer la división sintética del polinomio $f_5(x) =  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25$, entre el polinomio $f_2(x) = x^2+x-2$

De acuerdo con la sucesión 

\begin{eqnarray*}
b_5 &=& a_5 = 1.0 \\
b_4 &=& a_4 - r b_5 = -3.5 - (1) (1) = -4.5\\
b_3 &= &a_3 - r b_4 - s b_5 = 2.75 - (1) (-4.5) - (-2)(1) = 2.75 + 6.5 = 9.25\\
b_2 &=& a_2 - r b_3 - s b_4 = 2.125 - (1) (9.25) - (-2)(-4.5) = -16.125\\
b_1 &=& a_1 - r b_2 - s b_3 = - 3.875 - (1) (-16.125) - (-2)(9.25)  = 30.75\\
b_0 &=& a_0 - r b_1 - s b_2 = 1.25 - (1) (30.75) - (-2)(-16.125) = -61.75
\end{eqnarray*}

Nuestro residuo queda

\[
R = b_1(x+r) +b_0 = 30.75\times(x+1) +(- 61.75) = 30.75x - 31.00
\]

La siguiente tabla muestra el resumen de los coeficientes calculados. Note que el residuo no es cero.

\begin{center}
\begin{tabular}{c|cccccc}
                  &  1  &  -3.5  & 2.75 & 2.125 & -3.875 & 1.25 \\
$x^2+x-2$ &     &  -1   & 6.5 & -18.25 & 34.625 & -63.00\\
\hline
                  &       1 &   -4.5 &  9.25 &  -16.125 &  30.75 & -61.75\\
\end{tabular}
\end{center}



La ejecución en Matlab sería

\begin{verbatim}
[b, r] = Division_Sintetica2([1.25, -3.875, 2.125, 2.75, -3.5, 1], 1, -2)
Residuo 30.7500 x - 31.0000

b =

  -16.1250
    9.2500
   -4.5000
    1.0000


r =

  -31.0000
   30.7500
\end{verbatim}

Lo cual significa que 
\[
 \frac{x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25}{x^2+x-2} = x^3 -4.5x^2 +9.25x -16.125
\]


\subsection{Ejemplo}

Repetir el problema anterior pero suponiendo que $f_2(x) = x^2 - x - 2$

Pare este ejemplo tenemos la sucesión 

\begin{eqnarray*}
b_5 &=& a_5 = 1.0 \\
b_4 &=& a_4 - r b_5 = -3.5 - (-1) (1) = -2.5\\
b_3 &= &a_3 - r b_4 - s b_5 = 2.75 - (-1) (-2.5) - (-2)(1) =  2.25\\
b_2 &=& a_2 - r b_3 - s b_4 = 2.125 - (-1) (2.25) - (-2)(-2.5) = -0.625\\
b_1 &=& a_1 - r b_2 - s b_3 = - 3.875 - (-1) (-0.625) - (-2)(2.25)  = 0\\
b_0 &=& a_0 - r b_1 - s b_2 = 1.25 - (-1) (0) - (-2)(-0.625) = 0
\end{eqnarray*}
\begin{center}
\begin{tabular}{c|cccccc}
                  &  1  &  -3.5  & 2.75 & 2.125 & -3.875 & 1.25 \\
$x^2-x-2$ &     &   1   & -0.5 & -2.75& 3.875 & -1.25\\
\hline
                  &           1  &  -2.5    & 2.25   & -0.625   &      0 &  0\\
\end{tabular}
\end{center}

con un residuo

\[
R =  0 \times(x-1) +0 = 0
\]
La corrida en Matlab es

\begin{verbatim}
[b, r] = Division_Sintetica2([1.25, -3.875, 2.125, 2.75, -3.5, 1], -1, -2)
Residuo 0.0000 x + 0.0000

b =

   -0.6250
    2.2500
   -2.5000
    1.0000


r =

     0
     0
\end{verbatim}

Note que los residuos son iguales a cero, por lo tanto
\[
 \frac{x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25}{x^2-x-2} = x^3 -2.5x^2 +2.25x -0.625
\]

\section{Solución de un polinomio de orden 2 en la forma general}

Un polinomio de segundo orden lo podemos expresar como

\begin{equation*}
f_2(n) = a_0 + a_1 x + a_2 x^2
\end{equation*}

La solución utilizando la formula general tenemos que

\begin{equation*}
x=\frac{-a_1 \pm  \sqrt{a_1^2 - 4a_2 a_0} } {2 a_2}
\end{equation*}

Para su implementación en Matlab hacemos

\begin{verbatim}
function [x1, x2] = Formula_General(a)
x1 = (-a(2)+ sqrt(a(2)^2 - 4*a(3)*a(1)))/(2*a(3));
x2 = (-a(2)- sqrt(a(2)^2 - 4*a(3)*a(1)))/(2*a(3));
\end{verbatim}

\section{Método de Bairstow para la solución de polinomios }

El método de  Bairstow es un método iterativo, basado en el método de Newton Raphson. Dado un polinonio $f_n(x)$ se encuentran dos factores, un polinomio cuadrático $f_2(x) = x^2 + rx + s$ y un polinomio $f_{n-2}(x)$. El procedimiento general para el método de Bairstow es:

 
\begin{enumerate}
\item Dado $f_n(x)$ y $r_0$ y $s_0$
\item Utilizando el método de NR calculamos $f_2(x) = x^2 + r_0x + s_0$ y $f_{n-2}(x)$, tal que, el residuo de $fn(x)/ f_2(x)$ sea igual a cero.
\item Se determinan la raíces $f_2(x)$, utilizando la formula general.
\item Se calcula  $f_{n-2}(x)= f_n(x)/ f_2(x)$
\item Hacemos $f_n(x)= f_{n-2}(x)$
\item Si el grado del polinomio es mayor que tres regresamos al paso 2
\item Si no terminamos
\end{enumerate}

La principal diferencia de este método, respecto a otros, es que permite calcular todas las raíces de un polinomio (reales e imaginarias).

Para calcular la división de polinomios, hacemos uso de la división sintética. Así dado

\begin{equation*}
f_n(x) = a_nx^n + a_{n-1}x^{n-1} + … + a_2x^2 + a_1x + a_0
\end{equation*}
 
Al dividir entre $f_2(x) = x^2 – rx – s$, tenemos como resultado el siguiente polinomio

\begin{equation*}
f_{n-2}(x) = b_nx^{n-2} + b_{n-1}x^{n-3} + … + b_3x + b_2
\end{equation*}
 
con un residuo $ R = b_1(x+r) + b_0$, el residuo será cero solo si $b_1$ y $b_0$ lo son.

Los términos $b$, los calculamos utilizamos división sintética, la cual puede resolverse utilizando la siguiente relación de recurrencia

\begin{eqnarray*}
b_n  &=& a_n\\
b_{n-1} &=& a_{n-1} - rb_n\\
b_i &=& a_i - rb_{i+1} - sb_{i+2}
\end{eqnarray*}

Una manera de determinar los valores de $r$ y $s$ que hacen cero el residuo es utilizar el Método de Newton-Raphson. Para ello necesitamos una aproximación lineal de $b_1$ y $b_0$ respecto a $r$ y $s$ la cual calculamos utilizando la serie de Taylor

 \begin{eqnarray*}
b_1(r_{k+1}, s_{k+1}) &=& b_1(r_k, s_k) + \frac{\partial b_1(r_k, s_k)}{\partial r} (r_k -r_{k+1} ) + \frac{\partial b_1(r_k,s_k)}{\partial s} (s_k - s_{k+1})\\
b_0(r_{k+1}, s_{k+1}) &=& b_0(r_k, s_k) + \frac{\partial b_0(r_k, s_k)}{\partial r} (r_k -r_{k+1} ) + \frac{\partial b_0(r_k, s_k)}{\partial s} (s_k - s_{k+1} )\\
\end{eqnarray*}

donde los valores de $r_k$ y $s_k$ están dados y calculamos los nuevos valores $r_{k+1}$ y $s_{k+1}$ que hacen a $b_1(r_{k+1}, s_{k+1})$ y $b_0(r_{k+1}, s_{k+1})$ igual a cero. El sistema de ecuaciones que  tenemos que resolver es:

 \begin{eqnarray*}
0 &=& b_1(r_k, s_k) + \frac{\partial b_1(r_k, s_k)}{\partial r} (r_k -r_{k+1} ) + \frac{\partial b_1(r_k,s_k)}{\partial s} (s_k -s_{k+1})\\
0 &=& b_0(r_k, s_k) + \frac{\partial b_0(r_k, s_k)}{\partial r} (r_k -r_{k+1} ) + \frac{\partial b_0(r_k, s_k)}{\partial s} (s_k -s_{k+1})\\
\end{eqnarray*}

En forma matricial tenemos

\[
\left[
\begin{array}{cc}
\frac{\partial b_1(r_k, s_k)}{\partial r}  & \frac{\partial b_1(r_k,s_k)}{\partial s}   \\
\frac{\partial b_0(r_k, s_k)}{\partial r} & \frac{\partial b_0(r_k, s_k)}{\partial s}\\
\end{array}
\right]
\left[
\begin{array}{c}
dr_k  \\
ds_k\\
\end{array}
\right] =
\left[
\begin{array}{c}
-b_1(r_k, s_k)  \\
-b_0(r_k, s_k)\\
\end{array}
\right]
\]

donde $dr_k = r_k - r_{k+1} $ y $ds_k = s_k -s_{k+1}  $

Bairtow muestra que los coeficientes $c$ del polinomio de derivadas parciales $f'_n(x) = c_0 + c_1 x + c_2 x^2+ \cdots$ pueden obtenerse haciendo un procedimiento similar a la división sintética, así
 
\begin{eqnarray*}
c_n &=& b_n \\
c_{n-1} &=& b_{n-1} - rc_n\\
c_i &=& b_i - rc_{i+1} - sc_{i+2}
\end{eqnarray*}

donde
\begin{eqnarray*}
\frac{\partial b_0(r_k, s_k)}{\partial r} &=& c_1(r_k, s_k)\\
\frac{\partial b_1(r_k, s_k)}{\partial r} &=& c_2(r_k, s_k)\\
\frac{\partial b_0(r_k, s_k)}{\partial s} &=& c_2(r_k, s_k)\\
\frac{\partial b_1(r_k, s_k)}{\partial s} &=& c_3(r_k, s_k)\\
\end{eqnarray*}

El código en Matlab para implementar la división sintética y calculo de los coeficientes del polinomio de derivadas es

\begin{verbatim}
function [b, c] = Division_Sintetica_Derivadas(a, r, s)

n = length(a);

b=zeros(n,1);
c=zeros(n,1);
 
b(n) = a(n);
b(n-1) = a(n-1) - r*b(n);
c(n) = b(n);
c(n-1) = b(n-1) - r*c(n);
for k=n-2: -1:1
      b(k) = a(k) - r*b(k+1) - s*b(k+2);
      c(k) = b(k) - r*c(k+1) - s*c(k+2);
end;
\end{verbatim}

Finalmente el sistema de ecuaciones que debe resolverse es:

\[
\left[
\begin{array}{cc}
  c_2(r_k, s_k)   & c_3(r_k, s_k)   \\
  c_1(r_k, s_k)   & c_2(r_k, s_k)
\end{array}
\right]
\left[
\begin{array}{c}
  dr_k\\
  ds_k
\end{array}
\right]
=
\left[
\begin{array}{c}
  -b_1(r_k, s_k)\\
  -b_0(r_k, s_k) 
\end{array}
\right]
\]

y cada una de las iteraciones de Newton Raphson

\[
\left[
\begin{array}{c}
  r_{k+1}\\
  s_{k+1}
\end{array}
\right] = 
\left[
\begin{array}{c}
  r_{k}\\
  s_{k}
\end{array}
\right] 
-
\left[
\begin{array}{cc}
  c_2(r_k, s_k)   & c_3(r_k, s_k)   \\
  c_1(r_k, s_k)   & c_2(r_k, s_k)
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
 - b_1(r_k, s_k)\\
 - b_0(r_k, s_k) 
\end{array}
\right]
\]

La codificación en Matlab para llevar a cabo este procedimiento es:

\begin{verbatim}
function [c, r, s] = Itera_Bairstow(a, r0, s0)

N = length(a);
res =[r0, s0]';

for k=1:100
[b, c] = Division_Sintetica_Derivadas(a, res(1), res(2));

A = [c(3) c(4); c(2) c(3)];
d = [-b(2); -b(1)];

res = res - inv(A)*d;
if sqrt(norm(d)) < 0.0001 break; end
end;

r = res(1);
s = res(2);
c = b(3:N);
\end{verbatim}

 
\subsection{Ejemplo 1}

Dado el polinomio $f_5(x) =  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25$, determinar los  valores de $r$ y $s$ que hacen el resido igual a cero. Considere $r_0 = 1$ y $s_0 = -2$.

Solución.

\textbf{Iteración 1}

La división sintética con el polinomio $f_2(x) = x^2 +x - 2.0$ da como resultado

\begin{equation*} 
f_3(x) = x^3 - 4.5x^2 + 9.25x - 16.125  
\end{equation*}

\begin{equation*}
Residuo = [30.75, -61.75]
\end{equation*}

El polinomio de derivadas es

\[
c = [ -257.6250, 108.1250, -43.8750, 16.7500, -5.5000, 1.0000]
\]
 
Aplicando el método de Newton tenemos

\[
\left[
\begin{array}{r}
r_1\\
s_1\\
\end{array}
\right] = 
\left[
\begin{array}{r}
1\\
-2\\
\end{array}
\right] 
-
\left[
\begin{array}{rr}
-43.875 &16.750 \\
108.125 & -43.875\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
-30.7500\\
61.7500\\
\end{array}
\right]
=\left[
\begin{array}{r}
   -1.7637\\
   -7.4034
\end{array}
\right] 
\]

\textbf{Iteración 2}

La división sintética con el polinomio $f_2(x) = x^2 -1.7637x - 7.4034$ da como resultado

\begin{equation*}
f_3(x) = x^3 -  1.7363x^2 + 7.0911x - 1.7768
\end{equation*}
 
\begin{equation*}
 Residuo = [51.7564, 105.6858]
\end{equation*}

Los coeficiente del polinomio de derivadas es
\[
c_2 = [  677.3337, 208.1484, 27.6280, 14.5427, 0.0274, 1.0000]
\]
Aplicando el método de Newton tenemos

\[
\left[
\begin{array}{r}
r_2\\
s_2\\
\end{array}
\right] = 
\left[
\begin{array}{r}
   -1.7637\\
   -7.4034
\end{array}
\right] 
-
\left[
\begin{array}{rr}
27.62800 & 14.5427\\
208.1484 & 27.6280\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
-51.75640\\
-105.68578\\
\end{array}
\right]
=\left[
\begin{array}{r}
   -1.7164\\
   -3.9343\\
\end{array}
\right] 
\]
 
\textbf{Iteración 3}

La división sintética con el polinomio $f_2(x)=  x^2 - 1.7164x - 3.9343$ da como resultado


\begin{equation*}
f_3(x) = x^3 - 1.7836x^2 + 3.6229x + 1.3262
\end{equation*}
 
\begin{equation*}
Residuo = [12.6547, 28.1881]
\end{equation*}

Los coeficientes del polinomio de derivadas son
\[
c_3 = [  195.3505, 65.6792, 13.8350, 7.4418, -0.0672, 1.0000]
\]

Aplicando el método de Newton tenemos

\[
\left[
\begin{array}{r}
r_3\\
s_3\\
\end{array}
\right] = 
\left[
\begin{array}{r}
   -1.7164\\
   -3.9343\\
\end{array}
\right] 
-
\left[
\begin{array}{rr}
   13.8350  & 7.4418\\
   65.6792  &  13.8350 \\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -12.6547\\
   -28.1881
\end{array}
\right]
=\left[
\begin{array}{r}
   -1.5997\\
   -2.4507\\
\end{array}
\right] 
\]

En resumen

\begin{center}
\begin{tabular}{|r|r|r|r|r|}
\hline
$k$ & $r$ & $s$ & $b_1$ & $b_0$	\\
\hline
    0 &    1.0000 & - 2.0000 &    30.7500 &  -61.7500 \\
    1 &   -1.7637 &  -7.4034 &    51.7564 &  105.6858  \\
    2 &   -1.7164 &  -3.9343 &    12.6547 &  28.1881 \\
    3 &   -1.5997 &  -2.4507 &      2.8996 &   8.1547\\
    4 &   -1.3335 &  -2.1867 &      0.7601 &   2.5222\\
    5 &   -1.1183 &  -2.1130 &      0.2719  &   0.6077  \\
    6 &   -1.0271 &  -2.0232 &     0.0431  &   0.1119\\
    7 &   -1.0017 &  -2.0015 &     0.0028  &   0.0063 \\
    8 &   -1.0000 &  -2.0000 &     0.0000  &   0.0000\\
    9 &   -1.0000 &  -2.0000 &     0.0000  &   0.0000\\
  10 &   -1.0000 &  -2.0000 &     0.0000  &   0.0000\\
\hline
\end{tabular}
\end{center}

La solución es:

\begin{eqnarray*}
f_3(x) &=& x^3 - 2.5x^2 + 2.25x - 0.625 \\
f_2(x) &=& x^2 - x - 2
\end{eqnarray*}

Las raíces de $f_2(x) = x^2 - x - 2$, son $x_1 = 2$ y $x_2 = -1$

Para replicar el ejemplo damos

\begin{verbatim}
[c r s] =Itera_Bairstow([1.25, -3.875, 2.125, 2.75, -3.5, 1], 1, -2)

c =

   -0.6250
    2.2500
   -2.5000
    1.0000


r =

   -1.0000


s =

   -2.0000
\end{verbatim}

Si repetimos el ejemplo pero ahora considerando el polinomio $f_3(x) = x^3 - 2.5x^2 + 2.25x - 0.625$, podemos calcular el total de las raíces del polinomio original.

La implementación del Método de Bairstow es:

\begin{verbatim}
function raices = Bairstow (b, r0, s0)

raices = [];

N = length(b);

while(N > 3) 
    [b r s] = Itera_Bairstow(b, r0, s0);
    [x1, x2] = Formula_General(1, r, s);
    raices = [raices, x1, x2];

    N = N-2;
end;

if length(b) == 3
    [x1, x2] = Formula_General(b(3), b(2), b(1));
    raices = [raices, x1, x2];
end;

if length(b) == 2
    raices = [raices, -b(1)/b(2)];
end;
\end{verbatim}

\subsection {Ejemplo 2} 

Dado el polinomio $f_5(x) =  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25$, calcular la raíces utilizando el Método de Bairstow con valores iniciales $r_0 = 1$ y $s_0 = -2$ 

\textbf{Iteración 1}

Calculamos la deflación Polinomial 
\begin{equation*}
f_5(x) =   (-0.625   + 2.25 x   - 2.5 x^2    +x^3)   ( x^2 -x   -2)   
\end{equation*}

Las raíces calculadas son

\begin{eqnarray*}
x_1 &=& 2.0000\\
x_2 &=& -1.0000
\end{eqnarray*}

\textbf{Iteración 2}
\begin{equation*}
f_3(x) = (-0.625   + 2.25 x   - 2.5 x^2    +x^3)  = (x -0.5)( x^2 - 2x + 1.25)
\end{equation*}

Las raíces son
\begin{eqnarray*}
x_3 &=& 1 + j 0.5\\
x_4 &=& 1 - j0.5
\end{eqnarray*}

\textbf{Iteración 3}

Finalmente

\begin{equation*}
f_1(x) = x-0.5
\end{equation*}

tiene una raiz
\begin{equation*}
x_5 = 0.5
\end{equation*}

Para correr en Matlab hacer

\begin{verbatim}
Bairstow([1.25, -3.875, 2.125, 2.75, -3.5, 1], 1, -2)

ans =

   2.0000            -1.0000             1.0000 + 0.5000i   1.0000 - 0.5000i   0.5000   
\end{verbatim}

\subsection{Ejemplo 3}

Dado el polinomio $f_5(x) =  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25$, determinar las raíces de este polinomio. Considere $r_0 = -1$ y $s_0 = -1$.
 
\textbf{Iteración 1}

\begin{eqnarray*}
f_5(x) &=&  x^5 - 3.5x^4 + 2.75x^3 + 2.125x^2 - 3.875x + 1.25 \\
f_5(x) &=&( x^3 - 4x^2 + 5.25x - 2.5)*( x^2 +0.5x - 0.5)\\
\end{eqnarray*}

Las raíces de  $x^2 +0.5x - 0.5=0$ son 

\begin{eqnarray*}
x_1  &=& 0.5\\
x_2  &=&-1.0\\
\end{eqnarray*}

\textbf{Iteración 2}

\begin{eqnarray*}
f_3(x) &=& x^3 - 4x^2 + 5.25x - 2.5\\
f_3(x) &=& (x - 2)*(x^2 - 2x  +1.25)
\end{eqnarray*}
 

Las raíces de  $x^2 - 2x  +1.25=0$ son

\begin{eqnarray*}
x_3  &=& 1.0 + j0.5\\
x_4  &=& -1.0 - j0.5\\
\end{eqnarray*}
 
\textbf{Iteración 3}

\begin{equation*}
f_1(x) =(x - 2)
\end{equation*}
 
La raíz de este polinomio es

\begin{equation*}
x_5 = 2;
\end{equation*}
 

Todas la raíces de $f_5(x)$ son  $x = [0.5, -1.0, (1.0 + j0.5), (1 - j0.5), 2]$. Para correr en Matlab dar

\begin{verbatim}
>> Bairstow([1.25, -3.875, 2.125, 2.75, -3.5, 1], -1, -1)

ans =

   2.0000     -1.0000        1.0000 + 0.5000i   1.0000 - 0.5000i   0.5000    
\end{verbatim}