1. Modelaje de Sistemas Lineales

Ahora consideremos el sistema de ecuaciones lineales
ax + by = p

cx + dy = q

Nosotros podemos escribir esto más sólidamente como

AX = B

donde la matriz A de los coeficiente es:

a b

c d

el vector X de variables desconocidas es

x

y

el vector B en el lado derecho es

p

q

Si A es invertible, X = (1/A)B, o, usando anotación de MATLAB, X = A\B. Se prueba esto resolviendo AX = B.

Hagamos un poco de programación.
Sea A sea la matriz
0.8 0.1
0.2 0.9
y sea X el vector de la columna

1

0

Consideramos X para representar (por ejemplo) el estado de la población de una isla. La primera entrada (1) da el fragmento de la población en la mitad oriental de la isla, la segunda entrada (0) dé el fragmento en la otra mitad oriental. El estado de la población en las unidades de tiempo T son dadas después por la regla Y = AX. Esto expresa el hecho de que un individuo en medio de las estancias orientales puestas con probabilidad 0.8 y que se mueve del este con probabilidad 0.2 (notar que 0.8 + 0.2 = 1), y el hecho que un individuo en las estancias orientales puestas con probabilidad 0.9 y oeste de los movimientos esté con probabilidad 0.1. Así, los estados de la población sucesivos pueden ser predicho/computado por la multiplicación de la matriz repetida. Esto puede ser realizado por el programa de MATLAB siguiente:

>> A = [ 0.8 0.1; 0.2 0.9 ]

>> x = [1; 0]

>> for i = 1:20, x = a*x, end,

Hasta aquí hemos aprendido a escribir un tipo de loop simple. Ésta es una manera fácil de ordenar a la máquina, en otras palabras, hacer un trabajo muy repetitivo.
 

1.1. Definiendo Matrices

Si queremos definir la siguiente matriz en MATLAB:

entonces escribimos:

»A=[1 2 3 4;5 6 7 8;9 10 11 12;13,14,15,16];

»x=4:-1:1

genera el vector fila x=[4,3,2,1]. La instrucción

»C=A(3:4,1:3);

se refiere a la submatriz

de A. También D=A([1,3],3:4) genera

 

1.2. Matrices Especiales
En MATLAB podemos generar matrices especiales con las siguientes instrucciones:

rand(n,m) - matriz n x m de entradas aleatorias entre 0 y uno.
eye(n) - matriz identidad n x n.
zeros(n,m) - matriz cero de tamaño n x m.
ones(n,m) - matriz n x m con todas las entradas uno.

Combinando estas instrucciones podemos generar matrices bastante complicadas. Por ejemplo, la instrucción

»E=[eye(2),ones(2,3);zeros(2),[1:3;3:-1:1]]

genera la matriz

La instrucción round(x) redondea "x" al entero más cercano a "x". Podemos combinar funciones en MATLAB. Por ejemplo, round(10*rand(4)) genera una matriz con entradas aleatorias entre 0 y 10.
1.3. Aritmética de Matrices

Considere las siguientes matrices:

Entonces las operaciones A*B (producto matricial de A con B), A+B (suma de A mas B), 3*A (multiplicación escalar de 3 por A) tienen los siguientes resultados:

»A*B

ans =

??? Error using ==> +
Matrix dimensions must agree.

»3*A

ans =

»A’

ans =

4 2
5 3
»inv(A)

ans =

ans =

174 235
  94 127

Si precedemos las operaciones matriciales "*", "^" con el punto ".", entonces estas se hacen termino a termino. Por ejemplo A.*C y A.^2 generan:

» A.*C

ans =

-4 10
  4 12

» A.^2

ans =

16 25
  4   9
 

2. Solución de Sistemas Lineales
Considere el sistema lineal

Definimos la matriz de coeficientes y el lado derecho por las instrucciones:

»A=[1 -2 3; 4 1 -2; 2 -1 4];
»b=[1 -1 2]’;

Note que la transpuesta en b se usa para hacerlo un vector columna. Vamos a resolver este sistema por tres métodos:

» flops( 0 )
» x=A\b

x =

lleva a cabo eliminación Gaussiana en el sistema de arriba y produce como resultado:

ans =

73

entonces, aquí se necesitaron aproximadamente 73 operaciones de punto flotantes (sumas, restas, multiplicaciones ó divisiones) para resolver el sistema con eliminación Gaussiana.

Para el método de Gauss-Jordan tenemos:

» flops(0)
» rref([A b])

ans =

ans =

483

el cual requiere 483 operaciones de punto flotante.

Finalmente el método de la inversa se realiza con la siguiente secuencia de instrucciones:

» flops(0)
» x=inv(A)*b

x =

ans =

108

el cual toma 108 operaciones.

Vemos pues que la eliminación Gaussiana es el mejor de los tres métodos lo cual es cierto en general.

Usando MATLAB podemos estudiar la relación entre la solubilidad del sistema Ax=b y la no singularidad de la matriz de coeficientes A. En clase vimos que el sistema Ax=b tiene solución única para cualquier lado derecho b si y solo si la matriz A es no singular. ¿Qué sucede si A es singular? ¿Entonces Ax=b no tiene solución? Si A es singular el sistema Ax=b puede tener solución para algunos b’s pero de seguro hay al menos un b^* para el cual Ax=b^* no tiene solución. Vamos a generar una matriz singular con MATLAB:

» A=round(10*rand(6));

» A(:,3)=A(:,1:2)*[4 3]’

A =

» b=round(20*(rand(6,1)-0.5))

b =

» rref([A b])

ans =

Como la última fila es de la forma el sistema es inconsistente, i.e., no tiene solución. ¡Recuerde que A es singular! Esto no quiere decir que Ax=b nunca tenga solución. Si definimos c=A*b, con el b de arriba digamos, el sistema Ax=c tiene solución x=b (¿por qué?). De hecho si calculamos la forma echelon reducida de [A c] tenemos:

» c=A*b;
» rref([A c])

ans =

1 0 4 0 0 0 30
0 1 3 0 0 0 19
0 0 0 1 0 0 3
0 0 0 0 1 0 -9
0 0 0 0 0 1 3
0 0 0 0 0 0 0

el cual denota un sistema consistente dependiente con soluciones:

donde x₃ es arbitrario.
Recordemos que MATLAB posee una gran cantidad de funciones matriciales. De las más comunes que tenemos que repasarlas son:

3. Interpolación Polinomial

La función de interpolación es aquella función que pasa realmente por todos los puntos dados o aquella que "mejor" ajuste a esos puntos. Cuando se ve gráficamente, la línea pasa por cada uno de los puntos dados.

Para conveniencia seleccionamos un polinomio de grado n para n+1 pares de ordenadas (x,y).

En general el polinomio puede escribirse:

y = a0 + a1x + a2x2 + a3x3 + ....... + anxn

Para evaluar la interpolación polinomial de grado ‘n’ para cualquier conjunto de datos, vamos a evaluar los ‘n+1’ coeficientes, a0, a1, a2, a3, .....an.

Dado los ‘n+1’ pares de datos, podemos formas las ‘n+1’ ecuaciones

y1 = a0 + a1x1 + a2x12 + a3x13 + ...... + anx1n

y2 = a0 + a1x2 + a2x22 + a3x23 + ...... + anx2n

y3 = a0 + a1x3 + a2x32 + a3x33 + ...... + anx3n

. . . . . .

. . . . . .

. . . . . .

yn+1 = a0 + a1xn+1 + a2xn+12 + a3xn+13 + ... + anxn+1n

Este es el conjunto soluble de ‘n+1’ ecuaciones desconocidas

Los coeficientes desconocidos son a0, a1,......an. Si solo existen dos o tres ecuaciones podemos seguir con la Regla de Cramer para hallar estos coeficientes. Pero, en general vamos a recurrir al álgebra de matrices para manipular estas ecuaciones usando un formato de n sustituciones hacia atrás para obtener los coeficientes desconocidos.

Observe que el conjunto de ecuaciones puede expresarse convenientemente por la ecuación matricial:

[y] = [X][a]

donde,

[y] = [X] = [a] =
 
 

y1 1 x1 x12.....x1n+1 a1

y2 1 x2 x22.....x2n+1 a2

y3 1 x3 x32.....x3n+1 a3

. . . . . .

. . . . . .

. . . . . .

yn+1 1 xn+1 xn+12... xn+1n+1 an+1

Notamos que el producto de las matrices del lado derecho es una matriz columna de orden ‘n+1’x 1.

La solución que usa la eliminación de Gauss usamos la división por la izquierda de ‘MATLAB’ o

[a] = [X]\[y]

Usando este método, podemos encontrar los coeficientes de a para el polinomio de grado n que pasa exactamente a través de los ‘n+1’ puntos.

Si tenemos un conjunto grande de datos, cada uno de los datos se aparean incluyendo un error experimental, el polinomio de grado n no es una buena opción.

Los polinomios de grado mayores que cinco o seis tienen a menudo un comportamiento poco realista aunque la curva polinómica pase por cada punto dado. Como ejemplo, consideremos el siguiente conjunto de puntos:

x y

2 4

3 3

4 5

5 4

6 7

7 5

8 7

9 10

10 9

En este conjunto hay nueve puntos. A manera que el valor de x se incrementa, los valores de y también varían. Sin embargo, podemos observar que la distribución de los puntos no es lineal. Con ‘MATLAB’ vamos a crear dos vectores para estos datos:

x9 = [2:1:10];

y9 = [ 4 3 5 4 7 5 7 10 9 ];

Para observar estos puntos trazados, ejecutamos:

plot(x9,y9,’o’)

Como existe nueve puntos, es posible construir el polinomio de interpolación de grado 8:

y = a0 + a1x + a2x2 + a3x3 + ....... + a8x8

Para encontrar los coeficientes desconocidos, definimos el vector columna y

y = y9’

y la matriz X

X = [ones(1,9);x9;x9.^2;x9.^3;x9.^4;x9.^5;x9.^6;x9.^7;x9.^8]’

Notamos que X es definida usando la transpuesta, la función one( ) y el operador de de array ‘ .^ ‘. Con X y y así definidos, estos satisfacen la ecuación

[X][a] = [y]

Resolvemos para los coeficientes de la matriz a indicada anteriormente, entrando el comando

a = X\y

dando resultados en

a = 1.0e+003*

3.8140

-6.6204

4.7831

-1.8859

0.4457

-0.0649

0.0057

-0.0003

0.0000

Note que el noveno coeficiente a(8) parece tener el valor cero, realmente su valor es un número diferente de cero, aparece esto debido a que MATLAB está considerando 4 cifras significativas a la izquierda después del punto decimal. Para observar los coeficientes con cifras más significativas, entrar el comando:

format long

a

y el formato se cambia con 15 dígitos significativos. Claramente vemos que a(8) no es cero, es un número pequeño elevado a la potencia 8.

Ahora que ya contamos con los coeficiente, vamos a generar un número suficiente de puntos para crear una curva continua. Para el eje x, formamos un escala en el rango de

2 <= x <= 10 con incrementos de 0.1.

x = [ 2:.1:10 ];

Para y, calculamos el valor en el polinomio de octavo grado para cada uno de los componentes de x:

y =a(1)+ a(2).*x + a(3).*x.^2 + a(4).*x.^3 + a(5).*x.^4...

• a(6).*x.^5 + a(7).*x.^6 + a(8).*x.^7 + a(9).*x.^8;

Ahora plot(x,y) y los puntos dados (x9,y9).

plot(x,y,x9,y9,’o’)

Los resultados polinómicos parecen pasar por cada dato exactamente pero el polinomio no sirve para representar cualquier otro punto en el rango de x.
 
 

Por esta razón, no se debe intentar de usar la interpolación polinomial para representar datos experimentales.

3.1. Función MATLAB para los mínimos Cuadrados

Hay en MATLAB una función para que la ecuación encaje con el método de los mínimos cuadrados. El comando es

polyfit(x,y,n)

donde x, y son los vectores con los datos y ‘n’ es el orden del polinomio para el método de los mínimos cuadrados que se desea. El comando ‘polifit’ devuelve un vector cuyos elementos son los coeficientes del polinomio.

Los elementos del vector están en orden inverso, de lo que podemos anticipar que el primer elemento es el coeficiente del orden más alto de x y el último elemento es el coeficiente del orden más bajo (el orden 0 ó x0).
 

4. Números Reales y Complejos
4.1. Asignación de valores a variables:
Los números complejos se trabajan igual que los reales en lo que se refiere a asignación, a operaciones matemáticas y a comandos. A continuación, unos pocos ejemplos para mostrar como se realiza la asignación:
a = 25.203147;
b = 3;
c = 1 + 2j; ó también:
d = 1.5476 + 2.8*i; (el uso de j ó i es indiferente, desde que se tenga en cuenta la nota importante sobre el uso de las variables i y j )
d = 5.2347;

e = 3j;
 
 

4.2. Operaciones matemáticas simples:

Las operaciones simples son las siguientes:

5.0000 + 3.0000i

d = a ^ b da como resultado:

d =

61.3022 + 15.3369i
Nota importante sobre el uso de las variables i y j:
Puede usarse indistintamente las dos variables incorporadas (i ó j) y MATLAB no pone problema si se usan las dos al tiempo. Pero si se asignan las variables i y/o j en algún lugar del programa, esta variable perderá su valor como raíz de -1. Para clarificar esto es útil un ejemplo:

% Observación para cuando se trabaja con complejos.

i = 8

j = 9

c = 2 + 3*j
 

La salida del programa anterior es:
i =

8

j =

9

c =

29

Como se puede ver si se intentaba representar un complejo con la variable c, no se logró debido a que se cambiaron las variables i y j. Por lo tanto recomiendo que si se va a trabajar con complejos en un programa: Deje libres las variables i y j (no las utilice en contadores ó en otros propósitos, que no sean representar raiz de -1.
4.3. Comandos matemáticos para números (complejos y reales):

Los comandos matemáticos más empleados con números son:

4.3.1. Comando ABS

Calcula la norma de un complejo, o el valor absoluto de un real.
La sintaxis de la orden es:

Valor = abs(Número);

Valor es la norma del complejo si (Número es complejo) o el valor absoluto de Número (si es real).

Número puede ser un real o un complejo:
Si es Real: calcula el valor absoluto.

Si es Complejo: calcula la norma del complejo.
El siguiente ejemplo ilustra el uso de abs: (ver orden de programación DISP)

Calcula la raiz cuadrada de un complejo o de un real.

La sintaxis de la orden es:

Valor = sqrt(Número);

En Valor se almacena la raíz cuadrada del número.
Número puede ser un real o un complejo (si es real negativo, el resultado es un complejo)

El siguiente ejemplo ilustra el uso de sqrt:

Al correr el programa se obtienen como salida los siguientes resultados:

R =

12.347

raiz =

0 + 3.5138i

C =

2.6000 – 7.3i

raiz =

2.2748 – 1.6046i
 

4.3.3. Comando ANGLE
Calcula el ángulo de fase (en radianes) de una matriz (podría querer leer sobre matrices) con elementos complejos. Si la matriz sólo tiene un elemento, calcula el ángulo de fase de ese complejo.

La sintaxis de la orden es:

Valor = angle(Matriz);

Valor es una matriz que almacena el valor del ángulo de fase del complejo (de 0 a 2*pi) que ocupa la misma posición en Matriz (el ángulo de fase del elemento 1,1 lo almacena en la posición 1,1).
Matriz es una matriz (puede tener un solo elemento) cualquiera con componentes complejas (los reales forman parte de los complejos).
El siguiente ejemplo ilustra el uso de angle: (ver orden de programación DISP)

Al correr el programa se obtienen como salida los siguientes resultados:
C =

2.0000 0 - 2.7000i

3.0000 + 5.0000i 0.7000 - 5.0000i
 

c =

6.3000 + 7.2300i

Para la matriz:

v =
0 -1.5708

1.0304 - 1.4317
 

Para un complejo: (matriz de un solo elemento)

v =

0.8540
 

5. Integrales Definidas
Se solucionan numéricamente por medio del comando TRAPZ.
5.1. Comando TRAPZ
Calcula la integral definida entre dos límites de una función (área bajo la curva) representada por uno o dos vectores, como se explica más adelante. El cálculo de la integral se realiza numéricamente, por medio de una aproximación de la función a trapecios (En ningún momento calcula la integral simbólica).

Debido a que el cálculo de la integral es numérico, se deben construir vectores "decentes" para calcular la integral. Por esta razón, es fundamental aclarar las características de los vectores, con el fin de tener un criterio para decidir como construir el vector de forma apropiada. (para mayor claridad en el ejemplo de esta orden, puede ser necesario leer las secciones sobre FOR y Vectores y matrices).

La sintaxis de la orden es:
 

Valor = trapz([Vector,] Matriz);

Los símbolos [ ] significan que Vector es opcional.

Matriz puede ser una matriz o un vector. Una matriz si se desea calcular la integral definida para varias funciones en el mismo rango (entre los mismos límites). Un vector si se desea calcular la integral para una sola función (su tamaño tiene relación con el tamaño de Vector, esta relación se muestra en detalle en la explicación de Vector).

Vector es el vector de los valores para los cuales se desea calcular la integral, tal que si Matriz es:

Valor es donde se almacena el valor de la integral (un real si sólo se calculó para una función, y un vector fila si se calculó para varias).
5.2. Creación de vectores "decentes"
La integral se realiza aproximando la curva (función) a una serie de rectas, con el fin de aproximar el área bajo la curva a una serie de trapecios contiguos. Por lo tanto la aproximación es buena ("decente") si efectivamente la función se comporta como una recta (aproximadamente) en cada sub-intervalo determinado por el paso y número de puntos que se tomen. Una forma empírica de verificar que los vectores están bien construidos es por medio de la orden plot, ya que esta función dibuja los vectores, aproximando la función de la misma forma que trapz. Por lo tanto, si al dibujar la curva con plot, esta se ve "suave", los vectores están bien definidos.

El siguiente ejemplo ilustra el uso de trapz:

Al correr el programa se obtiene la siguiente salida:

Esto lo podemos lograr con las instrucciones:

» x=-5:.1:5;
» y=x.^2.*exp(-x.^2);
» plot(x,y)

La primera instrucción divide el intervalo [-5,5] en subintervalos de largo 0.1, la segunda instrucción evalúa la función en los puntos de la partición, y finalmente graficamos los resultados con plot. La instrucción plot tiene opciones para cambiar patrones del trazado, poner titulos, etc.

Supongamos ahora que queremos dibujar la superficie:

Esto lo hacemos con la secuencia de instrucciones:

» x=-5:.4:5;
» y=x;
» [X,Y]=meshgrid(x,y);
» Z=X.^2.*exp(-Y.^2);
» surf(X,Y,Z)

 

Las primeras dos instrucciones dividen los ejes de "x" y "y" en subintervalos de largo 0.4; la tercera instrucción genera una rejilla en el conjunto [-5,5]x[-5,5] con cuadraditos de lados 0.4 como se ilustra en la siguiente figura:

La cuarta instrucción evalúa la función en los puntos de la rejilla, y finalmente trazamos la superficie con surf.

Los temas tratados hasta ahora son suficientes para realizar programas sencillos y útiles. Los comandos disponibles en MATLAB son muchos más, pero los tratados aquí son los más frecuentemente necesitados. En caso de ser necesario emplear otras ordenes, es posible "buscar" la solución por medio de HELP (todo está en inglés), la cual lista los temas matemáticos que se pueden emplear (separados en librerías llamadas toolbox). HELP [toolbox] lista los comandos en la librería y HELP [comando] explica su uso y sintaxis.