Apuntes HI-1

Herramientas Informáticas
Paco Arandiga
Índice general
1. Algoritmos básicos en MATLAB 3

1.1. Programación en MATLAB . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2. Introducción de matrices . . . . . . . . . . . . . . . . . . . . 4
1.1.3. Operaciones con matrices, operaciones a coordenadas . . . . 5
1.1.4. Declaraciones, expresiones y variables; almacenamiento de una
sesión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.5. Formato de salida . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.6. Funciones en MATLAB . . . . . . . . . . . . . . . . . . . . 9
1.1.7. For, while, if — y relaciones . . . . . . . . . . . . . . . . . . 11
1.1.8. Submatrices y notación de dos puntos . . . . . . . . . . . . . 13
1.1.9. Comparación de la eficiencia de algoritmos: etime, cputime,
tic, toc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.1.10. Cadenas de texto, mensajes de error, input . . . . . . . . . . . 15
1.1.11. Archivos .m . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.1.12. Gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.1.13. Consulta . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2. Cálculo Matricial Básico . . . . . . . . . . . . . . . . . . . . . . . . 32
1.2.1. Normas de vectores . . . . . . . . . . . . . . . . . . . . . . . 32
1.2.2. Tipos de matrices . . . . . . . . . . . . . . . . . . . . . . . 34
1.2.3. Matrices simétricas. . . . . . . . . . . . . . . . . . . . . . . 37
2. Solución de ecuaciones no lineales 38

2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2. Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2.1. Tipos de convergencia . . . . . . . . . . . . . . . . . . . . . 40
2.3. Método de la Bisección . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.4. Método de régula falsi . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.5. Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . 47
2.6. Método de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.7. Raíces múltiples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.8. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.9. Problemas propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . 57
1
ÍNDICE GENERAL 2
3. Cálculo simbólico 59
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.2. Definiendo variables simbólicas . . . . . . . . . . . . . . . . . . . . 61
3.3. Expresiones y ecuaciones simbólicas . . . . . . . . . . . . . . . . . . 62
3.3.1. Substituciones . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3.2. Qué variables simbólicas hay en F? . . . . . . . . . . . . . . 64
3.4. Resolución de ecuaciones . . . . . . . . . . . . . . . . . . . . . . . . 65
3.5. La función inline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.5.1. Operadores relacionales y lógicos . . . . . . . . . . . . . . . 68
3.5.2. La función inline . . . . . . . . . . . . . . . . . . . . . . . . 69
3.6. Dominios, límites y continuidad . . . . . . . . . . . . . . . . . . . . 72
3.6.1. Cálculo de límites: Límites laterales. . . . . . . . . . . . . . . 74
3.6.2. Continuidad de funciones . . . . . . . . . . . . . . . . . . . 76
3.7. Simplificadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.8. Cálculo de derivadas. . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.9. Aplicaciones de la derivada: Máximos-mínimos-Puntos de Inflexión.
Optimización. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.10. Cálculo integral e integración numérica . . . . . . . . . . . . . . . . 93
4. Sistemas numéricos y fuentes de error 96

4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2. Representaciones numéricas . . . . . . . . . . . . . . . . . . . . . . 97
4.3. Aritmética de punto flotante . . . . . . . . . . . . . . . . . . . . . . 101
4.4. Consecuencias de la aritmética de precisión finita . . . . . . . . . . . 108
4.5. Inestabilidad numérica . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.6. Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Capítulo 1
Algoritmos básicos en MATLAB
3
1.1. PROGRAMACIÓN EN MATLAB 4
1.1. Programación en MATLAB

1.1.1. Introducción
MATLAB es un sistema interactivo basado en matrices para cálculos científicos
y de ingeniería. Se pueden resolver problemas numéricos relativamente complejos sin
escribir un programa en realidad. El nombre MATLAB es una abreviatura para MATrix
LABoratory.
El propósito de estas notas es ayudar en la iniciación a MATLAB. La mejor forma
de utilizarlas es poner manos a la obra. Se aconseja, en general, trabajar en el ordenador
a la vez que se leen las notas, así como a experimentar libremente con ejemplos.
Se puede utilizar la ayuda de la instrucción help para una información más detal-
lada. Después de entrar en MATLAB, la instrucción help mostrará una lista de fun-
ciones para las que se puede obtener ayuda mientras se está trabajando; la instrucción
help nombre_de_función nos dará información sobre una función específica. Así, la
instrucción help eig, nos dará información sobre la función eig, que calcula los
autovalores de una matriz. Se pueden ver algunas de las capacidades de MATLAB
usando la instrucción demo.
El alcance y la potencia de MATLAB van más allá de lo que podemos ver en estas
notas. En algún momento puede desear una información más detallada. Es el momento
de consultar algún manual más avanzado.
1.1.2. Introducción de matrices

MATLAB trabaja esencialmente con un solo tipo de objetos: una matriz numérica
rectangular con entradas posiblemente complejas; todas las variables representan ma-
trices. A veces, las matrices 1 × 1 se consideran escalares, y las matrices con una sola
fila o columna se consideran como vectores.
Hay varias formas diferentes para introducir una matriz en MATLAB. A saber:
Introduciendo una lista explícita de elementos,
Generándola mediante funciones y declaraciones,
Creándola en un archivo .m (ver secciones 1.1.8 y 1.1.10)
Cargándola de un archivo de datos externo (lo veremos más adelante).
Por ejemplo, cualquiera de las declaraciones
A = [1 2 3; 4 5 6; 7 8 9]
y
A = [
1 2 3
4 5 6
7 8 9 ]
crea la matriz 3×3 que se espera y la asigna a una variable A. Inténtelo. Los elementos
en una fila de una matriz pueden separarse tanto por comas como por espacios en
blanco.
Cuando alguno de los números se escribe en forma exponencial (por ejemplo
2.34e-9), deben evitarse los espacios en blanco. La escritura de una matriz grande
debe hacerse preferentemente en un archivo .m, donde es más sencillo corregir errores
(ver secciones 1.1.8 y 1.1.10)).
Las funciones internas rand, magic, y hilb, por ejemplo, proporcionan una for-
ma sencilla para crear matrices con las que experimentar. La instrucción rand(n),
resp. rand(m,n), creará una matriz n × n, resp. m × n, con entradas aleatoria-
mente generadas, distribuidas uniformemente entre 0 y 1. magic(n) creará una ma-
triz cuadrada mágica (las filas y las columnas suman la misma cantidad) con entradas
enteras; hilb(n) creará la matriz de Hilbert de orden n, la reina de las matrices mal
condicionadas. m y n, por supuesto, denotan enteros positivos. También se pueden
crear matrices utilizando bucles for (ver sección 1.1.6). Inténtelo.
Las entradas individuales de una matriz o de un vector se pueden obtener poniendo
los índices entre paréntesis de la forma usual. Por ejemplo, A(2, 3) denota la entrada en
la segunda fila y tercera columna de la matriz A y x(3) denota la tercera coordenada
del vector x. Inténtelo. Sólo se pueden usar como índices de vectores y de matrices
enteros positivos.
1.1.3. Operaciones con matrices, operaciones a coordenadas

Disponemos en MATLAB de las siguientes operaciones con matrices:
+ adición
− sustracción
∗ multiplicación
ˆ potenciación
0
traspuesta
\ división izquierda
/ división derecha
Estas operaciones para matrices se aplican también a escalares (matrices 1×1). Si los
tamaños de las matrices son incompatibles para la operación matricial se obtiene un
mensaje de error, exceptuando el caso en que uno de los operandos sea un escalar
y el otro una matriz (para la adición, sustracción, división y multiplicación). En esta
situación se opera el escalar con cada término de la matriz.
La división matricial merece un comentario especial. Si A es una matriz invertible

y b es una columna, resp. fila, compatible, entonces
>> x=A\b
es la solución de A ∗ x = b y, resp.,
>> x=A/b
es la solución de x ∗ A = b.
En la división izquierda, si A es cuadrada, se factoriza utilizando eliminación gaus-
siana. Con los factores se resuelve A ∗ x = b. Si la matriz A no es cuadrada, se fac-
toriza utilizando la ortogonalización de Householder con pivoteo de columnas. Con
los factores se resuelve el sistema indeterminado o sobredeterminado en el sentido
de los mínimos cuadrados. La división derecha se define a partir de la izquierda por
b/A = (A0 \b0 )0 .
Operaciones a coordenadas. Las operaciones de adición y sustracción operan
intrínsecamente a coordenadas pero las otras operaciones matriciales dadas antes no:
Son operaciones matriciales. Es importante observar que para estas otras operaciones,
∗, ˆ, \, y /, puede hacerse que operen a coordenadas precediendolas de un punto. Por
ejemplo, tanto
>> [1,2,3,4].*[1,2,3,4]
como
>> [1,2,3,4].^2
darán
>> [1,4,9,16]
Inténtelo. Esto es particularmente útil cuando se utilizan los gráficos de MATLAB.
1.1.4. Declaraciones, expresiones y variables; almacenamiento de

una sesión
MATLAB es un lenguaje de expresiones; las expresiones que se escriben son inter-
pretadas y evaluadas. Las instrucciones de MATLAB son, normalmente, de la forma
variable = expresión, o simplemente
expresión
Las expresiones se componen, usualmente, a partir de operadores, funciones y
nombres de variables. La evaluación de una expresión produce una matriz, que se
muestra en pantalla, y se asigna a la variable para su posterior uso. Si se omiten la

variable y el signo =, se crea una variable llamada ans (por answer) a la que se asigna
el resultado de la expresión.
Una instrucción termina, normalmente, con el retorno de carro (comando return).
Si se desea continuar una expresión en la línea siguiente, basta escribir tres (o más)
puntos antes del retorno de carro. Si por el contrario, deseamos escribir varias in-
strucciones en una misma línea, podemos hacerlo separandolas por comas o puntos y
comas.
Si el último carácter de una instrucción es un punto y coma el resultado no se
mostrará en pantalla, aunque por supuesto se realizará la asignación. Esto es esencial
para evitar pérdidas de tiempo al mostrar los resultados intermedios.
MATLAB distingue las letras mayúsculas de las minúsculas en los nombres de in-
strucciones, funciones y variables. Así, resolvente no es lo mismo que ReSoLvEnTe.
La instrucción who muestra las variables que se encuentran en el espacio de traba-
jo. Por ejemplo, Con who se lista las variables que tenemos. Per ejemplo,
>> who
Your variables are:
A ans
Para ver el tamaño de las variables hemos de escribir whos y aparecerá:
Name Size Bytes Class
A 2x2 32 double
ans 1x1 8 double
Cada elemento de una matriz real requiere 8 bytes de memoria.
Escribiendo what aparecerá
M-files in the current directory /home/paco/manual/matlab
obj plotf1 pn1 prod1 prod3 tra

afun peli plotf2 pn2 prod2 prod4 trape
esto es, todas *.m files (que ya explicaremos más adelante en que consisten) que ten-
emos en el directorio en el que estamos trabajando. Hay otros comandos que tambien
se pueden ejecutar directamente desde MATLAB como cd, ls, mkdir, etc.
Para eliminar una variable de la memoria se utiliza la instrucción clear nom-
bre_variable. Si se escribe sólo clear se borran todas las variables no permanentes.
La variable permanente eps (épsilon) da la precisión de la máquina—alrededor de
−16
10 en la mayoría de ellas. Es útil para determinar la tolerancia en procesos itera-
tivos.
Cualquier tipo de cálculo, gráfico, o impresión puede detenerse sin salir del pro-
grama con CTRL-C (CTRL-BREAK en PC).
Almacenamiento de sesiones. Cuando salimos de MATLAB se pierden todas las
variables. Para evitarlo se puede utilizar la instrucción save antes de salir. Esto hace
que las variables se almacenen en el archivo de disco matlab.mat. Al acceder de
nuevo a MATLAB, se pueden recuperar todas las variables con la instrucción load.
Con diary se pueden salvar los datos en un fichero que despues se puede editar.
Habria que escribir
>> diary nom
>>
.
.
.
>>
>>diary off
y todo lo que hace entre diary nom y diary off queda archivado en un fichero
ASCII llamado nom.
Al terminar, se puede editar el archivo como se desee.
Con save name podemos salvar datos en un fichero binario llamado name.dat.
Los datos almacenados aquí pueden ser recuperados con load. Si queremos salvar los
datos en formato ASCII lo podemos hacer de la forma siguiente: manera:
>> A=rand(4,2);
>> save enter.dat A -ascii
que crea un fichero llamado enter.dat que contiene
3.4119357e-01 7.2711322e-01 8.3849604e-01 3.7041356e-01
5.3407902e-01 3.0929016e-01 5.6807246e-01 7.0273991e-01
1.1.5. Formato de salida

Aunque todos los cálculos en MATLAB se efectúan en doble precisión, el formato
de la salida en pantalla puede ser controlado con las siguientes instrucciones.
format short coma fija con 4 decimales (el defecto)
format long coma fija con 14 decimales
format short e notación científica con 4 decimales
format long e notación científica con 15 decimales
Una vez que se ordena un formato, se mantiene hasta que se ordena un cambio.
La orden format compact evitará la mayor parte de las líneas en blanco, con
lo que se puede mostrar más información en pantalla. Es independiente de las demás
instrucciones de formato.
Ejercicio 1.1.1. Calcular A=rand(4,2) y ver como queda en cada uno de los for-
matos anteriores.
1.1.6. Funciones en MATLAB

Funciones para la construcción de matrices
Las siguientes funciones están disponibles en MATLAB:
eye matriz identidad

zeros matriz de ceros
ones matriz de unos
diag ver más adelante
triu parte triangular superior de una matriz
tril parte triangular inferior de una matriz
rand matriz generada aleatoriamente
hilb matriz de Hilbert
magic matriz mágica
toeplitz ver help toeplitz
Por ejemplo, zeros(m,n) produce una matriz nula m×n, y zeros(n) produce
otra cuadrada de orden n; si A es una matriz, entonces zeros(A) produce una matriz
de ceros del mismo orden que A.
Si x es un vector, diag(x) es la matriz diagonal con x en su diagonal; si A es
una matriz cuadrada, diag(A) es un vector formado por la diagonal de A. ¿Qué será
entonces diag(diag(A))? Inténtelo.
Las matrices se pueden construir por bloques. Por ejemplo, si A es 3 × 3, entonces
>> B = [A, zeros(3,2); zeros(2,3), eye(2)]
dará una cierta matriz 5 × 5. Inténtelo.
Funciones escalares
Algunas funciones de MATLAB operan esencialmente sobre escalares, aunque lo

hacen también sobre matrices (elemento a elemento). Las funciones más comunes en-
tre estas son:
sin asin exp abs round

cos acos log(natural) sqrt f loor
tan atan rem(resto) sign ceil
Funciones vectoriales
Otras funciones de MATLAB operan fundamentalmente sobre vectores (fila o colum-
na), aunque también pueden operar sobre matrices m × n (m ≥ 2) haciendolo en este
caso columna a columna, produciendo, por tanto, un vector fila que contiene el resul-
tado de su aplicación a cada columna. Para conseguir que actúen por filas basta usar la
traspuesta; por ejemplo, mean(A’)’. Veamos algunas de estas funciones:
max sum median any

min prod mean all
sort std
Por ejemplo, la entrada máxima de un matriz A se obtiene con max(max(A)) en vez

de max(A). Inténtelo.
Funciones matriciales
Las funciones matriciales más útiles de MATLAB son las siguientes:
size tamaño
eig autovalores y autovectores
chol factorización de Cholesky
svd descomposición en valores singulares
inv inversa
lu factorización LU
qr factorización QR
hess forma de Hessenberg
schur descomposición de Schur
rref forma escalonada reducida por filas
expm matriz exponencial
sqrtm matriz raíz cuadrada
poly polinomio característico
det determinante
norm norma 1, norma 2, norma de Frobenius, norma ∞
cond número de condición en la norma 2
rank rango
Las funciones de MATLAB admiten argumentos de salida simples o múltiples. Por

ejemplo,
y = eig(A) , o simplemente eig(A)
genera un vector columna conteniendo los autovalores de A mientras que
[U,D] = eig(A)
produce una matriz U cuyas columnas son los autovectores de A y una matriz diagonal
D con los autovalores de A en su diagonal. Pruebe.
1.1.7. For, while, if — y relaciones

Básicamente, las instrucciones para el control de flujo de MATLAB operan como
en la mayor parte de los lenguajes usuales.
for. Por ejemplo, las instrucciones
>> x = []; for i = 1:n, x=[x,i^2], end
>> x = [];
>> for i = 1:n
>> x = [x, i^2]
>> end
darán como resultado un cierto vector mientras que
>> x = []; for i = n:-1:1, x=[x,i^2], end
dará el mismo vector en orden inverso. Pruebe con ellas. Las instrucciones
>> for i = 1:m

>> for j = 1:n
>> H(i, j) = 1/(i+j-1);
>> end
>> end
>> H
producirán e imprimirán en pantalla la matriz de Hilbert m × n. El punto y coma de

la instrucción interior suprime la impresión no deseada de los resultados intermedios
mientras que el último H muestra el resultado final.
while. La forma general de un bucle while es
while relación
instrucciones
end
Las instrucciones se repetirán mientras la relación sea cierta. Por ejemplo, dado un
número a, las instrucciones siguientes calculan y muestran el menor entero no negativo
n tal que 2n ≥ a:
>> n = 0;
>> while 2^n < a
>> n = n + 1;
>> end
>> n
if. La forma general de un bucle if simple es

if relación
instrucciones
end
Las instrucciones se ejecutarán sólo si la relación es cierta. También son posibles las
ramificaciones múltiples, como se ilustra con
>> if n < 0
>> paridad = 0;
>> elseif rem(n,2) == 0
>> paridad = 2;
>> else
>> paridad = 1;
>> end
Si sólo tenemos dos ramificaciones podemos omitir, desde luego, la porción correspon-
diente a elseif.
Relaciones. Los operadores relacionales en MATLAB son
< menor que
> mayor que
<= menor o igual que
>= mayor o igual que
== igual
∼= no igual.
Hagamos notar que se usa “=” en las asignaciones mientras que para las relaciones se
usa “==”. Las relaciones pueden conectarse o cuantificarse por los operadores lógicos
& y
| o
∼ no.
Cuando se aplican a escalares, una relación es realmente el escalar 1 o 0 dependi-

endo de si la relación es verdadera o falsa:
Ejercicio 1.1.2. Ver los resultados obtenidos con
>> 3 < 5, 3 > 5, 3 == 5, 3 == 3
Cuando se aplica a matrices del mismo orden, una relación entre ellas da lugar a
una matriz de ceros y unos, dando el valor de la relación entre las correspondientes
entradas.
Ejercicio 1.1.3. Ver el resultado obtenido con

>> a = rand(5), b = triu(a), a == b.
Cuando se utiliza una relación entre matrices en un bucle while o if, la relación
se entiende verdadera si cada una de las entradas de la matriz de relación es no nula.
Por tanto, si se quiere ejecutar algo cuando las matrices A y B sean iguales, se puede
escribir:
if A == B
algo
end
pero si se desea ejecutar la misma instrucción si A y B son distintas, hay que recurrir
a:
if any(any(A ~= B))
algo
end
o, más simplemente,
if A == B else
algo
end
Recalcamos que la aparentemente obvia
if A ~= B, algo, end
no hará lo que deseamos ya que la instrucción sólo se ejecutará si todas las entradas
de A son distintas de las de B. Las funciones any y all pueden utilizarse de forma
creativa para reducir relaciones entre matrices a relaciones entre vectores y escalares.
Se requieren dos anys en el ejemplo anterior ya que any es un operador vectorial (ver
sección 1.1.6).
1.1.8. Submatrices y notación de dos puntos

Los vectores y submatrices son utilizados a menudo en MATLAB para conseguir
efectos de manipulación bastante complejos. La notación de dos puntos (que se utiliza
para generar vectores y submatrices), y la indexación por vectores son las llaves para
una manipulación eficiente de estos objetos. Su uso de forma creativa permite mini-
mizar el número de bucles (que enlentecen a MATLAB) y hacen que las instrucciones
sean más simples y legibles. Debe hacerse un esfuerzo especial para familiarizarse
con esta notación.
La expresión 1:5 (que ya encontramos en los bucles for) es realmente un vector
fila: el [1 2 3 4 5]. Los números no tienen que ser enteros ni el incremento uno.
Por ejemplo,
>> 0.2:0.2:1.2
da como resultado
>> [0.2 0.4 0.6 0.8 1.0 1.2]
mientras que con
>> 5:-1:1
se obtiene el vector
>> [5 4 3 2 1]
Las siguientes instrucciones, por ejemplo, generarán una tabla de senos. Pruebe.
>> x = [0.0:0.1:2.0]’;
>> y = sin(x);
>> [x y]
Hagamos notar que al operar sin a coordenadas, produce un vector y a partir de x.
La notación de dos puntos permite acceder a submatrices. Por ejemplo, A(1:4,3)
es el vector columna con las cuatro primeras entradas de la tercera columna de A.
Dos puntos sin más especificación denotan una fila o columna completa: A(:,3)
es la tercera columna de A, y A(1:4,:) son las cuatro primeras filas.
Se pueden usar como índices vectores enteros arbitrarios: A(:,[2 4]) está
formada por las columnas segunda y cuarta de A.
Estos índices se pueden usar a ambos lados de una instrucción de asignación:
A(:,[2 4 5]) = B(:,1:3) reemplaza las columnas 2, 4 y 5 de A por las
tres primeras de B. Se muestra y asigna la matriz A alterada completa. Pruebe.
Las columnas 2 y 4 de A pueden multiplicarse por la derecha por una matriz 2 × 2:
>> A(:,[2,4]) = A(:,[2,4])*[1 2;3 4]
Ejercicio 1.1.4. Dado
>> A=round(10*rand(4,5))
calcular
>> A(1:4,3)
>> A(:,3)
>> A(1:4,:)
>> A(:,[2 4])
>> B; A(:,[2 4 5]) = B(:,1:3)}
>> A(:,[2,4]) = A(:,[2,4])*[1 2;3 4]
Si denotamos por x un vector con n componentes, ¿cuál es el efecto de la instruc-
ción x = x(n:-1:1)? Haga la prueba.
Para comprobar la utilidad de esta notación, comparar estas instrucciones de MAT-
LAB con una rutina de Pascal, FORTRAN, o C que dé los mismos resultados.
1.1.9. Comparación de la eficiencia de algoritmos: etime, cputime,

tic, toc
La función clock da la hora actual aproximada hasta la centésima de segundo
(ver help clock). Dados dos tiempos t1 y t2, etime(t2,t1) proporciona el
tiempo transcurrido de t1 a t2. Se puede, por ejemplo, medir el tiempo que requiere
la resolución de un sistema de ecuaciones dado Ax = b usando eliminación gaussiana
como sigue:
>> tiempo = clock; x = A\b; tiempo = etime(clock,tiempo)
Hagamos notar que, en máquinas que operan a tiempo compartido, etime no es
una medida fiable de la eficiencia de un algoritmo ya que la velocidad de ejecución
depende de lo ocupada que esté la máquina en un momento determinado.
La función cputime da el tiempo de CPU en segundos que ha sido utilizado por
MATLAB desde que se inicio la sesión. Así
>> t1 = cputime; x = A\b; cputime-t
nos da el tiempo CPU que se ha utilizado para calcular x=A\b. Puede ser utilizado
para calcular tiempos parciales entre distintas operaciones.
Con los comandos tic, toc calculamos el tiempo trancurrido entre uno y el otro.
Por ejemplo,
>> tic, x = A\b; toc
1.1.10. Cadenas de texto, mensajes de error, input

Las cadenas de texto se introducen en MATLAB entre comillas simples. Por ejem-
plo,
>> s = ’Esto es una prueba’
asigna la cadena de texto dada a la variable s.
Las cadenas de texto pueden mostrarse con la función disp. Por ejemplo:
>> disp(’Este mensaje se está mostrando aquí’)
Los mensajes de error se muestran mejor con la función error
>> error(’Lo siento, la matriz debe ser sim\’etrica’)
ya que ésta hace que la ejecución salga del archivo .m.
En un archivo .m el usuario puede ser avisado para introducir datos interactiva-
mente con la función input. Si MATLAB se encuentra, por ejemplo, con la instruc-
ción
>> iter = input(’Introduzca el número de iteraciones: ’)
la cadena entre comillas se muestra y la ejecución se detiene mientras el usuario
introduce los datos. Tras pulsar el retorno de carro los datos se asignan a la variable
iter y continúa la ejecución.
1.1.11. Archivos .m
MATLAB puede ejecutar una sucesión de instrucciones almacenadas en archivos
de disco. Estos archivos se denominan “archivos .m", debido a que su sufijo debe ser
“m". Gran parte del trabajo con MATLAB será el de crear y refinar archivos .m.
Hay dos tipos de archivos .m: archivos de instrucciones y archivos de funciones.
Archivos de instrucciones. Un archivo de instrucciones consiste en una suce-
sión de instrucciones normales de MATLAB. Si tuviéramos un archivo denominado
nombre.m, las instrucciones del archivo pueden ser ejecutadas sin más que escribir
la instrucción nombre. Las variables en un archivo de instrucciones son globales y,
por tanto, cambiarán los valores del espacio de trabajo.
Los archivos de instrucciones son utilizados a menudo para introducir datos en una
matriz grande. En un archivo de este tipo es bastante sencillo corregir los errores sin
tener que repetir todo el trabajo. Si, por ejemplo, se escribe en el archivo datos.m
A = [
1 2 3 4
5 6 7 8
];
entonces la instrucción de MATLAB datos hará que se efectúe la asignación especi-

ficada en datos.m.
Un archivo .m puede hacer referencia a otros, incluyendo a él mismo.
Archivos de funciones. Los archivos de funciones hacen que MATLAB tenga ca-
pacidad de crecimiento. Se pueden crear funciones específicas para un problema con-
creto, y, a partir de su introducción, tendrán el mismo rango que las demás funciones
del sistema. Las variables en las funciones son locales aunque se pueden declarar las
variables para que sean globales.
Veremos, en primer lugar, un ejemplo sencillo de archivo de función:
function a = ental(m,n)
%ENTAL Matriz entera generada aleatoriamente.
% ental(m,n) produce una matriz mxn con entradas
% enteras entre 0 y 9
a = floor(10*rand(m,n));
Una versión más general de esta función es la siguiente:
function a = ental(m,n,a,b)
%ENTAL Matriz entera generada aleatoriamente.
% ental(m,n) produce una matriz mxn con entradas
% enteras entre 0 y 9
% ental(m,n,a,b) produce las entradas de la matriz entre a y b.
if nargin < 3, a = 0; b = 9; end
a = floor((b-a+1)*rand(m,n))+a;
Esto debe escribirse en el archivo ental.m (correspondiente al nombre de la fun-

ción). La primera línea declara el nombre de la función, argumentos de entrada, y
argumentos de salida; sin esta línea el archivo sería uno de instrucciones. La instruc-
ción
>> z = ental(4,5),
por ejemplo, hará que los números 4 y 5 pasen a las variables m y n en el archivo de
función y el resultado se asigna a la variable z. Como las variables en un archivo de
función son locales, sus nombres son independientes de los que se encuentren en el
espacio de trabajo.
Hagamos notar que el uso de nargin (“número de argumentos de entrada”) per-
mite asignar un valor por defecto a una variable que se omita—como a y b en el
ejemplo.
Una función puede tener también argumentos de salida múltiples. Por ejemplo:
function [media, desv] = estad(x)

%ESTAD Media y desviación típica. Para un vector x,
% estad(x) da la media y la desviación típica de x.
% Para una matriz x, estad(x) da dos vectores fila conteniendo,
% resp., la media y la desviación típica de cada columna.
[m n] = size(x);
if m == 1
m = n; % caso de un vector fila
end
media = sum(x)/m;
desv = sqrt(sum(x.^2)/m - media.^2)
Una vez situado en el archivo de disco estad.m, la instrución de MATLAB

>> [xm, xd] = estad(x),
por ejemplo, asignará la media y la desviación típica de x a las variables xm y xd,
respectivamente. Cuando se dispone de una función con argumento de salida múltiple,
se pueden efectuar asignaciones simples. Por ejemplo,
>> xm = estad(x)
(no son necesarios los corchetes alrededor de xm) asignará la media de x a xm.
El símbolo % indica que el resto de la línea es un comentario; MATLAB ignorará
el resto de la línea. Las primeras líneas de comentario, que documentan el archivo,
son accesibles con la instrucción help. Así, para que se muestren en pantalla basta
escribir help estad. Dicha documentación debe incluirse siempre en un archivo de
función.
Esta función ilustra algunas de las formas en que MATLAB puede usarse para
obtener un código eficiente. Hagamos notar, por ejemplo, que x.ˆ2 es la matriz de
los cuadrados de las entradas de x, que sum es una función vectorial (sección 1.1.6),
que sqrt es una función escalar (sección 1.1.7), y que la división en sum(x)/m
opera una matriz con un escalar.
La siguiente función, que da el máximo común divisor de dos enteros vía el algo-
ritmo de Euclides, ilustra el uso de un mensaje de error (ver sección siguiente).
function a = mcd(a,b)
% MCD Máximo común divisor
% mcd(a,b) es el máximo común divisor de a y b no nulos a la vez.
a = round(abs(a)); b = round(abs(b));
if a == 0 & b == 0
error(’El mcd no está definido si ambos números son nulos’)
else
while b ~= 0
r = rem(a,b);
a = b; b = r;
end
end
MATLAB esta preparado pera hacer operacions con matrices y vectores más rápido
que si las programaramos como lo haríamos si no existieran esas operaciones directas.
Así, es conveniente vectorizar el algoritmeo al máximo en las M-filas pera obtener
algoritmos más rápidos. Siempre que sea posible hay que convertir secuencias for y
while en operaciones matriciales o vectoriales. Por ejemplo, es equivalente calcular
>> i=0;
>> for t=0:0.001:10;
>> i=i+1;
>> y(i)=sin(t);
>> end
que
>> t=0:0.001:10;
>> y=sin(t);
A pesar de que en MATLAB no es necesario predimensionar los vectores a veces
puede ser conveniente hacerlo. Si no se pre-dimensionan, MATLAB redimensiona los
vectores cada paso de la iteración. Esto se elimina predimensionando y, de esta forma,
lo calcula mas ràpidamente. Para predimensionar (hay que saber la dimensión que
tendrà el vector y) hay que hacer y=zeros(1,10001).
Ejercicio 1.1.5. Comparar la diferencia, en tiempo computacional, de ejecutar los dos
programas anteriores para evaluar y=sin(t) sin tener los vectores predimension-
ados (hacer clear t y previamente) y ejecutarlos despues de haberlos predimen-
sionado previamente (t=zeros(1,10001);y=t;.
En MATLAB, dados dos vectores v i u de dimensión n × 1 es equivalente, calcular

su producto interior escribiendo p=v’*u que haciendo
function c=prodi(a,b)
%producto interior
for i=1:n
p(i)=v(i)*u(i);
end
Gràcies a esta capacidad tenemos que si queremos multiplicar dos matricess, las sigu-
ientes funciones nos dan el mismo resultado:
function C=prod1(A,B)
%Producto de la matrices A y B
[n,m]=size(A);
[p,q]=size(B);
for i=1:n
for j=1:m
C(i,j)=0;
for k=1:q
C(i,j)=C(i,j)+A(i,k)*B(k,j);
end
end
end
[n,m]=size(A);
[p,q]=size(B);
for i=1:n
for j=1:q
C(i,j)=A(i,:)*B(:,j);
end
end
[n,m]=size(A);
for i=1:n
C(i,:)=A(i,:)*B(:,:);
end
C(:,:)=A(:,:)*B(:,:); %Aco es equivalent a C=A*B;
Si la matriz tiene una estructura especial podemos aprovecharla para hacer menos
operaciones. Por ejemplo, si queremos evaluar C = A ∗ B con

A11 0 B11 0
A= ; B= (1.1)
0 A22 0 B22
siendo A11, A22, B11 y B22 matrices de dimensión n × n, con cualquiera de las
funciones prod anteriores estaríamos muchas veces multiplicando por cero y, por
tanto, haciendo operaciones no necesarias. Obtendríamos el mismo resuldo haciendo
[n,m]=size(A11);
C=[A11*B11, zeros(n);zeros(n),A22*B22];
pero con muchas menos operaciones.
Ejercicio 1.1.6. Dadas las matrices
>> A11=hilb(20); A22=magic(20); B11=pascal(20); B22=rand(20);
(utilizar el comando help para ver que matrices hemos obtenido), queremos evaluar
(1.1), hacerlo utilizando las distintas alternativas. Calcular el tiempo necesario para
cada una de las alternativas (Utilizar las funciones de MATLAB ETIME, TIC, TOC,
CLOCK i/o CPUTIME) y compararlos.
1.1.12. Gráficos
MATLAB puede producir gráficos planos y gráficos de malla de superficies tridi-
mensionales. Para ver algunas de sus posibilidades escriba plot demo.
Gráficos planos. La instrucción plot crea gráficos en el plano XY; si x e y son
vectores de la misma longitud, la orden plot(x,y) accede a la pantalla gráfica y
realiza un gráfico plano de los elementos de x contra los elementos de y. Por ejem-
plo, podemos dibujar la gráfica de la función seno sobre el intervalo [−4, 4] con las
instrucciones siguientes:
x = -4:.01:4; y = sin(x); plot(x,y)
Inténtelo. El vector x es una partición del dominio con paso 0.01 mientras que y
es un vector (sin es vectorial) con los valores que toma el seno en los nodos de esta
partición.
Para volver a la pantalla alfanumérica desde la gráfica, se pulsa cualquier tecla. Por
el contrario, para acceder a la pantalla gráfica, se usa la orden shg (show graph). Si
su máquina soporta ventanas múltiples con una ventana gráfica aparte, puede desear
mantener la ventana gráfica expuesta —aunque a un lado— y la ventana alfanumérica
activa.
2
Como un segundo ejemplo, puede dibujar la gráfica de y = e−x sobre el intervalo
[−1,5, 1,5] como sigue:
x = -1.5:.01:1.5; y = exp(-x.^2); plot(x,y)
Hagamos notar que ˆ está precedido por un punto para asegurarnos que opera a
coordenadas (ver sección 1.1.2).
Pueden hacerse también gráficos de curvas definidas paramétricamente. Por ejem-
plo,
t=0:.001:2*pi; x=cos(3*t); y=sin(2*t); plot(x,y)
La instrucción grid hará un cuadriculado en el gráfico actual.

Pueden ponerse títulos, comentarios en los ejes o en cualquier otra parte con los
siguientes comandos que tienen una cadena como argumento:
title título del gráfico
xlabel comentario en el eje x
ylabel comentario en el eje y
gtext texto posicionado interactivamente
text texto posicionado mediante coordenadas Por ejemplo, la instrucción
title(’La funci\’on m\’as bella’)
proporciona un título al gráfico. El comando gtext(’La mancha’) permite posi-

cionar una cruz en el gráfico con las flechas o el ratón, donde se situará el texto cuando
se pulse cualquier tecla.
Por defecto, los ejes se autoescalan. Para evitarlo se usa el comando axis. Si c =
[xmı́n , xmáx , ymı́n , ymáx ] es un vector con 4 elementos, entonces axis(c) establece el
escalado de ejes a los límites prescritos. axis, por sí mismo congela el escalado actual
para gráficos subsecuentes; Escribiendo axis de nuevo volvemos al autoescalado. El
comando axis(’square’) asegura que se use la misma escala en ambos ejes.
Dos formas de obtener dibujos múltiples se ilustran con
>> x=0:.01:2*pi;y1=sin(x);y2=sin(2*x);y3=sin(4*x);
>> plot(x,y1,x,y2,x,y3)
y formando una matriz Y conteniendo los valores funcionales como columnas
>> x=0:.01:2*pi; Y=[sin(x)’, sin(2*x)’, sin(4*x)’];

>> plot(x,Y)}
Otra forma es con hold. El comando hold congela la pantalla gráfica actual de forma
que los gráficos posteriores se sobreimponen en ella. Escribiendo hold de nuevo se
libera el hold. Los comandos hold on y hold off también se pueden utilizar.
Se pueden evitar los tipos de línea y de punto por defecto. Por ejemplo,
>> x=0:.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x);

>> plot(x,y1,’--’,x,y2,’:’,x,y3,’+’)}
produce líneas a trazos y de puntos para las dos primeras, mientras que para la tercera
se obtiene el símbolo + en cada nodo. Los tipos de líneas y de puntos son:
Tipos de línea: sólido (-), a trazos (-). puntos (:), punto y trazo (-.)
Tipos de puntos: punto (.), más (+), estrella (*), círculo (o), equis (x) Ver help
plot para los colores de las líneas y puntos.
El comando subplot se usa para dividir la pantalla de forma que puedan verse
hasta cuatro gráficos a la vez. Ver help subplot.
Gráficos de malla de superficies tridimensionales. Los gráficos de malla de su-
perficies tridimensionales se hacen con la función mesh. La instrucción mesh(z)
crea un gráfico tridimensional en perspectiva de la matriz z. La superficie de malla
está definida por las coordenadas z de los puntos sobre un cuadriculado rectangular en
el plano XY. Por ejemplo, pruebe con mesh(eye(10)).
Para dibujar la gráfica de una función z = f (x, y) sobre un rectángulo, se definen
en primer lugar los vectores xx e yy que dan particiones de los lados del rectángulo.
Con la función meshdom (mesh domain) se crea una matriz x, en la que cada fila es
igual a xx, y de igual forma una matriz y, con todas sus columnas iguales a yy, como
sigue: [x,y] = meshdom(xx,yy); Entonces se computa la matriz z, obtenida
evaluando f entrada a entrada sobre las matrices x e y, para aplicarle la función mesh.
2 2
Se puede, por ejemplo, dibujar la gráfica de z = e−x −y sobre el cuadrado [−2, 2]×
[−2, 2] como sigue (inténtelo):
>> xx = -2:.1:2;
>> yy = xx;
>> [x,y] = meshdom(xx,yy);
>> z = exp(-x.^2 - y.^2);
>> mesh(z)
Se podría, desde luego, cambiar las tres primeras líneas en lo anterior por
>> [x,y] = meshdom(-2:.1:2, -2:.1:2);
Para más detalles sobre mesh, ver la Guía del usuario.
Consulte la ayuda para plot3, mesh, y surf.
1.1.13. Consulta
Hay muchas características de MATLAB que no pueden incluirse en estas notas
introductorias. Mostraremos a continuación algunas de las funciones disponibles en
MATLAB, agrupadas por áreas (Fuente: MATLAB User’s Guide, versión 3.5.). Use la
instrucción help o consulte la Guía del usuario para una información más detallada
sobre las funciones.
Hay muchas funciones aparte de estas. Existen, en particular, varias “cajas de he-
rramientas” de funciones para áreas específicas; entre ellas, las de proceso de señales,
teoría de control, control robusto, identificación de sistemas, optimización, splines,
quimiometría, µ-análisis y síntesis, identificación, y redes neurales Las cajas de her-

ramientas, que son opcionales, pueden no estar instaladas en su sistema.. Para explorar
siempre se puede usar la instrucción help.
G ENERAL
help ayuda
demo demostraciones
who muestra variables en memoria
what muestra archivos .m en el disco
size número de filas y columnas
length longitud de un vector
clear limpia el espacio de trabajo
computer tipo de computadora
ˆC interrupción local
exit salida de MATLAB
quit lo mismo que exit
O PERADORES MATRICIALES O PERADORES PUNTUALES

+ suma + suma
− resta − resta
∗ multiplicación .∗ multiplicación
/ división derecha ./ división derecha
\ división izquierda .\ división izquierda
ˆ potenciación .ˆ potenciación
’ conjugada traspuesta .’ traspuesta
O PERADORES L ÓGICOS Y R ELACIONALES

< menor que & y
<= menor o igual que | o
> mayor que ∼ no
>= mayor o igual que
== igual
∼= no igual
C ARACTERES ESPECIALES
= instrucción de asignación
[ usado para formar vectores y matrices
] ver [
( precedencia aritmética
) ver (
. punto decimal
... la instrucción continúa en la siguiente línea
, separa índices y argumentos de función
; acaba filas, suprime la impresión
% comentarios
: indexación, generación de vectores
! ejecuta instrucción del sistema operativo
VALORES E SPECIALES
ans respuesta cuando no se asigna la expresión
eps precisión
pi π
√
i, j −1
inf ∞
NaN No Número (Not-a-Number)
clock reloj de pared
date fecha
nargin número de argumentos de entrada de una función
nargout número de argumentos de salida de una función
A RCHIVOS DE DISCO
chdir cambiar de directorio
delete borrar archivo
diary diario de la sesión
dir directorio de archivos en el disco
load cargar variables de un archivo
save guardar variables en un archivo
type mostrar función o archivo
what mostrar archivos .m en el disco
fprintf escribir en un archivo
pack compactar memoria vía save
M ATRICES E SPECIALES
compan compañera
diag diagonal
eye identidad
gallery esotérica
hadamard Hadamard
hankel Hankel
hilb Hilbert
invhilb inversa de Hilbert
linspace vectores igualmente espaciados
logspace vectores logarítmicamente espaciados
magic mágica cuadrada
meshdom dominio para puntos de malla
ones constante
pascal Pascal
rand elementos aleatorios
toeplitz Toeplitz
vander Vandermonde
zeros cero
M ANIPULACIÓN DE MATRICES
rot90 rotación
fliplr invierte el orden de las columnas
flipud invierte el orden de las filas
diag diagonal
tril parte triangular inferior
triu parte triangular superior
reshape reordena una matriz en otra
.’ traspuesta
: convierte una matriz en una columna simple; A(:)
F UNCIONES L ÓGICAS Y R ELACIONALES

any condiciones lógicas
all condiciones lógicas
find encuentra índices de valores lógicos
isnan detecta NaNs
finite detecta infinitos
isempty detecta matrices vacías
isstr detecta variables de cadena
strcomp compara variables de cadena
C ONTROL DE F LUJO
if ejecuta instrucciones condicionalmente
elseif usado con if
else usado con if
end termina if, for, while
for repite instrucciones un número de veces
while repite instrucciones mientras una sentencia lógica sea verdadera
break sale de los bucles for y while
return salida desde funciones
pause pausa hasta que se pulse una tecla
P ROGRAMACIÓN Y ARCHIVOS . M
input obtiene números desde el teclado
keyboard llamada al teclado como si fuera un archivo .m
error muestra mensaje de error
function define función
eval evalúa un texto en variables
feval evalúa función dada por una cadena
echo permite mostrar las instrucciones en pantalla
exist comprueba si las variables existen
casesen sensibilidad a las mayúsculas
global define variables globales
startup archivo de inicialización
getenv accede a una variable de entorno
menu genera un menú
etime tiempo gastado
T EXTO Y C ADENAS
abs convierte cadena en valores ASCII
eval evalúa texto como instrucciones
num2str convierte números en cadenas
int2str convierte enteros en cadenas
setstr indicador de cadenas
sprintf convierte números en cadenas
isstr detecta variables de cadena
strcomp compara variables de cadena
hex2num convierte cadenas hexadecimales en números
V ENTANA ALFANUMÉRICA
clc limpia pantalla
home mueve cursor al comienzo
format establece el formato de salida
disp muestra matriz o texto
fprintf imprime número formateado
echo permite la muestra de las instrucciones
G RÁFICOS
plot gráfico lineal en el plano XY
loglog gráfico logarítmico en el plano XY
semilogx gráfico semilogarítmico
semilogy gráfico semilogarítmico
polar gráfico polar
mesh superficie de malla tridimensional
contour plano de contornos
meshdom dominio para gráficos de superficie
bar gráficos de barras
stairs gráficos de escaleras
errorbar añade barras de errores
A NOTACIÓN G RÁFICA
title título
xlabel anotación en el eje x
ylabel anotación en el eje y
grid dibuja cuadriculado
text posiciona un texto arbitrariamente
gtext posiciona un texto con el ratón
ginput input gráfico
C ONTROL DE LA V ENTANA G RÁFICA

axis escalado manual de ejes
hold mantiene gráfico en pantalla
shg muestra la pantalla gráfica
clg limpia la pantalla gráfica
subplot divide la pantalla gráfica
I MPRESIÓN DE GRÁFICOS
print envía gráfico a impresora
prtsc volcado de pantalla
meta archivo de gráficos
F UNCIONES ELEMENTALES
abs módulo complejo
angle argumento complejo
sqrt raíz cuadrada
real parte real
imag parte imaginaria
conj conjugado complejo
round redondeo al entero más cercano
fix redondeo hacia cero
floor redondeo hacia −∞
ceil redondeo hacia ∞
sign función signo
rem resto
exp exponencial base e
log logaritmo natural
log10 logaritmo base 10
F UNCIONES T RIGONOMÉTRICAS
sin seno
cos coseno
tan tangente
asin arcoseno
acos arcocoseno
atan arcotangente
atan2 arcotangente de x/y
sinh seno hiperbólico
cosh coseno hiperbólico
tanh tangente hiperbólica
asinh arcoseno hiperbólico
acosh arcocoseno hiperbólico
atanh arcotangente hiperbólica
F UNCIONES E SPECIALES
bessel función de Bessel
gamma función Gamma
rat aproximación racional
erf función de error
inverf inversa de la función de error
ellipk integral completa elíptica de primera especie
ellipj integral elíptica de Jacobi
D ESCOMPOSICIONES Y FACTORIZACIONES
balance forma equilibrada
backsub sustitución regresiva
cdf2rdf convierte diagonales complejas en diagonales reales
chol factorización de Cholesky
eig autovalores y autovectores
hess forma de Hessenberg
inv inversa
lu factores de la eliminación gaussiana
nnls mínimos cuadrados con restricciones
null base ortonormal del núcleo
orth base ortonormal de la imagen
pinv pseudoinversa
qr factorización QR
qz algoritmo QZ
rref forma escalonada reducida por filas
schur descomposición de Schur
svd descomposición en valores singulares
C ONDICIONAMIENTO DE MATRICES
cond número de condición en la norma 2
norm norma 1, norma 2, norma de Frobenius, norma ∞
rank rango
rcond estimación de la condición (inverso)
F UNCIONES MATRICIALES ELEMENTALES

expm matriz exponencial
logm matriz logaritmo
sqrtm matriz raíz cuadrada
funm función arbitraria de matriz
det determinante
trace traza
kron producto tensorial de Kronecker
P OLINOMIOS
roots raíces de polinomios—método de la matriz compañera
roots1 raíces de polinomios—método de Laguerre
polyval evaluación de polinomios
polyvalm evaluación polinomio matricial
conv multiplicación
deconv división
residue desarrollo en fracciones parciales
polyfit ajuste por un polinomio
A NÁLISIS DE DATOS POR COLUMNAS

max valor máximo
min valor mínimo
mean valor medio
median mediana
std desviación típica
sort ordenación
sum suma de elementos
prod producto de elementos
cumsum suma acumulativa de elementos
cumprod producto acumulativo de elementos
diff derivadas aproximadas
hist histogramas
corrcoef coeficientes de correlación
cov matriz de covarianza
cplxpair reordena en pares complejos
P ROCESO DE S EÑALES
abs módulo complejo
angle argumento complejo
conv convolución
corrcoef coeficientes de correlación
cov covarianza
deconv deconvolución
fft transformada rápida de Fourier
fft2 FFT 2-dimensional
ifft FFT inversa
ifft2 FFT inversa 2-dimensinal
fftshift cambia las dos mitades de un vector
I NTEGRACIÓN N UMÉRICA
quad función de integración numérica
quad8 función de integración numérica
S OLUCIÓN DE ECUACIONES DIFERENCIALES

ode23 método de Runge-Kutta de orden 2/3
ode45 método de Runge-Kutta-Fehlberg de orden 4/5
E CUACIONES NO LINEALES Y O PTIMIZACIÓN

fmin mínimo de una función de una variable
fmins mínimo de una función de varias variables
fsolve solución de un sistema de ecuaciones no lineales
(ceros de una función de varias variables)
fzero cero de una función de una variable
I NTERPOLACIÓN
spline spline cúbico
table1 genera tablas 1-D
table2 genera tablas 2-D
1.2. CÁLCULO MATRICIAL BÁSICO 32
1.2. Cálculo Matricial Básico

Los escalares se representarán habitualmente por letras griegas minúsculas: α, β,
γ, . . .. Interpretaremos los vectores como vectores columna cuando no se especifica el
contrario y los representaremos habitualmente por letras romanas minúsculas: x, y, a,
b, . . .. Las matrices las representaremos por letras romanas mayúsculas: A, B, X, . . ..
Los vectores de la base canónica de Rn los representaremos por ei , donde ei será
el vector que tiene un 1 en la componente i-ésima y ceros en el resto. El conjunto de
las matrices m × n lo representamos por Rm×n , es decir, implícitamente identificamos
una matriz con el vector formado por sus elementos, en algún orden (por filas o por
columnas). El elemento (i, j) de una matriz A se representa por Aij o A(i, j) o aij , si
hemos especificado previamente que A = (aij ).
La notación que utilizaremos para referirnos a las filas, las columnas y las subma-
trices de una matriz es de tipo matlab, por ejemplo
A(i, :) ≡ fila i-ésima de A.
A(:, i) ≡ columna i-ésima de A.
Así mismo, si v1 , . . . , vn ∈ Rm , entonces representaremos la matriz, las columnas

de la cual son v1 , . . . , vn , por [v1 . . . vn ]. De la misma forma, si u1, . . . , um ∈ Rn son
vectores fila, entonces representaremos la matriz, las filas de la cual son u1, . . . , um ,
por  
u1
 .. 
 . .
um
Mirad la sección 1.1.8 para más detalles.
1.2.1. Normas de vectores

Definición 1.2.1. Una norma vectorial es una aplicación f : Rn → R+ ∪ {0} que
cumple:
1. f (x) = 0 ⇔ x = 0.
2. f (αx) = |α|f (x).
3. f (x + y) ≤ f (x) + f (y) (desigualdad triangular).
Ejemplo 1.2.2.
1. kxk∞ = máx1≤i≤n |xi |.

P
2. kxk1 = ni=1 |xi |.
Pn 1
3. kxkp = ( i=1 |xi |p ) p .
Pn 1
4. kxkp,w = ( |xi |p wi ) p , wi > 0.
i=1
P 1
5. La norma euclídea kxk2 = ( ni=1 x2i ) 2 es un caso particular de la norma k · kp ,
para p = 2.
Ejercicio 1.2.3. Hacer cuatro programs de matlab, que utilicen for y que se
llamen nori.m, nor1.m,norp.m y norw.m que calculen las normas 1-4
del ejemplo anterior. Hacer-lo de forma que las funciones sean llamadas de la
siguiente forma:
>> sol=nori(x);
>> sol=nor1(x);
>> sol=norp(x,p);
>> sol=norw(x,p,w);
Comprobar numéricamente que ∀x, norp(x,1)-nor1(x)=0.

Comprobar numéricament que si w=x*01+, entonces ∀x, norw(x,p,w)-norp(x,p)=0.
Hacer un único programa, llamado, nor.m, que calcule lo mismo que calculan
las cuatro funciones anteriores. El programa nuevo, para obtener los mismos
resultados que antes, deberia de utilizarse de la forma siguente:
>> sol=nor(x,0);
>> sol=nor(x,1);
>> sol=nor(x,p);
>> sol=nor(x,p,w);
Si haceis
>> help norm
vereis:
NORM vector norm.
NORM(V,P) = sum(abs(V).^P)^(1/P).
NORM(V) = norm(V,2).
NORM(V,inf) = max(abs(V)).
NORM(V,-inf) = min(abs(V)).
Intentar reproduir, utilizando la función norm.m, los resultados obtenidos con

les funcions nori.m, nor1.m,norp.m i norw.m.
Propiedades 1.2.4.
Diremos que dos normas f1 , f2 son equivalentes si ∃M, m > 0 tal que
mf1 (x) ≤ f2 (x) ≤ Mf1 (x).
No es difícil probar que

√
• kxk2 ≤ kxk1 ≤ nkxk2 .
√
• kxk∞ ≤ kxk2 ≤ nkxk∞ .
• kxk∞ ≤ kxk1 ≤ nkxk∞ .
Esto prueba que las normas 1, 2, ∞ son equivalentes en Rn . De hecho, se puede

probar también que en Rn todas las normas son equivalentes.
Toda norma en Rn es una función continua.
Desigualdad de Hölder:
1 1
|xT y| ≤ kxkp kykq , ∀x, y ∈ Rn , ∀p, q, + = 1. (1.2)
p q
Desigualdad de Cauchy-Schwartz: caso particular de la de Hölder, para p = 2:
|xT y| ≤ kxk2 kyk2, ∀x, y ∈ Rn . (1.3)
Ejercicio 1.2.5. Hacer un programa tal que dados p, q = 0, 1, 2, p = q, Calcule

C > 0 tal que
nor(x, p)
≤ C(≡ nor(x, p) ≤ Cnor(x, q)) ∀x ∈ Rn . (1.4)
nor(x, q)
Utilizar este programa para deducir, que no probar, que las normas kxk2 , kxk1 i kxk∞
son equivalentes. Deducir les constantes m y M en cada caso.
1.2.2. Tipos de matrices

Definición 1.2.6. La diagonal k-ésima de una matriz A está formada por los elemen-
tos Aij con j − i = k. La diagonal 0-ésima recibe el nombre de diagonal principal y
está formada por los elementos Aii i = 1, . . . , n. Ver la figura 1.1 para una visión es-
quemática de la localización de las diagonales. Una matriz se llama diagonal cuando
Aij = 0 para i 6= j, es decir, los elementos que no están en la diagonal principal son
nulos.
0 1 diagonals positives
-1
diagonals
negatives
Figura 1.1: Diagonales de una matriz
Definición 1.2.7. La matriz A ∈ Rn×n se llama triangular superior si aij = 0 para

i>j  
a11 a12 . . . a1n
 0 a22 . . . a2n 
A= . . . . . . . . . . . . . . . . . .  .
 (1.5)
0 0 . . . ann
Con el comando triu podemos calcular la parte triangular superior de una matriz.
Por ejemplo
>> A=rand(4); B=triu(A),
Notemos que son equivalentes:
>> triu(A),
>> triu(A,0),
Definición 1.2.8. La matriz A ∈ Rn×n se llama triangular inferior si aij = 0 para

i<j  
a11 0 . . . 0
 a21 a22 . . . 0 
A= . . . . . . . . . . . . . . . . . .  .
 (1.6)
an1 an2 . . . ann
Con el comando tril podemos calcular la parte triangular inferior de una matriz.
Por ejemplo
>> A=rand(4); B=tril(A),
Notad que son equivalentes:
>> tril(A),
>> tril(A,0),
Definición 1.2.9. Diremos que una matriz A = (aij ) m × n tiene un ancho de banda
inferior l ≤ m cuando aij = 0, si i > j + l, es decir, cuando las diagonales −m +
1, . . . , −l −1 son nulas (notad que la diagonal −m no existe). Diremos que una matriz
A m×n tiene un ancho de banda superior u ≤ n cuando aij = 0, si j > i+ u, es decir,

cuando las diagonales u + 1, . . . , n − 1 son nulas (notad también que la diagonal n
tampoco existe). Diremos que una matriz A m×n tiene un ancho de banda b (b ≤ m, n)
cuando aij = 0, si |i − j| > b, es decir, cuando las diagonales −m + 1, . . . , −b − 1,
b + 1, . . . , n − 1 son nulas.
Mirad el resultado que obtenemos evaluando:
>> A=rand(4); B=tril(triu(A,-1),2)
Ejercicio 1.2.10. Realizad un programa el cual, dada la matriz A y los enteros l y

u, devuelva una matriz que tenga un ancho de banda inferior l y un ancho de banda
superior u y que los valores diferentes de cero de esta matriz coincidan con los de A.
En el caso en que solo le introducimos un matriz A y un escalar b, el programa debe
devolver una matriz con un ancho de banda b.
P Una matriz A se llama estrictamente diagonalmente dominante por

Definición 1.2.11.
filas cuando i6=j |aij | < |aii |. Análogamente definimos la dominancia diagonal por
columnas.
Ejercicio 1.2.12. Haced un programa que dada una matriz A compruebe si la matriz
es diagonalmente dominante por filas.
Definición 1.2.13. La transposición de una matriz m × n A es una matriz n × m AT

que satisface
(AT )(i, j) = A(j, i) 1 ≤ i ≤ n, 1 ≤ j ≤ m. (1.7)
Esta operación corresponde a la expresión coloquial de “cambiar filas por columnas”.
Resumimos seguidamente las propiedades fundamentales del operador transposición:
(AT )T = A.
(A + B)T = AT + B T .
(αA)T = αAT .
(AB)T = B T AT .
Ejercicio 1.2.14. En matlab, para calcular la transpuesta de una matriz haremos A’.
Comprueba con dos matrices, que las propiedades anteriores son ciertas.
Definición 1.2.15. Una matriz A ∈ Rn×n es invertible si existe una matriz B ∈ Rn×n
tal que AB = BA = I; en este caso, la matriu que satisface esta condición es única,
se llama la inversa de A y se representa habitualmente por A−1 .
1.2.3. Matrices simétricas.

Definición 1.2.16. Una matriz A se llama simétrica si coincide con su transpuesta, es
decir: AT = A.
Definición 1.2.17. Una matriz se llama definida positiva cuando xT Ax > 0 ∀x 6= 0.
Las matrices simétricas y definidas positivas tienen un papel muy importante en

muchos campos de la matemática.
Ejercicio 1.2.18. Comprobad, utilizando matlab, si una matriz no es definida positiva.

Capítulo 2
Solución de ecuaciones no lineales
38
2.1. INTRODUCCIÓN 39
2.1. Introducción
Dada una función real de una variable f (x), tratamos de hallar los valores de x en
los que la función se anula, es decir, las soluciones o raíces de la ecuación f (x) = 0.
Si tenemos un polinomio de grado 2 (o uno), existe una formula para calcular la
solución. Si el polinomio es de grado 3 o 4 también se puede resolver. No existen
formulas generales para obtener la solución de polinomios de grado ≥ 5.
Si no es posible resolver analíticamente la ecuación f (x) = 0 despejando x, nos
vemos abocados a utilizar métodos iterativos, consistentes en obtener una sucesión de
valores x1 , x2 , . . . , xn , . . . cada vez más próximos a la raíz buscada.
En teoría, nos planteamos si tal raíz existe y si es o no única. En la práctica, lo más
importante es determinarla con la precisión necesaria, con el menor número posible de
operaciones. La respuesta a la cuestión de la existencia la proporciona el
Teorema 2.1.1. (Teorema de Bolzano) Si f : [a, b] → R es continua y f (a) y f (b)

tienen signos diferentes (f (a) ∗ f (b) < 0), entonces f se anula en algún punto entre a
y b.
En este caso, la unicidad de la solución puede garantizarse si la función es estric-
tamente creciente o estrictamente decreciente.
Ambos resultados proporcionan condiciones suficientes, pero no necesarias. Por

ejemplo, |x| tiene una única raíz en [−1, 1] y, sin embargo, ni cambia el signo ni es
monótona en dicho intervalo.
Veamos ahora cómo analizar la convergencia de los métodos iterativos y cómo
estimar la precisión del resultado.
2.2. Métodos Iterativos

En general, los métodos iterativos consisten en:
Obtener una aproximación inicial x1 de la solución. Esta aproximación es la

información previa de que disponemos. En algunos métodos, la elección de la
aproximación inicial determina la convergencia del método. La representación
gráfica de la función suele proporcionar aproximaciones iniciales aceptables.
Refinar la aproximación inicial mediante una fórmula iterativa que genera nuevos
valores (llamados también iterados), x2 , x3 , x4 , . . ., que idealmente, convergerán
a la solución buscada x∗ .
Establecer un criterio de parada o test de finalización, satisfecho el cual, se

detiene el proceso de obtención de iterados. En caso de convergencia, se toma el
último iterado como raíz aproximada.
2.2. MÉTODOS ITERATIVOS 40
Los métodos iterativos también pueden consistir en dado un intervalo. [a0 , b0 ], que
contenga a la solución, vamos calculando una sucesión de intervalos, [ak , bk ] cada vez
más pequeños y que también contengan a la solución.
A xk lo llamamos iterado k-ésimo (xk seria (ak + bk )/2 en el caso de intervalos) y
el error de esta aproximación viene determinado por
εk = |xk − x∗ |.
Que la sucesión de iterados xk converja a x∗ es equivalente a que εk tienda a cero

cuando k tiende a infinito. En la práctica, εk no es calculable, ya que no conocemos
x∗ . En cambio, obtenidos los iterados x1 , x2 , x3 , . . . , xk , xk+1 , . . . , podemos calcular
ek = |xk+1 − xk |, k = 1, 2, . . . .
Observemos que ek aproxima εk siempre que εk+1 εk . Por ello, podemos utilizar
en nuestros algoritmos como criterio de parada que la cantidad ek sea suficientemente
pequeña,
|xk+1 − xk | < tol o 0,5 ∗ |xk+1 − xk | < tol.
Debemos detener en cuenta que si el método convergente lentamente, estamos

subestimando el error, al contrario, si converge rápidamente, el error será seguramente
inferior a tol. Otro posible criterio es la satisfacción aproximada de la ecuación, es
decir, que
|f (xk )| < tol.
Conviene también limitar el número máximo de iteraciones, en previsión de que el

algoritmo diverja, oscile indefinidamente o converja muy lentamente.
2.2.1. Tipos de convergencia

Para comparar la eficiencia de los distintos métodos iterativos, vamos a medir la
rapidez con la que los iterados x1 , x2 , . . . , xk , . . ., convergen o no a la solución.
Al desconocer la solución, medimos la velocidad de convergencia de la sucesión
de incrementos ek , k = 1, 2, 3, . . .
Si los cocientes ek+1 /ek tienden a K, 0 < K < 1, diremos que el método con-
verge linealmente.
Si los cocientes ek+1 /ek tienden a cero y los ek+1 /e2k tienden a estabilizarse,
diremos que el método converge cuadráticamente.
Si los cocientes ek+1 /ek tienden a cero y los ek+1 /epk , para cierto 1 << p < 2,
tienden a estabilizarse, se dice que el método tiene convergencia superlineal.
2.2. MÉTODOS ITERATIVOS 41
En la práctica, el tipo de convergencia nos indica el trabajo necesario para alcanzar

determinada precisión. Con un método que converge linealmente, cada tantas itera-
ciones, por ejemplo cada 3, se obtiene un dígito adicional exacto en la solución. En
cambio, si hay convergencia cuadrática, cada tantas iteraciones se duplica el número
de decimales exactos.
Los métodos de bisección y regula falsi convergen linealmente. El método de New-
ton converge cuadráticamente, si se dan las condiciones adecuadas. El método de la
secante tiene generalmente convergencia superlineal.
Ejemplo de convergencia lineal
• La sucesión (xn ) presenta convergencia lineal con ratio de convergencia

R, 0 < R < 1, cuando
|xn+1 − x∗ |
lı́m =R
n |xn − x∗ |
• Ejemplo xn = 1 + 0,2n + (n + 1)−n converge a x∗ = 1.

|xn+1 −x∗ |
xn − x∗ |xn −x∗ |
2,000000e − 00
7,000000e − 01 3,500000e − 01
1,511111e − 01 2,158730e − 01
2,362500e − 02 1,563419e − 01
3,200000e − 03 1,354497e − 01
4,486008e − 04 1,401878e − 01
7,249986e − 05 1,616133e − 01
1,327684e − 05 1,831291e − 01
2,583231e − 06 1,945667e − 01
5,130000e − 07 1,985885e − 01
1,024386e − 07 1,996853e − 01
2,048135e − 08 1,999379e − 01
4,096043e − 09 1,999890e − 01
8,192013e − 10 1,999982e − 01
1,638400e − 10 1,999997e − 01
3,276800e − 11 2,000000e − 01
• En la tabla anterior vemos que la segunda columna, la que corresponde
−x∗ |
a los cocientes |x|xn+1
n −x |
∗ se estabiliza en 0,2, lo que quiere decir que esta
sucesión converge a cero linealmente con radio de convergencia 0.2.
• Observemos que se gana una cifra decimal en cada cierto numero de itera-
ciones.
Ejemplo convergencia cuadrática

2.3. MÉTODO DE LA BISECCIÓN 42
• (xn ) tiene convergencia cuadrática cuando
|xn+1 − x∗ |
lı́m 6= 0, ∞.
n (xn − x∗ )2
n n
• Ejemplo: xn = 1 + 0,22 + (n + 1)−2 converge a x∗ = 1.
|xn+1 −x∗ | |xn+1 −x∗ |

xn − x∗ |xn −x∗ | (xn −x∗ )2
1,2000e + 00
2,900000e − 01 2,416667e − 01 2,013889e − 01
1,394568e − 02 4,808855e − 02 1,658226e − 01
1,781879e − 05 1,277728e − 03 9,162181e − 02
1,310720e − 11 7,355831e − 07 4,128132e − 02
4,307532e − 23 3,286386e − 12 2,507314e − 01
1,844674e − 45 4,282439e − 23 9,941746e − 01
3,402824e − 90 1,844674e − 45 1,000000e + 00
1,157921e − 179 3,402824e − 90 1,000000e + 00
• En la tabla anterior vemos que la columna segunda, la cual corresponde a

−x∗ |
los cocientes |x|xn+1
n −x |
∗ no se estabiliza. En cambio, la tercera columna, la
∗
cual corresponde a |x(xn+1 −x |
∗ 2 si que se estabiliza, a 1 concretamente, lo que
n −x )
quiere decir que esta sucesión converge a cero cuadráticamente.
• Observemos que en cada paso se obtiene el doble de cifras decimales ex-

actas que en el anterior.
2.3. Método de la Bisección

Si f es una función continua en el intervalo [a, b] y f (a) y f (b) tienen distinto
signo, sabemos por el teorema de Bolzano que la ecuación f(x)=0 tiene una solución
a+b
en dicho intervalo. Tomamos el punto medio c = del intervalo y elegimos el
2
subintervalo [a, c] o [c, b] en el que f cambia de signo. Así localizamos la raíz en
el intervalo de longitud mitad de la inicial. Repitiendo varias veces este proceso se
determina la solución con la precisión deseada.
Algoritme bisecció
a0 m1 b0
a1 m2 b1
Figura 2.1: Mètodo de la bisección
Input: a, b, f , max_its, tol

Output: m
if ( signe(f (a)) = signe(f (b)) )
stop
end
for i = 1, . . . , max_its
m = 0,5 ∗ (a + b) %punt mig [a, b]
if ( 0,5 ∗ |a − b| < tol )
break
end
if ( signe(f (a)) ˜= signe(f (m)) )
b=m
else
a=m
end
end
Observemos que el intervalo [a, b] contiene la raíz. La cota del error esta dividida
por 2 a cada iteración
Función bisección.m
function [m i]=biseccio( a, b, def_func, max_its, tol)

% [m i]=biseccio( a, b, max_its, tol)
% Biseccio per a trobar solucio de f(x)=0
% La funció seria def_func
%
% INPUT:
% a: a_0
% b: b_0
% max_its: maxim número d’iteracions %’

% tol: tolerancia per a error relatiu iteracio
%
% OUTPUT:
% m: solucio
%
% LOCALS
% a: a_i % fa: f(a_i)
% b: b_{i} % fb: f(b_i)
% m: m_i % fm: f(m_i)
% i: index iteracions
% EXEMPLES:
% [m,i]=biseccio(1,2,’x^3-3*x+1’,100000,.0001)
% [m,i]=biseccio(0,1,’x^3-3*x+1’,100000,.0001)
func=inline(def_func, ’x’);
fa=func(a);
fb=func(b);
if ( fa*fb > 0 )
fprintf(’No n’’hi ha canvi de signe en [%e, %e]\n’, a, b);
return;
end
i=0;
while ( 0.5*abs(b-a) > tol & i < max_its )
m=0.5*(a+b);
fm=func(m);
fprintf(’%d\t% .16e\t% .16e\n’, i,m, fm);
if ( fm*fa < 0 )
b=m;
fb=fm;
else
a=m;
fa=fm;
end
i=i+1;
end
if ( i < max_its & 0.5*abs(b-a) < tol )

fprintf(’Convergit en %d iteracions,\nsolucio en interval...
2.4. MÉTODO DE RÉGULA FALSI 45
determinat per: [% .16e, % .16e]\n’, i, a, b);

end
2.4. Método de régula falsi

El método de bisección converge linealmente puesto que en cada paso se divide por
2 la amplitud del intervalo que contiene la raíz. Se puede acelerar la convergencia si
no sólo se tiene en cuenta el signo de la función en los extremos del intervalo, sino que
se elige un punto intermedio más próximo al extremo en el que la función sea menor
en valor absoluto.
El método de regula falsi procede como bisección, pero en lugar de tomar el punto
medio, calcula la intersección con el eje de abscisas de la recta que pasa por los puntos
(a, f (a)) y (b, f (b)). Es decir, y − f (b) = f (b)−f
b−a
(a)
(x − b). El punto de corte de la
secante con el eje de les abscisas: y = 0 ⇒ x = b − f (b) ∗ f (b)−f b−a
(a)
= ff(b)a−f (a)b
(b)−f (a)
.
quedando,
af (b) − bf (a)
c=
f (b) − f (a)
Las iteraciones convergen a la solución más rápidamente que las de bisección, pero la
a0 m1 b0
a1 m2 b1
Figura 2.2: Mètodo de regula falsi
amplitud del intervalo de búsqueda, b − a, no tiende a 0. Por ello, hay que modificar el
criterio de convergencia, terminando cuando dos iteraciones sucesivas difieran poco.
Algoritme regulafalsi
Régula falsi
2.4. MÉTODO DE RÉGULA FALSI 46
Input: a, b, f , max_its, tol

Output: m
if ( signe(f (a)) = signe(f (b)) )
stop
end
mold = b
for i = 1, . . . , max_its
m = (f (b)a − f (a)b)/(f (b) − f (a))
if ( 0,5 ∗ |a − b| < tol | 0,5 ∗ |m − mold | < tol )
break
end
if ( signe(f (a)) ˜= signe(f (m)) )
b=m
else
a=m
end
mold = m
end
Función regulafalsi.m
function [m i]=regulafalsi( a, b, def_func, max_its, tol)
% [m i]=regulafalsi( a, b, max_its, tol)
%
% INPUT:
% a: a_0
% b: b_0
%
% OUTPUT:
% m: solucio
%
% LOCALS
% a: a_i % fa: f(a_i)
% b: b_{i} % fb: f(b_i)
% m: m_i
% fm: f(m_i)
% EXEMPLES:
% [m,i]=regulafalsi(1,2,’x^3-3*x+1’,100000,.0001)
2.5. MÉTODO DE NEWTON-RAPHSON 47
% [m,i]=regulafalsi(0,1,’x^3-3*x+1’,100000,.0001)
fa=func(a);
fb=func(b);
mo=b;
if ( fa*fb > 0 )
fprintf(’No n’’hi ha canvi de signe en [%e, %e]\n’, a, b);
return;
end
i=0;
while ( 0.5*abs(b-a) > tol& 0.5*abs(m-mo) > tol & i < max_its )
m=(fb*a-fa*b)/(fb-fa);
fm=func(m);
if ( fm*fa < 0 )
b=m;
fb=fm;
else
a=m;
fa=fm;
end
i=i+1;
mo=m;
end
if ( i < max_its & 0.5*abs(b-a) < tol )

determinat per: [% .16e, % .16e]\n’, i, a, b);
end
2.5. Método de Newton-Raphson

La idea del método de Newton para resolver la ecuación f (x) = 0 es aproximar
la función f por su tangente l en una aproximación de la raíz y resolver la ecuación
l(x) = 0. Tomamos esta solución como nueva estimación de la raíz de f (x) = 0 y
repetimos el proceso hasta obtener la precisión deseada.
Sea x1 la estimación inicial de la solución. Hallamos la tangente a la gráfica de
f en el punto (x1 , f (x1 )) y tomamos la intersección x2 de la tangente con el eje OX
como nueva estimación de la solución.
x0 x2 x1
Figura 2.3: Mètodo de Newton
La ecuación de la recta tangente es
y = f (x1 ) + f 0 (x1 )(x − x1 )
y la intersección con el eje OX se obtiene haciendo y = 0:
f (x1 )
x2 ≡ x = x1 −
f 0 (x1 )
siempre que la derivada f 0 no se anule en x1 .

Procediendo de forma iterativa, en el paso k determinamos
f (xk )
xk+1 = xk − , k = 1, 2, . . . (2.1)
f 0 (xk )
Algoritme de Newton-Raphson
Newton-Raphson
Input: x0 , f , f 0 , max_its, tol
Output: x1
for i = 1, . . . , max_its
x1 = x0 − f (x0 )/f 0 (x0 )
if ( |x1 − x0 | < tol )
break
end
x0 = x1
end
Función newton.m
function [x1 i]=newton( x0, def_func, def_derfunc, max_its, tol)

% [x1 i]=newton( x0, max_its, tol)
% Newton per a trobar solucio de f(x)=0
% la funcio deu de estar definida en def_func,
% la derivada en def_derfunc
% INPUT:
% x0: x_0
% max_its: maxim número d’iteracions
%
% OUTPUT:
% x2: solucio
%
% LOCALS
% xin: x_0
% x0: x_n % x1: x_{n+1}
% fx0: f(x0) % fpx0: f’(x0)
% EXEMPLES:
% [m,i]=newton(2,’x^3-3*x+1’,’3*x^2-3’,100000,.0001)
% [m,i]=newton(1,’x^3-3*x+1’,’3*x^2-3’,100000,.0001)
derfunc=inline(def_derfunc, ’x’);
i=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bucle,
% condicio d’entrada: # iteracions < maxim # iteracions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ( i < max_its )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calcular fx0 i fpx0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fx0= func(x0);
fpx0= derfunc(x0);
if ( fpx0 ~= 0 )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Si fpx0 != 0 calcular x1, si no eixir
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1=x0-fx0/fpx0 ;
else
error(’tangent horizontal\n’);
end
fprintf(’%d\t% .16e\t% .16e\n’, i,x0, fx0);
if ( abs ( x1-x0 ) < tol )

fprintf(’Aproximacio solucio: %.15e\n valor funcio:...
%e, #its: %d\n’,x1, func(x1), i);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% renovar i, x0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0= x1 ;
i=i+1;
end
fprintf(’Aproximacio solucio: %.15e\nfuncio: %e, #its: %d\n’, ...

x1, func(x1), i);
El error cometido en la iteración k es proporcional al cuadrado del error cometido

en la iteración anterior. En la práctica, esto supone que el número de cifras decimales
exactas se duplica cada cierto número de iteraciones.
Las condiciones que garantizan la convergencia cuadrática del método de Newton
son las siguientes:
Teorema 2.5.1. Sea f ∈ C 2 ([a, b]). Si x∗ ∈ [a, b] es una raíz de la ecuación f (x) = 0
tal que f 0 (x∗ ) 6= 0 (es decir, es una raíz simple), entonces la sucesión de iterados
xk obtenida por el método de Newton converge a x∗ , para cualquier aproximación
inicial x1 elegida suficientemente próxima a la raíz. Además, la citada convergencia
es cuadrática.
Notemos que si la derivada de la función f se anula en un iterado, el siguiente paso

de Newton no está definido. En este caso, puede hacerse inestable, aunque haya una
raíz próxima. Sin embargo, si la derivada se anula en la raíz, el método de Newton
converge linealmente, en lugar de cuadráticamente.
2.6. MÉTODO DE LA SECANTE 51
2.6. Método de la secante

El método de Newton converge muy rápidamente, pero necesita evaluar la derivada
de la función en cada paso y es muy sensible a la estimación inicial. El método de la
secante tiene la ventaja de no utilizar la derivada y ser más robusto que el de Newton.
En cambio, es más lento que el de Newton, pero más rápido que el de bisección o
regula falsi.
Geométricamente, a partir de dos estimaciones x1 y x2 de la raíz, el método de
la secante aproxima la función f por la recta que pasa por los puntos (x1 , f (x1 )) y
(x2 , f (x2 )) y toma su intersección con el eje OX, x3 , como siguiente estimación. Ex-
presamos la iteración genérica de la secante en términos análogos a la de Newton como
f (xk )
xk+1 = xk − f (xk )−f (xk−1 )
, k = 1, 2, . . . . (2.2)
xk −xk−1
El método de la secante tambien puede ser decucido del de Newton si se substituye

f (xk ) en (2.1) por f (xxkk)−f
0 (xk−1 )
−xk−1
.
x0 x2 x3 x4 x1
Figura 2.4: Mètodo de la Secante
Algoritme de la secant
Secant
Input: x0 , x1 , f , max_its, tol
Output: x2
for i = 1, . . . , max_its
x2 = (f (x1 )x0 − f (x0 )x1 )/(f (x1 ) − f (x0 ))
if ( |x2 − x1 | < tol )
break
end
x0 = x1
x1 = x2
end
2.6. MÉTODO DE LA SECANTE 52
Función secante.m
function [m i]=secant( x0, x1, def_func, max_its, tol)

% [m i]=secant( x0,x1, max_its, tol)
%
% INPUT:
% x0: x0_0
% x1: x1_0
%
% OUTPUT:
% m: solucio
%
% LOCALS
% x0: x0_i % fx0: f(x0_i)
% x1: x1_{i} % fx1: f(x1_i)
% m: m_i; % fm: f(m_i)
% EXEMPLES:
% [m,i]=secant(1,2,’x^3-3*x+1’,100000,.0001)
% [m,i]=secant(0,1,’x^3-3*x+1’,100000,.0001)
t=linspace(x0,x1,100);hold off;plot(t,func(t));grid;hold on
i=0;
while ( 0.5*abs(x1-x0) > tol & i < max_its )
fx0=func(x0);
fx1=func(x1);
m=(fx1*x0-fx0*x1)/(fx1-fx0);
fm=func(m);
plot(x0,fx0,’*’,x1,fx1,’*’,m,fm,’o’);pause(1)
x0=x1;
x1=m;
2.7. RAÍCES MÚLTIPLES 53
i=i+1;
end
if ( i < max_its & 0.5*abs(x1-x0) < tol )

determinat per: [% .16e, % .16e]\n’, i, x0, x1);
end
√
El método de la secante tiene convergencia supralineal con p = (1 + 2)/2.
2.7. Raíces múltiples

Definición: x∗ es una raíz de orden p de f (x) = 0 cuando
0 = f (x∗ ) = f 0 (x∗ ) = . . . = f p−1(x∗ ).
• orden 1≡ raíz simple.

• orden 2≡ raíz doble.
Convergencia método de Newton raíces múltiples: Si x∗ es raíz de orden p, el
método de Newton converge linealmente con ratio de convergencia R = 1 − 1p .
En particular, para una raíz doble, velocidad de Newton ≡ velocidad método de
bisección: necesitamos otra métodos.
f (x)
Si x∗ es raíz de orden q entonces x∗ es raíz simple de u(x) = f 0 (x)
, entonces
Newton converge cuadráticamente para u.
Newton para u(x) = 0:
u(xn ) f (xn )f 0 (xn )
xn+1 = xn − = xn − .
u0 (xn ) f 0 (xn )2 − f (xn )f 00 (xn )
2.8. Problemas resueltos

Ejercicio 2.8.1. Determinar las raíces del polinomio p(x) = 2x3 − 3x − 1 mediante
el método de bisección mostrando el proceso paso a paso.
Solución. Representamos gráficamente el polinomio en el intervalo [−2, 2]:

x = linspace(-2,2);
p = [2 0 -3 -1];
y=polyval(p,x);
plot(x,y)
grid on
2.8. PROBLEMAS RESUELTOS 54
Gráficamente parece que x = −1 es una raíz del polinomio p(x), en efecto p(−1) =
0. Observando la Figura 2.5 partiremos del intervalo [a1 , b1 ] = [0, 2] para hallar la raíz
positiva. Procediendo de forma análoga, se obtendría la otra raíz negativa.
10
−5
−10
−15
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
Figura 2.5: p(x) = 2x3 − 3x − 1
Verificamos que el polinomio tiene distinto signo en los extremos del intervalo
p(0) = −1 < 0 y p(2) = 9 > 0.
Calculamos el punto medio del intervalo y evaluamos en éste el polinomio
c1 = (a1 + b1 )/2 = 1 y p(1) = −2 < 0.
A continuación, comprobamos que p(b1 )· (c1 ) < 0 por lo que el intervalo [c1 , b1 ] =
[1, 2] contiene a la raíz, nos centramos en él haciendo a2 = c1 y b2 = b1 , realizamos
una nueva iteración en [a2 , b2 ] = [1, 2]:
c2 = (a2 + b2 )/2 = 1,5, p(c2 ) = 1,25.
Ahora p(a2 ) · p(c2 ) < 0 entonces b3 = c2 y a3 = a2 nos quedamos con el subinter-

valo izquierdo ya que la raíz está en [a3 , b3 ] = [1, 1,5] cuyo punto medio
c3 = (a3 + b3 )/2 = 1,25, p(1,25) = −0,843750.
por lo que haremos a4 = c3 y b4 = b3 , nos centraremos en [a4 , b4 ] = [1,25, 1,5] y

continuaremos el proceso hasta que la longitud del intervalo [ak , bk ] sea tan pequeña
como queramos, en este caso tras 25 iteraciones se alcanza la solución 1.366025, para
tolerancia de 10−7 .
En la figura 2.6 mostramos los primeros iterados.
La ventaja del algoritmo de bisección de bisección es que siempre converge, si
partimos de un intervalo que contiene a la raíz. sin embargo, la convergencia es lenta
y puede darse el caso de que el punto medio calculado en un paso esté más alejado de
la raíz que el iterado del paso anterior.
2
a
a
a
a
0
b
b
b
b
−2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Figura 2.6: Método de la bisección
2
Ejercicio 2.8.2. Halla la menor raíz positiva de la ecuación e−x − cos x = 0 con una
tolerancia de 10−5 , utilizando el método de regula falsi y mostrando los pasos en las
primeras iteraciones.
Solución. En primer lugar, representamos gráficamente la función para elegir el inter-

valo de partida.
x = -3:0.1:3;
y = exp(-x.^2)-cos(x);
plot(x,y)
grid on
Observando la gráfica de la función (Figura 2.7) partimos del intervalo [a1 , b1 ] = [1, 2]
para hallar la raíz.
Evaluamos la función en los extremos:
2
f (a1 ) = e−a1 cos(a1 ) = −0,172423
2
f (b1 ) = e−b1 cos(b1 ) = 0,434462
Calculamos el punto c1 donde la recta que une los puntos (a1 , f (a1 )) y (b1 , f (b1 ))
corta al eje x:
c1 = (a1 · f (b1 ) − b1 · f (a1 ))/(f (b1 ) − f (a1 )) = 1,284111

f (c1 ) = −0,090521
A continuación, comprobamos en cuál de los subintervalos [a1 , c1 ] y [c1 , b1 ] se
encuentra la raíz y nos centramos en él. Como f (b1 ) · f (c1 ) < 0 entonces a2 = c1 y
b2 = b1 , es decir, la raíz se encuentra en [a2 , b2 ] = [1,284111, 2], realizamos una nueva
iteración:
1.2
0.8
0.6
0.4
0.2
−0.2
−3 −2 −1 0 1 2 3
Figura 2.7: f (x) = exp(−x2 ) − cos(x)
c2 = (a2 · f (b2 ) − b2 · f (a2 ))/(f (b2 ) − f (a2 )) = 1,407549

f (c2 ) = −0,024619
De nuevo f (b2 ) · f (c2 ) < 0 y por tanto, hacemos a3 = c2 y b3 = b2 , nos quedamos
con el subintervalo derecho [a3 , b3 ] = [1,407549, 4] y la siguiente iteración será:
c3 = (a3 · f (b3 ) − b3 · f (a3 ))/(f (b3 ) − f (a3 )) = 1,439320

f (c3 ) = −0,005119
0.5
0.4
0.3
0.2
0.1
−0.1
−0.2
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
Figura 2.8: Método regula falsi
Por lo que haremos a4 = c3 y b4 = b3 , nos centraremos en [a4 , b4 ] = [1,439320, 2].

En este caso vemos que el intervalo no se hace pequeño. Hemos de parar cuando
abs(ck − ck−1 ) es pequeño.
Ejercicio 2.8.3. Utiliza el método de Newton para aproximar la raíz quinta de 2 con
tres cifras decimales exactas, ilustrando el método paso a paso.
Comprueba la convergencia del método de Newton.
2.9. PROBLEMAS PROPUESTOS 57
√
Solución. Notemos que 5 2 es la solución de la ecuación f (x) = x5 − 2 = 0, cuya
derivada f 0 (x) = 5x4 utilizamos en el método de Newton.
Tomando una estimación inicial x1 = 1 y aplicando el método de Newton calcu-
lamos las siguientes iteraciones con
x51 − 2 −1 6
x2 = x1 − 4
=1− = = 1,2
5x1 5 5
5 5
x −2 1,2 − 2
x3 = x2 − 2 4 = 1,2 − = 1,152901
5x2 5 · 1,24
siguiendo este proceso con dos iteraciones más tenemos x5 = 1,148698, valor que
1
podemos comparar con el obtenido directamente al hacer 2 5 y observar que se tienen
más de tres cifras exactas.
Utilizamos el fichero newton.m para realizar los cálculos. Haciendo
>> [sol, x]=newton( 1,’x.^5-2’, ’5*x.^4’, 30, 1e-4)
obtenemos la solución sol = 1,148698, el vector de iterados x, el error 3,052991×10−5
y el número de iteraciones k = 5. Obsérvese que con una iteración más se obtiene la
solución con un error menor que 10−8, lo cual, en este caso, nos asegura que tenemos
al menos 7 decimales exactos.
Para estudiar la velocidad de convergencia, calculamos las tasas de convergencia a
partir del vector de iterados x:
>> n=max(size(x)); %Tambien podemos hacer n=length(x);

>> er=abs(x-sol);
>> tcl = er(2:n)./er(1:n-1); % Tasa de convergencia lineal.
>> tcl = er(2:n)./er(1:n-1).^2;%Tasa de convergencia cuadrática.
Observamos que las componentes del vector tcl tienden a cero y las tcc se estabilizan
por lo que la convergencia es cuadrática. Conclusión que también puede obtenerse a
partir de las gráficas de estos vectores. √
Un planteamiento más general consistiría en obtener una aproximación de n a con
a > 0, resolviendo la ecuación f (x) = xn − a = 0.
2.9. Problemas propuestos

Ejercicio 2.9.1. Dada la ecuaciones no lineales
1. 2 ∗ x − log(x) − 4 = 0;
2. sin(x) − x + 0,5 = 0;
3. exp(x) + exp(−3 ∗ x) − 4 = 0;
2.9. PROBLEMAS PROPUESTOS 58
√
4. x − 2 ∗ sin(x) − π/3 + x = 0;
5. x10 − 1 = 0;
Haz una representación gráfica para averiguar el número de soluciones y obten-

er una estimación inicial de las mismas.
Resuélvela aplicando los distintos métodos expuestos: bisección, regula falsi,

secante y Newton con tolerancia 10−5 .
Estudiar numéricamente la velocidad de convergencia de cada método.
Confeccionar una tabla indicando para cada método cuántas iteraciones son
necesarias para alcanzar la solución con la precisión requerida.
Capítulo 3
Cálculo simbólico
59
3.1. Introducción
Como su propio nombre indica, el cálculo simbólico consiste en manipular expre-
siones matemáticas haciendo abstracción de los valores numéricos concretos de las
variables involucradas, es decir tratándolas como símbolos. Según la escuela favorita
del autor de un artículo, este tipo de cálculo puede recibir también los nombres de
cálculo formal, cálculo algebraico, álgebra computacional, etc.
En un principio, las computadoras fueron diseñadas para realizar cálculos numéri-
cos (sólo con números) como las calculadoras de mano o el famoso ábaco. La razón era
muy simple; resultaba relativamente fácil reducir operaciones complejas a sucesiones
de sumas y restas o aproximar funciones complicadas por un polinomio determinado
(interpolación) o una función racional de polinomios (aproximantes de Padé). Toda la
teoría matemática de sucesiones y series permitió mas tarde elaborar algoritmos com-
plejos para obtener resultados aproximados de casi cualquier función y la rapidez de
los cálculos hizo el resto.
Pero de la misma forma que se puede reducir un producto de dos números a un
conjunto repetitivo (un algoritmo) de sumas ordenadas mediante una tabla de multi-
plicar de acceso rápido (y la división mediante otro algoritmo que use una tabla de
dividir), hay un conjunto finito de derivadas elementales que se pueden considerar
Reglas de Derivación como la derivada de un producto, la de un cociente, la de una
función de función, etc.. También es finito el número de funciones elementales que
sabemos derivar directamente, es decir, la Tabla de Derivar £No se podrán construir
algoritmos formales que usen esta Tabla de Derivar y apliquen repetitivamente las Re-
glas de Derivación para automatizar estos cálculos tediosos y permitir que las mentes
matemáticas dediquen mas tiempo a los problemas no resueltos aún?
En julio de 1836, Charles Babbage escribe en su libro de notas: “Hoy tuve por
primera vez una concepción general, pero muy confusa, sobre la posibilidad de hacer
una máquina que realice desarrollos algebraicos; es decir, sin ninguna referencia a los
valores de las letras...”. También Leibniz soñaba con construir un lenguaje formalizado
con el que una máquina de cálculo pudiera deducir todos los teoremas y resultados
apetecidos. Son estas intuiciones y sueños las que motivaron a tantos y tantos seres
humanos en la ardua tarea de construir máquinas con memoria suficiente para albergar
largas extensiones de desarrollos algebraicos y lo bastante rápidas para que el tiempo
necesario para realizarlos estuviera a escala humana.
Con cierto retraso, se pusieron a la tarea los creadores de software para generar
aplicaciones propietarias (de pago) o de libre uso (gratuitas) que poco a poco fueron
dando lugar a los actuales sistemas simbólicos. Entre ellos destacan el legendario Mac-
syma y el no menos potente Maple. Una rápida visión de la actual disponibilidad de
herramientas simbólicas la podéis ver en:
http://en.wikipedia.org/wiki/Comparison_of_computer_algebra_systems
Donde Matlab aparece como Symbolic Math Toolbox o con las características de
Maple que es el motor simbólico que usa. A su vez, Maxima, es un derivado gratuito
del Macsyma del MIT (Instituto Tecnològico de Massachusetts) y Scilab es lo mas
3.2. DEFINIENDO VARIABLES SIMBÓLICAS 61
parecido a la versión libre de Matlab en este campo.

En lo que respecta al cálculo simbólico, cada alumno puede optar por el sistema en
que se sienta mas cómodo ya que unos se basan en Algol, otros en Lisp, etc. Unos son
fáciles de obtener y otros cuestan caros. Nosotros vamos a hablar de las capacidades
simbólicas de Matlab por las razones siguientes:
1. Para no añadir el aprendizaje de una nueva aplicación restando tiempo a la posi-

bilidad de disfrutar de ella.
2. Porque muchas funciones simbólicas tienen el mismo nombre que la función

numérica similar, lo que ayuda a identificarla.
3. Porque ya conocemos una forma de interrelacionar Matlab y Latex. Aprender

una interrelación distinta con otra aplicación sería confuso como poco.
3.2. Definiendo variables simbólicas

En la ventana Workspace de Matlab van apareciendo los nombres de todas las
variables definidas junto con sus dimensiones (recordad que para Matlab toda vari-
able es una matriz) y el tipo de variable (en principio será double que es la precisión
numérica por defecto).
Al definir una variable simbólica el tipo cambia automáticamente a sym object
fácilmente convertible en una serie de caracteres si es una única expresión o un vector
columna de series de caracteres donde cada fila puede tener longitudes distintas. Tal
conversión es la que realiza el comando latex(VarSym) o el propio Matlab cuando
nos presenta en pantalla el resultado. Pero atención; una variable simbólica no es una
serie de caracteres. En determinados casos necesitaremos convertirla en una serie de
caracteres para aprovechar todas las capacidades de Matlab.
función Salida
syms Crea variables simbólicas
sym Convierte a variable simbólica
double Convierte a variable numérica
Ejemplos:
>syms x y b t; % Crea 4 variablas simbólicas x, y, b, t
>a=sym(pi) % Guarda en a la constante pi dato simbólico
a =
pi
>n=double(a) % Guarda en n el valor numérico de a
n =
3.1416
3.3. EXPRESIONES Y ECUACIONES SIMBÓLICAS 62
>A=[1 1/3;1 sqrt(2)] % A es matriz numérica 2 por 2

A =
1.0000 0.3333
1.0000 1.4142
>B=sym(A) % B es la expresión simbólica de A
B = [1, 1/3]
[1, sqrt(2)]
>C=double(B) % Recupera A en este caso
Para obtener este resultado de B, Matlab busca expresiones racionales de grandes
enteros (o irracionales si no encuentra lo anterior) que queden en el rango [num −
δ, num + δ] donde num es el dato que va a convertir y δ = num ∗ eps siendo eps la
precisión del dato. Si encuentra mas de una expresión simbólica elige la mas próxima
a num.
Como se supone que num es un resultado numérico obtenido por cálculos ante-
riores, no siempre se gana precisión pasando de numérico a simbólico. Otra cosa es
que se obtenga el valor de num desde un proceso enteramente simbólico. Mas adelante
veremos ejemplos de este hecho al estudiar la función vpa.
La función sym también admite la forma S=sym(A,flag) que permite con-
trolar el tipo de representación simbólica que asignará a S. Para más detalles hacer
>help sym
3.3. Expresiones y ecuaciones simbólicas

Una vez definida una variable simbólica, cualquier expresión válida en Matlab que
incluya dicha variable simbólica y sea del tipo f=expresion hace que la variable
f sea también una variable simbólica. Si la expresión no se asigna a ninguna variable,
entonces la variable por defecto ans se convierte en simbólica por el momento.
También se puede definir una expresión simbólica sin haber declarado previamente
como simbólicas las variables que contiene, por ejemplo:
>clear, syms x, f=x^2+2*x+1 % Crea expresión simbólica f

f =
x^2+2*x+1
>g=’z^2+1’
g =
z^2+1
Pese a que z no ha sido declarada explícitamente como variable simbólica, se obser-

va en el Workspace que g sí tiene el tipo sym object. La forma de declarar g consiste
en encerrar entre comillas simples la expresión, pero como ahora z no es sym object,
cualquier referencia a esta variable deberá hacerse también entre comillas simples (co-
mo si fuera una subexpresión simbólica)
En el ejemplo anterior hemos encadenado en la misma línea 3 comandos Matlab

separados por una , (en vez de ;) porque nos interesa ver los resultados parciales y para
ahorrar espacio en estos apuntes. Observad que la orden syms x no genera ninguna
salida porque sólo es una declaración de variable simbólica, no una asignación a otra
variable. Si ponemos >x=sym(’x’) sí se produce una respuesta al asignar el carácter
x a la variable simbólica del mismo nombre. Este hecho es totalmente coherente con
la declaración de g.
Nota: Si el carácter = aparece en una expresión simbólica, dejamos de hablar de
expresión para hablar de ecuación simbólica. En este caso el signo = no se considera
una asignación de valor a la variable de su izquierda, sino otro carácter de la expresión
simbólica.
3.3.1. Substituciones
Si en una expresión simbólica queremos sustituir una variable por otra o por una
constante para obtener su valor en un punto utilizamos la orden
subs(f,antiguas,nuevas)
Si hay que sustituir más de una variable se ponen entre llaves separadas por comas.
Por ejemplo, si queremos calcular f(3) y g(1) haríamos:
>subs(f,x,3), subs(g,’z’,1)
ans =
15
ans =
2
Nótese la referencia a la variable z no declarada como sym. Si no ponemos las

comillas simples el Matlab nos dirá que la variable z no está definida ya que no ha
recibido ningún valor ni está declarada como simbólica..
Ejemplo 3.3.1. Construir f = ax2 + bx + c. Sustituir la variable x por t. Para a=2,

b=1, c=0, obtener el valor de f cuando t=2. Para los mismos valores anteriores, obten-
er el valor de f para t=[1:4].
Solución
>syms x a b c t
>f=a*x^2+b*x+c;
>g=subs1(f,x,t)
g =
a*t^2+b*t+c
>h=subs(g,{a,b,c},{2,1,0}) % 3 variables y 3 valores ==> llaves
h =
2*t^2+t
>u=subs(h,t,2)
u =
10
>v=subs(h,t,[1:4]) % 1 variable con 4 valores
v =
3 10 21 36
Ejemplo 3.3.2. Realizar los cambios que se indican en las siguientes fórmulas:
a) En x3 + 2x + 5 ln(x) + 1 sustituir x por e.
b) En x + cosx realizar el cambio de variable x = y 2 + 1
a) >> syms x
>> subs(x^3+2*x+5*log(x)+1,exp(1))
>> ans =
31.5221
b) >> syms y
>> subs(x+cos(x),y^2+1)
>> ans =
y^2+1+cos(y^2+1)
Notad que la substitución para obtener ’v’ no precisa llaves { y } pues los 4 valores
reemplazan a la misma variable ’t’ y dan 4 respuestas. Por lo tanto se puede sustituir
una variable por todo un vector o matriz de constantes definido previamente.
3.3.2. Qué variables simbólicas hay en F?

Una expresión simbólica puede contener varias variables simbólicas y parámetros
(sean o no simbólicos). Por ejemplo F=a*t^2+b*t+c es la forma estándar de definir
un polinomio de grado 2 en la variable t, siendo a, b y c 3 parámetros. Si queremos
encontrar la forma simbólica de las raíces de F=0 vamos a tener que definir como sym
tanto t como a, b y c.
Para saber qué variables simbólicas hay definidas en una expresión, se usa la
función findsym, que busca todas las variables sym de una expresión (o todas las
definidas si no le damos ninguna expresión simbólica y la versión de Matlab lo per-
mite), exceptuando la constante pi y las variables i, j que Matlab usa para representar
la parte imaginaria de un complejo. La salida de esta función es una lista de variables
sym separadas por comas y ordenadas alfabéticamente. A la que esté más cerca de la
variable por defecto que es ’x’, la llamamos en adelante find1.Si hay dos a igual dis-
tancia de ’x’ se toma la primera posterior a ’x’. La lista no es cíclica, por lo que ’a’ es
siempre la más lejana a ’x’.
Es importante destacar que el propio Matlab llama a esta función cada vez que
vamos a trabajar con una expresión simbólica y no le indicamos explícitamente qué
variable debe considerar como variable principal. En tal caso, Matlab considera como
3.4. RESOLUCIÓN DE ECUACIONES 65
variable principal find1. Si queremos saber cual de todas las variables simbólicas es
find1 en la expresión A, basta hacer findsym(A,1)
Ejemplos:
función Salida find1

findsym(alpha+a+b) a, alpha, b b
findsym(cos(alpha)*b*x1 + 14*y,2) x1,y x1
findsym(y*(4+3*i) + 6*j) y y
3.4. Resolución de ecuaciones

Para resolver ecuaciones Matlab utiliza el comando solve con alguna de las sigu-
ientes sintaxis:
solve(’ecuación’, variable)
Variable es el nombre de la incógnita que se quiere averiguar. Si la ecuación tiene

una única variable, no es necesario escribirla.
Si la ecuación tiene la estructura f(x)=0 y la variable x está definida como sim-
bólica, el comando admite la siguiente estructura:
solve(f(x),x)
Ejemplo 3.4.1. Resolver las ecuaciones:
1. x2 + 2x + 1 = 0
2. x3 + 8 = 0
Solución: si no lo habías hecho ya en otro ejercicio, deberás definir la variable x

como simbólica. En caso contrario bastará con que escribas la ecuación entre comillas
simples.
1. >> solve(x^2+2*x+1)
ans =
-1
-1
2. >> solve(x^3+8)
ans =
-2
1+i*3^(1/2)
1-i*3^(1/2)
3.4. RESOLUCIÓN DE ECUACIONES 66
Pasamos las tres soluciones a su valor numérico.
>> double(ans)
ans =
-2.0000
1.0000 + 1.7321i
1.0000 - 1.7321i
Si tenemos sistemas de ecuaciones hemos de resolverlos de la siguiente forma:

>> [x,y] = solve(’x^2 + x*y + y = 3’,’x^2 - 4*x + 3 = 0’)
x =
1
3
y =
1
-3/2
Las soluciones serian (1, 1) y (3, −3/2).
Si la ecuación que vamos a resolver esta igualada a 0 (= 0) estos se puede omitir.
Obtendríamos las mismas soluciones haciendo:
>> [x,y] = solve(’x^2 + x*y + y - 3’,’x^2 - 4*x + 3 ’)
x =
1
3
y =
1
-3/2
Si hacemos
solve(’p*x = r’)
decide que la incógnita es x y devuelve
ans =
r/p
Si ahora hacemos
>> S = solve(’x^2+y^2 - 2*x - 1 = 0’,’x^2 - y^2 - 1 = 0’)
S =
x: [4x1 sym]
y: [4x1 sym]
que serian
3.5. LA FUNCIÓN INLINE 67
>> S.x
ans =
5^(1/2)/2 + 1/2
1/2 - 5^(1/2)/2
5^(1/2)/2 + 1/2
1/2 - 5^(1/2)/2
>> S.y
ans =
(5^(1/2)/2 + 1/2)^(1/2)
(1/2 - 5^(1/2)/2)^(1/2)
-(1/2*5^(1/2) + 1/2)^(1/2)
-(1/2 - 1/2*5^(1/2))^(1/2)
vemos que si tomamos
>> x=5^(1/2)/2 + 1/2

x =
1.6180
>> y=(5^(1/2)/2 + 1/2)^(1/2)
y =
1.2720
obtenemos
>> [x^2+y^2 - 2*x - 1,x^2 - y^2 - 1 ]

ans =
0 0
Si tenemos dos ecuaciones con tres incógnitas podemos calcular las soluciones de
dos en función de la tercera. Esto es:
>> [u,v] = solve(’a*u^2 + v^2 = 0’,’u - v = 1’)

u =
1 - (a + (-a)^(1/2))/(a + 1)
1 - (a - (-a)^(1/2))/(a + 1)
v =
-(a + (-a)^(1/2))/(a + 1)
-(a - (-a)^(1/2))/(a + 1)
3.5. La función inline

Esta sirve para definir funciones definidas a trozos. En primer lugar veamos
3.5.1. Operadores relacionales y lógicos

Para dibujar gráficas de funciones definidas a trozos, hay que hacer uso de las
llamadas variables lógicas. Lo mejor es ver un ejemplo. Creemos un vector del 1 al 8
equiespaciado
>> x=1:8
x =
1 2 3 4 5 6 7 8
Ahora escribimos
>> x>5
ans =
0 0 0 0 0 1 1 1
Lo que ha ocurrido es que donde no se cumple la condición (x>5), aparece un 0 y

donde sí se cumple un 1. La nueva variable creada es una variable lógica.
Para crear estas variables lógicas pueden utilizarse los siguientes operadores rela-
cionales
O PERADORES RELACIONALES
< Menor que
> Mayor que
<= Menor o igual
>= Mayor o igual
== Igual
∼= distinto
Además estos operadores relacionales se pueden combinar con los operadores lógi-
cos
O PERADORES LÓGICOS
& y
| o
∼ no
Observar la notación del punto delante de la operación del producto (.*) para
operar elemento a elemento.
Ahora la dibujaríamos con
>> plot(x,y,’.’);grid on
(conviene dibujarla con un trazado de puntos o cruces para que se vean las discon-
tinuidades si las hay)
3.5.2. La función inline

Hay otra forma de evaluar una función en un punto incluso cuando ésta es una
función a trozos, consiste en utilizar la función de matlab inline.
Supongamos que queremos evaluar la función en el punto x=2, x=3,x=-1, hacemos
f=inline(’(x^2+2*x-8)/(x+1)’)
f =
Inline function:
f(x) = (x^2+2*x-8)/(x+1)
>> f(2)
ans =
0
>> f(3)
ans =
1.7500
>> f(-1)
Warning: Divide by zero. % Aqui nos avisa de que dividimos
% por cero, luego cuidado con el
% punto x=-1, anula el denominador.
> In inlineeval at 13
In inline.subsref at 25
ans =
-Inf
Hay otra forma de sustituir una variable por un valor, con el comando subs,
>> subs(’x+1’,3)
ans =
4
pero éste no funciona sobre funciones a trozos, pues requiere que la variable sea
simbólica.
Veamos un ejemplo de inline sobre una función a trozos:
>> f=inline(’(x.^2+1).*(x<0)+1.*((0<=x)&(x<1))+(-x+3).*(1<=x)’)
f =
Inline function:
f(x) = (x.^2+1).*(x<0)+1.*((0<=x)&(x<1))+(-x+3).*(1<=x)
>> f(-2)
ans =5
>> f(0)
ans =1
>> f(2)
ans =1
Hay dos nociones distintas aunque relacionadas que son importantes para el Cál-
culo:
Una son las expresiones simbólicas como sen x o x2.
La otra son las “reglas” o algoritmos que permiten obtener una salida numérica a
partir de unas entradas también numéricas.
La segunda “definición” es más general; sirve, por ejemplo, para definir una fun-
ción f(x) como x2 si x es negativa o 0 y sen x si x es positiva.
Por otro lado, cualquier expresión simbólica implica una regla de evaluación. Es
decir, si sabemos que f (x) = x2 entonces sabemos que f (4) = 42 = 16.
En MATLAB, la diferencia fundamental entre una función y una expresión sim-
bólica radica en que una función puede ser invocada con argumentos y una expresión
simbólica no.
Por otra parte, una expresión simbólica puede ser derivada mientras que una fun-
ción no puede.
En MATLAB las funciones son creadas de dos formas: como ficheros .m y como
funciones inline.
La forma típica de definir una expresión simbólica es la siguiente:
syms x
f=x^2-sin(x)
f =
x^2-sin(x)
Las dos líneas siguientes muestra que podemos derivar f, pero que no podemos eval-
uarla, al menos de una forma obvia. Nótese que MATLAB reconoce la variable. En
caso de que haya varias variables simbólicas, podemos especificar aquella respecto a
la cual queremos derivar
diff(f)
ans =
2*x-cos(x)
f(4)
??? Index exceeds matrix dimensions.
Podemos evaluar f(4) sustituyendo x por 4 de la siguiente forma:
subs(f,x,4)
ans =
16.7568
Podemos también convertir f en una función inline con el comando:
fin=inline(char(f))
fin =
Inline function:
fin(x) = x^2-sin(x)
Lo que está sucediendo aquí es que el comando inline requiere una cadena como en-
trada y char convierte f, expresión simbólica, en la cadena ’x^2-sin(x)’. (Si sim-
plemente hubiéramos escrito fin=inline(f) obtendríamos un mensaje de error
puesto que f no es una cadena). La función inline fin acepta ahora argumentos:
fin(4)
ans =
16.7568
De forma similar podemos construir una función a partir de la derivada de f:
fxin=inline(char(diff(f)))
fxin =
Inline function:
fxin(x) = 2*x-cos(x)
La función MATLAB char reemplaza el argumento que recibe por la cadena que lo rep-
resenta, haciéndolo así accesible a funciones que requieren cadenas como argumento,
como por ejemplo inline.
Sin embargo, cuando char es aplicado a una expresión simbólica, el resultado aún
es una expresión simbólica y puede ser derivada:
diff(char(f))
ans =
2*x-cos(x)
La otra forma de crear una función evaluable es escribiendo una función en fichero .m
Esta es la forma principal de definir funciones en la mayor parte de aplicaciones de
MATLAB. Recordemos que un fichero .m se crea de forma separada.
Supongamos que hemos escrito el fichero .m fun1.m que contiene lo siguiente:
function out=fun1(x)
out=x^2-sin(x);
Ahora la función fun1 puede ser invocada con un argumento:
fun1(4)
ans =
16.7568
Una de las cosas que podemos querer hacer con una función es representar su
gráfica. La operación más elemental en MATLAB es dibujar un punto con unas coor-
denadas específicas: plot(4,4)
La salida de este comando es el punto azul centrado en la figura. Para dibujar una
curva MATLAB dibuja una secuencia de puntos conectados mediante segmentos de
recta.
3.6. DOMINIOS, LÍMITES Y CONTINUIDAD 72
La entrada para tales dibujos consiste en dos vectores (listas de números). El primer
argumento es el vector de coordenadas x y el segundo el vector de coordenadas y.
MATLAB conecta los puntos cuyas coordenadas aparecen en posiciones consecu-
tivas de los vectores de entrada.
Dibujemos la función que definimos en la transparencia anterior. En primer lugar
se necesita un vector de coordenadas x:
X1=-2:.5:2
X1 =
Columns 1 through 7
-2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000
Columns 8 through 9
1.5000 2.0000
X1 es un vector de nueve componentes, comenzando en -2 y avanzando hasta 2 con

incrementos de 0,5. Ahora se debe preparar un vector de coordenadas y aplicando
nuestra función sobre las coordenadas x. Para ello debemos modificar nuestra función
de tal forma que pueda operar sobre las componentes individuales de un vector.
La función MATLAB vectorize reemplaza los operadores *, ^ y / por
.*, .^ y ./ respectivamente.
Recordemos que MATLAB trabaja fundamentalmente sobre vectores y matrices,
y su interpretación por defecto de la multiplicación, división y exponenciazión es que
son operaciones sobre matrices.
El punto antes del operador indica que se debe aplicar componente a componente.
fin=inline(vectorize(f))
fin =
Inline function:
fin(x) = x.^2-sin(x)
Y1=fin(X1)
Y1 =
Columns 1 through 7
4.9093 3.2475 1.8415 0.7294 0 -0.2294
0.1585
Columns 8 through 9
1.2525 3.0907
plot(X1,Y1)
3.6. Dominios, límites y continuidad

Tras conocer como representar funciones y algunos aspectos de cálculo simbólico
en matlab, vamos a utilizarlos para trabajar en conceptos como cálculo de dominios,
límites y continuidad de funciones.
−2
−4
−6
−8
−4 −3 −2 −1 0 1 2 3 4
√
Figura 3.1: función f (x) = x2 − x − 6
√
Ejemplo 3.6.1. Determinar el dominio de la función f (x) = x2 − x − 6.
>> sym x
>> solve(x^2-x-6)
>> ans =
3
-2
Los dos valores que anulan son x = 3 y x = −2.

Dibujando el polinomio (Figura 3.1), detectamos cuando es positivo y cuando neg-
ativo.
>> fplot(’x^2-x-6’,[-4 4 -8 5]);

>> hold on; fplot(’x*0’,[-4 4 -8 5]) % Dibujo el eje x para
% situar los cambios de
% signo del polinomio.
Esto es suficiente para darse cuenta que ] −∞, −2[ es positiva, ] −2, 3[ es negativa
y ]3, +∞[ es positiva.
Luego el dominio de f (x) es ] − ∞, −2[∪]3, +ı́nf].
Si intentamos dibujar la gráfica de la función inicial (Figure 3.2), sin pensar en el
dominio de la función,
>> fplot(’sqrt(x^2-x-6)’,[-4 4 -1 5])

Warning: Imaginary parts of complex X and/or Y arguments ignored.
> In fplot at 158
Nos avisa de que estamos tomando valores de x que producen raíces de números
negativos, es decir, complejos, aún así dibuja la gráfica, aunque podríamos pensar que
entre [−2, 3] la función vale cero y eso no es cierto, pues en ese intervalo la función
no existe.
−1
−4 −3 −2 −1 0 1 2 3 4
√
Figura 3.2: función f (x) = x2 − x − 6
Ésto nos permite darnos cuenta que hay imagen para los x que caen dentro de
] − ∞, −2[∪]3, +∞], pero cuidado cómo ya he dicho antes, podríamos pensar que
entre [−2, 3] la función vale cero y eso no es cierto, pues en ese intervalo la función
no existe.
3.6.1. Cálculo de límites: Límites laterales.

.
Matlab permite calcular límites mediante el comando limit, que admite las sigu-
ientes sintaxis:
limit(f,x,a) calcula el límite de una expresión simbólica f cuando x tiende

al punto a. Si sólo hay una variable se puede omitir.
limit(f,a) calcula el límite de una expresión simbólica f , con respecto a la

variable independiente de la función, cuando ésta tiende al punto a, si hay más
de una variable toma la primera por orden alfabético.
limit(f) calcula el límite anterior cuando a = 0.
limit(f,x,a,’right’) calcula el límite por la derecha de a.
limit(f,x,a,’left’) calcula el límite por la izquierda de a.

función Salida
limit(f,x,a) lim. cuando x → a
limit(f,a) Como no está x usa find1 → a
limit(f) Si tampoco está a, usa a=0
limit(f,x,a,’left’) lim. por la izquierda.
limit(f,x,a,’right’) lim. por la derecha.
Como se observa en el cuadro anterior disponemos de todas las opciones. Veamos

algunos ejemplos tras definir >clear, syms x a t h;
función Salida
limit(sin(x)/x) 1
limit((x-2)/(x^2-4),2) 1/4
limit((1+2*t/x)^(3*x),x,inf) exp(6*t)
limit(1/x,x,0,’right’) inf
limit(1/x,x,0,’left’) -inf
limit((sin(x+h)-sin(x))/h,h,0) cos(x)
Otro ejemplo. Si hemos hecho v = [(1 + a/x)^x, exp(-x)]; entonces:

limit(v,x,inf,’left’) % debe dar [exp(a), 0]
con lo que hemos obtenido dos límites simultáneos ya que v tiene 2 expresiones sim-
bólicas (vector fila de expresiones).
2 ∗ x + 3 (x+1)
Ejemplo 3.6.2. Para calcular lı́m ( ) hemos de hacer
x→∞ 2 ∗ x + 5
>> limit(((2*x+3)/(2*x+5))^(x+1))
>> ans =
3/5
>> limit(((2*x+3)/(2*x+5))^(x+1),x,inf)
>> ans =
exp(-1)
y resuelve la indeterminación 1∞ .
1
Ejemplo 3.6.3. Calcular lı́m .
x→0 x
Tal y como sabemos por teoría, hay una división por cero y por tanto hay que
estudiar lateralidad.
Si calculamos con matlab el límite total, se obtiene
>> limit(1/x)
>> ans =
NaN
lo que nos obliga a estudiar lateralidad para poder concluir.
>> limit(1/x,x,0,’right’)
>> ans =
Inf
>> limit(1/x,x,0,’left’)
>> ans =
-Inf
Concluyendo como ya debemos saber que el límite no existe, pues los laterales existen
pero son distintos.
0
Veamos que Matlab resuelve perfectamente las indeterminaciones .
0
x2 − 4
Ejemplo 3.6.4. Calcular lı́m .
x→2 x − 2
>> limit((x^2-4)/(x-2),x,2)
>> ans =
4
3.6.2. Continuidad de funciones

Para estudiar la continuidad de una función utilizando Matlab, lo más conveniente
es compaginar toda la información gráfica y analítica. La información gráfica la ob-
tendremos representado la función convenientemente en su dominio y la analítica es-
tudiando el valor de la función en el punto de estudio y los límites laterales.
√
x2 + 3 − 1
Ejemplo 3.6.5. Estudiar la continuidad de la función f (x) = .
x2 − 1
Represento la función,
>> fplot(’(sqrt(x^2+3)-1)/(x^2-1)’,[-6 6 -6 6]); grid on
El estudio gráfico (Figura 3.3) nos permite afirmar que la función está definida en
todos los reales menos en x = −1, x = 1, que analíticamente se puede deducir pues
son los únicos puntos que anulan el denominador y el numerador está definido para
todo número real. Observar como Matlab dibuja el comportamiento asintótico en esos
puntos.
Gráficamente también puedo afirmar que los límites totales en −1 y en 1 no existen,
pues para ambos puntos los límites laterales dan diferente. Compruebo esta afirmación
analíticamente:
Hay una forma de evaluar la función en un punto concreto: comando inline
>> f=inline(’(sqrt(x^2+3)-1)/(x^2-1)’)
f =
Inline function:
−2
−4
−6
−6 −4 −2 0 2 4 6
√
x2 + 3 − 1
Figura 3.3: función f (x) =
x2 − 1
f(x) = (sqrt(x^2+3)-1)/(x^2-1)
>> f(1)
Warning: Divide by zero.
>> ans =
Inf
>> f(-1)
>> ans =
Inf
Notemos que la función no está definida en x = 1 y tampoco en x = −1.
>> syms x % Necesario para utilizar limit.

>> limit((sqrt(x^2+3)-1)/(x^2-1),x,1)
>> ans =
NaN
Estudio lateralidad,
>> limit((sqrt(x^2+3)-1)/(x^2-1),x,1,’right’)
>> ans =
Inf
>> limit((sqrt(x^2+3)-1)/(x^2-1),x,1,’left’)
>> ans =
-Inf
Este estudio me permite concluir que en x = 1 hay una discontinuidad asintótica.

Hacemos lo mismo en el punto x = −1.
>>limit((sqrt(x^2+3)-1)/(x^2-1),x,-1)
>> ans =
NaN
Estudio lateralidad,
>> limit((sqrt(x^2+3)-1)/(x^2-1),x,-1,’right’)
>> ans =
-Inf
>> limit((sqrt(x^2+3)-1)/(x^2-1),x,-1,’left’)
>> ans =
Inf
Este estudio me permite concluir que en x = −1 hay una discontinuidad asintóti-

ca.
Ejemplo 3.6.6. Estudiar la continuidad de la función

1
x+1
si x < 2
f (x) = (3.1)
(2x + 2) si x ≥ 2
Estudiamos el dominio de la función:
>>f=inline(’1./(x+1).*(x<2)+(2*x+2).*(x>=2)’)
f =
Inline function:
f(x) = 1./(x+1).*(x<2)+(2*x+2).*(x>=2)
>> f(-1)
>> ans =
Inf
>> f(2)
>> ans =
6
Si x = −1, la función no está definida,
>> x=-3:0.01:3;
>> y=1./(x+1).*(x<2)+(2*x+2).*(x>=2);
>> plot(x,y,’*’);grid on
100
80
60
40
20
−20
−40
−60
−80
−100
−3 −2 −1 0 1 2 3
Figura 3.4: función 3.1
Aparentemente se ve una discontinuidad (Figure 3.4) asintótica y otra de salto

finito.
Estudio analítico:
Lo primero es definir la variable independiente x como simbólica para que el co-
mando limit funcione adecuadamente.
>> syms x
Como no disponemos de una función simbólica que almacene la función a trozos

f (x), para calcular los límites laterales, debemos de elegir manualmente la expre-
sión analítica (rama) que le corresponde a la función según nos acercamos a 2 por la
derecha o por la izquierda, y lo mismo con el −1.
>> limit(1/(x+1),x,2,’left’)
>> ans =
1/3
>> limit(2*x+2,x,2,’right’)
>> ans =
6
Discontinuidad de salto finito en x = 2.
>> limit(1/(x+1),x,-1,’left’)
>> ans =
-Inf
>> limit(1/(x+1),x,-1,’right’)
>> ans =
Inf
Discontinuidad asintótica en x = 1.
0.8
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
−6 −4 −2 0 2 4 6
Figura 3.5: función (3.2)

sin(x) si x < π
f (x) = (3.2)
cos(x) si x > π
Si la función presenta alguna discontinuidad ésta ocurrirá en l punto de ruptura
x = π.
Represento la función en [−6, 6], definiéndola como función a trozos,
>> x=-6:0.01:6;
>> y=sin(x).*(x<pi)+cos(x).*(x>=pi);
>> plot(x,y,’*’)
Muestra claramente (Figura 3.5) la discontinuidad de salto finito en x = π.
Para estudiar la discontinuidad analíticamente, necesitamos estudiar el valor de
la función en el punto y los límites laterales.
>> f=inline(’sin(x).*(x<pi)+cos(x).*(x>=pi)’)
f =
Inline function:
f(x) = sin(x).*(x<pi)+cos(x).*(x>=pi)
>> f(pi)
>> ans =
-1
Lo primero es definir la variable independiente x como simbólica para que el comando
limit funcione adecuadamente.
>> syms x
Ahora no disponemos de una función simbólica que almacene la función a trozos f (x).
Para calcular los límites laterales, del mismo modo, debes de elegir tú la expresión
analítica (rama) que le corresponde a la función según te acerques a π por la derecha
o por la izquierda,
40
30
20
10
−10
−20
−30
−6 −4 −2 0 2 4 6
Figura 3.6: función (3.3)
>> limit(cos(x),x,pi,’right’)
>> ans =
-1
Si vamos por la derecha de π, la función que debo coger es cos(x).
>> limit(sin(x),x,pi,’left’)
>> ans =
0
Si vamos por la izquierda de π, la función que debo coger es sin(x).

Por tanto, se concluye que la función es continua sólo a la izquierda de π, presen-
tando en x = π una discontinuidad de salto finito.

−3(x+1)
2x2 −x−6
si x < 1
2 (3.3)
2x − 8x + 7 si x ≥ 1
Representamos la función a trozos,
>> x=-6:0.01:6;
>> y=(-3.*(x+1)./(2*x.^2-x-6)).*(x<1)+(2*x.^2-8*x+7).*(x>=1);
>> plot(x,y,’*’);grid on
Observamos (Figura 3.6) un punto de discontinuidad asintótica, que estará en un

valor que anula un denominador, y luego habrá que llevar cuidado también con el
punto de ruptura x = 1.
· Estudiamos si dentro del conjunto x < 1 hay algún valor que anule el denomi-
nador 2 ∗ x2 − x − 6.
3.7. SIMPLIFICADORES 82
>> solve(’2*x^2-x-6’)
>> ans =
2
-3/2
Efectivamente el punto x = −3/2 cae dentro de los menores que 1 y de x = 2 me
olvido.
Estudio analítico entorno al punto x = −3/2:
>> syms x
>> limit((-3*(x+1)/(2*x^2-x-6)),x,-3/2,’left’)
>> ans =
Inf
>> limit((-3*(x+1)/(2*x^2-x-6)),x,-3/2,’right’)
>> ans =
-Inf
Discontinuidad asintótica en x = −3/2.
Estudio del punto de ruptura, x = 1:
>> limit((-3*(x+1)/(2*x^2-x-6)),x,1,’left’)>>
ans =
6/5
>> limit(2*x^2-8*x+7,x,1,’right’)
>> ans =
1
Discontinuidad de salto finito en x = 1.
3.7. Simplificadores
Supongamos que F es una expresión simbólica con la forma F=P*exp(sqrt(3*x^2/a^4)),
donde P representa un polinomio en la variables ’x’ que puede tener grado elevado. Si
estamos buscando extremales (máximos y mínimos relativos) de F debemos anular su
derivada. Pero en ella va a aparecer la exponencial multiplicando a cada término del
polinomio resultante y sabemos que esta exponencial nunca podrá anularse en el rango
de los reales.
Resulta pues inútil derivar F y endosarle todo el mogollón que salga a cualquier
resolvedor simbólico de ecuaciones, bajo la excusa de que si lo puede hacer él, para
qué tengo que hacer yo el trabajo?. Nos podemos encontrar sorpresas inesperadas y
no serán aceptables quejas del tipo aes
˛ que a mí se me ha colgado! o esta otra aes˛
que mi ordenador va muy lento!. Hay que buscar una forma de definir esa exponencial
como una subexpresión simbólica, luego pedir que simplifique o factorize con esa
subexpresión y, finalmente, buscar los ceros del otro factor.
El Matlab dispone de varias funciones simplificadoras que deben aplicarse en los
procesos de derivación e integración (sobretodo en el primero). Estas son:
Uso del comando simplify.

Matlab permite hacer cálculos algebraicos como simplificar, factorizar o expandir.
Algunas simplificaciones sencillas las efectúa directamente el programa.
El comando que permite simplificar una expresión es simplif y, su sintaxis es:
simplif y(expresion)
donde las variables deberán definirse previamente como simbólicas.
Ejemplo 3.7.1. Simplificar las siguientes expresiones:
a) x3 + 5xy − 5x3 − 11xy − y 3
(x2 − 4)(x + 5)
b)
x+2
x+4
c)
x2 + 7x + 12
d) sin2 x + cos2 x
En el primer apartado bastará con escribir la fórmula, tras definir las dos vari-
ables x e y como simbólicas
>> syms x y
>> x^3+5*x*y-5*x^3-11*x*y-y^3
>> ans =
-4*x^3-6*x*y-y^3
En el apartado b) necesito que recurra a una identidad notable, la suma por diferencia,
diferencia de cuadrados, luego requiere del comando simplif y,
>> simplify((x^2-4)*(x+5)/(x+2)) % prueba sin usar simplify y verás
% cómo no realiza el cálculo.
>> ans =
(x+5)*(x-2)
En el c) se requiere el comando simplif y, porque el polinomio del denominador no
está factorizado.
>> simplify((x+4)/(x^2+7*x+12))
>> ans =
1/(x+3)
En el d) también, porque la fórmula contiene funciones trigonométricas:
>> simplify(sin(x)^2+cos(x)^2)>>
ans =
1
Prueba sin usar simplify y verás cómo no efectúa el cálculo.
Uso del comando pretty.

Si se quiere que la fórmula simplificada aparezca de una forma similar a como se
escribe con el lenguaje matemático ordinario se utiliza el comando pretty, de sintaxis:
pretty(expresion)
3
Ejemplo 3.7.2. Simplificar la expresión (x + 3)(x + 5) + 1 + expresando el
x+2
resultado como una función racional (cociente de dos polinomios).
>> pretty((x+3)*(x+5)+1+3/(x+2))
3
(x + 3) (x + 5) + ----- + 1
x + 2
Uso del comando factor

El comando factor permite expresar un polinomio con coeficientes racionales como
producto de polinomios de menor grado con coeficientes racionales.
Su sintaxis es
f actor(expresion)
En el caso de que el polinomio tenga raíces complejas, aparecerá como factor el pro-
ducto de cada raíz por su conjugada (es decir, el polinomio con coeficientes reales).
Ejemplo 3.7.3. Factorizar x6 + 1 y ((x + 3)(x + 4) − (x + 2)2 )2
>> factor(x^6+1)
>> ans =
(x^2+1)*(x^4-x^2+1)
>> factor(((x+3)*(x+4)-(x+2)^2)^2)
>> ans =
(3*x+8)^2
Uso del comando expand.

expand(S) escribe cada elemento de una expresión simbólica S como una pro-
ducto de sus factores. se suele utilizar en los polinomios, aunque tambiém espande las
funciones trigonométricas, exponenciales y logarítmicas. El comando que realiza una
serie de operaciones cuando su ejecución no se corresponde necesariamente con una
simplificación de la fórmula o cuando hay otra fórmula que permite escribir la expre-
sión de otro modo es expand.
Su sintaxis es
expand(expresion)
3.8. CÁLCULO DE DERIVADAS. 85
>> expand((3*x+8)^2)
>> ans =
9*x^2+48*x+64
>> expand(cos(x+y))
>> ans =
cos(x)*cos(y)-sin(x)*sin(y)
>> expand(exp(x+y))
>> ans =
exp(x)*exp(y)
Uso del comando simple.

El comando simple sirve para que el sistema combine fórmulas matemáticas y
exprese el resultado de la forma más simple. Sintaxis
simple(expresion)
>> syms x
>> p=simple(sin(x)^2+cos(x)^2);
>> p
p = 1
Tal y como hemos visto con el paquete simbólico podemos calcular límites, ahora
vamos a utilizarlo para derivar.
3.8. Cálculo de derivadas.

Para calcular la derivada de orden n de una función de una variable se utiliza el
comando diff cuya sintaxis es,
diff(f,variable,n)
que calcula la derivada de orden n de la función (expresión simbólica) f , respecto de
la variable escogida. Si no se especifica el orden de la derivada, calcula la derivada
primera. Si no se especifica la variable y en la fórmula de la función aparece la x,
derivará respecto de esa variable. Si la función depende de una única variable, no es
necesario especificar ésta para calcular la derivada.
d 2 d2
Ejemplo 3.8.1. Calcular (x + cos(x + 1) · sin(x2 )) y 2 (x ln(x))
dx dx
>> syms x
>> diff(x^2+cos(x+1)*sin(x^2),x)
>> ans =
2*x-sin(x+1)*sin(x^2)+2*cos(x+1)*cos(x^2)*x
Para que salga en formato matemático basta con enlazar con el comando pretty
>> pretty(diff(x^2+cos(x+1)*sin(x^2),x))
2 2
2 x - sin(x + 1) sin(x ) + 2 cos(x + 1) cos(x ) x
Para la segunda derivada
>> pretty(diff(x*log(x),x,2))
1/x
Estudio de la derivabilidad de una función.

Recordemos que si f (x) es una función definida en un intervalo [a, b[, si existe el
f (x0 + h) − f (x0 )
lı́m+ y es finito para todo x0 ∈ [a, b[, a dicho límite se le llama
h→0 h
derivada por la derecha en el punto x0 , y se representa por:
f (x0 + h) − f (x0 )
f+0 (x0 ) = lı́m+ donde h > 0.
h→0 h
Análogamente si f (x) está definida en el intervalo ]a, b], se define la derivada por
la izquierda en el punto x0 , ∈]a, b] cómo:
f (x0 + h) − f (x0 )
f−0 (x0 ) = lı́m− donde h < 0
h→0 h
Teniendo definidas la derivada por la derecha y por la izquierda, concluimos que,
f (x) es derivable en un punto x0 si su derivada por la derecha y por la izquierda existen
y coinciden, es decir cuando existe el límite total, y entonces la derivada de f (x) en el
punto x0 es:
f (x0 + h) − f (x0 ) 0
f 0 (x0 ) = lı́m , f− (x0 ) = f+0 (x0 ) = f 0 (x0 )
h→0 h
Vamos a utilizar lo que sabemos de Matlab para estudiar si una función es derivable
en un punto.
Ejemplo 3.8.2. Calcula la derivada de la función f (x) = 4x2 − 5x + 7 haciendo uso
de la definición y del comando diff.
>> syms x h
>> f=4*x^2-5*x+7 % Guardamos en f la función.
f =
4*x^2-5*x+7
>> g=(subs(f,x+h)-f)/h % Guardamos en g el cociente incremental
% que hay que calcular en la definición
% de derivada en un punto.
g =
(4*(x+h)^2-5*h-4*x^2)/h
>> g=expand(g) % Desarrollamos con el comando expand.
g =
8*x+4*h-5
>> der_f=limit(g,h,0) % Calculo el límite de g cuando h
% tiende a cero, es decir, la
% derivada, la guardo en der_f.
der_f =
8*x-5
>> diff(f) % Ahora hago uso del comando diff, el resultado
% debe ser el mismo.
>> ans =
8*x-5
Ejemplo 3.8.3. Estudiemos la derivabilidad de la función valor absoluto f (x) = |x|

en el punto x = 0.
>> syms x h % Necesitamos dos variables simbólicas, x y h.

>> f=abs(x) % Guardamos en f la función valor absoluto de x.
f =abs(x)
>> g=(subs(f,0+h)-subs(f,0))/h % Guardamos en g el cociente
% incremental que hay que
% calcular en la definición
% de derivada en un punto.
g =abs(h)/h
>> limit(g,h,0,’right’) % Derivada por la derecha.
>> ans =
1
>> limit(g,h,0,’left’) % Derivada por la izquierda.
>> ans =
-1
Al ser distintas la función no es derivable en 0. Si directamente realizamos el cálculo

de la derivada del valor absoluto con diff
>> diff(abs(x))
>> ans =
sign(x)
Si sustituyo la derivada en x = 0
>> subs( diff(abs(x)),0)

>> ans
=0
8000
7000
6000
5000
4000
3000
2000
1000
−1000
−2000
−5 −4 −3 −2 −1 0 1 2 3 4 5
Figura 3.7: función 3.4
nos dice que no es derivable pero no aporta información de las derivadas a la derecha
y a la izquierda.
Ejemplo 3.8.4. Estudiar la derivada de la función

x−3
x2 −2
, x≤2
f (x) = 2x−6 (3.4)
x2
, x>2
en x = 2.
Empiezo dibujando la función a ver la información que obtengo,
>> x=-5:0.001:5;
>> y=(x-3)./(x.^2-2).*(x<=2)+(2*x-6)./x.^2.*(x>2);
>> plot(x,y,’.’)
Observamos en la Figura 3.7 dos discontinuidades asintóticas. En el punto x = 2,

es aparentemente continua, lo estudio analíticamente.
Hallo el valor de f (2), recordad que como es una función a trozos puedo usar
inline, o bien puedo usar el comando subs, eligiendo manualmente la rama que cor-
responde a x = 2. Como me interesan cálculos simbólicos, voy a definir la x como
simbólica y por tanto debo de usar, el comando subs para hallar f (2).
Si elijo la rama manualmente.
>> syms x
>> subs((x-3)/(x^2-2),x,2)
>> ans =
-0.5000
Con el comando inline,

>> f=inline(’(x-3)/(x.^2).*(x<=2)+(2*x-6)/x.^2.*(x>2)’)
f =
Inline function:
f(x) = (x-3)/(x.^22).*(x<=2)+(2*x-6)/x.^2.*(x>2)
>> f(2)
>> ans =
-0.5000
>> limit((2*x-6)/x^2,x,2,’left’)
>> ans =
-1/2
>> limit((x-3)/(x^2-2),x,2,’right’)
>> ans =
-1/2
Por tanto, función es continua en x = 2, puede ser derivable o no serlo.

Estudio la derivabilidad en x = 2:
Debo de calcular f (2 + h) para h > 0 y para h < 0, encontrándonos con el problema
de siempre, al no disponer de una función simbólica que almacene f (x) de forma
simbólica, debemos calcular f (2 + h) eligiendo manualmente la rama que toca, según
vamos por la derecha o por la izquierda de 2.
>> syms x h
>> rama_izq=(x-3)/(x^2-2);
>> g=(subs(rama_izq,x,2+h)-f(2))/h
>> g =
((-1+h)/((2+h)^2-2)+1/2)/h
>> limit(g,h,0,’right’)
>> ans =
3/2
La derivada por la izquierda en x = 2 vale 3/2 (f (2) sólo sabrá lo que es si

definiste f con inline)
>> rama_dcha=(2*x-6)/x^2;
>> g=(subs(rama_dcha,x,2+h)-f(2))/h
>> g =
((-2+2*h)/(2+h)^2+1/2)/h
>> limit(g,h,0,’left’)
>> ans =
1
La derivada por la derecha de x = 2 vale 1.

Por tanto, la función es continua en x = 2 pero no es derivable en x = 2.
3.9. APLICACIONES DE LA DERIVADA: MÁXIMOS-MÍNIMOS-PUNTOS DE INFLEXIÓN. OPTIMIZAC
3.9. Aplicaciones de la derivada: Máximos-mínimos-Puntos

de Inflexión. Optimización.
Ejemplo 3.9.1. Determinar los máximos-mínimos y puntos de Inflexión de la función
1
f (x) = e x .
Debemos saber que el dominio de la función es R\{0}. Puedo empezar dibujando
la función, a ver la información que me aporta la gráfica, para después corroborar
dicha información analíticamente.
>> fplot(’exp(1/x)’,[-10 10 -1 8])
−1
−10 −8 −6 −4 −2 0 2 4 6 8 10
1
Figura 3.8: función f (x) = e x .
La gráfica (Figura 3.8) me indica que no hay ni máximos ni mínimos relativos;

parece que existe un punto de inflexión a la izquierda del cero?
>> syms x
>> f=exp(1/x);
>> der_f=diff(f,x)
>> der_f =
-1/x^2*exp(1/x)
>> solve(der_f,x)
Warning: Explicit solution could not be found.
> In solve at 140
In sym.solve at 49
>> ans =
[ empty sym ]
1
La derivada es − x12 e x que no se anula en el dominio de la función. De ahí que no
posea ni máximos ni mínimos relativos.
Calculo la derivada segunda, para determinar los puntos de Inflexión:
>> der2_f=diff(f,x,2)
>> der2_f =
2/x^3*exp(1/x)+1/x^4*exp(1/x)
>> solve(der2_f,x)
>> ans =
-1/2
Al igualar la segunda derivada a cero, obtenemos un posible punto de inflexión, que
admitiremos como tal o descartaremos según su valor en la tercera derivada,
>> der3_f=diff(f,x,3)
>> der3_f =
-6/x^4*exp(1/x)-6/x^5*exp(1/x)-1/x^6*exp(1/x)
>> subs(der3_f,x,-1/2)
>> ans =
4.3307 % Al ser distinto de cero, afirmamos que
% la función tiene un punto de
% inflexión que se alcanza en x=-1/2.
% El punto de inflexión es (-1/2, f(-1/2)).
>> subs(f,-1/2)
>> ans =
0.1353 % el punto de inflexión es (-0.5, 0.1353).
35
P
30
25
20
15
10
0
A C B
−5
−5 0 5 10 15 20
Figura 3.9: Dibujo exemple 3.9.2
Ejemplo 3.9.2. Un hombre conduce un vehículo y debe atravesar un terreno cenagoso.

Se encuentra en un puto P a 32 Km del punto más próximo (línea recta) de una
carretera (punto A) y desea llegar en el menor tiempo posible a un punto B que se
encuentra a 16 Km de A, en la línea de la carretera. A cuántos Km de A debe llegar a
la carretera (punto C) si conduce a una velocidad a 48 Km/h por el terreno cenagoso
y a 80 Km/h por la carretera? (ver Figura 3.9)
Debemos de obtener la función a minimizar en términos de una única variable, y a
partir de ahí utilizar matlab como herramienta para determinar su mínimo absoluto.
>> syms x
>> tiempo=sqrt(x^2+32^2)/48+(16-x)/80;
% función a minimizar en el intervalo [0,16].
>> der_tiempo=diff(tiempo,x)
% hallo su derivada, la guardo en der_tiempo.
>> der_tiempo =
1/48/(x^2+1024)^(1/2)*x-1/80
>> pretty(der_tiempo) % la muestro en formato matemático.
x
1/48 -------------- - 1/80
2 1/2
(x + 1024)
>> solve(der_tiempo,x) % Igualo la derivada a cero, y resuelvo.
>> ans =
24
El mínimo relativo de la función tiempo se alcanza en x = 24, que queda fuera del in-
tervalo [0, 16], por tanto lo desecho. Por el teorema de Weiertrass dado que la función
tiempo es continua en [0, 16], existe un mínimo absoluto, que se alcanzará en uno de
los extremos del intervalo.
>> subs(tiempo,x,0)
>> ans =
0.8667
>> subs(tiempo,x,16)
>> ans =
0.7454
El menor valor de la función tiempo se alcanza en x = 16. Para llegar en el menor

tiempo posible debe de ir directamente por el terreno cenagoso hasta el punto B, sin
tocar la carretera.
Ejercicios propuestos.
Ejercicio 3.9.3. Estudiar la derivabilidad de la función f (x) = |x − 3|.
Ejercicio 3.9.4. Dada la función

x+3
x−2
,
x<1
f (x) = 2
x − 5x, x ≥ 1
Estudiar su derivabilidad en x = 1.
Ejercicio 3.9.5. Determinar los máximos-mínimos y puntos de inflexión de la función

f (x) = 3x2 + 5x + 2.
3.10. CÁLCULO INTEGRAL E INTEGRACIÓN NUMÉRICA 93
Ejercicio 3.9.6. Dos postes con longitudes 6 y 8 m. se colocan verticalmente sobre el

suelo con sus bases separadas por una distancia de 10 m. Calcular la longitud mínima
de un cable que vaya desde la punta de uno de los postes hasta un punto en el suelo
entre los dos postes y luego hasta la punta del otro poste.
3.10. Cálculo integral e integración numérica

Aún no hemos hablado de la integración. MATLAB dispone de un “integrador”
simbólico, denominado int, que puede integrar fácilmente f=x^2-sin(x):
sol=int(f,0,2)
sol =
5/3+cos(2)
Sin embargo, si reemplazamos f por la función h definida de la siguiente forma:
h=sqrt(x^2-sin(x^4))
h =
(x^2-sin(x^4))^(1/2)
int será incapaz de evaluar la integral:
int(h,0,2)
Warning: Explicit integral could not be found.
> In C:\MATLABR11\toolbox\symbolic\@sym\int.m at line 58
ans =
int((x^2-sin(x^4))^(1/2),x = 0 .. 2)
Sin embargo, si escribimos double (indicando números de doble precisión) antes

de la expresión integral, MATLAB retornará el resultado de una integración numérica.
double(int(h,0,2))
Warning: Explicit integral could not be
found.
> In
C:\MATLABR11\toolbox\symbolic\@sym\int.m
at line 58
ans =
1.7196
Podemos comprobar la plausibilidad de esta respuesta dibujando h entre 0 y 2 y

estimando el área bajo la curva:
ezplot(h,[0,2]);grid
(x2 − sin(x4))1/2
1.5
0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x
p
Figura 3.10: Dibujo h(x) = x2 − sin(x4 )
El valor numérico retornado por MATLAB (Figura 3.10) es algo menos que la
mitad del área de un cuadrado de 2 unidades de lado lo cual es consistente con nuestro
gráfico.
La integración numérica invocada por la combinación de double e int no es nativa,
esto es no es propio sino de MAPLE del cual han sido tomadas las rutinas de cálculo
simbólico de MATLAB.
MATLAB también dispone de un integrador numérico denominado quadl. Las ruti-
nas double(int(...)) y quadl(...) proporcionan respuestas ligeramente distintas.
quadl(inline(vectorize(h)),0,2)
ans =
1.7196
Para derivar e integrar una expresion simbolica f, disponemos de las instrucciones

diff e int, que actuan como sigue:
int(f) Calcula una primitiva de f respecto de la variable simbolica preferente.
int(f,s) Calcula una primitiva de f respecto de la variable s.
int(f,a,b) Calcula la integral deŕnida de f respecto de la variable simbolica

preferente.
int(f,s,a,b) Calcula la integral deŕnida de f respecto de la variable s.
Nota: Por defecto, la variable preferente en una expresion simbolica es la letra x. Si

esta no interviene en la expresion, se toma la letra minuscula mas proxima a ella segun
el orden alfabetico y que no sea ni la i ni la j. En caso de que haya dos (una anterior y
otra posterior), se considera variable preferente el caracter posterior.
Ejemplo 3.10.1. Construyanse las expresiones simbolicas f = ax + b y g = y2 + z y

calculense:
R
1. f dx
R
2. f da
R
3. gdy
R1
4. 0
gdz
Solucion
>> syms a b x;
>> f=a*x+b;g=’y^2+z’;% practicamos los dos modos de definir expresiones
>> I1=int(f) % integra respecto de la variable preferente x
I1 =
1/2*a*x^2+b*x
>> I2=int(f,a) % integra respecto de a

I2 =
1/2*a^2*x+b*a
>> I3=int(g) % integra respecto de la variable preferente y.

% Equivale a int(g,’y’).
I3 =
1/3*y^3+z*y
>> I4=int(g,’z’,0,1) % integra respecto de la variable z.

I4 =
y^2+1/2
Capítulo 4
Sistemas numéricos y fuentes de error
96
4.1. Introducción
El proceso de cálculo lleva associada, inevitablemente, una serie de complicaciones
debidas a la propagación de los errores en los cálculos. Hay tres clases de errores
principales en cálculo numérico:
Incertidumbre en los datos: la incertidumbre en los datos es habitual cuando se

trata de resolver problemas prácticos. Su aparición puede ser debida a errores
en las medidas físicas de las cantidades que definen nuestro problema, a errores
de almacenamiento de los datos en el ordenador (en este caso los consideramos
errores de redondeo) o porque los datos de nuestro problema son, a la vez, la
solución de otro problema y, por lo tanto, vienen afectados de error.
Redondeo: los errores de redondeo, en cierta manera relacionados con la incer-

tidumbre de los datos, son una consecuencia directa de la aritmética de precisión
finita que usan los ordenadores. Dedicaremos este capítulo a estudiar con un
poco de detalle esta fuente de errores.
Truncamiento: el análisis del error de truncamiento asociado a un método nu-

mérico concreto es otra de les tareas más importantes del analista numérico.
Muchos métodos numéricos, por ejemplo el método de Newton para calcular una
raiz de una ecuación no lineal, se deriva usando la serie de Taylor de la función
considerada. Los términos que se omiten constituyen el error de truncamiento,
y el análisis de la magnitud de este error nos permite evaluar la capacidad de
aproximación del algoritmo en cuestión.
La naturaleza introductoria de estos apuntes no nos permite profundizar en las

técnicas del análisis de la propagación del error, sin embargo es absolutamente nece-
sario que el estudiante que se encuentra por primera vez con el análisis numérico tome
conciencia de la importancia que la propagación de los errores puede tener en los
cálculos prácticos. En este capítulo introduciremos la aritmética de punto flotante y
mostraremos, con algunos ejemplos, los peligros que pueden aparecer a la hora de
calcular con aritmética de precisión finita, cuáles pueden ser las causas de su apari-
ción y cómo se pueden evitar algunas de sus consecuencias adversas en determinadas
situaciones.
4.2. Representaciones numéricas

Los ordenadores no pueden almacenar los números con una cantidad infinita de
dígitos, sino que usan una aproximación que solo contiene un número fijo de bits
(dígitos binarios) o bytes (paquetes de 8 bits). A menudo, es posible elegir entre difer-
entes tipos de representaciones en un mismo ordenador. Las distintas aproximaciones
difieren en el número de bits que utilizan en las representaciones, es decir, difieren en
4.2. REPRESENTACIONES NUMÉRICAS 98
la longitud de palabra usada, pero también pueden diferir en otro aspecto muy impor-
tante: el tipo de representación, de punto fijo, o de punto flotante.
Para familiarizarnos con los diversos tipos de representaciones numéricas, hay que
recordar que nosotros trabajamos habitualmente en base 10, es decir, entendemos un
número cualquiera, por exemple 2543, como
2543 = 3 · 1 + 4 · 10 + 5 · 102 + 2 · 103 .
Así, cualquier entero en base 10 se puede asociar con un polinomio, en el caso anterior
p(x) = 3 + 4x + 5x2 + 2x3 , de manera que la evaluación de este polinomio para x =
10, la base de nuestro sistema numérico, nos da el entero mencionado. Así podemos
utilizar la notación siguiente:
N = (an an−1 . . . , a0 )10 = an 10n + an−1 10n−1 + . . . + a0 100 ,
donde 0 ≤ ai < 10, para representar cualquier número entero. En realidad no hay
ningún motivo especial para utilizar el número 10 como base de las representaciones
numéricas. Como los ordenadores leen impulsos eléctricos, y el estado natural es sí o
no, es más conveniente representar los números en los ordenadores mediante el sistema
binario, es decir, en base 2. Así un entero no negativo N se representa en el sistema
binario como
N = (an an−1 . . . , a0 )2 = an 2n + an−1 2n−1 + . . . + a0 20 ,
con 0 ≤ ai < 2. El número N tiene, pues, asociado un polinomio p(x), de manera que
p(2) = N, y los coeficientes de este polinomio constituyen su representación numérica
en base 2. Prácticament todos los computadores modernos utilizan el sistema binario
internamente. No obstante, los usuarios prefieren, por razones obvias, visualizar re-
sultados en el sistema decimal. Las conversiones de un sistema numérico a otro son
muy sencillas. La conversión de un número binario a decimal puede llevarse a cabo
directament evaluando p(2), por ejemplo:
(101)2 = 1 · 22 + 1 · 20 = 5,
(1100101001)2 = 1 · 29 + 1 · 28 + 1 · 25 + 1 · 23 + 1 = 809.
Esta conversión se puede implementar de manera muy eficiente usando el algoritmo

de la figura 4.1 para evaluar un polinomio
p(x) = an xn + an−1 xn−1 + . . . a1 x + a0 (4.1)
en un punto concreto, x = β en nuestro caso.
Ejercicio 4.2.1. Comparad el número de operaciones necesarias para evaluar p(x)

en (4.1) directamente y con el algoritmo de la figura 4.1.
Input: an , . . . , a0 , β
bn−1 = an
for i = n − 1, . . . , 1
bi−1 = ai + bi β
end
r = a0 + b0 β
Output: p(β) = r
Figura 4.1: Algoritmo para calcular el valor de un polinomio p(x) en un punto x = β.
Input: N, β
Calcular b1 i a0 tal que N = b1 · β + a0 .
n=1
while bn ≥ β
Calcular bn+1 i an tal que bn = bn+1 · β + an
n=n+1
end
an = bn
Output: (an an−1 . . . a1 a0 )
Figura 4.2: Algoritmo para pasar de un número entero N a base β. N =

(an an−1 . . . a1 a0 )β .
Ejercicio 4.2.2. Convertid el número binario (1100101001)2 en decimal utilizando el

algoritmo de la figura 4.1.
La conversión inversa, de decimal a binario, se hace usando la ley fundamental de

la división, D = d · c + r con 0 ≤ r < c, para confeccionar el algoritmo de la figura
4.2.
La conversión siguiente sirve como ejemplo de aplicación:
187 = 93 · 2 + 1 = (46 · 2 + 1) · 2 + 1 = 46 · 22 + 1 · 2 + 1
= (23 · 2 + 0) · 22 + 1 · 2 + 1 = 23 · 23 + 0 · 22 + 1 · 2 + 1
= (11 · 2 + 1) · 23 + 0 · 22 + 1 · 2 + 1
= 11 · 24 + 1 · 23 + 0 · 22 + 1 · 2 + 1 = · · ·
= (10111011)2.
Si x ∈ R+ , definiremos la parte entera de x como el mayor entero tal que es menor

o igual que x. Lo denotaremos xI = [x]. La parte fraccionaria se define, pues, como
xF = x − xI . Hay que tener en cuenta que la parte entera de un número real siempre
tiene una cantidad finita de dígitos significativos, pero la parte fraccionaria puede tener
infinitos dígitos significativos, por ejemplo x = 1/3 = ,33333 . . ..
En general, la representación decimal de x se puede especificar dando:
n
X ∞
X
xI = ak 10 , k
xF = bk 10−k ,
k=0 k=1
donde 0 ≤ ak < 10, k = 0, . . . , n i 0 ≤ bk < 10, ∀k. Si bk = 0, ∀k ≥ n0 , diremos que

la fracción se acaba.
Como x = xI + xF , utilitzaremos también la notación compacta siguiente:
x = (an an−1 . . . a0 .b1 b2 . . .)10 .
Cuando se utiliza el sistema binario, las expresiones formales son del mismo tipo,
es decir, si
x = (an an−1 . . . a0 .b1 b2 . . .)2 ,
entonces
n
X ∞
X
xI = k
ak 2 , 0 ≤ ak < 2; xF = bk 2−k , 0 ≤ bk < 2.
k=0 k=1
Hemos visto antes cómo llevar a cabo la conversión de enteros positivos decimales
a binarios y viceversa. La conversión de fracciones decimales a binarias y viceversa se
puede hacer utilizando un mecanismo también muy sencillo. Sea un xF , con 0 < xF <
1, podemos calcular su representación binaria (.b1 b2 . . .)2 de la manera seguiente:
P
Si xF = ∞ −k
k=1 bk 2 , tenemos
∞
X ∞
X
−k+1
2xF = bk 2 = b1 + bk+1 2−k ,
k=1 k=1
por tanto, bP
1 es la parte entera de 2xF , es decir, b1 = [2xF ]. Su parte fraccionaria es
(2xF )F = ∞ k=1 bk+1 2
−k
= b2 2−1 + b3 2−2 + . . .
Si repetimos este proceso, tenemos que b2 = [2(2xF )F ], etc. La implementación
de este proceso se puede hacer con el algoritmo de la figura 4.3
Por ejemplo, si xF = 0,625, tendremos
2(0,625) = 1,25 ⇒ b1 = 1, c1 = ,25

2(0,25) = 0,5 ⇒ b2 = 0, c2 = ,5
2(0,5) = 1,0 ⇒ b3 = 1, c3 = 0
y, por tanto, xF = (,101)2 .

4.3. ARITMÉTICA DE PUNTO FLOTANTE 101
Input: x, β (0 < x < 1)

c0 = x
i=0
while ci 6= 0
bi+1 = [βci ]
ci+1 = βci − bi+1
i=i+1
end
Output: (b1 b2 . . . bi )
Figura 4.3: Algoritmo para pasar un número decimal 0 < x < 1, a base β, x =
P i −k
k=1 bk β . Si la expresión tiene una cantidad infinita de dígitos (es decir, bi 6= 0
para todo i), el algoritmo no acaba nunca.
Este es un ejemplo especialment sencillo porque la expresión binaria de xF tiene

una cantidad finita de dígitos, pero está claro que esto no es la regla general como se
puede apreciar en el ejemplo siguiente: Sea xF = 0,1,
2(0,1) = 0,2 ⇒ b1 =0
2(0,2) = 0,4 ⇒ b2 =0
2(0,4) = 0,8 ⇒ b3 =0
2(0,8) = 1,6 ⇒ b4 =1
2(0,6) = 1,2 ⇒ b5 =1
2(0,2) = 0,4 ⇒ b6 =0
2(0,4) = 0,8 ⇒ b7 =0
..
.
y, por tanto, xF = (,0001100110011 . . .)2 .
4.3. Aritmética de punto flotante

Está claro que el almacenaje y la manipulación de los números enteros o reales
en un ordenador ha de involucrar sólo una cantidad finita de dígitos significativos.
Por esta razón los ordenadores utilizan sistemas numéricos de punto flotante para las
representaciones numéricas.
Definición 4.3.1. Un sistema numérico de punto flotanto F ⊂ R es un subconjunto de
los números reales formado por números de la forma
y = ±(.d1 d2 . . . dt )β × β e . (4.2)
El sistema F está caracterizado por cuatro parámetros enteros:
la base β,
la precisión o longitud de palabra t,
el rango de los exponentes emin ≤ e ≤ emax .
La mantisa m = (.d1 d2 . . . dt )β cumple 0 ≤ m < 1 y 0 ≤ di < β. para asegurar una

representación única, d1 , que llamaremos primer dígito significativo, ha de cumplir
d1 6= 0. En este caso, la representación (4.2) se dice normalizada, y t es la longitud de
palabra de la representación.
Es importante darse cuenta de que sólo hay una cantidad finita de números en un
sistema numérico de punto flotante. Todos los números estrictamente positivos de un
sistema numérico de punto flotante como en la definición 4.3.1 cumplen
β emin −1 ≤ y ≤ (1 − β −t ) · β emax .
Por ejemplo, si β = 2, t = 3, emin = −1 y emax = 3, entonces los números en punto

flotante positivos son
0,25, 0,3125, 0,3750, 0,4375, 0,5, 0,625, 0,750, 0,875,
1,0, 1,25, 1,50, 1,75, 2,0, 2,5, 3,0, 3,5, 4,0, 5,0, 6,0, 7,0,
y gráficamente el sistema numérico de punto flotante es:
0 1 2 3 4 5 6 7
En la Tabla 4.1 podemos ver las mantisas y los exponentes correspondientes a cada
uno de esos valores:
El almacenaje de un número real cualquiera en el ordenador supone, por lo tan-
to, elegir un elemento del sistema numérico del ordenador que será, pues, la repre-
sentación de punto flotante del número en cuestión. La obtención de la representació
de punto flotante puede ser, pues, considerada como una aplicación de R → F . Hay
dos formas de traducir un número real x a su representación de punto flotante de lon-
gitud t y base β: representación de punto flotante por corte y representación de punto
flotante por redondeo.
Si escribimos primero x = ±(.d1 d2 d3 . . . dt dt+1 dt+2 . . .)β · β e donde d1 6= 0 y
0 ≤ di < β, su representación en punto flotante, al cortar, se define como
f lT (x) = ±(.d1 d2 . . . dt )β · β e . (4.3)

m .100 .101 .110 .111 .100 .101 .110 .111

e -1 -1 -1 -1 0 0 0 0
y 0.25 0.3125 0.375 0.4375 0.5 0.625 0.75 0.875
m .100 .101 .110 .111 .100 .101 .110 .111

e 1 1 1 1 2 2 2 2
y 1 1.25 1.5 1.75 2 2.5 3 3.5
m .100 .101 .110 .111

e 3 3 3 3
y 4 5 6 7
Cuadro 4.1: Números que se pueden representar cuando emin = −1, emax = 3, β = 2
it=3
En cambio, la representación por redodeo se define como

f lR (x) = ±(d1 β −1 + d2 β −2 + . . . + dt−1 β −(t−1) + d˜t β −t ) · β e , (4.4)
donde
dt , 0 ≤ (.dt+1 dt+2 . . .)β < 1/2
d˜t =
dt + 1, 1/2 ≤ (.dt+1 dt+2 . . .)β < 1.
Ejemplo 4.3.2. Si suponemos que t = 5 i β = 10, es fácil comprobar que
f lT (285,007) = (,28500) · 103 , f lR (285,007) = (,28501) · 103 .
Hay que observar que los primeros t − 1 dígitos de la mantisa de f lR (x) pueden
no coincidir con los dígitos d1 , d2 , . . . , dt−1 , según sea el valor de dt+1 . El ejemplo
siguiente sirve de ilustración.
Ejemplo 4.3.3. Con t = 5 y β = 10 se comprueba que
f lR (285,997) = (,28600) · 103 .
Este resultado se obtiene a partir de (4.4) de la manera siguiente:
f lR (285,997) = f lR ((2 · 10−1 + 8 · 10−2 + 5 · 10−3 + 9 · 10−4 +
9 · 10−5 + 7 · 10−6 ) · 103 )
= (2 · 10−1 + 8 · 10−2 + 5 · 10−3 + 9 · 10−4 +
(9 + 1) · 10−5 ) · 103
= (2 · 10−1 + 8 · 10−2 + 5 · 10−3 + (9 + 1) · 10−4 ) · 103
= (2 · 10−1 + 8 · 10−2 + 6 · 10−3) · 103
= (,28600) · 103 .
El hecho de que los sistemas numéricos que utilizan los ordenadores tengan longi-
tud de palabra finita tiene como consecuencia que en el ordenador siempre se utilizan
aproximaciones a los números reales que en realidad quisiéramos usar. Hay, por lo
tanto, que estimar y controlar, en la medida en que eso sea posible, los errores que
se cometen cuando manipulamos estas aproximaciones y, lo que es mas importante,
cómo se propagan a través de una cadena de cálculos.
Si x ∈ R i x̃ es una aproximación a x, definimos:
El error absoluto de la aproximación: Ea (x̃) = |x − x̃|.
|x−x̃|
El error relativo de la aproximación: Er (x̃) = |x|
.
El error absoluto no da, por si mismo, demasiada información sobre la validez de
la aproximación. Por ejemplo, si tenemos x = 1000002 y x̃ = 1000000, tenemos un
error absoluto de 2 unidades, lo mismo que para y = 2,5 y ỹ = 0,5. No obstante, en
el primer caso tenemos 6 dígitos coincidentes entre x y x̃, mientras que en el segundo
caso no tenemos ninguno.
El error relativo, en cambio, da información sobre el número de dígitos significa-
tivos coincidentes. Comparar, por ejemplo Er (x̃) ≈ 2 · 10−6 , con Er (ỹ) = 2/2,5 ≈
8 · 10−1 . En el primer caso, el error relativo es mucho mas pequeño que en el segundo
caso, lo que indica que hay muchos mas dígitos correctos en el primer caso que en el
segundo. Esta relación entre el error relativo y el número de cifras significativas correc-
tas se puede cuantificar. Ese es el resultado fundamental que establece la proposición
siguiente y el corolario subsiguiente, y que justifica que los sistemas numéricos en
todos los ordenadores actuales trabajen con redondeo.
Proposició 4.3.4. Si x ∈ R y f lT (x) y f lR (x) son las representaciones de punto
flotante (normalizadas) con t dígitos y base β, entonces:
Er (f lT (x)) ≤ β 1−t , (4.5)
1 1−t
Er (f lR (x)) ≤ ·β . (4.6)
2
Prueba. Sea x = (.d1 d2 . . . dt dt+1 dt+2 . . .)β · β e la representación de punto flotante
normalizada de x, entonces f lT (x) = (.d1 d2 . . . dt )β · β e . Por tanto,
T |x − f lT (x)| (.dt+1 dt+2 . . .)β · β e−t

Er (f l (x)) = = .
|x| (.d1 d2 . . .)β · β e
Como (.dt+1 dt+2 . . .)β ≤ 1 y (.d1 d2 . . .)β ≥ β −1 (ya que d1 ≥ 1), tenemos
β e−t
Er (f lT (x)) ≤ −1 e
= β 1−t .
β ·β
Veamos ahora el caso del redondeo. Si 0 ≤ (.dt+1 dt+2 . . .)β < 1/2, tenemos
1 e−t
|x − f lR (x)| = (.dt+1 dt+2 . . .)β · β e−t < ·β ,
2
por tanto,
|x − f lR (x)| 1/2β e−t 1
Er (f lR (x)) = ≤ −1 e = · β 1−t .
|x| β ·β 2
En cambio, si 1/2 ≤ (.dt+1 dt+2 . . .)β ≤ 1,
f lR (x) = (d1 β −1 + . . . + (dt + 1)β −t ) · β e ,
y, por tanto,
+∞ t−1
!
X X
R −k −k −t
|x − f l (x)| = dk β − dk β + (dt + 1)β · βe
k=1 k=1
+∞
X
= −β −t + dk β −k · β e
k=t+1
1 e−t
= |1 − (.dt+1 dt+2 . . .)β | · β e−t ≤ ·β ,
2
de donde se deduce de nuevo
1/2β e−t 1
Er (f lR (x)) ≤ −1 e
= · β 1−t .
β ·β 2
Corolario 4.3.5. Si x ∈ R y f l(x) es una representación de punto flotante (normal-
izada) con t dígitos y base β, entonces
f l(x) = x(1 + δ) con |δ| ≤ u, (4.7)
donde u se llama unidad de redondeo, o también precisión máquina. Si uT y uR son

las unidades de redondeo correspondientes a las representaciones de punto flotante
por corte y redondeo respectivamente, tenemos que
1 1−t
uT = β 1−t uR = ·β .
2
Prueba. Si x 6= 0, tenemos
f l(x) − x
f l(x) = f l(x) − x + x = x + x = x(1 + δ),
x
donde δ = (f l(x) − x)/x. Si x = 0, f l(x) = 0 y no hay que probar nada.
Hace treinta años, la situación era mucho mas complicada que ahora. Cada orde-
nador tenia su propio sistema de aritmética de punto flotante. Algunos eran binarios;
otros eran decimales. Incluso había un ordenador Ruso que usaba aritmética ternaria
(β = 3). Entre los ordenadores binarios, algunos usaban 2 como base y otros usa-
ban 8 o 16. Y cada uno tenía diferente precisión. En 1985, el IEEE Standars Board y el
American National Standards Institute adoptaron el ANSI/IEEE Standard 754-1985 de

aritmética de punto flotante binaria. Fue la culminación de casi una década de trabajo
hecho por un grupo de 92 matemáticos e ingenieros que trabajaban en la Universidad
o en compañias de ordenadores.
Todos los ordenadores construidos desde 1985 utilizan la aritmetica de punto flotante
IEEE.
El modelo de aritmética del IEEE es cerrado, en el sentido de que cada operación
aritmética tiene un resultado, con señales especiales para operaciones especiales. Por la
definición de F tenemos situaciones en las que x > máx{|y| : y ∈ F } o x < mı́n{|y| :
y ∈ F }. En el primer caso se habla de overflow y en el segundo de underflow. Las
situaciones y operaciones especiales y sus resultados pueden consultarse en la tabla
4.2.
Tipo de excepciones Ejemplo √ Resultado por defecto

Operación inválida 0/0, 0 × ∞, −1 NaN (Not a Number)
Overflow ±∞
Divide por zero No zero finito/0 ±∞
Underflow números sin sentido
Inexacto Si f l(x w y) 6= x w y resultado redondeado correctamente
Cuadro 4.2: Excepciones de la aritmética IEEE i resultados por defecto.
El estándar del IEEE especifica dos formatos básicos para representaciones de pun-
to flotante, formatos en simple precisión y formatos en doble precisión. Las especifi-
caciones se dan en la tabla siguiente:
Tipo Longitud Mantisa Exponente Unidad redondeo Rango

−24 −8
Simple 32 bits 23 + 1 bits 8 bits 2 ≈ 5,96 × 10 10±38
−53 −16
Doble 64 bits 52 + 1 bits 11 bits 2 ≈ 1,11 × 10 10±308
MATLAB ha utilitzado tradicionalmente el formato IEEE en doble precisión. Con

la precisión simple no se es mucho más rápido en los ordenadores modernos.
Como β = 2, los dígitos solo pueden ser 0 o 1. Teniendo esto en cuenta, idearon
una forma de ganar un bit de información. Así, en vez de escribir (4.3) o (4.4), expresan
los números en la forma normalizada
f l(x) = ±(1.d1 d2 . . . dt )2 · 2e = ±(1 + m) · 2e (4.8)
con t = 23 o t = 52 dependiendo de la precisión (simple o doble). La cantitad m es la

mantisa y e es el exponente. La mantisa satisface:
0≤m<1
y debe poder ser representable en binario utilizando como mucho t = 23 o t = 52 bits.

En otras palabras, 2t m es un entero en el intervalo:
0 ≤ 2t m < 2t .
El exponente e es un entero en el intervalo
−126 ≤ e ≤ 127 (o − 1022 ≤ e ≤ 1023).
La finitud de m es una limitación en la precisión. La finitud de e es una limitación
en el rango. Cualquier número que no se pueda expresar así ha de ser aproximado por
números que si lo puedan ser.
Los números de aritmética en punto flotante y doble precisión están almacenados
en 64 bits, con 52 bits para m, 11 bits para e, y un bit para el signo del número. El
signo de e se recupera almacenando e + 1023, que está entre 1 y 211 − 2. Los dos
extremos, 0 y 211 − 1, se reservan para números excepcionales en aritmética de punto
flotante (NaN, 0, Inf , etc.).
La parte fraccionaria completa de un número en punto flotante no es, m, es 1 + m,
que tiene 53 bits. Como hemos dicho antes, el primer 1 no necesita ser almacenado.
Concretamente, el formato IEEE tiene 65 bits de información solo almacenando 64
bits.
El número más grande que tendríamos en doble precisión sería:
52
X 1
· · 1} ·21023 = (1 +
x = 1. |1 ·{z ( )i ) · 21023
i=1
2
52
1
(1− 2152 )
= (1 + 2
) · 21023
1 − 12
1
= (2 − ) · 21023
252
≈ 1,797693134862316 · 10308
El número positivo normalizado más pequeño sería:
· · 0} ·2−1022 = 1 · 2−1022
x = 1. |0 ·{z
52
≈ 2,225073858507201 · 10−308 .
El hecho de no haber utilizado en los números normalizados los exponentes −1023
y 1024 (o 0 y 211 − 1 según se guardarían) se utiliza para guardar números excep-
cionales (NaN, 0, Inf ) y también nos permite guardar números desnormalizados (que
no son de la forma 1 + m).
Así, si el exponente son todo 0s y la mantisa no es zero (en otro caso sería la
representación del número zero) entonces el valor es un número desnormalizado. Rep-
resentaría un número de la forma (−1)s × (0 + m) × 2−1022 , donde s es el bit de signo
(0 o 1) y m la mantisa.
4.4. CONSECUENCIAS DE LA ARITMÉTICA DE PRECISIÓN FINITA 108
De esta forma, el número positivo más pequeño que podemos representar en el

ordenador sería el 2−52 × 2−1022 = 2−1074 = 4,940656458412465 · 10−324 .
Los números positivos más pequeños que 21074 no pueden ser representados (un-
derflow positivo). Los números positivos más grandes que (2 − 2−52 ) · 101023 tampoco
(overflow positivo).
El problema del underflow es menor ya que significa una pérdida de precisión pero
tenemos garantizado que el número está muy cerca de zero.
El estándar del IEEE también especifica un modelo de aritmética para las repre-
sentaciones de punto flotante: todas las operaciones aritméticas se han de hacer como
si se efectuasen primero con precisión infinita y después redondeásemos el resultado,
es decir, si x, y son dos números de punto flotante, w es cualquiera de las operaciones
elementales (+, −, ·, /) y w ∗ denota la operación correspondiente en punto flotante,
entonces
x w ∗ y = f l(x w y),
por tanto, teniendo en cuenta el corolario 4.3.5,
x w ∗ y = f l(x w y) = (x w y)(1 + δ), |δ| ≤ u, (4.9)
donde u es la unidad de redondeo. La expresión f l(x w y) se utiliza habitualmente para

hacer referencia al resultado de la operación en precisión finita. El modelo estándar
dice que el valor calculado para x w y es tan bueno como el valor exacto, redondeado,
en el sentido que la cota del error relativo es la misma. No obstante, el modelo no
obliga a que δ = 0 cuando x w y ∈ F .
En MATLAB tenemos
>> (2-1/2^52)*2^1023+1
ans =
1.797693134862316e+308
que es exactamente (2-1/2^52)*2^1023.
Las consecuencias del modelo de aritmética que utiliza el ordenador en el rendimien-
to de los algoritmos numéricos son muchas y de naturaleza muy variada. Aquí solo
harblaremos de dos de las consecuencias de la aritmética de precisión finita: la can-
celación catastrófica y la inestabilidad numérica.
4.4. Consecuencias de la aritmética de precisión finita

Hay que observar que la aritmética de punto flotante no es exacta. Aunque la suma
de dos números de punto flotante es un número de punto flotante, puede darse el caso
que la longitud de palabra de los términos de la suma y del resultado (ambas cosas con
precisión infinita) no sea la misma.
Algunas consecuencias de la aritmética de precisión finita que hay que tener en
cuenta son las seguientes:
1. La suma o la resta de números dispares pude no tener efecto sobre el mayor. Así,
con t = 3 y β = 10:
(,518) · 103 +∗ (,437) = 518 +∗ ,437 = f l(518,437) = ,518 · 103 .
Vemos, pues, que x +∗ y = x, con y 6= 0, es decir, todos los dígitos de y se han

perdido en la suma final, a causa de la longitud de palabra utilizada.
2. Con t = 3, si a = 3, 1/a = 0.3̂ y f l(1/a) = 0,333. De esta forma, f l(a) .∗ f l(1/a) =

,999 6= 1.
3. En general, (a +∗ b) +∗ c 6= a +∗ (b +∗ c). Por ejemplo, si a = 63,1,

b = 4,24 y c = ,247, mientras que (a +∗ b) +∗ c = (,675) · 102 , tenemos que
a +∗ (b +∗ c) = (,676) · 102 .
Los errores generados por los cálculos en aritmética de precisión finita, pueden
adquirir una importancia relativa considerable. Por poner un ejemplo sencillo, supong-
amos ahora que queremos calcular z = x − y con x = 7,6545421, y = 7,6544200
en un ordenador que trabaja con t = 6, por tanto, disponemos de aproximaciones
x̃ = (,765454) · 101 , ỹ = (,765442) · 101 para las cuales
1
Er (x̃) ≤= · 10−5 = 5 · 10−6 , Er (ỹ) ≤ 5 · 10−6 .
2
El valor
z̃ = x̃ − ỹ = (0,000012) · 101 = (0,120000) · 10−3
es una aproximación a z = ,00001221 para la que se cumple
Ea (z̃) = |z − z̃| = |x − x̃ + y − ỹ| ≤ Ea (x̃) + Ea (ỹ).
No obstante, si calculamos el error relativo, tenemos que
21
Er (z̃) = ≈ 1,7 · 10−2 .
1221
Este incremento de casi cuatro órdenes de magnitud en el error relativo indica una
disminución considerable del número de dígitos significativos correctos en el resul-
tado de la diferencia de las dos aproximaciones. Este fenómeno, que aparece cuando
se restan números muy parecidos, se suele llamar cancelación catastrófica. La can-
celación catastrófica solo puede ser evitada si anticipamos su posible aparición.
Ejemplo 4.4.1. Queremos calcular las raices de la ecuación ax2 + bx + c = 0 uti-

lizando la fórmula √
−b ± b2 − 4ac
x= .
2a
Supongamos que b > 0, la raiz

√
−b + b2 − 4ac
x1 = (4.10)
2a
puede padecer el efecto de la cancelación catastrófica cuando b2 −√4ac ≈ b2 , ya
que nos veremos obligados a restar dos números muy parecidos (b y b2 − 4ac). La
alternativa para evitarlo es bastante sencilla:
√ √ √
−b + b2 − 4ac 1 (−b + b2 − 4ac)(−b − b2 − 4ac)
x1 = = √
2a 2a −b − b2 − 4ac
1 (−b)2 − (b2 − 4ac) 2c
= √ = √ . (4.11)
2a −b − b2 − 4ac −b − b2 − 4ac
En esta última fórmula, que analíticamente es equivalente a la anterior, (4.10), no hay
diferencias de números parecidos y, por lo tanto, su utilización para el cálculo de x1 ,
en el caso de que b2 − 4ac ≈ b2 , es más recomendable.
Ejemplo 4.4.2. Calculad las raices de ax2 + bx + c = 0 trabajando con aritmética de

punto flotante con β = 10 y t = 5 si a = 1, b = 5,4321 · 104 y c = 0,1. Comparad las
soluciones que obtenemos con la solución exacta x1 = −1,8409 · 10−6 .
Primero utilizamos la fórmula (4.10):
b2 → (5,4321 · 104 )2 = 2,9508 · 109

b2 − 4ac → 2,9508 · 109 − 4 ∗ c = 2,9508 · 109
√ p
b2 − 4ac → 2,9508 · 109 = 5,4321 · 104
√
−b + b2 − 4ac → (−5,4321 + 5,4321) · 104 = 0
√
(−b + b2 − 4ac)/2a → 0/(2) = 0.
Y ahora utilizando la fórmula (4.11):
b2 → (5,4321 · 104)2 = 2,9508 · 109

b2 − 4ac → 2,9508 · 109 − 4 ∗ c = 2,9508 · 109
√ p
b2 − 4ac → 2,9508 · 109 = 5,4321 · 104
√
−b − b2 − 4ac → (−5,4321 − 5,4321) · 104 = −,10864 · 106
√
2c/(−b − b2 − 4ac) → 2 ∗ 0,1/(−,10864 · 106 ) = −1,8409 · 10−6 .
Ejercicio 4.4.3. ¿Qué expresión utilizaríamos cuando b < 0?
El efecto de un error individual sobre el resultado de evaluar una expresión f (x) en

un punto concreto x se puede estimar a priori si la expresión es suficientemente suave
y conocemos el error en la aproximación a x, x̃. Si utilizamos desarrollos de Taylor,

podemos escribir:
1
f (x) = f (x̃) + f 0 (x̃)(x − x̃) + f 00 (x̃)(x − x̃)2 + . . .
2
00
Así, si (x − x̃) es suficientemente pequeño y f no es muy grande, tenemos
f (x) − f (x̃) ≈ f 0 (x̃)(x − x̃)
y, por tanto,
|f (x) − f (x̃)| |f 0 (x̃)| |x − x̃|
≈ |x| ,
|f (x)| |f (x)| |x|
de donde
|f 0 (x̃)|
Er (f (x̃)) ≈ |x| Er (x̃). (4.12)
|f (x)|
El error relativo en la evaluación de la expresión f (x) en un valor x̃ ≈ x está deter-
minado por el error relativo de la aproximación, Er (x̃) y el llamado factor de amplifi-
cación de la expresión f (x) en el punto x,
|f 0(x̃)|
Af = |x|.
|f (x)|
La relación (4.12) es muy útil a la hora de decidir si una expresión concreta es

adecuada para los cálculos con precisión finita, o si, al contrario, hay que manipularla
analíticamente y transformarla en otra
√ expresión más adecuada para los cálculos prác-
ticos. Por ejemplo, para f (x) = x y x̃ ≈ x, el factor de amplificación del error
es:
|x| 1
√ √ ≈ ,
2 x x̃ 2
es decir, el error relativo se reduce a la mitad y obtener raices cuadradas es una op-
eración computacionalmente segura.
Así mismo, la expresión
10
f (x) =
1 − x2
presenta problemas a causa de la cancelación catastrófica cerca de x = 1. Si calculam-
os su factor de amplificación, tenemos
20x
|f 0(x)||x| | (1−x 2 )2 x| 2x2
= 10 = 1 cuando x ≈ 1.
|f (x)| |1−x2 |
|1 − x2 |
Finalmente, mostraremos ahora cómo utilizar (4.12) para decidir entre diversas
expresiones alternativas.
4.5. INESTABILIDAD NUMÉRICA 112
√
Ejemplo 4.4.4. Para calcular ( 2−1)6 se puede utilizar cualquiera de las expresiones
siguientes (se deja como ejercicio probar que todas son equivalentes):
1 √ 1 √ 1
√ ; (3 − 2 2)3 ; √ ; 99 − 70 2 i √ .
( 2 + 1)6 (3 + 2 2)3 99 + 70 2
Queremos calcular el resultado de esta operación con aritmética de precisión finita
con β = 10 y t = 5. Se trata de decidir cuál de todas las fórmulas anteriores es
computacionalmente más adecuada. Definimos
1
f0 (x) = (x − 1)6 ; f1 (x) = ; f2 (x) = (3 − 2x)3 ;
(x + 1)6
1 1
f3 (x) = ; f4 (x) = 99 − 70x; f5 (x) = ,
(3 + 2x)3 99 + 70x
y evaluamos el factor de amplificación que corresponde a cada fi ,
|fi0 (x̃)| √
Ai = |x|, x, x̃ ≈ 2.
|fi (x)|
√
Como 2 = 1,41421356237310, si trabajamos con t = 5 (y redondeo), tenemos
x̃ = 1,4142 con |Er (x̃)| ≈ 5 · 10−5 y, por tanto, |Er (f (x̃))| ≈ Ai · 5 · 10−5 . Los
cálculos con todas las expresiones se recojen en la tabla siguiente:
i Ai fi (1,4142) Error absoluto Error relativo
0 2,05 · 101 5,0496 · 10−3 1,03 · 10−6 2,05 · 10−4
1 3,51 · 100 5,0508 · 10−3 1,66 · 10−7 3,29 · 10−5
2 4,95 · 101 5,0530 · 10−3 2,37 · 10−6 4,68 · 10−4
3 1,46 · 100 5,0508 · 10−3 1,66 · 10−7 3,29 · 10−5
4 1,96 · 104 6,0000 · 10−3 9,49 · 10−4 1,88 · 10−1
5 5,0 · 10−1 5,0508 · 10−3 1,66 · 10−7 3,29 · 10−5
Los factores de amplificación (segunda columna de la tabla) predicen que el resul-

tado menos fiable, por ser A4 excesivamente grande, será f4 (x̃), mientras que el más
fiable ha de ser, f5 (x̃), seguido de f3 (x̃) y f1 (x̃). Todas las predicciones se comprueban
examinando las columnas 4 i 5 de la tabla.
4.5. Inestabilidad numérica

La propagación descontrolada de los errores puede llegar a invalidar los resultados
de una secuencia de cálculos. Si llamamos Rn () al crecimiento de un error inicial de-
spués de n operaciones y éste se comporta de la forma Rn () ≈ C ·n·, donde C es una
constante independiente de n, diremos que el crecimiento es lineal. Este crecimiento
es normal, usualmente inevitable, y no peligroso. Por otra parte, si Rn () ≈ k n · ,
donde k > 1 es una constante independiente de n, diremos que el crecimiento es ex-

ponencial. Este tipo de crecimiento puede llegar a ser desastroso e incluso producir
resultados numéricos absurdos. Los algoritmos que presentan este tipo de crecimiento
se dicen numéricamente inestables y se han de sustituir por otros que no presenten este
problema.
Ejemplo 4.5.1. Queremos calcular los términos de la successión
Z 1
xn
yn = , n = 0, 1, 2, . . . (4.13)
0 x+5
No es difícil probar que yn ≥ 0, ya que xn /(x + 5) ≥ 0 para 0 < x < 1. Además

se trata de una sucesión decreciente (yn+1 < yn para todo n). Por tanto, la sucesión
formada por las yn tiene límite cuando n → ∞ y no es difícil probar que, de hecho,
yn → 0.
Para calcular los términos yn sin hacer ninguna integral, se puede utilizar la fór-
mula recurrente siguiente (la prueba se deja como exercicio).
1
yn+1 + 5yn = . (4.14)
n+1
El primer término es necesario para iniciar el proceso recurrente,
Z 1
1
y0 = dx = log(x + 5)]10 = log(6) − log(5) = log(6/5). (4.15)
0 x+5
i simple (4.14) doble (4.14) simple (4.16) sol. exacta yi

0 1,8232 · 10−1 1,8232 · 10−1 1,8232 · 10−1 1,8232 · 10−1
1 8,8392 · 10−2 8,8392 · 10−2 8,8392 · 10−2 8,8392 · 10−2
2 5,8040 · 10−2 5,8039 · 10−2 5,8039 · 10−2 5,8039 · 10−2
3 4,3134 · 10−2 4,3139 · 10−2 4,3139 · 10−2 4,3139 · 10−2
4 3,4329 · 10−2 3,4306 · 10−2 3,4306 · 10−2 3,4306 · 10−2
5 2,8355 · 10−2 2,8468 · 10−2 2,8468 · 10−2 2,8468 · 10−2
6 2,4893 · 10−2 2,4325 · 10−2 2,4325 · 10−2 2,4325 · 10−2
7 1,8390 · 10−2 2,1233 · 10−2 2,1233 · 10−2 2,1233 · 10−2
8 3,3048 · 10−2 1,8837 · 10−2 1,8837 · 10−2 1,8837 · 10−2
9 −5,4131 · 10−2 1,6926 · 10−2 1,6926 · 10−2 1,6926 · 10−2
10 3,7066 · 10−1 1,5368 · 10−2 1,5368 · 10−2 1,5368 · 10−2
11 −1,7624 · 100 1,4071 · 10−2 1,4071 · 10−2 1,4071 · 10−2
12 8,8952 · 100 1,2977 · 10−2 1,2977 · 10−2 1,2977 · 10−2
13 −4,4399 · 101 1,2040 · 10−2 1,2040 · 10−2 1,2040 · 10−2
14 2,2207 · 102 1,1229 · 10−2 1,1229 · 10−2 1,1229 · 10−2
15 −1,1103 · 103 1,0522 · 10−2 1,0521 · 10−2 1,0521 · 10−2
16 5,5514 · 103 9,8911 · 10−3 9,8963 · 10−3 9,8963 · 10−3
17 −2,7757 · 104 9,3678 · 10−3 9,3419 · 10−3 9,3419 · 10−3
18 1,3878 · 105 8,7166 · 10−3 8,8462 · 10−3 8,8462 · 10−3
19 −6,9392 · 105 9,0487 · 10−3 8,4005 · 10−3 8,4005 · 10−3
20 3,4696 · 106 4,7565 · 10−3 7,9975 · 10−3 7,9975 · 10−3
21 −1,7348 · 107 2,3836 · 10−2 7,6314 · 10−3 7,6314 · 10−3
22 8,6740 · 107 −7,3728 · 10−2 7,2974 · 10−3 7,2974 · 10−3
23 −4,3370 · 108 4,1212 · 10−1 6,9913 · 10−3 6,9913 · 10−3
24 2,1685 · 109 −2,0189 · 100 6,7099 · 10−3 6,7099 · 10−3
25 −1,0843 · 1010 1,0135 · 101 6,4503 · 10−3 6,4503 · 10−3
Cuadro 4.3: Cálculo de yn en (4.13) con precisión finita. Segunda columna: con (4.14)
y simple precisión. Tercera columna: (4.14) y doble precisión. Cuarta columna: con
(4.16) e y33 = 0 como valor inicial, con simple precisión. La última columna es la
solución exacta.
Las columnas 2 y 3, en la tabla 4.3, muestran los resultados de la aplicación de

un algoritmo recurrente basado en (4.14) obtenidos con programas (en Fortran77)
que utilizan precisión simple (columna 2) y doble (columna 3). Cuando trabajamos
con precisión simple, se observa que el comportamiento decreciente de la sucesión de
valores numéricos se pierde en y8 (y8 > y7 ), lo que indica una acumulación anormal
de los errores. Además, y9 < 0 es un resultado numérico absurdo ya que todos los yn
han de ser positivos.
En los cálculos con doble precisión la acumulación anormal de errores produce
inconsistencias a partir de y19 (y19 > y18 ) con resultados absurdos a partir de y22 .
Este proceso recursivo es un ejemplo típico de algoritmo numéricament inestable:
la propagación de los errores en los cálculos acaba por arruinar completamente los
resultados numéricos y obtenemos resultados absurdos.
Este efecto puede ser analizado teóricamente estudiando el efecto de un error con-
creto, por ejemplo el error de redondeo cometido al almacenar el dato inicial y0 , en el
resto del proceso recursivo.
Inicialmente comenzamos con y0∗ = f l(y0 ). Así, si suponemos que las operaciones
siguientes no están afectadas por errores (claramente una idealización), tenemos
y0∗ = y0 +
1
y1∗ = − 5y0∗ = 1 − 5(y0 + ) = y1 − 5
1
1
y2∗ = − 5y1∗ = y2 + 52 ,
2
y es una tarea sencilla probar (por inducción) que
yn∗ = yn + (−1)n 5n .
El efecto del proceso recursivo hace que |Rn ()| = 5n , es decir, produce un crec-
imiento exponencial del error que acaba por invalidar los resultados obtenidos. Es
interesante notar que podemos predecir bastante aproximadamente en qué momento
los resultados numéricos comienzan a ser absurdos. En efecto, si trebajamos con pre-
cisión simple, uR ≈ 6 · 10−8 (mirad la tabla de estándards de l’IEEE), por tanto, el
error inicial verifica
|| = |y0∗ − y0 | = |f l(y0 ) − y0 | ≈ uR |y0 | ≈ 10−8.
Así, el error arrastrado después de k aplicaciones del proceso recursivo, 5k ||, tiene
el valor siguiente:
k 5 6 7 8 9 10
.
k
5 || 3 · 10−5 1,6 · 10−4 8 · 10−4 4 · 10−3 2 · 10−2 10−1
Observad que el error arrastrado después de 9 aplicaciones, ≈ 0,02, es ya mayor

que la solución exacta y9 y, por tanto, el valor calculado y9∗ no tiene ningún dígito
correcto. Con estas observaciones vemos que es posible predecir el índice a partir del
cual los resultados numéricos comienzan a notar el efecto de la inestabilidad numérica
del algoritmo.
Ejercicio 4.5.2. Utilizad el mismo mecanismo para predecir a partir de qué índtce se
comenzarán a notar los efectos del crecimiento exponencial del error arrastrado por
el algoritmo recursivo cuando se utiliza doble precisión. Comprobadlo con la tabla
4.3.
Estas observaciones nos permiten concluir que el proceso recursivo (4.14) no es

válido pera los cálculos con aritmética de precisión finita: es un ejemplo de proceso
numéricamente inestable. Para calcular yn sin hacer integrales, hay que diseñar otro
proceso que sea numéricamente estable. En este caso el análisis de la propagación del
error que hemos hecho antes nos indica que una buena estrategia, por lo que se refiere
a la propagación del error, es considerar la relación de recurrencia al revés, es decir,

1 1
yn = − yn+1 . (4.16)
5 n+1
Supongamos que se quiere calcular yn para un n fijado, por ejemplo y25 . Co-
mo la relación de recurrencia va hacia atrás, necesitamos comenzar por yn+k , con k
suficientemente grande. La determinación de un k apropiado se puede hacer con el
análisis del arrastre del error. Sabemos que yn → 0, por tanto, para k suficientemente
∗
grande yn+k ≈ 0, y es razonable considerar yn+k = 0 como comienzo del nuevo pro-
∗
ceso iterativo. Con esta aproximación tendremos un error yn+k = yn+k + que, con
la estrategia (4.16), se propaga de la manera siguiente:

∗ 1 1 ∗
yn+k−1 = − yn+k = yn+k−1 −
5 n+k 5
∗
yn+k−2 = yn+k−2 + 2
5
..
.

yn∗ = yn + (−1)k .
5k
Como |yn+k | = || y |yn+k | < |yn |, tenemos
|yn∗ − yn | /5k |yn+k | 1 1

Er (yn∗ ) = = = k
< k.
|yn | |yn | |yn | 5 5
De acuerdo con estas estimaciones, si queremos una aproximación a yn de mane-

ra que |Er (yn∗ )| ≤ δ, podemos utilizar el proceso recurrente (4.16) comenzando con
yn+k = 0, donde k verifica 1/5k ≤ δ, es decir, hay que elegir
k > [− log δ/ log 5].
Por ejemplo, si δ = 10−5, k = 8 es el menor entero tal que 1/5k < 10−5 . Así, si
comenzamos con y25+8 = 0, calcularemos y25 con un error relativo menor que 10−5 .
Está claro que todos los otros valores yi , i = 24, . . . , 0 también tendrán, por lo menos,
esta exactitud. Los resultados del proceso recurrente con esta estrategia de comienzo
se recojen en la columna 4 de la tabla 4.3. Los resultados con precisión simple o con
doble coinciden con los dígitos que aparecen en la tabla.
4.6. PROBLEMAS 117
4.6. Problemas
Problema 1. Discutid la estabilidad numérica de las expresiones siguientes y dad una
alternativa, si es preciso:
1 1−x
1. − , x ≈ 0.
1 + 2x 1 + x
1 − cos x
2. , x ≈ 0.
x
3. ex − e−x , x ≈ 0.
Problema 2. El año 250 AC, Arquimedes estimó el número π de la manera siguiente:
tomando el círculo de diámetro 1, por tanto, de longitud de la circunferencia π, in-
scribió dentro un cuadrado. El perímetro del cuadrado es menor que la longitud de la
circunferencia y, por tanto, una cota inferior para π. Arquimedes consideró después un
octágono, un 16-ágono, etc., y dobló cada vez el número de lados del polígono inscrito
produciendo cada vez mejores aproximaciones a π. Utilizando polígonos inscritos y
circunscritos de 96 lados pudo demostrar que:
223 22
<π< .
71 7
Si pn es el perímetro del polígono inscrito con 2n lados, hay una fórmula de recur-
rencia para estas estimaciones:
q p
n
pn+1 = 2 2(1 − (1 − (pn /2n )2 ). (4.17)
1. ¿Es éste un procedimiento estable para calcular π? ¿Porqué?
2. ¿Podeis diseñar un procedimiento alternativo mejor?
3. Comprobad numéricamente las afirmaciones que haceis.
Problema 3. Queremos calcular las integrales de la forma:

Z 1
yn = ex xn dx, n = 0, 1, 2, . . . ,
0
utilizando la fórmula recurrente:
y0 = e − 1, (4.18)
yn+1 + (n + 1)yn = e. (4.19)
1. Deducid la fórmula recurrente (4.19).
2. Probad que lı́mn→∞ yn = 0 .

4.6. PROBLEMAS 118
3. Confeccionad un programa de ordenador, que calcule y0 , y1, . . . , y30 utilizando

la fórmula recurrent (4.19). Comentad los resultados obtenidos con el orde-
nador.
4. Poned yn en función de yn+1 en la fórmula de recurrencia (4.19). Utilizando la

nueva fórmula recurrente calculad y5 , y10 , y15, y20 , comenzando la recurrencia
con ym = 0 con m suficientemente grande para que el error inicial no tenga
ningún efecto sobre el resultado que se busca.
5. Justificad los resultados numéricos.
Problema 4. La sucesión bn se define de manera recurrente como sigue:
b1 = 1, (4.20)
p
1 + b2n − 1
bn+1 = . (4.21)
bn
1. Provad que lı́mn→∞ bn = 0.
2. ¿Es éste un procedimiento numéricament estable? Justificad la respuesta com-

putacionalmente calculando b30 .
√
5−1
Problema 5. Llamamos ψ = 2
y queremos calcular las sucesivas potencias de ψ.
1. Provad que ψ n cumple la relación de recurrencia:
ψ n+1 = ψ n−1 − ψ n . (4.22)
2. Confeccionad un programa de ordenador para calcular las potencies de ψ a

partir de los valores iniciales ψ 0 = 1, ψ 1 = 0,61803398, mediante la relación
de recurrencia (4.22).
3. A partir de los resultados que se obtienen con el código numérico, ¿dirias que
éste es un algoritmo estable?

Apuntes HI-1

Cargado por

Copyright:

Formatos disponibles

Apuntes HI-1

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Apuntes HI-1

Cargado por

Copyright:

Formatos disponibles

Herramientas Informáticas

1. Algoritmos básicos en MATLAB 3

2. Solución de ecuaciones no lineales 38

4. Sistemas numéricos y fuentes de error 96

Algoritmos básicos en MATLAB

1.1. Programación en MATLAB

1.1.2. Introducción de matrices

Introduciendo una lista explícita de elementos,

Generándola mediante funciones y declaraciones,

Creándola en un archivo .m (ver secciones 1.1.8 y 1.1.10)

Cargándola de un archivo de datos externo (lo veremos más adelante).

Por ejemplo, cualquiera de las declaraciones

1.1.3. Operaciones con matrices, operaciones a coordenadas

La división matricial merece un comentario especial. Si A es una matriz invertible

Inténtelo. Esto es particularmente útil cuando se utilizan los gráficos de MATLAB.

1.1.4. Declaraciones, expresiones y variables; almacenamiento de

muestra en pantalla, y se asigna a la variable para su posterior uso. Si se omiten la

Your variables are:

obj plotf1 pn1 prod1 prod3 tra

1.1.5. Formato de salida

1.1.6. Funciones en MATLAB

Las siguientes funciones están disponibles en MATLAB:

eye matriz identidad

>> B = [A, zeros(3,2); zeros(2,3), eye(2)]

dará una cierta matriz 5 × 5. Inténtelo.

Algunas funciones de MATLAB operan esencialmente sobre escalares, aunque lo

sin asin exp abs round

max sum median any

Por ejemplo, la entrada máxima de un matriz A se obtiene con max(max(A)) en vez

Las funciones de MATLAB admiten argumentos de salida simples o múltiples. Por

1.1.7. For, while, if — y relaciones

>> x = []; for i = 1:n, x=[x,i^2], end

darán como resultado un cierto vector mientras que

>> x = []; for i = n:-1:1, x=[x,i^2], end

>> for i = 1:m

producirán e imprimirán en pantalla la matriz de Hilbert m × n. El punto y coma de

if. La forma general de un bucle if simple es

Cuando se aplican a escalares, una relación es realmente el escalar 1 o 0 dependi-

Ejercicio 1.1.2. Ver los resultados obtenidos con

>> 3 < 5, 3 > 5, 3 == 5, 3 == 3

Ejercicio 1.1.3. Ver el resultado obtenido con

Recalcamos que la aparentemente obvia

1.1.8. Submatrices y notación de dos puntos

1.1.9. Comparación de la eficiencia de algoritmos: etime, cputime,

1.1.10. Cadenas de texto, mensajes de error, input

entonces la instrucción de MATLAB datos hará que se efectúe la asignación especi-

Esto debe escribirse en el archivo ental.m (correspondiente al nombre de la fun-

function [media, desv] = estad(x)

Una vez situado en el archivo de disco estad.m, la instrución de MATLAB

En MATLAB, dados dos vectores v i u de dimensión n × 1 es equivalente, calcular

x = -4:.01:4; y = sin(x); plot(x,y)

x = -1.5:.01:1.5; y = exp(-x.^2); plot(x,y)

t=0:.001:2*pi; x=cos(3*t); y=sin(2*t); plot(x,y)

La instrucción grid hará un cuadriculado en el gráfico actual.

title(’La funci\’on m\’as bella’)

proporciona un título al gráfico. El comando gtext(’La mancha’) permite posi-

y formando una matriz Y conteniendo los valores funcionales como columnas

>> x=0:.01:2*pi; Y=[sin(x)’, sin(2*x)’, sin(4*x)’];

>> x=0:.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x);

t=0:.001:2pi; x=cos(3t); y=sin(2*t); plot(x,y)

>> x=0:.01:2pi; Y=[sin(x)’, sin(2x)’, sin(4*x)’];

>> x=0:.01:2pi; y1=sin(x); y2=sin(2x); y3=sin(4*x);