viernes, 2 de diciembre de 2011

OPERACIONES DE POLINOMIOS EN MATLAB


POLINOMIOS EN MATLAB


En MATLAB un polinomio se representa mediante un vector fila que contiene los coeficientes de las potencias en orden decreciente: empezando por el coeficiente principal y terminando por el término independiente.

Por ejemplo, el polinomio p(x)=3x2-2x-1 se representa con
        p = [3  -2  -1];
MATLAB contempla las siguientes operaciones básicas con polinomios:
Cálculo de las raíces a partir de la lista coeficientes, por medio del comando roots( ), por ejemplo
        r=roots(p)
nos devuelve
        r =    1.0000
              -0.3333
El resultado es un vector columna de ceros.
Cálculo de los coeficientes a partir del vector columna de ceros, por medio del comando poly( ), por ejemplo
        poly(r)
nos devuelve
        ans =    1.0000   -0.6667   -0.3333
Observe que el polinomio devuelto siempre es mónico.
Multiplicación de dos polinomios dados por la lista de sus coeficientes, por medio del comando conv( , ). Por ejemplo, para comprobar que (x-5) (x+1)=x2-4x-5 basta ejecutar
        conv([1  -5], [1  1])
obteniendo
        ans =    1    -4    -5
La división se realiza por medio del comando deconv( , ): si p(x)=s(x) q(x) + r(x), se puede usar el formato
        [s, r] = deconv(p, q])
Evaluación de un polinomio dado por la lista de sus coeficientes p en un valor x, por medio del comando polyval(p, x). Por ejemplo, para comprobar que p(1)0= basta realizar
        polyval(p, 1)
obteniendo
        ans =    0
polyval( ) realiza la evaluación siguiendo el algoritmo de Horner o de multiplicación anidada.

Si x es un vector o una matriz, MATLAB devuelve la matriz con el polinomio evaluado en cada elemento.
Recordemos que dados unos nodos de interpolación x=[x0, x1, ..., xn], los polinomios básicos de Lagrange se definen por la fórmula
li(t)=
(t-x0)... (t-xi-1)(t-xi+1)... (t-xn)

(xi-x0)...(xi-xi-1)(xi-xi+1)... (xi-xn)
 ,    i=0, 1, ..., n .
Entonces el polinomio pn de menor grado que interpola la nube de puntos (x0, f0), ..., (xn,fn) está dado por
pn(t)=f0 l0(t) + ... + fn ln(t) .

Problema 1   Vamos a construir los polinomios básicos de Lagrange para los nodos x=[ -3, -1, 1], usando los comandos vistos hasta ahora. Un procedimiento, aunque poco eficiente, podría ser:
  1. construir el polinomio nodal (el polinomio mónico que se anula en todos los nodos de interpolación) por medio de la función poly;
  2. para cada uno de los nodos eliminar el factor que contiene ese nodo, dividiendo con deconv el polinomio nodal entre el factor lineal que se anula en dicho nodo;
  3. normalizar cada polinomio, dividiéndolo entre su valor en el nodo omitido; utilizar polyval.
Dibuje las gráficas de cada uno de los polinomios obtenidos en el intervalo [-4, 2] y compruebe visualmente que en efecto los polinomios obtenidos satisfacen la propiedad li(xj)=dij. Para ello tabule los valores de cada polinomio en 200 puntos equiespaciados del intervalo, creados con
 
        linspace(-4, 2, 200);
Problema 2   Usando los resultados del problema anterior construya el polinomio de menor grado que interpole la función f(t)=exp (t) en los nodos -3, -1, 1. Dibuje su gráfica junto con la de la función f en el intervalo [-4,2] y compruebe visualmente que el polinomio obtenido interpola a f en los nodos indicados.

Para diferenciar las dos curvas puede después de plot usar el comando
 
        legend('Funcion', 'Polinomio');
Problema 3  [Opcional]   Escribir una función lagrange con el formato
 
        function [c]=lagrange(x, f)
que dados los nodos (vector x) y los valores correspondientes de ordenada (vector f) devuelva los coeficientes del polinomio interpolador de Lagrange.

Polinomios
Se define un polinomio de grado n como
(2)p_n(x) = a_0 + a_1 x + a_2 x^2 + \ldots + a_{n-1}x^{n-1} + a_n x^n
No es más que una función en la que el valor de la variable se eleva sucesivamente a una potencia hasta n y se multiplica por una constante. Utilizando el símbolo del sumatorio la expresión anterior puede compactarse a:
p_n(x) = \sum_{i=0}^n a_i x^i
Si nos fijamos un momento en la expresión (2) observaremos que un polinomio puede expresarse fácilmente en forma de vector utilizando sus coeficientes. El orden puede deducirse fácilmente con el número de coeficientes. Matlab utiliza vectores para expresar los polinomios con la única salvedad que los almacena del modo inverso al que hemos escrito (2). El polinomio x^3-x+1sería en Matlab
>> p = [1, 0, -1, 1];
La operación más común con un polinomio es evaluarlo en un punto dado, para ello utilizaremos la función polyval.
polyval(p, x)
Evalúa el polinomio p en el punto x
Ejemplo
>> p = [1, 0, -1, 1];
>> polyval(p,3)
ans =  25
La importancia de los polinomios es que, siendo una función, todas las operaciones elementales (suma, resta, multiplicación y división) pueden reducirse sólo a operaciones con sus coeficientes. De esta manera podemos convertir operaciones simbólicas en operaciones puramente numéricas. Tomemos por ejemplo estas dos funciones: p(x) = 4x^3-xy q(x) = x^2 + 6. Sumar y restar estas dos funciones es trivial, pero no multiplicarlas. Como se trata de una operación con coeficientes Matlab la hará sin inmutarse
>> p = [4, 0, -1, 0];
>> q = [1, 0, 6];
>> conv(p,q)
ans =
    4    0   23    0   -6    0
conv(u, v)
Calcula la convolución de dos vectores de coeficientes. En el caso de vectores, la convolución es la misma operación que el producto.
Efectivamente p(x)*q(x) = 4x^5+23x^3-6x.
Dividir dos polinomios nos servirá para aprender cómo tratar las funciones con dos argumentos de salida. De define la división de dos polinomios como
p(x) = q(x)*c(x) + r(x)
Entonces la división entre p(x)y q(x)tiene como resultado dos polinomios más, el cociente c(x)y el residuo r(x). Si a la salida de deconv se le asigna sólo una variable obtendremos el cociente
deconv(u, v)
Calcula la deconvolución de dos vectores de coeficientes. En el caso de polinomios esta operación es equivalente al cociente del primero por el segundo.
Devuelve dos argumentos, el cociente y el residuo.
>> c = deconv(p,q)
c =
    4  0
Si necesitamos también el residuo tendremos que hacer lo siguiente
>> [c,r] = deconv(p,q)
c =
    4  0

r =
    0  0  -25  0
Hay otras operaciones que son operadores lineales aplicados a los polinomios como por ejemplo la derivada y la integral.
polyderiv(p)
Calcula los coeficientes de la derivada del polinomio p. Si le proporcionamos un segundo argumento q calculará la derivada del producto de polinomios.
polyinteg(p)
Calcula los coeficientes de la primitiva del polinomio p.
Sabemos que los polinomios de orden n tienen el mismo número de raíces. Esto nos permite descomponer cualquier polinomio de la siguiente manera:
p_n(x) = \sum_{i=0}^n a_i x^i = \prod_{i=0}^n (x-r_i)
roots(p)
Calcula las raíces del polinomio p.
Las raíces no son sólo importantes como solución de la ecuación p_n(x) = 0sino que sirven, por ejemplo, para buscar factores comunes entre dos polinomios. Otra función bastante útil para los que utilizan Matlab para el análisis de sistemas dinámicos lineales es la función residue que calcula la descomposición en fracciones parciales del cociente de dos polinomios
residue(p, q)
Calcula la descomposición en fracciones parciales del cociente de dos polinomios p y q donde el primero es el numerador y el segundo el denominador.
Por ejemplo
>> b = [1, 1, 1];
>> a = [1, -5, 8, -4];
>> help residue
>> [r,p,k,e] = residue(b,a)
r =

  -2.0000
   7.0000
   3.0000

p =

   2.00000
   2.00000
   1.00000

k = [](0x0)
e =

   1
   2
   1
El vector r es el numerador de cada uno de los términos, el vector p son los polos del sistema y el vector e es la multiplicidad de cada uno de los polos. Entonces la descomposición en fracciones parciales será:
\frac{s^2+s+1}{s^3-5s^2+8s-4} = \frac{-2}{s-2} +
\frac{7}{(s-2)^2} + \frac{3}{s-1}
Ejercicio de síntesis
Existe una manera de representar la forma de una función cualesquiera en un punto dado mediante un polinomio. Dicho polinomio converge con mayor orden en los alrededores del punto a medida que se van añadiendo términos. Se trata del desarrollo de Taylor.
La única información que necesitamos de la función es su valor y el de sus derivadas en el punto dado x_0. La expresión general es
p_n(x-x_0) = f(x_0) + \sum_{i = 1}^n f^{(i)}(x_0)\frac{(x-x_0)^i}{i!}
Para entender mejor cómo este polinomio se ajusta a la función podemos utilizar el desarrollo de la función exponencial en x=0.
e^x = 1 + x + \frac{1}{2} x^{2} + \frac{1}{6} x^{3} + \frac{1}{24}
x^{4} + \frac{1}{120} x^{5} +
\operatorname{\mathcal{O}}\left(x^{6}\right)
Este polinomio puede crearse de muchas maneras pero posiblemente la más sencilla sea utilizar los polinomios en Matlab para tener que generar sólo los coeficientes.
>> exp_serie = @(x,n) polyval(1./[factorial(linspace(n,1,n)),1],x)
exp_serie =

@(x, n) polyval (1 ./ [factorial(linspace (n, 1, n)), 1], x)
Nota
Esta línea de código sirve para aprender una regla muy importante sobre cómo debemos escribir un programa. Las líneas demasiado largas son difíciles de leer, por lo tanto son un peligro incluso para nosotros mismos. Es recomendable romperlas en algún punto donde romperíamos una operación matemática: después de un operador, justo después de abrir un paréntesis. Para hacerlo debemos escribir tres puntos ....
Podemos utilizar esta función para entender de un modo mucho más visual el concepto de convergencia de una serie. Sabemos que a medida que añadamos términos el error que comete el desarrollo de Taylor cerca del punto se reduce. ¿Pero de qué forma? Una confusión habitual es pensar que al aumentar orden del desarrollo aumenta la región donde se acerca a la función pero esto sólo es cierto accidentalmente. Sólo existe una mayor convergencia cerca del punto.
Para verlo mejor calcularemos el error de la aproximación en los puntos 0.2 y 0.1 para distintos órdenes.
exp_serie = @(x,n) polyval(1./[factorial(linspace(n,1,n)),1],x)

x_01 = [exp_serie(0.1,1),
               exp_serie(0.1,2),
               exp_serie(0.1,3),
               exp_serie(0.1,4),
               exp_serie(0.1,5),
               exp_serie(0.1,6),
               exp_serie(0.1,7)];

x_02 = [exp_serie(0.2,1),
               exp_serie(0.2,2),
               exp_serie(0.2,3),
               exp_serie(0.2,4),
               exp_serie(0.2,5),
               exp_serie(0.2,6),
               exp_serie(0.2,7)];


disp('error en 0.1')
err_01 = abs(exp(0.1)-x_01)
disp('error en 0.2')
err_02 = abs(exp(0.2)-x_02)

disp('logaritmo del error en 0.1')
logerr_01 = log(err_01)
disp('logaritmo del error en 0.2')
logerr_02 = log(err_02)
El programa anterior tiene la siguiente salida:
error en 0.1
err_01 =

   5.1709e-03
   1.7092e-04
   4.2514e-06
   8.4742e-08
   1.4090e-09
   2.0092e-11
   2.5091e-13

error en 0.2
err_02 =

   2.1403e-02
   1.4028e-03
   6.9425e-05
   2.7582e-06
   9.1494e-08
   2.6046e-09
   6.4932e-11

logaritmo del error en 0.1
logerr_01 =

   -5.2647
   -8.6743
  -12.3683
  -16.2837
  -20.3804
  -24.6307
  -29.0137

logaritmo del error en 0.2
logerr_02 =

   -3.8442
   -6.5693
   -9.5753
  -12.8009
  -16.2070
  -19.7660
  -23.4577
Podemos ver que si tomamos logaritmos la diferencia entre los valores permanece aproximadamente constante.

POLINOMIOS


Con MATLAB se puede trabajar con polinomios de forma sencilla. Es suficiente tener en cuenta que un polinomio no es nada más que un vector, en que el orden de los coeficientes va de mayor
a menor grado.
Ejemplos. 

EDU>> p=[3 5 2 8 6]   % 3*x^4+5*x^3+2*x^2+8*x+6

p =
3 5 2 8 6


EDU>> q=[6 2 1 7 8]   % 6*x^4+2*x^3+x^2+7*x+8

q =
6 2 1 7 8

Además, MATLAB incluye funciones específicas para operar con polinomios. Por ejemplo,
si queremos evaluar lo que vale un polinomio en un punto.
Ejemplo.

EDU>>
polyval(p,-1)  % Evaluación de 3*x^4+5*x^3+2*x^2+8*x+6 en x=-1
ans =
-2

También es posible multiplicar dos polinomios.
Ejemplo.

EDU>>
conv(p,q)  % producto de p por q

ans =
1 36 25 78 113 74 78 106 48


O obtener el cociente que se obtiene al dividirlos.

EDU>>
deconv(p,q) % cociente resultado de dividir p entre q

ans =
0.5000

EDU>> roots(p)    % Raices del polinomio p


ans =
-1.7793         
0.4292 + 1.1502i
0.4292 - 1.1502i
-0.7458    
Polinomios
En Matlab los polinomios se representan por vectores cuyas componentes son los coeficientes del polinomio.
Sea
http://www.fisica.unav.es/%7Eangel/matlab/poli.gif
Este polinomio se representa por un vector p
» p=[1 -3 +2]
p =
1   -3   2
Para hallar las raíces del polinomio, se hace
» roots(p)
ans =
2
1
y si se quiere hallar el valor de P(x) para un determinado valor de x (por ejemplo, para x=0)
» polyval(p,0)
ans =
2

1 comentario: