DEPARTAMENTO DE MATEMÁTICA APLICADA A LA INGENIERÍA INDUSTRIAL

GUION PARA LA PRÁCTICA 1 DE ÁLGEBRA

Nociones sobre resolución numérica de sistemas de ecuaciones lineales

Antes de trabajar con este guion es necesario haber leı́do con detalle los tres archivos siguientes que ya
habrá descargado de Moodle:
• 1_Instrucciones_para_Practicas_MatInd.pdf,
• 2_Inicio_Sesion_MATLAB.pdf,
• 3_Introduccion_MATLAB.pdf.

Índice
1. Errores de redondeo y aritmética de precisión finita 2

2. Existencia y unicidad de la solución de un sistema lineal 4

3. Uso de rutinas o funciones en MATLAB: la función gauss 6

4. La barra invertida de MATLAB 8

5. Número de condición de una matriz 9

6. Cálculo simbólico 11

INSTRUCCIONES. Este guion pretende ser un documento autocontenido y está pensado para que lo trabaje
con anterioridad a la realización de la práctica en el aula de ordenadores. Debe, por tanto, leerlo con el detalle
y la atención necesaria para entender todo lo que en él se desarrolla. Siga por orden las instrucciones que se
indican a continuación. Es muy importante asegurarse de que se llega al resultado final en todos los pasos.
1. Antes de empezar, abra una sesión de MATLAB según se describe en el archivo
2_Inicio_Sesion_MATLAB.pdf
que se ha bajado de Moodle y que ya ha debido leer con suficiente detalle. Como en él se indica, al finalizar
el proceso debe haber creado un directorio de trabajo llamado NMat, donde NMat debe ser su número de
matrı́cula con cinco dı́gitos. Asegúrese de que este es el directorio activo cuyo contenido, de momento vacı́o,
es el que debe aparecer en la ventana Current folder.
2. Descargue de Moodle el fichero gauss.m y guárdelo en el directorio de trabajo creado en el paso anterior
—que debe ser el directorio activo—. El fichero gauss.m debe entonces mostrarse en la ventana Current
folder.
3. En la sesión de trabajo que acaba de abrir, deberá ir ejecutando todas las instrucciones de MATLAB que
contiene este guion. Al ejecutar estas instrucciones tenga en cuenta que nunca será necesario que teclee los
caracteres ✭✭ %✮✮ que aparecen ni tampoco el texto que figura a su derecha. Además, recuerde incluir siempre
el carácter ✭✭;✮✮ al final de lı́nea cuando se le indique.
4. Teclee en su sesión de MATLAB la instrucción
>> format short e
para seleccionar el formato en que MATLAB va a mostrar los números por pantalla a lo largo de todo este
guion. No cambie este formato a no ser que se le indique. Es importante tener presente que los diferentes
formatos solo muestran diferentes presentaciones del mismo número. Ası́, el cambio de formato no cambia
la manera en que MATLAB efectúa las operaciones ni su nivel de precisión.

1
1. Errores de redondeo y aritmética de precisión finita
Se trata en esta práctica de hacer un recorrido introductorio por las cuestiones y dificultades que surgen al
resolver numéricamente un sistema de ecuaciones lineales mediante algoritmos que deben programarse en un
ordenador.
A la hora de utilizar un ordenador para efectuar cálculos numéricos, la primera propiedad importante a
tener en cuenta es que el conjunto de números que el ordenador es capaz de manejar es finito: se denominan
✭✭números de máquina✮✮ y constituyen un conjunto F finito de números racionales. Como consecuencia, todos los
números reales con los que se trabaja en un ordenador se aproximan por números de máquina, que son los que el
ordenador reconoce. Esta aproximación se llama redondeo y el error que se comete al utilizarla es el denominado
error de redondeo.
En otras palabras, lo que ocurre al usar un ordenador es que todo número real x ∈ R se sustituye por un
número máquina x ! ∈ F cercano que se obtiene truncando el desarrollo binario (en base 2) de x. Además, como
F es finito, infinitos números reales próximos entre sı́ se sustituyen por el mismo número máquina. Los números
x∈Ryx ! ∈ F cumplen la relación
|x − x!| 1
≤ eps, (1)
|x| 2
donde eps se especifica de antemano y es el denominado épsilon de la máquina. Se aprecia de la desigualdad
(1) que el número eps mide la precisión relativa con la que MATLAB trabaja.
El conjunto F y eps dependen del ordenador que se utilice y de como este codifique los números. Sin entrar
en mucho detalle, diremos que MATLAB utiliza el denominado estándar IEEE de aritmética en coma flotante de
doble precisión, que fija el conjunto F y también eps, cuyo valor se obtiene tecleando en su sesión de MATLAB

>> eps
ans =
2.2204e-16
Además, en el conjunto F de MATLAB

>> realmin
ans =
2.2251e-308
da el número positivo más pequeño, y

>> realmax
ans =
1.7977e+308
es el número positivo más grande.
Las siguientes instrucciones de MATLAB ponen de manifiesto el significado de eps:

>> a=1+eps
a =
1.0000e+00
>> a-1
ans =
2.2204e-16

Es decir, para MATLAB, los dos números reales distintos 1 y 1 + eps son también números máquina distintos.
En particular, observe que la expresión 1.0000e+00 no representa al número entero 1. Sin embargo,

>> b=1+eps/2
b =
1
>> b-1
ans =
0

2
Es decir, para MATLAB, los dos números reales distintos 1 y 1 + eps/2 son el mismo número máquina y, por
tanto, su diferencia es cero para la aritmética del ordenador.
MATLAB incorpora una función eps(x) que mide la distancia de |x| al número máquina más próximo y
mayor que |x|. Obsérvese que la precisión con la que MATLAB maneja el número x —que es lo que indica
eps(x) si x no es un número máquina— depende de la magnitud de x, ya que de (1) se deduce la desigualdad

|x|
|x − x
!| ≤ eps.
2
En concreto, a mayor tamaño de x menor precisión, es decir, como ya hemos mencionado, MATLAB mantiene
la precisión relativa. Por ejemplo,

>> eps(1)
ans =
2.2204e-16
>> eps(100)
ans =
1.4211e-14
>> eps(1.0e+07)
ans =
1.8626e-09

Por otro lado, no hay que confundir la precisión con el formato de presentación. Analice el ejemplo siguiente,
para el cual vamos a cambiar temporalmente el formato, con el fin de que se aprecie mejor la idea subyacente.
>> format long
>> f=7+eps(7)
f =
7.000000000000001
>> f-7
ans =
8.881784197001252e-16
>> g=7+eps(7)/2
g =
7
>> g-7
ans =
0
Vuelva a ejecutar las instrucciones anteriores después de teclear nuevamente

>> format short e

y aprecie la diferencia en presentación que muestra MATLAB en los números anteriores. Recordamos una vez
más que el cambio de formato solo afecta a la presentación que MATLAB hace de los números y no a los
números en sı́ ni a la precisión de la operaciones que se realizan con ellos.
Es necesario también ser consciente del siguiente tipo de situaciones:

>> 0.25+0.25-0.5
ans =
0
>> 0.2+0.1-0.3
ans =
5.5511e-17

Esta aparente paradoja ocurre porque, a diferencia de 0.25, el número 0.1 no es un número máquina —no es
un elemento del conjunto F de MATLAB— debido a que su representación en base 2 es infinita. De hecho, se
podrı́a poner una colección infinita de ejemplos de este tipo de errores de redondeo que se producen por el hecho
de utilizar aritmética de precisión finita. Uno más es

3
>> 1-3*(4/3-1)
ans =
2.2204e-16
que no da cero como deberı́a ser. Averigüe por qué se produce este último resultado e intente construir algún
otro ejemplo en el que aparezcan errores de redondeo.
Este tipo de errores debidos al redondeo pueden propagarse cuando el ordenador hace sucesivas operaciones
con aritmética de precisión finita —también llamada aritmética en coma flotante—, dando lugar a resultados
inexactos. Esto es, por ejemplo, lo que puede ocurrir al calcular la matriz inversa de una matriz dada. Al hacerlo,
MATLAB realiza una serie de operaciones —sumas, restas, multiplicaciones y divisiones— que pueden en algún
caso acumular errores de redondeo, lo que se conoce como inestabilidad numérica del algoritmo utilizado. Esto
se ilustra en el siguiente cálculo donde se utiliza una matriz de Hilbert. Recuerde que una matriz de Hilbert es
una matriz de la forma
1
A = (aij )ni,j=1 , aij = , i, j = 1, . . . , n.
1+i+j
y tiene siempre determinante no nulo.
>> A=hilb(5); % Matriz de Hilbert 5x5.
>> det(A) % A tiene determinante distinto de cero, si bien es peque~ no (cercano a eps).
ans =
3.7493e-12
>> Ai=inv(A); % Observe que el sı́mbolo ; hace que no se muestre el resultado por pantalla.
>> A*Ai % El resultado de esta operación deberı́a ser la matriz identidad 5x5.
ans =
1.0000e+00 -5.0369e-13 1.3571e-12 -6.8000e-12 5.2190e-12
-9.9585e-13 1.0000e+00 1.2944e-12 -2.8220e-12 2.9268e-12
-8.1704e-13 -6.2490e-14 1.0000e+00 -4.2378e-12 3.1583e-12
-7.1054e-13 -2.2737e-13 1.8190e-12 1.0000e+00 1.8190e-12
-6.3074e-13 2.7188e-14 1.4693e-12 -2.4876e-12 1.0000e+00
Se observa, sin embargo, que al multiplicar la matriz por su inversa, calculada con la aritmética de precisión
finita del ordenador, el resultado no coincide con la matriz identidad.
Por tanto, a la hora de programar algoritmos en un ordenador, dos de los problemas importantes a considerar
son, por un lado, disponer de criterios que informen a priori de la posible inestabilidad numérica del algoritmo
y, por otro, diseñar estos algoritmos de forma que se evite —si es posible— que estos errores de redondeo se
propaguen o acumulen de forma catastrófica.

2. Existencia y unicidad de la solución de un sistema lineal

Considere a partir de ahora sistemas de ecuaciones lineales de la forma

Ax = b, (2)

en los que A es una matriz real cuadrada n × n y el término independiente b es un vector de Rn .

A la hora de resolver el sistema lo primero es saber si admite solución y, en tal caso, si la solución es única.
Se disponen de varios criterios para resolver este problema que requieren del cálculo de ciertas magnitudes. Dos
condiciones equivalentes, que son necesarias y suficientes para la existencia y unicidad de solución del sistema
de ecuaciones (2), son:
El determinante de la matriz A es distinto de cero.
La matriz A es de rango completo; es decir, r(A) = n.
Recuerde que en el caso rectangular más general, es decir, si A ∈ Rm×n , entonces 0 ≤ r(A) ≤ mı́n(n, m). El
cálculo de determinantes siempre requiere muchas operaciones. Por tanto, no es adecuado utilizar este recurso
si se dispone de alguna alternativa, máxime si la matriz es grande y las operaciones se van a realizar con la
aritmética de precisión finita de una máquina. No obstante, MATLAB incorpora la función det que calcula
determinantes de matrices cuadradas utilizando la denominada factorización LU de la matriz, que se estudiará

4
más adelante en la asignatura —teclee help det en la lı́nea de comandos para obtener información sobre
esta función—. Se muestran a continuación algunos ejemplos del uso de det que debe teclear en su sesión de
MATLAB. En estos ejemplos interviene el comando rand, que sirve para generar números pseudoaleatorios
uniformemente distribuidos entre 0 y 1.

>> rng(10) % Establece una semilla para el generador de números aleatorios.

>> A=5*rand(10); % Genera una matriz 10x10 que se asigna a la variable A.
>> det(A) % Calcula el determinante de A que es distinto de cero.
ans =
1.6533e+05
>> B = hilb(30); % Genera una matriz de Hilbert 30x30 que se sabe es invertible.
>> det(B) % Calcula el determinante de B y da cero (inestabilidad numérica).
ans =
0

El último resultado pone de manifiesto la posible aparición de fenómenos de inestabilidad numérica provocados
por la acumulación de errores de redondeo que, también aquı́, producen un resultado que es erróneo. Conviene
mencionar que la función det produce un mensaje de error si se trata de aplicar a una matriz que no es cuadrada:

>> A=rand(3,4); % matriz 3x4

>> det(A)
Error using det
Matrix must be square.

Es importante destacar en este momento que MATLAB dispone de un conjunto muy completo y preciso de
mensajes de error para las funciones que contiene. A la hora de trabajar en este entorno de programación es
fundamental prestar atención y entender con detalle los mensajes que MATLAB produce.
A la vista de los problemas que presenta el cálculo de determinantes con aritmética de precisión finita,
para estudiar el problema de existencia y unicidad que nos ocupa resulta más conveniente utilizar el rango de
la matriz A. Ello es debido a que, en general, el número de operaciones necesario para calcular el rango es
menor que si se utiliza el determinante. MATLAB también incorpora una función que lo calcula para matrices
rectangulares cualesquiera. Se trata de la función rank —teclee help rank en la linea de comandos para obtener
información sobre esta función—. Se muestran a continuación algunos ejemplos del uso de rank que debe teclear
en su sesión de MATLAB.

>> A=rand(40,8); % Genera una matriz 40x8 y la asigna a la variable A.

>> rank(A)
ans =
8

¿Qué significa el valor del rango obtenido? Recuerde lo estudiado en clase de Álgebra.
>> B=rand(40); % Genera una matriz cuadrada 40x40 y la asigna a la variable B.
>> rank(A)
ans =
40

En este caso el rango de la matriz cuadrada B es máximo, por tanto, todo sistema de ecuaciones lineales de la
forma Bx = b, con b ∈ R40 , tiene solución única. Por otra parte, se puede comprobar la propiedad según la
cual el rango del producto de dos matrices es menor o igual que el mı́nimo de los rangos de cada una de ellas.
>> rank(B*A)
ans =
8

Sin embargo, también al calcular el rango numéricamente se pueden producir inestabilidades que den lugar a
resultados erróneos. Por ejemplo:

5
>> C=hilb(12); % Genera una matriz de Hilbert que (se sabe) es invertible.
>> rank(C) % El rango de la matriz no es máximo (inestabilidad numérica).
ans =
11
Si el resultado fuese correcto ¿Qué significarı́a el valor del rango obtenido en relación con las filas y columnas
de la matriz C?
Para terminar, insistir en una idea que ya se ha mencionado. Debido a la influencia de los errores de redondeo
en los resultados, es necesario disponer de criterios que nos permitan decidir si los resultados obtenidos de un
programa ejecutado en un ordenador son fidedignos o no.

3. Uso de rutinas o funciones en MATLAB: la función gauss

Hasta el momento nos hemos ocupado del trabajo con MATLAB en modo interactivo, utilizando la ventana
de comandos Command Window. Sin embargo, también es posible trabajar en modo programado utilizando para
ello los llamados M-files. Estos son ficheros de texto con extensión .m, que contienen conjuntos de instrucciones de
MATLAB y que se ejecutan desde la lı́nea de comandos. Aunque sobre estos ficheros y su manejo se profundizará
en la Práctica 2 de esta asignatura, daremos aquı́ alguna noción básica sobre su uso que nos permita trabajar
con la función gauss.
Hay dos tipos elementales de m-ficheros:
• los scripts, que están formados simplemente por una serie de instrucciones sucesivas que se ejecutan en
MATLAB una tras otra como si estuviéramos en modo interactivo;
• las functions (o m-funciones), que son el equivalente en MATLAB a las subrutinas o procedimientos de los
lenguajes de programación tradicionales.
Se trata aquı́ de ilustrar cómo se usa el segundo tipo de M-files de MATLAB trabajando con uno de ellos que
resuelve sistemas de ecuaciones lineales mediante eliminación gaussiana con pivotación parcial por columnas.
Para ello, se ha debido descargar de Moodle el correspondiente fichero gauss.m, que debe estar visible en su
directorio activo.
Haga doble click sobre este fichero. Se abrirá una ventana encima de la ventana de comandos de MATLAB
mostrando en el editor el contenido del fichero. Observe la primera lı́nea de texto de este fichero. Tiene la
estructura
function [variables de salida]=nombre_de_la_funcion(variables de entrada separadas por comas)
que es la que deben tener las funciones de usuario de MATLAB. En ella se especifica el nombre de la función,
la(s) variable(s) de entrada que necesitará el programa para realizar los cálculos y la(s) variable(s) de salida en
la que estarán los resultados. En el caso de la función gauss que vamos a manejar, esa primera lı́nea es:
function [x]=gauss(A,b)
indicando que la función se llama gauss y que, por tanto, el fichero en el que está contenida es gauss.m.
También indica que hay dos variables de entrada: la matriz A y el vector b, que definen el sistema de ecuaciones
Ax = b; y que la salida que produce el programa es el vector x, que contiene la solución del sistema. El resto
de instrucciones constituyen una de las posibles maneras que hay de programar en MATLAB el algoritmo de
Gauss con pivotación parcial por columnas. Su uso se ilustra a continuación mediante varios ejemplos que debe
reproducir en su sesión de MATLAB.
>> rng(10);
>> A=rand(5); % Genera una matriz 5x5 y la asigna a la variable A.
>> b=rand(5,1); % Genera un vector columna 5x1 y lo asigna a la variable b.
>> sol=gauss(A,b) % Resuelve el sistema Ax=b y asigna la solución a la variable sol.
sol =
-6.0635e+00
1.7016e+00
1.3320e+00
-3.3197e+00
1.0781e+01

6
Existe una forma obvia de comprobar que la solución obtenida es correcta: ver si el residuo Ax − b es el vector
nulo; en el ejemplo considerado:
>> A*sol - b
ans =
-3.8858e-16
-3.3307e-16
8.8818e-16
2.2204e-16
1.1102e-16
Observe que, aunque no es el vector nulo, todas las entradas aparecen multiplicadas por un factor 10−16 , por
lo que todas ellas tienen 15 posiciones decimales nulas. La solución truncada a 15 cifras decimales es entonces
numéricamente correcta.
Otra manera de medir la proximidad a cero de la diferencia Ax − b es obtener su norma 2 que MATLAB es
capaz de calcular mediante la función norm. El resultado es:
>> norm(A*sol-b)
ans =
1.0547e-15
Es decir, la diferencia en norma es del orden de 10−15 . Este valor también nos permite concluir que la solución
es numéricamente aceptable, pese a que no es cero.
Considere ahora otro ejemplo:
>> A=hilb(10);
>> rank(A)
ans =
10
>> rng(10);
>> b=rand(10,1);
>> sol=gauss(A,b);
>> A*sol-b
ans =
1.9147e-05
2.1675e-05
7.0020e-06
3.5040e-07
-1.0270e-05
-9.9493e-06
1.3383e-05
-9.6849e-07
-4.2427e-06
-2.1627e-06
>> norm(A*sol-b)
ans =
3.5954e-05
Los valores obtenidos indican que la solución es mucho menos precisa que en el ejemplo anterior, debido a la
inestabilidad numérica, y ello a pesar de que el rango de la matriz A calculado con la función rank es 10.
En cualquier caso, habrá que ser muy cuidadosos si utilizamos esta solución para sacar conclusiones en algún
problema concreto.
Finalmente queremos comentar que no se ha hecho referencia a la regla de Cramer como método fiable para
resolver sistemas de ecuaciones lineales debido al elevadı́simo número de operaciones aritméticas que requiere
su aplicación. Baste decir, a tı́tulo meramente informativo, que para resolver un sistema de 10 ecuaciones con 10
incógnitas, la regla de Cramer requiere aproximadamente unos 500 millones de operaciones elementales (suma,
restas, multiplicaciones y divisiones), mientras que la reducción gaussiana que hemos usado aquı́ solamente
necesita 740.

7
4. La barra invertida de MATLAB
Por supuesto MATLAB incorpora una variada colección de funciones que permiten resolver sistemas de
ecuaciones lineales utilizando algoritmos diversos. La que se utiliza más a menudo es la barra invertida ✭✭\✮✮
que, en el caso más general cuando la matriz de coeficientes no tiene ninguna propiedad especial, hace uso de
la factorización LU , método este último que, en general, es más eficiente que la reducción gaussiana empleada
por la función gauss introducida en la sección anterior.
Una manera de medir la eficiencia de un algoritmo es medir el tiempo de computación que le lleva a MATLAB
ejecutarlo. La pareja de comandos tic, toc mide el tiempo de CPU empleado por MATLAB entre uno y otro
comando. Por ejemplo, las instrucciones

>> rng(1)
>> A= rand(100);
>> tic; det(A); toc;
Elapsed time is 0.007566 seconds.

indican que MATLAB ha empleado 0,007566 segundos en calcular el determinante de la matriz de entradas
aleatorias A, que tiene dimensiones 100 × 100. Compare con el siguiente ejemplo análogo
>> rng(1)
>> A= rand(10000);
>> tic; det(A); toc;
Elapsed time is 10.734310 seconds.

Estas consideraciones de eficiencia son esenciales cuando el cálculo que se aborda es muy complejo o se manejan
matrices grandes.
Es importante escribir los comandos tic, toc en la misma lı́nea —y separados por punto y coma— que
la instrucción cuyo tiempo de computación se quiere calcular. De otra manera, el resultado tendrá en cuenta
el tiempo que se ha tardado en realizar todas las tareas intermedias, como, por ejemplo, pulsar las teclas
para escribir los comandos. Por otro lado, es claro que los valores concretos que tome la pareja de funciones
tic, toc dependen en gran medida de la máquina que se esté utilizando e, incluso, de si la computadora está
realizando en ese preciso momento otras tareas. Por eso no se extrañe si los valores que obtiene en su sesión de
MATLAB difieren de los aquı́ presentados. Teclee help tic y help toc para obtener más información sobre
estas funciones.
Podemos ahora utilizar esta herramienta para comparar la eficiencia de la barra invertida de MATLAB y la
función gauss.
>> rng(1)
>> A=rand(1000);
>> b=rand(1000,1);
>> tic; gauss(A,b); toc;
Elapsed time is 1.622426 seconds.
>> tic; A\b; toc;
Elapsed time is 0.019318 seconds.
En este ejemplo la barra invertida de MATLAB es del orden de 100 veces más rápida que la función gauss. La
diferencia tiende a crecer cuando las dimensiones de la matriz son todavı́a mayores. Experimente con diferentes
dimensiones de la matriz A. ¿Puede encontrar alguna correlación entre el orden de A y la diferencia de tiempo
existente entre la barra invertida de MATLAB y gauss?
Para finalizar esta sección, considere el siguiente ejemplo:

>> C=[1 -6 7 -9; 1 -5 0 0; 0 1 -5 0; 0 0 1 -5]

C =
1 -6 7 -9
1 -5 0 0
0 1 -5 0
0 0 1 -5
>> det(C)

8
ans =
-1
>> A=C^7;
>> rng(1)
>> b=rand(4,1);
>> sol=A\b;
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.151102e-18.

>> A*sol-b
ans =
-1.1671e+00
5.9360e-01
-2.5116e-01
-8.7702e-01
>> norm(A*sol-b)
ans =
1.5958e+00
Claramente, la solución obtenida no es aceptable. Observe el warning que da MATLAB en el que avisa de que
la matriz A presenta una estructura problemática desde el punto de vista numérico. De hecho, si se calcula
su rango el valor que se obtiene es 3 —compruébelo en su sesión de MATLAB—. Esto significa que, siendo
todas las filas de esta matriz linealmente independientes —ya que su determinante es igual a −1—, hay una
de ellas que aproximadamente es combinación lineal de las restantes, lo que induce errores de redondeo y
cancelaciones que hacen inaceptable la solución numérica que se obtiene. Como ya se indicó al final del documento
3_Introduccion_MATLAB.pdf, en estos casos en que los resultados son totalmente erróneos, es posible que los
valores numéricos que obtenga en su sesión de MATLAB sean diferentes de los que aquı́ se presentan.
Matrices con este comportamiento se dice que están mal condicionadas. Observe que en el warning que da
MATLAB aparece un número, rcond, que permite dar una estimación de este mal condicionamiento numérico.
En la siguiente sección se dan unas breves pinceladas sobre esta manera de medir el comportamiento de una
matriz frente a cálculos numéricos. Pero antes de pasar a ella, practique con la barra invertida en su sesión de
MATLAB usando los ejemplos que se han considerado en la sección anterior.

5. Número de condición de una matriz

Considere atentamente los dos sistemas de ecuaciones
" "
x+y =2 x+y =2
x + 1,001y = 2 x + 1,001y = 2,001

El sistema de la izquierda tiene solución x = 2, y = 0, mientras que la solución del de la derecha es x = y = 1.

Ası́, un cambio pequeño en los coeficientes del término independiente provoca un cambio mucho mayor en los
resultados del sistema de ecuaciones. La matriz de coeficientes
# $
1 1
A=
1 1,001

resulta tener entonces una sensibilidad excepcional a pequeñas perturbaciones, que es lo que la hace difı́cil de
manejar numéricamente, pues no son otra cosa que pequeñas perturbaciones los redondeos que realiza MATLAB.
Si A es una matriz de columnas independientes, se define su número de condición espectral —hay otros—
como
σmax
c2 (A) = ,
σmin
donde σmax , σmin representan, respectivamente, el mayor y el menor valor singular de A. Los valores singulares
de A son las raı́ces cuadradas de los autovalores de la matriz AT A, que es simétrica y definida positiva. Si A
posee columnas dependientes, se dice que su número de condición es infinito.
Una matriz se dice que está bien condicionada si su número de condición está cerca de uno y diremos que está
mal condicionada cuando ese número sea significativamente mayor que uno. La función cond es la que calcula el

9
número de condición en MATLAB —que también posee una forma aproximada de hacerlo mediante la función
rcond, la cual da una estimación del inverso del número de condición—. Observe las siguientes instrucciones:
>> B=[1 1;1 1.001]
B =
1.0000e+00 1.0000e+00
1.0000e+00 1.0010e+00
>> cond(B)
ans =
4.0020e+03
>> rcond(B)
ans =
2.4975e-04
Teclee help cond y help rcond para obtener más información sobre estas funciones.
La importancia del número de condición c2 (A) en la resolución numérica del sistema lineal Ax = b proviene
de la desigualdad fundamental
%x − x′ % %b − b′ %
≤ c2 (A) , (3)
%x% %b%
donde b′ representa una perturbación del dato b y x′ es la solución del sistema lineal correspondiente. La
desigualdad (3) establece que el error relativo en la solución es menor o igual que el error relativo en los datos
multiplicado por el número de condición de la matriz A. Dicho en otras palabras, el número de condición mide
la máxima amplificación que puede sufrir el error relativo en la solución respecto al error relativo en los datos.
Más adelante en la asignatura de Álgebra se tratarán con detalle estas cuestiones.
Los razonamientos anteriores se trasladan igualmente al problema de averiguar si un vector de residuos
Ax′ − b = b′ − b pequeño implica que x′ es una solución aproximada precisa del sistema de ecuaciones lineales
Ax = b. La respuesta viene dada por la desigualdad (3). Por ejemplo,
>> C=[17 -864 716 -799; 1 -50 0 0; 0 1 -50 0; 0 0 1 -50];
>> A=C^3
A =
19453 -989027 836017 -858925
1075 -52572 -60227 66317
-83 6636 -124284 -799
1 -150 7500 -125000
>> b=sum(A’)’
b =
-992482
-45407
-118530
-117649
>> sol=A\b
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.791266e-16.
sol =
4.3945e+01
1.8589e+00
1.0172e+00
1.0003e+00
>> A*sol-b
ans =
1.1642e-10
-7.2760e-12
-1.4552e-11
1.4552e-11
>> cond(A)
ans =
4.0521e+18

10
Observe que el vector b es la suma de las columnas de la matriz A, de manera que la solución exacta x del
sistema Ax = b es el vector (1, 1, 1, 1), que tiene un tamaño del orden de 100 . El residuo Ax′ − b es del orden de
10−10 y el orden de b es aproximadamente 107 ; sin embargo, la solución numérica x′ difiere del verdadero valor
x en una magnitud del orden de 101 . Esta discrepancia tan grande entre los errores relativos de la solución y
de los datos es posible debido al elevado número de condición de la matriz A.
El mal condicionamiento debe analizarse desde el punto de vista de la precisión de la aritmética empleada
en el cálculo. Si se utiliza una aritmética de t dı́gitos (16 en MATLAB), los errores de redondeo serán del orden
de 10−t . En esta situación, si el número de condición es del orden de 10s , entonces la solución proporciona,
grosso modo, t − s dı́gitos fiables, es decir, se pierden s cifras significativas. Esto es compatible con el hecho de
que en algún cálculo concreto el número de dı́gitos correctos puede ser mayor. Esto se debe principalmente a dos
razones: por un lado, el número de condición representa una cota superior del error relativo, pudiendo de hecho
ser menor; y por otro, del mismo modo que los errores de redondeo pueden sumarse unos a otros y propagarse,
es posible que se cancelen unos con otros dando lugar a una solución más precisa de lo que se esperaba.
Para terminar, use la función cond en los ejemplos considerados en las secciones 3 y 4 y, a la vista de lo que
se acaba de comentar, decida en cada caso el número de dı́gitos de la solución que se puede considerar fiable.

6. Cálculo simbólico
En ocasiones puede ser útil realizar cálculos simbólicos en MATLAB, es decir, cálculos que se realizan de
manera exacta o con precisión infinita. La mayor parte de las funciones simbólicas de MATLAB tienen el mismo
nombre que sus análogas numéricas; MATLAB selecciona en cada caso la función adecuada dependiendo de cómo
se hayan definido las variables involucradas en el cálculo: si la variable es numérica —lo cual aparece detallado
en la ventana Workspace bajo la columna Value—, entonces el cálculo será numérico; si, por el contrario, la
variable es simbólica el cálculo se efectuará de modo simbólico. Por defecto, las variables en MATLAB son
numéricas; una variable se define como simbólica mediante el comando sym. Por ejemplo,
>> a=log(2) % Variable numérica.
a =
6.9315e-01
>> b=sym(a) % Variable simbólica definida a partir del número de máquina anterior.
b =
6243314768165359/9007199254740992
>> c=log(sym(2)) % Variable simbólica.
c =
log(2)
>> d=double(c). % Transforma la variable simbólica en numérica.
d =
6.9315e-01

La primera variable es numérica —observe la ventana Workspace—; las otras dos siguientes son simbólicas, si
bien la segunda variable se ha formado a partir de un número en aritmética finita y no puede representar el
valor exacto de log 2, mientras que la tercera variable es puramente simbólica. La función double transforma
una variable simbólica en una numérica escrita en doble precisión. Teclee help sym y help double para mayor
información.
El cálculo simbólico permite en ocasiones descubrir los errores numéricos. Retomemos el ejemplo de la sección
anterior:

>> C=[17 -864 716 -799; 1 -50 0 0; 0 1 -50 0; 0 0 1 -50];

>> A=C^3; % Matriz numérica. Las entradas de A son números enteros.
>> rank(A) % Cálculo numérico del rango de A.
ans =
3
>> det(A) % Cálculo numérico del determinante de A.
ans =
1.7001e+03 % Resultados contradictorios entre rank(A) y det(A).
>> As=sym(A); % Matriz simbólica, no se pierde precisión al ser A entera.

11
>> rank(As) % Cálculo simbólico del rango de A.
ans =
4
>> det(As) % Cálculo simbólico del determinante de A.
ans =
-1
Si el cálculo simbólico proporciona respuestas correctas, ¿por qué no recurrir siempre a él? La respuesta es que el
cálculo simbólico da resultados correctos siempre que sea capaz de realizar los cálculos. En el caso de algoritmos
complejos o matrices grandes el procedimiento simbólico es muy ineficiente o, simplemente, incapaz de llegar
al resultado, por lo que su uso en las aplicaciones industriales es tremendamente limitado. Tenga cuidado al
asignar variables simbólicas: si las variables son demasiado grandes, MATLAB puede ✭✭quedarse colgado✮✮. A
este respecto, teclee ctrl+c (sistema operativo Windows) o cmd+. (sistema operativo Mac Os) en el caso de
que necesite abortar una operación. Para finalizar la práctica, considere el siguiente ejemplo:

>> rng(1)
>> A=rand(100);
>> tic; d1=det(A); toc;
Elapsed time is 0.012221 seconds.
>> As=sym(A);
>> tic; d2=det(As); toc;
Elapsed time is 22.171493 seconds.
>> d1
d1 =
4.3544e+24
>> double(d2)
ans =
4.3544e+24

Observe la gran diferencia que existe en el tiempo de ejecución entre el procedimiento simbólico y el numérico.
Experimente —con cuidado— qué ocurre si, en los cálculos anteriores, se cambia el tamaño de la matriz.

COMENTARIO FINAL. Las cuestiones que se han abordado en este guion relativas a la estabilidad numérica
no son exclusivas de MATLAB, sino que aparecen al usar cualquier lenguaje o procedimiento de cálculo numérico.
Por otra parte, los problemas de mal condicionamiento no ocurren siempre, de hecho, la mayorı́a de los sistemas
lineales que surgen en las aplicaciones industriales resultan ser bien condicionados. Sin embargo, es necesario
tener presente que algunos problemas importantes del cálculo numérico, tales como el cálculo de autovalores,
dan lugar —si se atacan directamente— a sistemas mal condicionados, por lo que en la práctica los algoritmos
han de ser diseñados de manera que sorteen esa dificultad. Esa es la razón por la que en este guion se ha querido
poner de manifiesto la necesidad de analizar la estabilidad de los algoritmos y de disponer de criterios para
evaluar lo precisas que son las soluciones cuando estas se obtienen utilizando la aritmética de precisión finita
de un ordenador.

12