Algunas funciones útiles para manipular polinomios en MATLAB

El día de hoy les traigo algunas funciones que podrían serles de utilidad para implementar un método numérico en MATLAB.

La siguiente es una lista de funciones incorporadas para manipular polinomios, a menos que se especifique lo contrario, tanto la entrada como la salida de cada función es un vector de coeficientes que definen un polinomio en  potencias descendentes de la variable independiente.

Función	Descripción
conv	Calcular el polinomio resultante del producto de dos polinomios
deconv	Calcular el polinomio resultante de la división de dos polinomios
poly	Crear un polinomio con raíces especificas
polyder	Calcular el polinomio resultante al derivar un polinomio.
polyval	Evaluar un polinomio  en un punto particular de x.
polyvalm	Evaluar una expresión polinómica de una matriz cuadrada A. Las potencias de A son evaluadas como productos matriz-matriz
Polyfit	Calcular los coeficientes de un polinomio , de grado n, que se ajusta a un conjunto de mínimos cuadrados de los datos de  (x , y)
residue	Calcular la expansión de las fracciones parciales del radio de dos polinomios, o, dados los vectores conteniendo los residuos, polos y coeficientes directos, calcula el radio correspondiente de los polinomios.
roots	Encuentra las “n” raíces de un polinomio de grado n.

Código en MATLAB para desarrollar el método de bisección

El programa responsable del proceso iterativo será el siguiente.

function z=Metbiseccion(a,b,err,fun)

while (abs(b-a)> err);

    fa=fun(a);

    fb=fun(b);

    z=a+((b-a)/2)

    fz=fun(z);

    if(fa*fz<0);

        b=z;

    else

        a=z;

    end

end

En un script distinto se implementará el siguiente código, el cual será el receptor de la ecuación problema.

function y=fun(x);

y=x.^3+4*x.^2-10;  %( aquí es donde debes escribir la ecuación problema)

end

No es estrictamente necesario hacer otro script para la función anterior, ya que se puede hacer uso de funciones simbólicas, por lo tanto sería posible escribir el siguiente código en el programa principal como sustitución para el código anteriormente descrito.

y=@(x.^3+4*x.^2-10)

Si realizas este procedimiento no olvides borrar el "@fun" y tampoco olvides escribir y=@(x.^3+4*x.^2-10) antes de llamar al programa Metbisección.

El programa principal, en un scrpt nuevo, será el siguiente.

err=0.001;

a=1; (Intervalo menor)

b=2; (Intervalo mayor)

z=Metbiseccion(a,b,err,@fun);

x=(1:0.1:2.5); ( este intervalo debe ser sustituido por uno en el cual se piensa que la raíz podría encontrarse)

hold on

plot (x,fun(x))

grid on

plot (z, fun(z), '*')

Recordar que "err" se refiere al margen de error deseado, esto debido a que el método de bisección , como todo método numérico, solo es capaz de aproximar el valor de la raíz con respecto del valor real, por lo tanto siempre habrá un error muy pequeño. Recordar también que  los script deben estar guardados en la misma carpeta, de lo contrario MARLAB marcará error. Sin duda hay muchas maneras para automatizar el código, el código aquí mostrado representa un esqueleto, puedes modificarlo y mejorarlo de la manera que mejor te convenga. 

Si te interesa aprender más sobre el método en lápiz y papel, visita nuestra entrada del método de bisección en la sección de "Métodos numéricos en papel".

Método de Gauss-Jordan en MATLAB, solución de sistemas de ecuaciones lineales

Este código te permitirá resolver un sistema de ecuaciones lineales de cualquier tamaño, para comenzar, debes crear un Script en el cual debes vaciar el siguiente código.

%Programa principal, este es el único script que debes correr
A=[1,1,0,3;2,1,-1,1;3,-1,-1,2;-1,2,3,-1]; %Los datos que contiene representan una matriz de 4x4, debes sustituirlos por tus datos propios
B=[4,1,-3,4]'; %Los datos que contiene representan una matriz de 4x1, debes sustituirlos por tus datos propios
x=EliminacionGaussiana(A,B)

Importante: 
A: representa la parte "izquierda" del sistema de ecuaciones, es decir, solo debes añadir los valores que se encuentran antes del signo = para cada ecuación, recuerda que este sistema trabaja en forma matricial. B representa una matriz vertical la cual contiene los resultados de cada una de las ecuaciones, recuerda que el apostrofe que se encuentra despues del ultimo corchete en B, le indica a MATLAB que la matriz es vertical, por lo tanto es muy importante no removerlo.

Posteriormente debes crear un nuevo Script, en el cual se deberá vaciar el siguiente código.

function x=EliminacionGaussiana(A,B)
% A es una matriz de orden NxN
% B es una matriz de orden Nx1
% x es una matraz de orden Nx1 que contiene la solucion de Ax=B
%Scienceprocedures.blogspot.mx
[N N]=size(A);
x=zeros(N,1);
C=zeros (1,N+1);
Aug=[A B];
for q= 1:(N-1)
    [Y,j]=max(abs(Aug(q:N,q)));
    C=Aug(q,:);
    Aug(q,:)=Aug(j+q-1,:);
    Aug(j+q-1,:)=C;
    if Aug(q,q)==0
        'El valor de A es irregular. No hay solucion o no es unica';
        break
    end 
    for k=q+1:N
        m=Aug(k,q)/Aug(q,q);
        Aug(k,q:N+1)=Aug(k,q:N+1)-m*(Aug(q,q:N+1));
    end
    
end
x=resource(Aug(1:N,1:N),Aug(1:N,N+1));
end 

Para continuar, debes crear un tercer Script con el código que se muestra a continuación.

function x=resource(A,B)
n=length(B);
x=zeros(n,1); x(n)=B(n)/A(n,n);
for k=n-1:-1:1
    x(k)=(B(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end 


Consideraciones: Recuerda que debes guardar los tres script en la misma carpeta, de lo contrario no correrán adecuadamente.
El programa principal contiene los datos para resolver un sistema de 4 ecuaciones con cuatro incógnitas a manera de ejemplo para el usuario.

Código en MATLAB para implementar el método de interpolación lineal.

Para implementar el método en MATLAB, se debe crear un script que contenga el siguiente código, este será el programa principal, en el cual se deben vaciar los valores de interés, este es el unico programa que se debe correr. El programa contiene los datos del problema que resolvimos en nuestra entrada "Métodos numéricos, ejemplo del método de aproximación lineal".

%Cienciaparacualquiera.blogspot.com
x=[1,4,5,7,9]
y=[2,4,6,7,10]
plot (x,y,'*r')
x0=0:0.1:10
P0=Metododeaproximacion(x,y,x0)
hold on
plot (x0,P0)
grid on
xlabel('x')
ylabel('y')

title('Plot using linear Interpolation method')

En un nuevo script, se debe vaciar el siguiente código, el cual contiene el método numérico, este programa no debe ser corrido, pero debes asegurarte de guardarlo en la misma carpeta que el script anterior.

function P0=Metododeaproximacion(x,y,x0)
s1=0;
s2=0;
s3=0;
p1=0;
p2=0;
p3=0;

for i=1:length(x)
    s1=(s1+(x(i))^2);
     s2=s2+(x(i));
     s3=s3+(x(i)*y(i));
end
s1=s1
s2=s2
s3=s3

for i=1:length(x)
    p1=p1+x(i);
    p2=(length(x));
    p3=p3+(y(i));
end
p1=p1;
p2=p2;
p3=p3;
a=(((p2*s3)-(p3*s2))/((p2*s1)-(s2*p1))) %"a" es el valor de la pendiente (a = slope value)
b=(((p1*s3)-(s1*p3))/((p1*s2)-(s1*p2)))
P0=((a.*x0)+(b))

PI=@(x0)(((a.*x0)+(b)))

Nota: El codigo aquí mostrado solo representa un esqueleto, si necesitas añadirle más funciones, eres libre de hacerlo.

Si quieres aprender más sobre el método, revisa la siguiente entrada
http://cienciaparacualquiera.blogspot.mx/2015/11/metodos-numericos-el-metodo-de.html

Implementación del método de Euler en MATLAB

Para implementar el método de Euler para ecuaciones diferenciales, se debe crear un script en MATLAB en el cual se deberá vaciar el siguiente código.

function Z=meteuler(n,a,b,y0,x0,P)
Z=0
h=(b-a)/n;
x=zeros(1,n+1)
y=zeros(1,n+1)
x(1)=x0;
y(1)=y0;
x=a:h:b
for i=1:n
y(i+1)=y(i)+h*((P(x(i),y(i))))
end;
Z=[x' y']
plot(x,y,'r')
grid on
hold on
plot (x, P(x,y),'b')
L=P(x,y)
hold on
%plot (1.1,P(1.1,y(1.1)),'*k')
%Cienciaparacualquiera.blogspot.mx
%plot (x,P(x,y))
title('grafica solucion, azul= derivada, rojo=Ec. original')
%Blue plot= Derivative.
%Red plot= Original equation.


Para continuar, el siguiente programa debe ser desarrollado en un nuevo script

%programa principal Euler
%Cienciaparacualquiera.blogspot.mx
%a=intervalo inferior (Lower interval)
%b=intervalo superior (Upper interval)
%n= numero de iteraciones ( number of iterations)
%x0= primer valor de x (es igual al valor de a) (First value of x, which is equal to "a")
%y0= valor de y cuando x=0 (Initial condition for y)
a=0
b=1
n=10
x0=a
h=0.1
x=(0:h:1)
y0=0
P=@(x,y)(cos(x));
Z=meteuler(n,a,b,y0,x0,P,h)
hold on

Importante

• El usuario debe vaciar los datos necesarios en el programa principal.
• El programa principal contiene los datos para graficar la ecuación diferencial y'=cos(x), la cual se muestra de color azul, mientras que la ecuación original (solución) se muestra de color rojo, los datos pueden ser sustituidos de la manera que el usuario disponga.
• Es importante tomar en cuenta que el método de Euler no ofrece una solución para una ecuación diferencial, es decir, no es posible encontrar la ecuación original a través de este método, sin embargo, este método nos ofrece una buena aproximación de la gráfica de la función original. 
• Recuerda guardar ambos scripts en la misma carpeta, de lo contrario MATLAB no podrá llevar a cabo el procedimiento y arrojara un error. 

Si quieres comprender más acerca de las ecuaciones diferenciales, revisa la siguiente entrada 

http://cienciaparacualquiera.blogspot.mx/2015/10/ecuaciones-diferenciales.html

.

Ejemplo de Integración por el método del trapecio.

El método del trapecio es un método numérico que es de gran utilidad para aproximar la integral de una función, es el método más sencillo para aproximar la integral de una función, sin embargo también es, comparativamente, el más inexacto. Ejemplificaremos la manera de aplicar el método en lápiz y papel a través de un ejemplo. 

A continuación se muestra lo formula general para la solución de una integral por el método del trapecio.

Para ejemplificar el uso de la ecuación anterior debemos plantear una ecuación para encontrar su integral. Utilizaremos la siguiente ecuación para nuestro ejemplo.

 El siguiente paso para implementar el método del trapecio es definir h. Si la integral está definida en el intervalo [a , b] entonces h se definirá de la siguiente manera. 

h=(b-a)/n

n= número de subintervalos. 

El valor de n es definido por la persona que realiza el procedimiento. 

Nota: Mientras más subintervalos se utilicen, el resultado obtenido del procedimiento será mas exacto. 

Para nuestro ejemplo utilizaremos un valor de n igual a 4, por lo tanto el valor de h quedará de la siguiente manera.

h=( 5-1)/(4)=1

Si reformulamos la ecuación del método del trapecio para utilizarlo en nuestra ecuación problema, la formula quedaría de la siguiente manera.

Nota: En cualquier ecuación que se quiera resolver por este método, tanto el primero como el ultimo termino de la ecuación serán los únicos que no estén multiplicados por 2, el resto lo estarán. 

Para efectos de lo anterior

 f(x1)= f(a)

f(x2)= f(a+h)

f(x3)=f(a+2h)

f(x4)=f(a+3h)

f(x5)=f(a+4h)

Sustituyendo términos la ecuación quedaría de la siguiente manera. 

f (a) =1

f (a+h) = f(2) = 1/2

f (a+2h) = f(3) = 1/3

f (a+3h) = f(4) = 1/4

f (a+4h) = f(5) = 1/5

Si lo revisamos, el valor real de la integral de nuestra ecuación problema es de 1.6094, por lo que al usar este método obtuvimos una muy buena estimación del valor de la integral, si hubiéramos utilizado un valor de "n" más grande hubiéramos obtenido una aproximación mucho más exacta. 

Métodos numéricos, el método de bisección.

El método de bisección en métodos numéricos es utilizado para obtener la raíz de una ecuación. El termino "raíz" se refiere al valor de x que tiene la capacidad de hacer que la ecuación planteada tenga un valor igual, o muy cercano, a cero. 

Para implementar el método de bisección en lápiz y papel primero debemos plantear una ecuación problema, nuestra ecuación problema será la siguiente.

f(x)=( x^3)+(4*x^2)-10

El siguiente paso para implementar el método es proponer un intervalo entre el cual creemos que se encuentra la raíz de la ecuación, al valor más pequeño del intervalo le llamaremos "a" y al mayor le llamaremos "b".

La única condición que debe cumplir el intervalo propuesto es que el producto de f(a) y f(b) sea menor a cero.

f(a) * f(b) < 0 

The product of f(a) and f(b) must be minor to cero for the ethod to work.

Para nuestra ecuación problema propondremos el intervalo abierto de 1 a 2.

f(a)= 1+4-10 = -5

f(b)= 8+16-10 = 14

Por lo anterior comprobamos que nuestro intervalo cumple con la condición mencionada en el punto anterior.

El siguiente paso será definir el valor de P.

P= a + (b-a)/2

Por lo tanto, para nuestro problema, P tendrá el siguiente valor.

P=1+ (2-1)/2 = 1+ (1/2)= 3/2

El siguiente paso consiste en sustituir el valor de P en nuestra ecuación.

f(P)=  (3.375)+9-10= 2.375

El siguiente paso del método será decisivo, requiere de la implementación de las siguientes 3 reglas.

- Si f (P)= 0 entonces P es la raíz de nuestra ecuación. 
If f(P) = 0 then P is the root of the equation

- Si f (P) y f(a) tiene el mismo signo, entonces la raíz de la ecuación se encontrará en el intervalo de [P , b]
If f(P) and f(a) have the same sign, then the root will be between the interval [P , b] 

- Si f(P) y f(a) tienen signo distinto, entonces la raíz de nuestra ecuación se encuentra en el intervalo de [a , P].
If f(P) and f(a) have different sign, then the root will be between the interval [a , P] 

Los anterior significa que , para nuestra ecuación problema el nuevo intervalo se encontrará en

[1 , 2.375]

f(a)= 1+4-10 = -5

f(b)= 8+16-10 = 14

f(P)=  (3.375)+9-10= 2.375

El procedimiento anteriormente señalado debe repetirse hasta que se cumpla la primera de las 3 reglas antes mencionadas, o bien, el resultado sea tan pequeño que pueda volverse cero. Por lo tanto, continuaremos con el procedimiento hasta encontrar la raíz de nuestra ecuación problema.

* Segunda iteración.

 [1 , 2.375]

f(a)= 1+4-10 = -5

f(b)= (13.39) + 22.56 - 10 = 25.95

 P= 1+ (2.375-1)/2 = 1+ 0.6875 = 1.6875

f(P)= 4.805 + 11.39 - 10 = 6.1956

*Tercera iteración [1 , 1.6875]

f(a)= 1+4-10 = -5

f(b)=  4.805 + 11.39 - 10 = 6.1956

P=1+(1.6875-1)/2 = 1+0.3437 = 1.3437

f(P) = 2.4260 + 7.222 -10 = -0.3518

*Tercera iteración [ 1.3437 , 1.6875]

f(a)= 2.4260 + 7.222 -10 = -0.3518

f(b)=  4.805 + 11.39 - 10 = 6.1956

P= 1.3437 +(1.6875-1.3437)/2= 1.5156

f(P)= 3.4813+ 9.1881 - 10 = 2.6694

*Cuarta iteración [1.3437 , 1.5156]

f(a)= 2.4260 + 7.222 -10 = -0.3518

f(b)=  3.4813+ 9.1881 - 10 = 2.6694

P= 1.3436+(1.5156 - 1.3437)/2 = 1.4296

f(P) = 2.8220 + 8.1750 - 10 = 0.9970

*Quinta iteración. [1.3437 , 1.4296]

f(a)= 2.4260 + 7.222 -10 = -0.3518

f(b) = 2.8220 + 8.1750 - 10 = 0.9970

P= 1.3437 + ( 1.4296-1.3437)/2 = 1.3866

f (P) = 2.666 + 7.69 - 10 = 0.35063

*Sexta iteración [1.3437 , 1.3866]

f(a)= 2.4260 + 7.222 -10 = -0.3518

f(b) = 2.666 + 7.69 - 10 = 0.35063

P= 1.3437+ (1.3866-1.3437)/2 = 1.36515

f(P)= 2.5441 + 7.4539 - 10 = - 0.00190

Se puede observar que en la sexta iteracion f (P) es aproximadamente cero, por lo tanto, la raíz de la ecuación es aproximadamente 1.3651.

Como se puede ver, el procedimiento es largo y tedioso, y, si te estabas preguntando para que sirve la tan mencionada raíz, pues resulta que la raíz representa un valor de x que satisface la ecuación, por lo tanto, si aplicamos este método a una ecuación cuadrática sería posible obtener el valor de x. A través del uso de este método es posible encontrar solución a ecuaciones cuadráticas  o, incluso, de niveles más altos siempre y cuando se cumplan todas las reglas antes mencionadas. Por supuesto no tendría sentido aplicar el método a mano debido a la cantidad de esfuerzo que ello representa y al desconocimiento de la cantidad de iteraciones necesarias para llegar a una solución precisa, por fortuna el uso de los lenguajes de programación hace mucho mas fácil está tarea. Si quieres saber más sobre el código en MATLAB para este procedimiento, revisa nuestra sección de Métodos numéricos en MATLAB