APUNTESMATLAB
APUNTESMATLAB
B
SE
UNIVER
VILLA
e an
Ecuaciones Diferenciales
y Análisis Numérico
Apuntes de MATLAB
orientados a métodos numéricos elementales
Rosa Echevarría
Dpto. de Ecuaciones Diferenciales y Análisis Numérico
Universidad de Sevilla
14 de febrero de 2022
DAD DE
SI
SE
UNIVER
VILLA
Índice e an
Ecuaciones Diferenciales
y Análisis Numérico
3
2.3. Operaciones aritméticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.3.1. Reglas de prioridad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.4. Forma en que se muestran los resultados . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.5. Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.5.1. Instrucciones de asignación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.5.2. Variables pre-definidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.6. Funciones matemáticas elementales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.7. Operaciones de comparación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.8. Operadores lógicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.9. Orden general de evaluación de expresiones . . . . . . . . . . . . . . . . . . . . . . . . 56
2.10. Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.10.1. Construcción de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.10.2. Acceso a los elementos de vectores y matrices . . . . . . . . . . . . . . . . . . . 59
2.10.3. Operaciones con vectores y matrices . . . . . . . . . . . . . . . . . . . . . . . . 62
2.10.4. Operaciones de comparación y lógicas con matrices . . . . . . . . . . . . . . . . 64
2.11. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3. Gráficos 71
3.1. Generalidades sobre la representación gráfica de funciones . . . . . . . . . . . . . . . . 71
3.1.1. Representación gráfica de funciones de una variable real . . . . . . . . . . . . . 72
3.1.2. Curvas planas definidas mediante ecuaciones paramétricas . . . . . . . . . . . . 73
3.1.3. Curvas planas en coordenadas polares . . . . . . . . . . . . . . . . . . . . . . . 73
3.1.4. Curvas en tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.1.5. Gráficas de funciones de dos variables: superficies . . . . . . . . . . . . . . . . . 76
3.1.6. Superficies definidas mediante ecuaciones paramétricas . . . . . . . . . . . . . . 77
3.1.7. Representación mediante curvas de nivel de una función de dos variables . . . . 78
3.2. Funciones gráficas de MATLAB fáciles de usar . . . . . . . . . . . . . . . . . . . . . . 80
3.2.1. Gráficas de funciones de una variable real . . . . . . . . . . . . . . . . . . . . . 80
3.2.2. Curvas planas definidas implícitamente . . . . . . . . . . . . . . . . . . . . . . . 81
3.2.3. Curvas planas definidas por ecuaciones paramétricas . . . . . . . . . . . . . . . 82
3.2.4. Curvas en tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.2.5. Superficies definidas mediante ecuaciones explícitas . . . . . . . . . . . . . . . . 84
3.2.6. Superficies definidas mediante ecuaciones paramétricas . . . . . . . . . . . . . . 86
3.2.7. Representación mediante curvas de nivel de una función de dos variables . . . . 87
3.3. Funciones de bajo nivel para el dibujo de curvas . . . . . . . . . . . . . . . . . . . . . . 89
3.3.1. La orden plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.3.2. Color y tipo de línea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.3.3. Personalizar las propiedades de las líneas . . . . . . . . . . . . . . . . . . . . . . 92
3.3.4. Personalizar los ejes: orden axis . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.3.5. Anotaciones: título, etiquetas y leyendas . . . . . . . . . . . . . . . . . . . . . . 95
3.3.6. Dibujar cosas encima de otras . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.3.7. Añadir texto a la gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.3.8. Gestión de la ventana gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.3.9. Otros tipos de gráficos de datos planos . . . . . . . . . . . . . . . . . . . . . . . 99
3.3.10. Varios ejes en la misma gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.4. Funciones para el dibujo de superficies . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.4.1. Las funciones mesh y meshgrid . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.4.2. La función surf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.4.3. La función contour y variantes . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.4.4. La función colorbar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.4.5. La función view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.4.6. La función colormap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.5. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
SE
UNIVER
VILLA
Índice de tablas e an
Ecuaciones Diferenciales
y Análisis Numérico
3.1. Símbolos para definir el color y el tipo de línea a dibujar en la orden plot . . . . . . . 92
3.2. Algunas propiedades de las líneas, controlables con la orden plot . . . . . . . . . . . . 93
3.3. Algunas propiedades de los títulos y las etiquetas de los ejes . . . . . . . . . . . . . . . 96
9
DAD DE
SI
SE
UNIVER
VILLA
Índice de figuras e an
Ecuaciones Diferenciales
y Análisis Numérico
11
4.2. Condicional doble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.3. Condicional doble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.4. Repetición condicional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4.5. Ruptura de bucle break . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.6. Ruptura de bucle continue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
1 B
SE
Introducción a los
UNIVER
VILLA
ordenadores e an
Ecuaciones Diferenciales
y Análisis Numérico
Las palabras españolas informática y ordenador provienen de las francesas informatique y or-
dinateur. La palabra francesa informatique se construyó como una contracción de las palabras
INFORmation y autoMATIQUE.
La palabra computador (o computadora, mayoritariamente utilizada en Sudamérica), proviene
del término inglés computer. La ciencia que se ocupa de los ordenadores se denomina, en inglés,
Computer Science, aunque también existe la palabra Informatics.
• Software del sistema o sistema operativo, el conjunto de programas necesarios para que el
ordenador tenga capacidad de trabajar (funcionamiento de la pantalla, del teclado, movi-
mientos del ratón, etc.).
• Software de aplicación que son los programas que maneja el usuario (tratamiento de textos,
hojas de cálculo, bases de datos, etc).
El modelo básico de los ordenadores actuales se atribuye a von Neumann1 . Consta principalmente de
la unidad central de proceso, la memoria y los dispositivos de entrada/salida, conectados entre sí como
se muestra de forma esquemática en el diagrama mostrado en la Figura 1.1.
Unidad central de proceso o CPU (Central Process Unit) es el “cerebro” del ordenador: ejecuta las
instrucciones de los programas y controla el funcionamiento de los distintos componentes del ordenador.
Suele estar integrada en un chip denominado microprocesador. Los dos componentes más importantes
de la CPU son:
1
John von Neumann (1903-1957), matemático húngaro-estadounidense, con aportaciones importantes en numerosas
áreas, y uno de los pioneros de los modernos ordenadores: a él se atribuye la “arquitectura de von Neumann”, base de los
ordenadores modernos. Construyó el EDVAC, uno de los primeros ordenadores programables con programa almacenado
en el propio ordenador, junto con los datos, poniendo en práctica la idea teórica debida a otro de los “padres” de la
Informática y de las Ciencias de la Computación, Alan Turing (1912-1954).
13
1. Introducción a los ordenadores 14
DISPOSITIVOS DE
ENTRADA
La unidad de control o CU (Control Unit), dirige y coordina la mayoría de las operaciones del
ordenador. Interpreta cada instrucción enviada por un programa e inicia la acción apropiada para
realizar esa instrucción. Para cada instrucción, la unidad de control repite un grupo de cuatro
operaciones básicas, que constituye lo que se llama ciclo de instrucción:
La unidad aritmético-lógica o ALU (Arithmetic Logic Unit), que realiza las operaciones
elementales que constituyen el programa (sumar, multiplicar, comparar, etc).
La velocidad de un procesador se suele medir en hercios (Hz): número de operaciones por segundo
que puede realizar.
1.1.2 Memoria
Son los componentes de hardware en los que se almacena la información procesada por el ordena-
dor. Generalmente, se distinguen entre la memoria central y la memoria secundaria (véase más
adelante). La principal característica de las memorias es su:
Capacidad: indica la cantidad de datos que puede almacenar. Se mide en bits (bit: binary unit),
o múltiplos de bit. Los bits suelen agruparse de ocho en ocho: 1 byte = 8 bits. Las unidades de
almacenamiento están recogidas en la Tabla 1.1.
bit 0ó1
byte B 8 bits
Kilobyte KB 210 B = 1024 B
Megabyte MB 210 KB = 1024 KB
Gigabyte GB 210 MB = 1024 MB
Terabyte TB 210 GB = 1024 GB
Memoria ROM (Read Only Memory), (memoria de sólo lectura). Es permanente, esto es no
se borra al apagar el ordenador y no se puede alterar (o se puede hacer difícilmente). Almacena
códigos de programa grabados en fábrica necesarios para el funcionamiento del ordenador, como
por ejemplo, la secuencia de instrucciones que hay que ejecutar cuando se enciende el ordenador
(BIOS: Basic Input/Output System; POST: Power On Self Test).
Memoria RAM (Random Acces Memory), (memoria de acceso aleatorio). Almacena las ins-
trucciones de los programas mientras que se ejecutan así como los datos que éstos necesitan
o generen. Su contenido es volátil, es decir, se borra al apagar el ordenador. La denominación
Acceso Aleatorio es antigua y se comenzó a utilizar para diferenciarla de otro tipo de memoria
(cintas magnéticas) que era de acceso secuencial: para acceder a un dato había que recorrer todos
los anteriores.
Con el objetivo de que el procesador pueda obtener los datos de la memoria central más rápidamente,
la mayoría de los procesadores usan un tipo especial de memoria llamada memoria caché ó RAM
caché (pequeña y de acceso muy rápido). En ella, además se almacenan los datos más recientes y las
instrucciones inminentes.
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1100
16 17 18 19 20
Contenido de la
dirección 15
La memoria central es un circuito que se puede imaginar como una enorme tabla que almacena infor-
mación en cada una de sus celdas o posiciones de memoria (véase Figura 1.2).
Cada celda es una agrupación de bytes que se denomina palabra. Dependiendo del ordenador, existen
palabras de 4 bytes (32 bits), 8 bytes (64 bits). Las palabras están numeradas correlativamente, comen-
zando desde 0. El número de una palabra es la dirección de esa palabra. La información almacenada
en una palabra es su contenido.
Dispositivos de entrada/salida son todos los componentes de hardware anexos que permiten la comu-
nicación entre el ordenador y el usuario. Los dispositivos de entrada son, por ejemplo, teclados, ratones,
lectores de discos, micrófonos, etc. Los dispositivos de salida son, por ejemplo, monitores, impresoras,
altavoces, etc.
1.2 Software
El software es el conjunto de instrucciones que hacen que el ordenador funcione y lleve a cabo tareas
concretas, como procesar textos o conectarse a internet. El software también controla la forma en que
se maneja el propio hardware.
Un programa es un conjunto detallado de instrucciones que indican al procesador la forma de llevar
a cabo una tarea concreta.
El software de sistema son los programas indispensables para el funcionamiento del ordenador. Son
el sistema operativo, los compiladores e intérpretes y los utilitarios.
• Los compiladores y traductores son programas capaces de traducir un programa escrito en algún
lenguaje de programación al lenguaje máquina (véase la Sección 1.3), que es el único que entiende
el procesador. Cada lenguaje necesita su propio compilador o traductor.
• Se consideran utilitarios a aquellos programas que realizan tareas complementarias como an-
tivirus, compresores de archivos, editores de texto, recuperación de archivos, visualizadores de
imágenes, etc.
1.3 Lenguajes
Un lenguaje de programación es un conjunto de símbolos y reglas sintácticas y semánticas que
sirven para describir las órdenes que controlan el comportamiento físico y lógico de una máquina.
Los procesadores de las máquinas sólo son capaces de entender y obedecer programas escritos en
lenguaje-máquina, cuyas instrucciones son cadenas binarias (formadas por 0 y 1) que se pueden
“cargar” directamente en la memoria central, sin necesidad de traducción. Sin embargo, el lenguaje-
máquina es específico de cada tipo de ordenador, por lo que los programas escritos en dicho lenguaje
no son portables en general de una máquina a otra. Además, debido a su representación totalmente
numérica, son muy difíciles de escribir, leer y corregir.
El lenguaje ensamblador facilita (sólo un poco) esas tareas, ya que permite la escritura de las ins-
trucciones básicas del lenguaje máquina en una forma más legible para el programador2 .
A modo de ejemplo, una instrucción básica ejecutable por la CPU podría ser:
01101100
| {z } 10100100
| {z } 10110101
| {z }
COPIAR 164 183
CP 164 183
Usualmente hay una correspondencia uno a uno entre las instrucciones simples de un código ensam-
blador y las de un código máquina. Por ello, el ensamblador sigue siendo un lenguaje de bajo nivel :
por un lado el programador necesita conocer en profundidad la arquitectura de la máquina, por otro
los programas escritos en ensamblador no son portables.
Los lenguajes de alto nivel, por el contrario, no obligan al programador a conocer los detalles del
ordenador que utiliza. Las instrucciones se escriben en un formato flexible y más “humano”. Además,
se pueden escribir, de forma sencilla, instrucciones mucho más complicadas: cada instrucción en un
lenguaje de alto nivel corresponde a varias (incluso muchas) de lenguaje-máquina. La instrucción del
los ejemplos anteriores, se podría escribir:
A = B
El ordenador sólo “comprende” los programas escritos en lenguaje-máquina. Cualquier otro debe ser
traducido.
La traducción de un programa escrito en ensamblador a lenguaje-máquina la realiza un programa
específico denominado también ensamblador.
Un programa escrito en un lenguaje de alto nivel se denomina programa fuente. Para traducirlo al
lenguaje-máquina será necesario usar un compilador o un intérprete.
Ejemplo 1.1
5 22 -1
Números reales: números pertenecientes a un subconjunto finito de los números reales (constan
de una parte entera y una parte fraccionaria).
Ejemplo 1.2
0.09 -31.423 3.0
Ejemplo 1.3
true false
Ejemplo 1.4
alfabéticos: A B C ... X Y Z a b c ... x y z
numéricos: 0 1 2 3 4 5 6 7 8 9
especiales: ; . = + - * & $ < > etc.
Números complejos: son datos formados por un par de datos reales y sirven para tratar
números complejos.
Ejemplo 1.5
2+3i -3+i (i es la unidad imaginaria)
Ejemplo 1.6
’Esto es una cadena de caracteres’ ’string’ ’123abc’
Matrices: son conjuntos de datos numéricos organizados para formar una matriz o un vector.
Ejemplo 1.7
[1, 0, 3, 4] (vector fila de dimensión 4)
Bit El bit (binary unit) es la unidad mínima de información que puede almacenarse en un ordenador
o en cualquier dispositivo digital. Un bit es un (diminuto) dispositivo electrónico capaz de tener dos
estados: “apagado”, que se asimila al cero y “encendido”, que se asimila al uno.
Debido a que cada bit sólo dispone de dos estados posibles, los ordenadores utilizan para representar
datos numéricos el sistema binario de numeración, en el que sólo hay dos dígitos: 0 y 1.
Así, con un bit, sólo pueden representarse dos valores.
Byte Para representar más cantidad de valores es necesario usar más bits. Si se usan 2 bits se pueden
representar 22 = 4 valores distintos: 00, 01, 10, 11. Si se usan 4 bits se pueden representar 24 = 16
valores distintos: 0000, 0001, 0010, 0011, . . . etc. En general, con n bits se pueden representar 2n valores
distintos3 .
Los bits se suelen agrupar en grupos de 8, formando un byte. En un byte sólo se pueden representar
28 = 256 valores distintos. Por ejemplo, se pueden representar los números enteros no negativos desde
0 hasta 255.
Palabra Para almacenar datos, a su vez, se suelen considerar agrupaciones de 4 u 8 bytes (32 ó 64
bits), a las que se llama palabra (word). Estas agrupaciones determinan el conjunto de valores que se
pueden almacenar en ellas: en 4 bytes (32 bits) se pueden almacenar 232 =4 294 967 296 valores distintos,
mientras que en 8 bytes (64 bits) se pueden almacenar 264 > 18 × 1018 valores distintos.
3
Variaciones con repetición de 2 elementos (0, 1) tomados de n en n = 2n .
Es claro que el cardinal (número de elementos) del conjunto de datos de un determinado tipo que se
pueden almacenar en un ordenador dependerá del número de bits que se dediquen a ello. Por ejemplo,
en una palabra de 32 bits, no se podrán almacenar más de 232 números enteros distintos. Si se dedican
más bits, el cardinal de dicho conjunto aumentará, pero siempre será finito. Es decir, el conjunto de
números con los que se puede trabajar en un ordenador es finito.
Este hecho tiene una importancia fundamental desde el punto de vista del cálculo numérico, ya que
en un ordenador, contrariamente a lo que sucede cuando trabajamos con el conjunto de los números
reales, el conjunto de números con los que podemos trabajar es discreto: entre dos números-máquina
consecutivos no hay ningún otro.
A los números (enteros o reales) que sí se pueden representar en el ordenador se les suele llamar
números-máquina.
N(10) = ak ak−1 . . . a1 a0
= ak · 10k + ak−1 · 10k−1 + · · · + a1 · 101 + a0 · 100 , ai ∈ {0, 1, 2, . . . , 9} ∀i
Ejemplo 1.8
R(10) = 0.d1 d2 d3 . . . dn . . .
= d1 · 10−1 + d2 · 10−2 + · · · + dn · 10−n + . . .
= 0.d1 + 0.0d2 + · · · + 0.00 . . . dn + . . . , di ∈ {0, 1, . . . , 9} ∀i
Ejemplo 1.9
3 2 5
R = 0.325 = 3 · 10−1 + 2 · 10−2 + 5 · 10−3 = + + = 0.3 + 0.02 + 0.005
10 100 1000
Sistema de numeración binario En este sistema de numeración la posición de cada dígito indica la
potencia de 2 por la que hay que multiplicar el dígito:
Para un número N entero se tiene
N(2) = bk bk−1 . . . b1 b0
= bk · 2k + bk−1 · 2k−1 + · · · + b1 · 21 + b0 · 20 , bi ∈ {0, 1} ∀i
Ejemplo 1.10
Ejemplo 1.11
Paso del sistema decimal al binario Para pasar de un número entero N en el sistema decimal a binario
se dividen sucesivamente dicho número y los cocientes sucesivos por 2, hasta llegar a un cociente que
sea igual a 1. Entonces se tiene N(2) = Cn Rn Rn−1 . . . R2 R1 , con Cn = 1, siendo Rk el resto de la
k-ésima división, como se muestra en la Figura 1.3.
:2 N
C1 R1
C2 R2
C3 R3
... ...
Cn Rn
Ejemplo 1.12
Pasar a binario el número N = 77(10)
:2 77
38 1
19 0
9 1
4 1
2 0
1 0
N = 77(10) = 1001101(2)
Para hacer la operación contraria basta con expresar el significado de la numeración posicional:
N = 1001101(2) = 1 + 22 + 23 + 26 = 1 + 4 + 8 + 64 = 77(10)
El procedimiento (iterativo) para obtener la representación binaria de un número con parte fraccionaria
en sistema decimal es:
Poner R0 = R; calcular 2R0 y tomar c1 = [2R0 ] (parte entera de 2R0 ).
Poner R1 = 2R0 − c1 ; calcular 2R1 y tomar c2 = [2R1 ] (parte entera de 2R1 ).
...
Poner Rk = 2Rk−1 − ck ; calcular 2Rk y tomar ck+1 = [2Rk ] (parte entera de 2Rk ).
...
Este algoritmo se detendrá si:
1. Alguno de los Rk = 0, en cuyo caso el resto de dígitos de la parte fraccionaria binaria serán
también ceros.
2. Algún Rk se repite, es decir, coincide con alguno anterior. En este caso hay un grupo de dígitos
de la parte fraccionaria que se repite de forma periódica.
Si no ocurre ninguno de los casos anteriores, el número en cuestión tiene una parte fraccionaria binaria
infinita.
Ejemplo 1.13
Pasar a binario el número R = 0.23(10)
Ejemplo 1.14
Pasar a binario el número R = 0.2(10)
Ejemplo 1.15
Pasar a binario el número R = 0.3125(10)
R = 0.3125(10) = 0.0101(2)
Obsérvese en los ejemplos 1.13 y 1.14 cómo un número decimal con parte fraccionaria finita en base
10 puede tener una parte fraccionaria infinita en base 2.
Lo contrario no puede suceder: si 0.c1 c2 c3 . . . cn es la expresión en base 2 de un número fraccionario,
se tiene:
c1 c2 cn
0.c1 c2 c3 . . . cn (2) = c1 · 2−1 + c2 · 2−2 + · · · + cn · 2−n = + 2 + ··· + n
2 2 2
donde los ci valen 0 ó 1. Por tanto la expresión anterior no es más que una suma finita de términos
1 1 k
de la forma k = = 0.5k , que tiene siempre un número finito de cifras decimales.
2 2
Ejemplo 1.16
Pasar a decimal el número R = 0.0011101(2)
1 1 1 1
0.0011101(2) = 3
+ 4+ 5+ 7
2 2 2 2
= 0.125 + 0.0625 + 0.03125 + 0.0078125
= 0.2265625(10)
R = 0.0011101(2) = 0.2265625(10)
Ejemplo 1.17
Pasar a binario el número N = 178.023(10)
Pasamos a numeración binaria la parte entera del número Ne = 178, por un lado, y por otro la
parte fraccionaria Nf = 0.023
Operaciones aritméticas en binario Las operaciones en binario se realizan de la misma forma que en
el sistema decimal:
Suma
0+0=0 1 1 1
0+1=1+0=1 1 1 0 1 0 1 = 5310
1 + 1 = 10 = 0 (acarreo 1) + 1 1 1 1 0 = + 3010
1 0 1 0 0 1 1 = 8310
Resta
0−0=0 1 1 0 1 0 1 = 5310
1−0=1 − 1 1 1 1 0 = − 3010
0 − 1 = 1 (acarreo 1) 1 1 1 1
0 1 0 1 1 1 = 2310
1−1=0
Multiplicación
1 1 0 1 0 1 = 5310
× 1 0 1 = × 510
1 1 0 1 0 1
1 1 0 1 0 1
1 0 0 0 0 1 0 0 1 = 26510
División
1 1 0 1 0 1 1 0 1 110101 11
= 1010 +
− 1 0 1 1 0 1 0 101 101
0 0 1 1 0
− 1 0 1 53 3
= 10 +
0 0 1 1 5 5
Notación hexadecimal La manipulación por humanos de números binarios (cadenas largas de ceros
y unos) es tediosa y propensa a la comisión de errores. Por ello, en muchas ocasiones se utilizan otras
notaciones. La más utilizada es la notación hexadecimal, basada en la numeración en base 16.
En ella se necesitan 16 dígitos para representar un número. Puesto que la numeración decimal sólo dis-
pone de 10 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9), se complementa con las primeras letras del abecedario,
hasta conseguir los 16 dígitos: A=10, B=11, C=12, D=13, E=14, F=15. Así, en la numeración hexadecimal
todos los números se representan mediante los dígitos
0 1 2 3 4 5 6 7 8 9 A B C D E F
El valor en base 10 de un número hexadecimal viene dado, como en otros casos, por la interpretación
posicional de los dígitos.
Ejemplo 1.18
Hallar la representación decimal del número H = 2AE1(16)
Para pasar un número decimal a hexadecimal se procede como para otras bases: se divide el número
sucesivamente por 16 hasta llegar a un cociente que sea menor que 16 y y se toman como dígitos el
último cociente seguido de los restos, en orden inverso.
Ejemplo 1.19
Hallar la representación hexadecimal del número D = 41486(10)
La ventaja que ofrece esta notación es que, al ser 16 una potencia de 2 (24 ), la conversión binario ↔
hexadecimal puede hacerse de forma muy sencilla.
Para obtener la representación hexadecimal de un patrón de bits (un número binario) se puede proceder
como sigue:
– Se agrupan los dígitos binarios de 4 en 4 comenzando por la derecha
– Se rellena con ceros a la izquierda si es necesario
– Se sustituye cada grupo de 4 bits por el correspondiente dígito hexadecimal.
En sentido contrario, para obtener la representación binaria de un número hexadecimal basta con
sustituir cada dígito hexadecimal por su codificación binaria.
Ejemplo 1.20
Escribir el patrón de bits P = 1010001000001110 en notación hexadecimal
1010001000001110 = 1010
| {z } 0010
| {z } 0000
| {z } 1110
| {z } = A20E
A 2 0 E
En un ejemplo anterior se ha visto que A20E es la representación en base 16 del número decimal
41486. Se deja al lector la verificación de que
Notación octal La notación octal es similar a la hexadecimal, pero con el sistema de numeración en
base 8 (8 = 23 ). En este sistema se utilizan 8 dígitos para representar cualquier número, 0, 1, 2, 3,
4, 5, 6, 7, y se necesitan 3 bits para representar cada uno de estos 8 dígitos (23 = 8).
La conversión binario ↔ octal se realiza de forma similar al caso hexadecimal, pero agrupando ahora
los dígitos binarios de 3 en 3.
Ejemplo 1.21
Escribir el patrón de bits P = 1010001000001110 en notación octal
En un ejemplo anterior se ha visto que P es la representación binaria del número decimal 41486.
Luego
41486(10) = A20E(16) = 1010001000001110(2) = 12016(8)
Observación Distintos lenguajes de programación tienen diversas formas de escribir los números bi-
narios, octales y hexadecimales. Una de las más habituales consiste en añadir un prefijo formado por
un cero (0) seguido de una letra que indica el sistema: 0b para indicar un número binario, 0o para
indicar un número octal y 0x para indicar un número hexadecimal:
0b100110 indica el número binario 100110
0o7103 indica el número octal 7103
0x21A1 indica el número hexadecimal 21A1
8 bits 0 . . . 255
16 bits 0 . . . 65535
32 bits 0 . . . 4294967295
Tabla 1.2: Rango de valores enteros sin signo representables,
en función del número de bits.
Para almacenar estos números, se pasan a binario y se completan con ceros a la izquierda hasta rellenar
todos los bits disponibles.
Los enteros sin signo se utilizan, por ejemplo, para contar y para designar direcciones de memoria.
Ejemplo 1.22
¿Cómo se almacena el entero sin signo N = 109(10) en un registro de 16 bits?
0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1
Signo-magnitud En este sistema –el más básico– se reserva un bit para indicar el signo del número.
Normalmente es el primero y se asocia el 0 con el signo + y el 1 con el signo − . En los bits restantes
se almacena la representación binaria de la magnitud del número.
Al hacer esto se dispone de un bit menos para almacenar la magnitud. Así, mientras que con 8 bits
se pueden almacenar, como entero sin signo, desde el 0 hasta el 255, con este sistema sólo se pueden
almacenar desde el −127 hasta el +127. En general, con N bits, se podrán representar desde el −(2N −1 −
1) hasta el +(2N −1 − 1).
Ejemplo 1.23
¿Cómo se almacena el entero N = −109(10) en un registro de 16 bits en el sistema signo-
magnitud?
1 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1
Para decodificar un patrón de bits en signo-magnitud con N bits, sólo habría que pasar a numeración
decimal el número binaro contenido en los últimos N − 1 bits, y luego asignarle signo + o − según que
el primer bit sea cero o uno.
Ejemplo 1.24
Decodificar el patrón de bits siguiente, en el sistema signo-magnitud con 16 bits:
1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1
Hay que observar que, en este sistema, hay dos representaciones distintas para el cero, lo cual no es
deseable (véase la tabla 1.4).
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = +0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = −0
Tabla 1.4: Las dos representaciones del cero en el sistema
signo-magnitud con 16 bits.
Este sistema no se utiliza en los ordenadores actuales para almacenar números enteros con signo, ya
que las operaciones de sumar y restar no son fáciles con ella y por la doble representación del cero.
Sin embargo, resulta útil para aplicaciones en las que no haya que realizar operaciones aritméticas con
los números, como por ejemplo, en la transmisión de señales.
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = +0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 = −0
Tabla 1.5: Las dos representaciones del cero en el sistema de
complemento a 1 con 16 bits.
Ejemplo 1.25
¿Cómo se almacena el entero N = 109(10) en un registro de 16 bits en el sistema complemento a
1? ¿Y el número N = −109(10) ?
0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1
El número −109(10) = −1101101(2) , negativo, se almacena con un 1 en el bit del signo y comple-
mentado a 1 en el resto, es decir:
1 1 1 1 1 1 1 1 1 0 0 1 0 0 1 0
Para decodificar un patrón de bits en complemento a 1 con N bits, se procede como sigue:
— Si el primer bit es un cero, se procede como en el sistema signo-magnitud: se pasa a numeración
decimal el número binario contenido en los últimos N − 1 bits. Se le pone signo + .
— Si el primer bit es un uno, se aplica el complemento a 1 al número binario contenido en los N − 1
últimos bits y el número resultante se pasa a numeración decimal. Se le pone signo − .
Ejemplo 1.26
Decodificar el siguiente patron de bits, en el sistema de complemento a 1, con un registro de 16
bits
1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1
Puesto que el primer bit es un uno, se aplica el complemento a 1 a la parte de la magnitud (15
bits):
C115 110001 = 111111111001110
1 0 1 1 1 0 0 1 0 0 0 Número binario
0 1 0 0 0 1 1 0 1 1 1 Complemento a 1
+ 1 Sumar 1
0 1 0 0 0 1 1 1 0 0 0 Complemento a 2
Otra forma de obtenerlo es la siguiente: comenzando por la derecha, se dejan los bits como están hasta
llegar al primer 1 (incluido). A partir de este, se aplica el complemento a 1 a todo el resto:
1 0 1 1 1 0 0 1 0 0 0 Número binario
↓ ↓ ↓ ↓ ↓ ↓ ↓ k k k k
0 1 0 0 0 1 1 1 0 0 0 Complemento a 2
Ejemplo 1.27
¿Cómo se almacena el entero N = 109(10) en un registro de 16 bits en el sistema complemento a
2? ¿Y el número N = −109(10) ?
0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = +0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = −0
Sin embargo, por convenio, se asigna la representación del −0 al número −2N −1 , es decir, al −128
con registros de 8 bits, al −32768 con registros de 16 bits, etc. Así el rango de números que se pueden
representar en el sistema de complemento a 2 queda aumentado en un número negativo más: es el
mostrado en la tabla 1.6.
Así, para decodificar un patrón de bits en complemento a 2 con 16 bits se procede como sigue:
— Si el primer bit es un cero, se procede como en el sistema signo-magnitud: se pasa a numeración
decimal el número binario contenido en los últimos 15 bits. Se le pone signo +.
— Si el primer bit es un uno y todos los demás son ceros, el número representado es −215 = −32768.
— Si el primer bit es un uno, se aplica el complemento a 2 al número binario contenido en los últimos
15 bits, y el número binario resultante se pasa a numeración decimal. Se le pone signo − .
Ejemplo 1.28
Decodificar el siguiente patron de bits, en el sistema de complemento a 2, con un registro de 16
bits
1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1
Puesto que el primer bit es un uno, se aplica el complemento a 2 a la parte de la magnitud (15
bits):
C215 000000000110001 = 111111111001111
última suma (la correspondiente al bit del signo) produzca acarreo de 1, se descarta el 1 acarreado. El
patrón de bits resultante es la representación en complemento a 2 con n bits de la suma.
Ejemplos 1.29
Realizar las siguientes sumas en la representación en complemento a 2 con 8 bits:
73 + 26 73 − 26 − 73 + 26 − 73 − 26
Escribimos la representación en complemento a 2 con 8 bits (C28 ) de todos los números que aparecen:
+ 73(10) = 0 1 0 0 1 0 0 1
+ 26(10) = 0 0 0 1 1 0 1 0
+ 99(10) 0 1 1 0 0 0 1 1 = 01100011C28 = +1100011(2) = +99(10)
+ 73(10) = 0 1 0 0 1 0 0 1
− 26(10) = 1 1 1 0 0 1 1 0
+ 47(10) 1 0 0 1 0 1 1 1 1 = 00101111C28 = +101111(2) = +47(10)
− 73(10) = 1 0 1 1 0 1 1 1
+ 26(10) = 0 0 0 1 1 0 1 0
− 47(10) 1 1 0 1 0 0 0 1 = 11010001C28 = -101111(2) = −47(10)
− 73(10) = 1 0 1 1 0 1 1 1
− 26(10) = 1 1 1 0 0 1 1 0
− 99(10) 1 1 0 0 1 1 1 0 1 = 10011101C28 = -1100011(2) = −99(10)
Representación sesgada Este sistema consiste en, una vez fijado un número entero positivo K, deno-
minado sesgo, sumar este número a todos los demas. Así, en lugar de almacenar el entero con signo
N , se almacena N + K como un entero sin signo. Este sistema también se denomina exceso a K .
Normalmente, se elige como sesgo K = 2n−1 o bien K = 2n−1 − 1, donde n es el número de bits que
se emplean para el almacenamiento.
• Caso K = 2n−1 .
Puesto que N +K se va a almacenar como un entero sin signo, sólo se pueden almacenar números N tales
que N + K ≥ 0. Así, el número más pequeño que se puede representar por este sistema es Nmin tal que
Nmin +2n−1 = 0, es decir Nmin = −2n−1 . El número más grande es Nmax tal que Nmax +2n−1 = 2n −1,
es decir Nmax = 2n−1 − 1. Ahora el cero tiene una única representación.
Por ejemplo, si se utilizan n = 8 bits, será K = 27 = 128, y se habla de exceso a 128. El número más
pequeño que se puede representar en exceso a 128 es Nmin = −128 y el más grande es Nmax = 127.
• Caso K = 2n−1 − 1.
Ahora el número más pequeño que se puede almacenar es Nmin tal que
Nmin +2n−1 −1 = 0, es decir Nmin = −2n−1 +1, y el más grande Nmax tal que Nmax +2n−1 −1 = 2n −1,
es decir Nmax = 2n−1
Por ejemplo, si se utilizan n = 8 bits, será K = 27 − 1 = 127, y se habla de exceso a 127. El número
más pequeño que se puede representar en exceso a 127 es Nmin = −127 y el más grande es Nmax = 128.
Ejemplo 1.30
¿Cómo se almacena el número −95(10) en un registro de 8 bits con exceso a 127? ¿Y el número
95(10) ?.
Para almacenar el número −95(10) se comienza por sumarle 127, −95 + 127 = 32(10) y se pasa
32(10) a binario, 33(10) = 100000(2) . El resultado se rellena con ceros a la izquierda hasta completar
los 8 bits. El valor binario resultante siempre comenzará por cero.
0 0 1 0 0 0 0 0
Para el número 95(10) se tiene: 95 + 127 = 222(10) = 11011110(2) . El valor binario resultante
comienza por uno, ya que el valor a almacenar es mayor que 128.
1 1 0 1 1 1 1 0
Para decodificar un patrón de 8 bits en exceso a 127, simplemente hay que decodificarlo como si fuera
un entero sin signo y luego restarle el sesgo.
El sistema exceso a 2n−1 − 1 se utiliza como parte de la codificación de números reales, como se verá
a continuación.
Decodificación de un patrón de bits La pregunta que quizá cabe hacerse es: ¿cómo sabe la memoria
del ordenador en qué sistema está almacenado un determinado patrón de bits en un registro de la
memoria? Es decir, ¿cómo sabe cómo hay que decodificarlo? La respuesta es que NO LO SABE.
Una vez almacenado un determinado número (o dato) en la memoria, con cualquier sistema, no queda
ningún rastro que indique cuál fue ese sistema. Es el programa que “lea” o el dispositivo que utilice ese
dato el que debe gestionar que se decodifique/interprete de la manera adecuada.
A modo de ejemplo, se presenta en la tabla 1.7 cómo sería la representación de los números en los
distintos sistemas explicados, en un registro imaginario de 4 bits. Obsérvese que, en todos los casos
(salvo el sesgado), el primer bit de la izquierda es un 1 para los números negativos.
Ejemplo 1.31
El patrón de bits
0 0 0 1 0 0 1 1 1 0 1 0 1 1 0 1
comienza por un 0, luego representa el mismo número entero en cualquiera de los sistemas: es el
entero positivo 5037.
1 20 1
0 21
1 22 4
1 23 8
0 24
1 25 32
0 26
1 27 128
1 28 256
1 29 512
0 210
0 211
1 212 4096
0 213
0 214
0 215
5037
Ejemplo 1.32
El patrón de bits siguiente comienza por un 1, luego puede representar distintos números.
1 0 0 1 0 0 1 1 1 0 1 0 1 1 0 1
Los números reales se representan según la norma denominada IEEE 7545 , establecida en 1985 para
estandarizar la forma de almacenarlos en todas las marcas de ordenadores.
Esta norma establece dos formatos básicos para representar los números reales: simple precisión
y doble precisión. En ambos casos lo que se representa es una expresión normalizada del número
denominada coma flotante.
Como flotante La coma flotante (en inglés, floating point) es un modo “normalizado” de representa-
ción de números con parte fraccionaria que se adapta al orden de magnitud de los mismos, permitiendo
así ampliar el rango de números representables.
Consiste en trasladar la coma (o punto) decimal hasta situarla detrás de la primera cifra significativa
(la primera, comenzando por la izquierda, que no es nula) a la vez que se multiplica el número por la
potencia adecuada de la base de numeración.
Ejemplo 1.33
5
IEEE son las siglas del Institute of Electrical and Electronics Engineers.
Sea, pues, un número R dado por la siguiente representación binaria en coma flotante:
R = (±)m1 .m2 m3 . . . mk · 2±e ,
donde m1 m2 m3 . . . mk se llama parte significativa o mantisa y el número ±e (entero, binario) es el
exponente, es decir, el número de lugares que hay que trasladar la coma a derecha o izquierda para
obtener el número inicial. Obsérvese que m1 tiene que ser siempre igual a 1, ya que es, por definición,
la primera cifra no nula empezando por la izquierda.
Ejemplo 1.34
Simple precisión Para almacenar R en simple precisión se utiliza una palabra de 32 bits, que se
distribuyen de la siguiente manera:
Tabla 1.8: Algunos casos especiales en la norma IEEE 754. Cuando el exponente es todo de ceros, no
1 1 1 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
existe el bit implícito.
Ejemplo 1.35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
Representar
1 1en0simple
0 0 0precisión
0 0 1 0el 1número
1 0 1R0= 0−5.625
0 0 (10) 0 0 0 0· 2010 0 0 0 0 0 0
0 0 =0−1.01101
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
1 1 0 0 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Se tiene:
El mayor exponente positivo que se puede almacenar en simple precisión es, como hemos visto,
128.
El menor exponente negativo: −127.
El número de mayor magnitud: ≈ 2128 ≈ 3 × 1038 .
El número de menor magnitud: ≈ 2−127 ≈ 5 × 10−39 .
R. Echevarría Dpto. EDAN - Universidad de Sevilla
1. Introducción a los ordenadores 38
Doble precisión Para la representación en doble precisión se utilizan 64 bits distribuidos de la siguiente
forma:
El mayor exponente positivo que se puede almacenar en doble precisión es, como hemos visto,
210 = 1024.
El menor exponente negativo: −1023.
El número de mayor magnitud: ≈ 10308 .
El número de menor magnitud: ≈ 10−307 .
Overflow y underflow Es claro que usando palabras con un número finito de bits (32 o 64) sólo
podemos representar una cantidad finita de números, en consecuencia podemos encontrar los siguientes
fenómenos:
Overflow fenómeno que se produce cuando el resultado de un cálculo es un número real con exponente
demasiado grande (en el caso de la simple precisión, > 128). Cuando se produce un overflow
durante la ejecución de un programa, los ordenadores modernos generan un “evento” de error y
como resultado devuelven el símbolo Inf.
Underflow fenómeno que se produce cuando el resultado de un cálculo es un número real con ex-
ponente demasiado pequeño (en los ejemplos arriba mencionados < −127). Cuando se produce
durante la ejecución de un programa, normalmente se sustituye el número por cero.
1.5.6 Errores
Al realizar cálculos en el ordenador, es inevitable cometer errores. Grosso modo, estos errores pueden
ser de tres tipos:
1. Errores en los datos de entrada: vienen causados por los errores al realizar mediciones de
magnitudes físicas.
3. Errores de redondeo: aparecen debido a la cantidad finita de números que podemos representar
en un ordenador. Por ejemplo, consideramos las dos siguientes mantisas-máquina consecutivas
que podemos almacenar en 23 bits (recordar que el primer dígito es siempre 1 y por tanto no se
guarda):
R1 = |1.00000000000000000000000
{z }
24 dígitos
R2 = |1.00000000000000000000001
{z }
24 dígitos
Ningún número real entre R1 y R2 es representable en un ordenador que utilice 23 bits para la
mantisa. Tendrá que ser aproximado bien por R1 , bien por R2 , cometiendo con ello un error E
que verifica:
De esta forma, en el almacenamiento de la mantisa se comete un error menor o igual que 10−7 o,
dicho de otra forma, sólo los 7 primeros dígitos significativos de un número decimal almacenado
en el ordenador en simple precisión se pueden considerar exactos.
Ejemplo 1.36
R = 0.1(10) = 0.00011
¯ (2)
= 0.00011001100110011 . . .
=4
z}|{
= = 1.10011001100110011001100 − 100
| {z } 1100 · · · × 2
24
R no es un número-máquina (su mantisa binaria tiene infinitos dígitos).
La aproximación de R por redondeo en simple precisión es:
−100
R = |1.10011001100110011001101
{z } ×2 .
24 dígitos
Aritmética en coma flotante Para llevar a cabo una operación aritmética en un ordenador, hay que
tener en cuenta que los números con los que trabajamos pueden no ser exactamente los de partida, sino
sus aproximaciones por redondeo. Incluso aunque los números sean números-máquina (representables
en palabras de 32 bits), los resultados que obtengamos puede que no lo sean. A continuación vamos a
ver algunos ejemplos de cálculos efectuados en la aritmética de coma flotante. Por facilidad, lo hacemos
en el caso (hipótetico) de un ordenador que trabaje con números en expresión decimal con 7 dígitos
para la mantisa, y con aproximación por redondeo.
Ejemplo 1.37
Para sumar dos números en coma flotante, primero se representan con el mismo exponente, luego
se suman y el resultado se redondea:
0.123456700000 ·103
+ 0.000003232323 ·103
0.123459932323 ·103 (sólo se almacenan los 7 primeros dígitos)
0.1234599 ·103 (resultado-máquina de x1 + x2 )
Ejemplo 1.38
Para multiplicar (la división se hace por algoritmos más complicados) dos números en coma
flotante se multiplican las mantisas y se suman los exponentes. El resultado se normaliza:
x1 = 0.4734612 · 103
x2 = 0.5417242 · 105
0.4734612 ·103
× 0.5417242 ·105
0.25648538980104 ·108
0.2564854 ·108 (resultado-máquina de x1 · x2 )
n = 00000000000 = 0
Ahora bien, en realidad, la norma IEEE 754 contiene una ampliación del caso en que se tiene un
exponente todo de ceros:
y el patrón
11 51
z }| { z }| {
0 00 . . . 0 00 . . . 0 1
correspondería al número
51
2 B
SE
Introducción a MATLAB.
UNIVER
VILLA
Operaciones elementales e an
Ecuaciones Diferenciales
y Análisis Numérico
matlab es un potente paquete de software para computación científica, orientado al cálculo numérico,
a las operaciones matriciales y especialmente a las aplicaciones científicas y de ingeniería. Ofrece lo que
se llama un entorno de desarrollo integrado (IDE), es decir, una herramienta que permite, en una sóla
aplicación, ejecutar órdenes sencillas, escribir programas utilizando un editor integrado, compilarlos (o
interpretarlos), depurarlos (buscar errores) y realizar gráficas.
Puede ser utilizado como simple calculadora matricial, pero su interés principal radica en los cientos
de funciones tanto de propósito general como especializadas que posee, así como en sus posibilidades
para la visualización gráfica.
matlab posee un lenguaje de programación propio, muy próximo a los habituales en cálculo numérico
(Fortran, C, . . . ), aunque mucho más tolerante en su sintaxis, que permite al usuario escribir sus
propios scripts (conjunto de comandos escritos en un fichero, que se pueden ejecutar con una única
orden) para resolver un problema concreto y también escribir nuevas funciones con, por ejemplo,
sus propios algoritmos, o para modularizar la resolución de un problema complejo. matlab dispone,
además, de numerosas Toolboxes, que le añaden funcionalidades especializadas.
Numerosas contribuciones de sus miles de usuarios en todo el mundo pueden encontrarse en la web de
The MathWorks: [Link]
2.1 Introducción
Al iniciar matlab nos aparecerá una ventana más o menos como la de la Figura 2.1 (dependiendo del
sistema operativo y de la versión)
Si la ubicación de las ventanas integradas es diferente, se puede volver a ésta mediante:
Menú Desktop → Desktop Layout → Default
Se puede experimentar con otras disposiciones. Si hay alguna que nos gusta, se puede salvaguardar con
Menú Desktop → Desktop Layout → Save Layout . . .
dándole un nombre, para usarla en otras ocasiones. De momento, sólo vamos a utilizar la ventana
principal de matlab: Command Window (ventana de comandos). A través de esta ventana nos comu-
nicaremos con matlab, escribiendo las órdenes en la línea de comandos. El símbolo >> al comienzo
de una línea de la ventana de comandos se denomina prompt e indica que matlab está desocupado,
disponible para ejecutar nuestras órdenes.
Las explicaciones sobre las funciones/comandos que se presentan en estas notas están muy resumidas
y sólo incluyen las funcionalidades que, según el parecer subjetivo de la autora, pueden despertar más
interés. La mayoría de las funciones tienen mas y/o distintas funcionalidades que las que se exponen
43
2. Introducción a MATLAB. Operaciones elementales 44
aquí. Para una descripción exacta y exhaustiva es preciso consultar la Ayuda on-line.
Los tipos básicos de datos que maneja matlab son números reales, booleanos (valores lógicos) y
cadenas de caracteres (string). También puede manipular distintos tipos de números enteros, aunque
sólo suele ser necesario en circunstancias específicas.
En matlab, por defecto, los números son codificados como números reales en coma flotante en doble
precisión. La precisión, esto es, el número de bits dedicados a representar la mantisa y el exponente,
depende de cada (tipo de) máquina.
matlab manipula también otros objetos, compuestos a partir de los anteriores: números complejos,
matrices, cells, estructuras definidas por el usuario, clases Java, etc.
El objeto básico de trabajo de matlab es una matriz bidimensional cuyos elementos son números
reales o complejos. Escalares y vectores son considerados casos particulares de matrices. También se
pueden manipular matrices de cadenas de caracteres, booleanas y enteras.
El lenguaje de matlab es interpretado, esto es, las instrucciones se traducen a lenguaje máquina
una a una y se ejecutan antes de pasar a la siguiente. Es posible escribir varias instrucciones en la
misma línea, separándolas por una coma o por punto y coma. Las intrucciones que terminan por punto
y coma no producen salida de resultados por pantalla.
Algunas constantes numéricas están predefinidas (ver sección 2.5.2)
matlab distingue entre mayúsculas y minúsculas: pi no es los mismo que Pi.
matlab conserva un historial de las instrucciones escritas en la línea de comandos. Se pueden recuperar
instrucciones anteriores, usando las teclas de flechas arriba y abajo. Con las flechas izquierda y derecha
nos podemos desplazar sobre la línea de comando y modificarlo.
Se pueden salvaguardar todas las instrucciones y la salida de resultados de una sesión de trabajo de
matlab a un fichero:
>> helpwin
Además, a través del navegador del Help se pueden descargar, desde The MathWorks, guías detalladas,
en formato pdf, de cada capítulo.
Las funciones deben guardarse en un fichero con el mismo nombre que la función y sufijo .m. Lo que
se escribe en cualquier línea detrás de % es considerado como comentario.
Ejemplo 2.1
El siguiente código debe guardarse en un fichero de nombre areaequi.m.
La primera línea de una M-función siempre debe comenzar con la claúsula (palabra reservada)
function. El fichero que contiene la función debe estar en un sitio en el que MATLAB lo pueda
encontrar, normalmente, en la carpeta de trabajo.
Se puede ejecutar la M-función en la misma forma que cualquier otra función de MATLAB:
Ejemplos 2.2
>> areaequi(1.5)
ans =
0.9743
>> rho = 4 * areaequi(2) +1;
>> sqrt(areaequi(rho))
ans =
5.2171
Los breves comentarios que se incluyen a continuación de la línea que contiene la claúsula function
deben explicar, brevemente, el funcionamiento y uso de la función. Además, constituyen la ayuda
on-line de la función:
Ejemplo 2.3
>> help areaequi
areaequi(long) devuelve el area del triangulo
equilatero de lado = long
Funciones anónimas Algunas funciones sencillas, que devuelvan el resultado de una expresión, se
pueden definir mediante una sóla instrucción, en mitad de un programa (script o función) o en la línea
de comandos. Se llaman funciones anónimas. La sintaxis para definirlas es:
1
También es posible definir funciones anidadas, esto es, funciones «insertadas» dentro del código de otras funciones.
(Se informa aquí para conocer su existencia. Su utilización es delicada.)
Las funciones anónimas pueden tener varios argumentos y hacer uso de variables previamente definidas:
Workspace (espacio de trabajo) es el conjunto de variables que en un momento dado están definidas en
la memoria del MATLAB.
Las variables creadas desde la línea de comandos de MATLAB pertenecen al Base Workspace (espacio
de trabajo base; es el que se puede «hojear» en la ventana Workspace). Los mismo sucede con las
variables creadas por un script que se ejecuta desde la línea de comandos. Estas variables permanecen
en el Base Workspace cuando se termina la ejecución del script y se mantienen allí durante toda la
sesión de trabajo o hasta que se borren.
Sin embargo, las variables creadas por una M-función pertenecen al espacio de trabajo de dicha función,
que es independiente del espacio de trabajo base. Es decir, las variables de las M-funciones son locales:
MATLAB reserva una zona de memoria cuando comienza a ejecutar una M-función, almacena en esa
zona las variables que se crean dentro de ella, y «borra» dicha zona cuando termina la ejecución de la
función.
Esta es una de las principales diferencias entre un script y una M-función: cuando finaliza la ejecución
de un script se puede «ver» y utilizar el valor de todas las variables que ha creado el script en el
Workspace; en cambio, cuando finaliza una función no hay rastro de sus variables en el Workspace.
Para hacer que una variable local pertenezca al espacio de trabajo de otra función, hay que declararla
global en cada uno de ellos: la orden
en una función hace que las variables a, suma y error sean almacenadas en un sitio especial y puedan
ser utilizadas por otras funciones que también las declaren como globales.
2.2 Números
Enteros Los números enteros se escriben sin punto decimal
Ejemplo 2.6
23 321 -34
Reales Los números reales se pueden escribir en notación fija decimal y en notación científica, utili-
zando el punto decimal (no la coma):
Ejemplo 2.7
23. -10.1 22.0765
-2
2.e-2 = 2 × 10 = 0.02
2.07e+5 = 2.07 × 105 = 207000
Ejemplo 2.8
2 + 3 i -4.07 - 2.3 i
Descripción Símbolo
Exponenciación ∧
Suma +
Resta −
Multiplicación ∗
División /
Las operaciones aritméticas NO se efectúan siempre en el orden en que están escritas. El orden viene
determinado por las reglas siguientes:
1. Exponenciaciones.
2. Multiplicaciones y divisiones.
3. Sumas y restas.
Ejemplos 2.9
2 + 3 * 4 = 2 + 12 = 14
(2 + 3) * 4 = 5 * 4 = 20
1 / 3 * 2 = 0.3333. . . * 2 = 0.6666
1 / (3 * 2) = 1 / 6 = 0.166666...
2 + 3 ∧ 4 / 2 = 2 + 81 / 2 = 2 + 40.5 = 42.5
4 ∧ 3 ∧ 2 = (4 ∧ 3) ∧ 2 = 64 ∧ 2 = 4096
b a+b
a+b/2 = a + , (a+b)/2 =
2 2
1 1
1/2*x = x, 1/(2*x) =
2 2 x
3
Ejercicio 2.10 Calcular 2 /2 mediante una única expresión.
1
Resultado : 1.0905
Resultado : 0.3882
Ejemplos 2.12
>> 0.002 será mostrado como 0.0020
>> 0.0000123 será mostrado como 1.2300e-05
>> -709.2121211 será mostrado como -709.2121
>> 1003.010101 será mostrado como 1.0030e+03
2.5 Variables
Cuando hay que hacer cálculos largos interesa guardar resultados intermedios para utilizarlos más
adelante. ¿Dónde se guardan estos resultados? En variables.
Una variable es un nombre simbólico que identifica una parte de la memoria, en la que se pueden
guardar números u otro tipo de datos. El contenido de una variable se puede recuperar y modificar
cuantas veces se desee, a lo largo de una sesión de trabajo o durante la ejecución de un programa.
En matlab, los nombres de las variables deben estar formados por letras y números, hasta un máximo
de 19, y comenzar por una letra. Se debe tener en cuenta que se distingue entre letras mayúsculas y
minúsculas.
En la mayoría de los lenguajes de programación es necesario especificar el tipo de dato que va a contener
una variable antes de usarla, declarándolo con las ordenes específicas. En el lenguaje de programación
de matlab, las variables no necesitan ningún tipo de declaración y pueden almacenar sucesivamente
distintos tipos de datos: enteros, reales, escalares, matriciales, caracteres, etc. Se crean, simplemente,
asignándoles un valor.
variable = expresión
que debe ser interpretada como: evaluar el resultado de la expresión y almacenarlo en la dirección de
memoria correspondiente a variable.
El signo = significa: almacenar el resultado de la derecha en la variable que está a la izquierda. Por
ello, el orden en esta instrucción un fundamental: NO se puede escribir expresión = variable.
La acción de almacenar un valor en una variable hace que se pierda el valor que, eventualmente, tuviera
dicha variable previamente.
Ejemplos 2.14
a = 2 :: guardar en la variable a el valor 2
b = -4 :: guardar en la variable b el valor -4
raiz = sqrt(2*b+8*a) :: guardar en la variable raiz el resultado
de la expresión de la derecha
a = a + 1 :: sumar 1 al contenido de a (guardar 3 en a)
Ejemplos 2.15
>> a = 1/3 :: sin punto y coma
ans =
0.3333
Ejemplos 2.16
>> pi % para ver el valor de pi con 15 cifras,
ans = % usar previamente la orden format long
3.1416
Ejemplos 2.17
>> sqrt(34*exp(2))/(cos(23.7)+12)
ans =
1.3058717
>> 7*exp(5/4)+3.54
ans =
27.97240
>> exp(1+3i)
ans =
- 2.6910786 + 0.3836040i
Descripción Símbolo
Igual a ==
No igual a ∼=
Menor que <
Mayor que >
Menor o igual que <=
Mayor o igual que >=
Ejemplos 2.18
>> 3<6
ans =
1 (logical) (true)
>> 0>=1
ans =
0 (logical) (false)
>> 'A'=='B'
ans =
0 (logical) (false)
Aunque matlab muestre los valores lógicos como 0 y 1, no se deben confundir con los correspondientes
valores numéricos. Así por ejemplo, el intento de calcular sin(3<4) dará error ya que la función sin
no admite valores lógicos como argumento de entrada.
Sin embargo, cuando aparecen valores lógicos en expresiones aritméticas junto con datos numéricos,
matlab transforma el valor lógico true en el valor numérico 1 y el valor lógico false en el valor
numérico 0 antes de realizar la operación.
Ejemplos 2.19
>> (3<6) + 2
ans =
3
>> (0>=1) * pi
ans =
0
Descripción Símbolo
Negación ∼
Conjunción &
Disyunción |
El primero de ellos, el operador de negación lógica ∼ , es un operador unario, es decir, actúa sobre un
solo operando lógico, dando como resultado el opuesto de su valor, como muestra la siguiente tabla:
A ∼A
true false
false true
Ejemplos 2.20
>> ∼ (6 > 10)
ans =
1 (logical)
>> ∼ ('b' > 'a')
ans =
0 (logical)
Los otros dos operadores lógicos actúan siempre entre dos operandos. Los resultados de su aplicación
se muestran en la tabla 2.6.
A B A& B A| B
true true true true
true false false true
false true false true
false false false false
Ejemplos 2.21
>> (1 > 5) & (5 < 10)
ans =
0 (logical)
>> (0 < 5) | (0 > 5)
ans =
1 (logical)
Cuando aparecen valores numéricos en una operación lógica, matlab transforma el valor numérico 0
en el valor lógico false y cualquier otro (distinto de 0) en true.
Ejemplos 2.22
>> 3 & true
ans =
1 (logical)
>> 0 & true
ans =
0 (logical)
En una expresión pueden aparecer varios operadores lógicos. En ese caso, el orden en que se evalúan,
es el siguiente:
1. La negación lógica ∼
2. La conjunción y disyunción & y | (de izquierda a derecha en caso de que haya varios).
Ejemplo 2.23
>> ∼ (5 > 0) & (5 < 4)
ans =
0 (logical)
Si en una expresión hay paréntesis, lo primero que se evalúa es su contenido. Si hay paréntesis
anidados, se comienza por los más internos. Si hay varios grupos de paréntesis disjuntos, se
comienza por el que esté más a la izquierda.
Ejemplos 2.24
>> ∼ 5 > 0 & 4 < 4 + 1
ans =
0 (logical)
>> -5 < -2 & -2 < 1
ans =
1 (logical)
>> -5 < -2 < 1
ans =
0 (logical)
2.10 Matrices
Las matrices bidimensionales de números reales o complejos son los objetos básicos con los que trabaja
matlab. Los vectores y escalares son casos particulares de matrices.
La forma más elemental de introducir matrices en matlab es describir sus elementos de forma exhaus-
tiva (por filas y entre corchetes rectos [ ]) : elementos de una fila se separan unos de otros por comas
y una fila de la siguiente por punto y coma.
Observación.
El hecho de que, al introducirlas, se escriban las matrices por filas no significa que internamente,
en la memoria del ordenador, estén así organizadas: en la memoria las matrices se almacenan
como un vector unidimensional ordenadas por columnas.
1 5 9 2 6 10 3 7 11 4 8 12
Los siguientes operadores y funciones sirven para generar vectores con elementos regularmente espa-
ciados (útiles en muchas ocasiones) y para trasponer matrices.
Se pueden también utilizar los vectores/matrices como objetos (bloques) para construir otras matrices.
>> v1 = 1:4;
Las siguientes funciones generan algunas matrices especiales que serán de utilidad.
MATLAB posee, además, decenas de funciones útiles para generar distintos tipos de matrices. Para
ver una lista exhaustiva consultar:
Help → MATLAB → Functions: By Category → Language Fundamentals → Arrays and Matrices
A =
1 2 3
6 3 1
-1 0 0
2 0 4
>> A(2, 3)
ans =
1
>> A(4, 2) = 11
A =
1 2 3
6 3 1
-1 0 0
2 11 4
Se puede acceder a una fila o columna completa sustituyendo el índice correspondiente por : (dos
puntos)
También se puede acceder a una “submatriz” especificando los valores de las filas y columnas que la
definen:
Otra forma de acceder a elementos específicos de una matriz es utilizando vectores de valores lógicos
en lugar de subíndices, como se muestra en el siguiente ejemplo (más adelante se verán ejemplos más
interesantes).
Además de estos operadores, matlab dispone de ciertos operadores aritméticos que operan elemento
a elemento. Son los operadores .* ./ y .ˆ, muy útiles para aprovechar las funcionalidades
vectoriales de matlab. Véase la tabla 2.10.
Por otra parte, la mayoría de las funciones matlab están hechas de forma que admiten matrices como
argumentos. Esto se aplica en particular a las funciones matemáticas elementales y su utilización debe
entenderse en el sentido de elemento a elemento. Por ejemplo, si A es una matriz de elementos aij ,
exp(A) es otra matriz con las mismas dimensiones que A, cuyos elementos son eaij .
2.11 Ejercicios
1. Comprobar el valor de las siguientes expresiones:
π
π tan2 √ 3
5
M= 4
2 = 1.3985 N= = 25.9782
e + sen2 (17◦ ) exp(−0.66)
1 2 3
Ç å
2. Definir la matriz A= 4 5 6 y crear, a partir de A, las siguientes matrices:
7 8 9
2 3 1 3
Ç å Ç å
4 5
B= 7 8 , C = (7, 8, 9) , D= 5 6 , E= 4 6
8 9 7 9
1 2 3 0 0 1 2 3 4 0
Ç å Ç å
F = (1, 4, 7, 2, 5, 8, 3, 6, 9), G= 4 5 6 0 1 , H= 4 5 6 4 1
7 8 9 0 0 7 8 9 4 0
1 0 0 0 0 0 1 1 1 1 2 3 4 5 6
Ç å Ç å Ç å
A= 0 1 0 0 0 , B= 1 0 1 1 1 , C= 3 4 5 6 7
0 0 1 0 0 1 1 0 1 1 4 5 6 7 8
4. Generar los vectores que se indican utilizando los operadores adecuados de MATLAB (esto es,
sin escribir expresamente todas sus componentes):
5. Generar los vectores que se indican, de longitud n variable. Esto es, escribir las órdenes adecuadas
para construir los vectores que se indican de forma que puedan servir para cualquier valor de n.
a) w1 = (1, 3, 5, 7, . . . ) e) w5 = (0, 2, 2, . . . , 2, 2, 0)
b) w2 = (π, 2π, 3π, 4π . . . )
c) w3 = (n, n − 1, n − 2, . . . , 2, 1) f ) w6 = (1, 3, . . . , n−1, n, n−2, . . . , 4, 2) (su-
d) w4 = (2n, 2n − 2, . . . , 4, 2) poniendo, aquí, que n es par)
cos(x/4) x
a) f (x) = c) f (x) =
ln(2 + x2 ) ln(2 + x)
1
b) f (x) = e−x
2 +2
sen(x/2) d ) f (x) = x + x2 − x
e
7. Utilizando las funciones de matlab de generación de matrices y vectores, crear las siguientes
matrices de dimensión N × N variable. Esto es, generar las matrices que se indican de forma que
las órdenes utilizadas puedan servir para cualquier valor de N :
à−N −N + 1 · · · −2 −1í
4 0 ··· 0 0 1 π ··· 0 0
á ë á ë
0 2 0 0 1 0 1 0 ··· 0 0
.. .. .. .. .. .. ..
A= . , B= . , C= 0 1 . 0 0
. . . .
0 2 0 0 1 π .. .. .. ..
. . . .
0 0 ··· 0 4 0 0 ··· 0 1 0 ··· ··· 1 0
8. La función rand(n, m) devuelve una matriz n×m de números aleatorios de distribución uniforme
entre 0 y 1 (esto es, que todos los números entre 0 y 1 tienen la misma probabilidad). Construir
un vector-fila v de números aleatorios de dimensión par N ≥ 10 . Extraer de él los siguientes
conjuntos de elementos:
a) w1 = (v1 , v3 , v5 , . . . , vN −1 )
b) w2 = (vN , vN −1 , . . . , v2 , v1 )
c) w3 = (v1 , v3 , v5 , . . . , vN −1 , v2 , v4 , . . . , vN )
d ) w4 = (v1 , v2 , v3 , . . . , v N , vN , vN −1 , . . . , v N +1 )
2 2
9. Utilizando la función rand, construir vectores-fila de números aleatorios con cada una de las
propiedades siguientes:
10. La función randi se puede utilizar para obtener matrices de números aleatorios enteros (ver en
el Help). Construir un vector-fila de números enteros aleatorios entre 1 y 6.
12. Generar un vector v de 20 números enteros aleatorios entre 1 y 15. Extraer las componentes de
v que sean iguales a 5 ó a 10.
13. Generar un vector v de números aleatorios de longitud 20. Buscando en la ayuda en línea de
MATLAB (help o helpwin) información sobre la función que se indica en cada caso, calcular
14. Siendo v un vector cualquiera de números enteros, escribir las órdenes adecuadas para obtener
lo que se pide en cada uno de los casos siguientes:
15. La función magic de MATLAB construye una matriz cuadrada conteniendo un cuadrado mágico
(ver Help). Dado un número n > 3 almacenado en una variable, construir una matriz M que
contenga un cuadrado mágico de dimensión n.
a) Comprobar mediante la función sum que todas las filas y columnas de M tienen la misma
suma.
b) Calcular el rango, el determinante y la traza de M (buscar en el Help las funciones adecua-
das).
c) Resolver el sistema lineal de ecuaciones M x = b, siendo b = (1, 2 . . . n)t . (Observación: la
matriz mágica construida es regular para n > 3 impar, y es singular para n par).
3 B
SE
UNIVER
VILLA
Gráficos e an
Ecuaciones Diferenciales
y Análisis Numérico
y
( x 7, y 7)
( x1 , y1 )
x1 x2 x3 x4 x5 x6 x7 x
La línea así obtenida tendrá mayor apariencia de “suave” cuanto más puntos se utilicen para construirla,
ya que los segmentos serán imperceptibles (véase la Figura 3.2).
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5
71
3. Gráficos 72
En esta sección se explican, muy brevemente, las formas más habituales de representación gráfica de
“objetos” matemáticos bi y tri-dimensionales.
Construir un conjunto de puntos (tantos como se quiera) en el intervalo [a, b], que serán las
abscisas de los puntos que determinan la poligonal a construir. Normalmente, dichos puntos se
toman regularmente espaciados y en número suficiente como para que la gráfica tenga aspecto
“suave”:
a = x1 < x2 < . . . < xn = b
Cuando una curva viene definida por una relación del tipo y = f (x) se dice que está definida de forma
explícita.
En ocasiones, una curva viene descrita por una relación, también explícita, pero del tipo:
0
−6 −4 −2 0 2 4 6 8
Una relación del tipo f (x, y) = 0 puede también representar, implícitamente, una curva: la formada
por los puntos (x, y) del plano sobre los cuales la función f toma el valor cero. Se puede dibujar esta
curva dibujando la curva de nivel k = 0 de la función f (ver Sección 3.1.7).
Otra forma de definir una curva plana es mediante sus ecuaciones paramétricas, en la cual los puntos
(x, y) que forman la curva vienen dados por dos funciones que dependen de una variable auxiliar:
t x y
10
Y 0 0 1
1 -1.5 2.4
8
4 3 2.6 7.0
2
4 6.3 6
t=0
X
5 7.9 3.1
0
6 6.8 1.1
−2
−2 0 2 4 6 8 10 12
7 5.0 1.7
Figura 3.4: Representación de la curva de ecuaciones paramé- 8 5.0 4.4
tricas x = t−3 sen(t), y = 4−3 cos(t) para t ∈ [0, 10]. Obsérvese
que no hay eje t. 9 7.8 6.7
10 11.6 6.5
Mediante ecuaciones paramétricas es posible describir muchas más curvas y más complicadas que
mediante una ecuación explícita. Algunas serían prácticamente imposibles de visualizar sin la ayuda
de herramientas gráficas informáticas (véanse Figuras 3.5 y Figuras 3.6).
2 30
1.5
20
10
0.5
0 0
−0.5
−10
−1
−20
−1.5
−2 −30
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −30 −20 −10 0 10 20 30
θ, que es el ángulo que forma el segmento OP con una semirrecta fija de origen O denominada
eje polar.
En tal sistema de coordenadas, se dice que el par (r, θ) son las coordenadas polares del punto P
(ver Figura 3.7).
p/2
y
(x ,y)
P (q, r) P
(q, r)
r r
y=rsen(q)
O q q x
0
Polo Eje polar x=rcos(q)
Figura 3.7: Sistema de coordenadas polares. Figura 3.8: Coordenadas cartesianas y pola-
res.
El paso de las coordenadas polares a cartesianas y viceversa se efectúa mediante las siguientes fórmulas,
tomando el polo como origen de coordenadas y el eje polar como semi-eje positivo de abscisas (ver la
Figura 3.8):
x = r cos(θ), y = r sen(θ);
p y
x2 + y 2 , θ = arctan
r=
x
Una relación del tipo r = f (θ) define de forma explícita una curva en coordenadas polares. Ver ejemplos
en las Figuras 3.9 y 3.10.
Los programas de que permiten realizar gráficas suelen disponer de las funciones adecuadas para
dibujar curvas utilizando directamente las coordenadas polares. En este caso habrá que proporcionar
90 90
15 2
120 60 120 60
1.5
10
150 30 150 1 30
5
0.5
180 0 180 0
270 270
Figura 3.9: Curva de ecuación, en coordena- Figura 3.10: Curva de ecuación, en coorde-
das polares, r = θ, θ ∈ [0, 9π/2] nadas polares, r = 2 sen(6θ), θ ∈ [0, 2π].
En caso de que no se disponga de dichas funciones, habrá que utilizar las fórmulas
xi = ri cos(θi ), yi = ri sen(θi )
La gráfica de una curva tridimensional se dibuja, igual que la bidimensional, uniendo mediante segmen-
tos rectos (en 3D) los puntos consecutivos de un conjunto discreto y ordenado. Mediante el software
adecuado, estos segmentos se “proyectan” sobre el plano del dibujo para obtener impresión tridimen-
sional.
La forma más sencilla de describir matemáticamente una curva tridimensional es mediante sus ecua-
ciones paramétricas. Estas ecuaciones describen los valores de las coordenadas (x, y, z) de cada punto
de la curva en función de una variable auxiliar, llamada parámetro:
x = f (t)
y = g(t) para t ∈ [a, b]
z = h(t)
Para dibujar su gráfica habrá, pués, que construir las coordenadas de un conjunto discreto y ordenado
de puntos de la curva. De forma similar a como se hizo en el caso bidimensional, el procedimiento es
el siguiente (véanse los ejemplos de las Figuras 3.11 y 3.12):
Construir un conjunto de valores del parámetro t ∈ [a, b]: {a = t1 < t2 < . . . < tn = b}
x = {x1 , x2 , . . . , xn }, xi = f (ti ), i = 1, 2, . . . , n
y = {y1 , y2 , . . . , yn }, yi = g(ti ), i = 1, 2, . . . , n
z = {z1 , z2 , . . . , zn }, zi = f (ti ), i = 1, 2, . . . , n
30
1
25
0.5
20
15 0
10
−0.5
5
0 −1
1 2
0.5 1 1.5 1
0.5 0.5
0 1
0 0
−0.5 0.5
−0.5 −0.5
−1 −1 0 −1
Figura 3.11: Gráfica de la curva 3D de ecua- Figura 3.12: Gráfica de la curva x(t) =
ciones paramétricas x(t) = cos(t), y(t) = cos(3t), y(t) = 2 cos2 (t), z(t) = sen(2t),
sen(t), z(t) = t, t ∈ [0, 8π]. t ∈ [−π, π].
La relación explícita
z = f (x, y)
con f : Ω ⊂ R2 7→ R, representa una superficie en el espacio R3 : a cada punto (x, y) del dominio Ω del
plano R2 , la función f le hace corresponder un valor z que representa la “altura” de la superficie en ese
punto.
Para dibujar la superficie es preciso disponer de una “discretización” del dominio Ω en el que está
definida la función, es decir un conjunto de polígonos (normalmente triángulos o rectángulos) cuya
unión sea Ω.
La forma de proporcionar los datos en uno y otro caso es diferente. Un mallado rectangular de un
dominio Ω = [a, b] × [c, d] queda definido mediante las particiones de los intervalos [a, b] y [c, d] cuyo
producto cartesiano produce los nodos de la malla: {x1 , x2 , . . . , xn } e {y1 , y2 , . . . , ym }.
Para definir un mallado mediante triángulos es preciso, por un lado numerar sus vértices y disponer de
sus coordenadas, (xi , yi ), 1 = 1, . . . , n y, por otro, numerar sus triángulos y describirlos enumerando,
para cada uno, sus tres vértices.
Elevando cada vértice del mallado según el valor de f en ese punto se consigue una representación de
la superficie como una red deformada, como en las Figuras 3.16 y 3.17.
Y
Y
X
X
Dar un color a cada arista dependiendo del valor de la función en sus extremos, como en la Figura 3.18,
puede resultar útil.
Rellenando de color cada retícula del mallado, la superficie se hace opaca. El color de las caras puede ser
constante en toda la superficie, como en la Figura 3.20, constante en cada cara, como en la Figura 3.21,
o interpolado, es decir, degradado en cada cara, en función de los valores en los vértices, como se hace
en la Figura 3.22.
Una superficie en el espacio de tres dimensiones pueden también venir definida mediante ecuaciones
paramétricas.
En este caso, para construir la gráfica de la superficie es preciso crear una discretización del dominio
donde varían los parámetros, [a, b]×[c, d], y utilizar las ecuaciones paramétricas para calcular los puntos
correspondientes sobre la superficie.
Figura 3.17: Red rectangular deformada. Figura 3.18: Red rectangular deformada. El
color de las aristas depende del valor de la
función.
Figura 3.19: Cara triangular rellena de color Figura 3.20: Mallado deformado con caras de
plano. color constante.
Figura 3.21: Mallado deformado con color Figura 3.22: Mallado deformado con caras de
plano en cada cara, dependiente de la altura. color interpolado a partir de los valores en los
vértices.
Una forma habitual de representar gráficamente los valores de una función de dos variables, f (x, y) = 0
es dibujando sus lineas o curvas de nivel.
Se llama curva de nivel de valor k de la función f (x, y) a la curva formada por los puntos del plano XY
sobre los cuales la función f toma el valor k, es decir la curva ímplicitamente definida por la ecuación
f (x, y) = k
El dibujo de las curvas de nivel correspondientes a un conjunto de valores k proporciona una buena
información del comportamiento de la función f .
30
0.32
25
0.3
20 0.28
0.26
15
0.24
10
0.22
5 0.2
0.18
5 10 15 20 25 30
Para dibujar una curva definida por una relación de la forma y = f (x) donde f : [a, b] 7→ R es una
función de una variable real dada por una expresión analítica, se puede usar la función fplot de
MATLAB, en alguna de estas versiones:
fplot(f)
fplot(f, [a,b])
Ejemplo 3.1
Representar, usando fplot, la gráfica de la función
x sen(x2 )
f (x) = en [−5, 5]
2
Por defecto, la función fplot hace variar la variable independiente en el intervalo [−5, 5]. Por lo
tanto, basta escribir en la línea de comandos
fplot(@(x) sin(x.^2).*x/2)
2
1.5
0.5
-0.5
-1
-1.5
-2
-5 -4 -3 -2 -1 0 1 2 3 4 5
Ejemplo 3.2
Representar, usando fplot, la gráfica de la función
f (x) = ln(x)
La función f (x) = ln(x) no está definida para valores de x ≤ 0. La función fplot se limitará a
dibujar la función en la parte del intervalo [−5, 5] en el que sí está definida: x ∈ (0, 5]. Escribe, en
la línea de comandos
0.5
-0.5
-1
-1.5
-2
Ejemplo 3.4
Representar, usando fplot, la gráfica de la función
x sen(x2 )
f (x) = para x ∈ [0, π]
2
f = @(x) x.*sin(x.^2)/2;
fplot(f, [0,pi])
0.5
-0.5
-1
fimplicit(f)
fimplicit(f, [a, b])
fimplicit(f, [a, b, c, d])
f es la función (de dos variables) que describe la ecuación, y debe ser descrita mediante un «handle»
(manejador) de una función de dos variables (de nuevo, anónima o M-función).
[a,b] (opcional) la curva se dibuja en el cuadrado del plano [a, b] × [a, b]. Si no se especifica,
fimplicit usa unos valores por defecto.
[a,b,c,d] (opcional) la curva se dibuja en el rectángulo del plano [a, b] × [c, d].
Ejemplo 3.5
Dibujar, usando fimplicit, la circunferencia x2 + y 2 = 7 en el cuadrado [−5, 5] × [−5, 5].
Se trata de una curva definida en forma implícita. Lo primero es escribir la ecuación en forma
homogénea, es decir, con segundo miembro nulo:
x2 + y 2 − 7 = 0.
De lo que se trata es de dibujar la curva formada por los puntos del plano sobre los cuales la
función de dos variables f (x) = x2 + y 2 − 7 toma el valor 0. Para ello basta con escribir en la
línea de comandos
-1
-2
-3
-4
-5
-5 -4 -3 -2 -1 0 1 2 3 4 5
fplot(f, g)
fplot(f, g, [a,b])
f, g son las dos funciones que definen la curva. Pueden ser especificadas, como antes, mediante
funciones anónimas o M-funciones.
Ejemplo 3.7
Dibujar, usando fplot, la curva definida por las ecuaciones (paramétricas)
se obtendrá la figura de la izquierda. Si t recorriera sólo el intervalo [0, π], con la orden
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
-0.2 -0.2
-0.4 -0.4
-0.6 -0.6
-0.8 -0.8
-1 -1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Ejercicio 3.8 Escribir un script que dibuje, usando ezplot y funciones anónimas, la curva definida
por las ecuaciones:
x = f (t) = cos(t) + 1/2 cos(7t) + 1/3 sen(17t)
, para t ∈ [0, 2π]
y = g(t) = sen(t) + 1/2 sen(7t) + 1/3 cos(17t)
x = f (t)
y = g(t) para t ∈ [a, b]
z = h(t)
fplot3(f, g, h)
fplot3(f, g, h, [a,b])
f,g,h son las tres funciones que definen la curva. Se proporcionan en la forma habitual.
[a, b] es el intervalo donde varía el parámetro t. Si no se especifican, fplot3 toma por defecto
[a, b] = [−5, 5].
Ejemplo 3.9
Dibujar, usando fplot3, la siguiente curva tridimensional:
√
x = 3 cos(t), y = t sen(t2 ), z= t, t ∈ [0, 2π]
2.5
1.5
0.5
5
3
2
0 1
0
-1
-5 -2
-3
fmesh(f)
fmesh(f, [a,b])
fmesh(f, [a, b, c, d])
[a, b] (opcional) establece el rectángulo para dibujar la superficie: [a, b] × [a, b]. Si no se especifica,
fmesh toma [a, b] = [−5, 5].
Ejemplo 3.11
2 2
Dibujar la superficie z = x e−(x +y )
0.2
0.1
-0.1
-0.2
-0.3
-0.4
3
2
3
1 2
0 1
-1 0
-1
-2 -2
-3 -3
La función fsurf se utiliza igual que fmesh, pero además, rellena de color sólido los paneles entre la
malla.
Ejemplo 3.12
2 2
Dibujar la superficie z = x e−(x +y )
0.2
0.1
-0.1
-0.2
-0.3
-0.4
3
2
3
1 2
0 1
-1 0
-1
-2 -2
-3 -3
fmesh(f, g, h)
fmesh(f, g, h, [a, b])
fmesh(f, g, h, [a, b, c, d])
f, g, h son las expresiones de las funciones de dos variables f (s, t), g(s, t) y h(s, t)
[a,b] (opcional) define el dominio en que varían s y t: s, t ∈ [a, b]. Si no se especifica, fmesh toma
[a, b] = [−5, 5].
[a,b,c,d] (opcional) define el dominio en que varían s y t: (s, t) ∈ [a, b] × [c, d], es decir, s ∈ [a, b]
y t ∈ [c, d].
Ejemplo 3.13
Dibujar la superficie cilíndrica definida por las ecuaciones paramétricas
x(t, ϕ) = (2 + cos(t)) cos(ϕ)
y(t, ϕ) = (2 + cos(t)) sen(ϕ) , t ∈ [0, 2π], ϕ ∈ [0, 2π]
z(t, ϕ) = t
zfun = @(t,s) t; 3
0
3
2
3
1 2
0 1
-1 0
-1
-2 -2
-3 -3
0
3
2
3
1 2
0 1
-1 0
-1
-2 -2
-3 -3
fcontour(f)
fcontour(f, [a,b])
2 +y 2 )
Ejercicio 3.15 Dibujar las lineas de nivel de la función z = x e−(x :
-1
-2
-3
-3 -2 -1 0 1 2 3
En el siguiente ejemplo se muestra cómo usar una opción de fcontour para rellenar de color los espacios
entre las curvas de nivel.
2 +y 2 )
Ejercicio 3.16 Dibujar las lineas de nivel de la función z = x e−(x :
-1
-2
-3
-3 -2 -1 0 1 2 3
Ejercicio 3.17
0.4
0.2
-0.1
-0.2
-0.3
-0.4
3
2
3
1 2
0 1
-1 0
-1
-2 -2
-3 -3
Ejercicio 3.18 Dibujar usando fmesh la superficie z = f (x, y) y las curvas de nivel de la función
f (x, y), siendo f (x, y) = sen(x/2) sen(y/2).
Observación
Las funciones f*** que se han presentado en esta sección han sido diseñadas para que el usuario pueda
obtener la gráfica de una función o dibujar una curva con suma facilidad y sin tener que preocuparse
en exceso de las exigencias del lenguaje.
Por ello, “corrigen” un buen número de fallos propios de los usuarios noveles, eventuales o poco cuida-
dosos.
Por ejemplo, no se produce un error al dar la orden fplot(@(x) sqrt(x),[-3,3]) a pesar de que la
función a dibujar no está definida para valores negativos de x. Como tampoco se produce al dar la
orden fplot(@(x) 1./x,[-10,10]), aunque esta función no está definida en x = 0.
Como en casi todas las cosas, lo que se gana en facilidad por un lado se paga en eficacia por otro.
En la sección siguiente se ven órdenes, más básicas que éstas, pero que permiten realizar gráficas más
personalizadas y menos estándar.
El uso de las funciones f*** debería limitarse a ocasiones en que se desea tener una idea rápida de la
forma de una curva o superficie, mediante una sola orden directamente en la línea de comandos.
plot(x,y)
siendo x e y dos vectores de las mismas dimensiones conteniendo, respectivamente, las abscisas y las
ordenadas de los puntos de la gráfica (los puntos (xi , yi ) de la Figura 3.1).
Ejemplo 3.19
x2 + 2
Usando la función plot, dibujar la curva y = para x ∈ [−2, 3]
x+5
Comenzamos por construir un soporte de nodos {xi } regularmente espaciados en el intervalo [−2, 3].
Para ello utilizamos la función linspace con la opción por defecto de 100 puntos (suficiente para
x2 + 2
la mayoría de las gráficas). Luego calculamos los valores de la función f (x) = en todos esos
x+5
puntos (usando operaciones vectoriales cuando haga falta) y dibujamos la curva.
y = (x.^2+2)./(x+5); 1.6
plot(x,y) 1.4
1.2
Es posible dibujar dos o más curvas a la vez, proporcionando a la función plot varios pares de vectores
(abscisas, ordenadas).
Ejemplo 3.20
1
Usando la función plot, dibujar las curvas y = x + x2 − x e y = ex + sen(x2 ) en el intervalo
e
[−π, π].
Puesto que hay que dibujar ambas curvas en el mismo intervalo [−π, π], usamos el mismo vector
de abscisas para las dos.
40
x = linspace(-pi, pi); 35
y2 = exp(x) + sin(x.^2); 25
15
10
-5
-4 -3 -2 -1 0 1 2 3 4
Ejemplo 3.21
Usando la función plot, dibujar los arcos de curva:
y = 1 + sen(x2 ) − x, x ∈ [−2, 6]
ex
y = ex + sen(x2 ) + x, x ∈ [−6, 2]
Las curvas a dibujar tienen distintos intervalos, luego tenemos que crear dos vectores. Vamos,
ahora, a definir las funciones que definen las curvas como funciones anónimas.
10
xf = linspace(-2, 6); 4
xg = linspace(-6, 2); 2
yf = f(xf); 0
yg = g(xg); -2
-6
-8
-6 -4 -2 0 2 4 6
cos(x/4)
f (x) = ,
ln(1 + x2 )
Ejemplo 3.23
Usando la función plot, dibujar la circunferencia de centro el punto (−1, 4) y radio 2.
Lo más cómodo para dibujar una circunferencia es utilizar sus ecuaciones paramátricas: los puntos
de la circunferencia de centro (a, b) y radio r vienen dados por
®
x = a + r cos(t)
t ∈ [0, 2π]
y = b + r sin(t)
x = -1 + 2*cos(t); 5
y = 4 + 2*sin(t); 4.5
plot(x, y, 'LineWidth', 2) 4
3.5
evitar esto.
Tabla 3.1: Símbolos para definir el color y el tipo de línea a dibujar en la orden plot
Ejemplo 3.26
Dibujar las parábolas y = x2 − 1 e y = −2x2 + x + 5 en el intervalo [−3, 3] usando distintos colores,
marcadores y tipo de líneas.
10
f = @(x) x.^2 - 1; 5
g = @(x) -2*x.^2 + x + 5;
x = linspace(-3,3, 40); 0
-10
-15
-20
-3 -2 -1 0 1 2 3
Algunas de las propiedades más utilizadas están descritas en la tabla 3.2. Este procedimiento se puede
combinar con el método abreviado descrito en la sección 3.3.2.
Tabla 3.2: Algunas propiedades de las líneas, controlables con la orden plot
Para información sobre el resto de propiedades de las líneas, se puede consultar, en el Help de matlab,
Line Properties.
Ejemplo 3.27
Dibujar la curva de ecuación y = ex + sen(x2 ) en el intervalo [−2, 2] usando distintos colores,
marcadores y tipo de líneas.
x = linspace(-2, 2, 30); 5
'LineWidth', 2) 0
-1
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
La orden
axis equal
establece el mismo factor de escala para ambos ejes. Esto haría, por ejemplo, que la gráfica de una
circunferencia se vea “redonda” en lugar de parecer una elipse.
La orden
axis off
axis on
Ejemplo 3.28
x2 + 2
Dibujar la curva y = para x ∈ [−2, 3], modificando los límites de los ejes.
x+5
y = (x.^2+2)./(x+5); 2
0.5
-0.5
-1
-3 -2 -1 0 1 2 3 4
Ejemplo 3.29
Dibujar la circunferencia de centro el punto (−1, 4) y radio 2, imponiendo el mismo factor de escala
en los ejes X e Y ..
t = linspace(0, 2*pi); 6
x = -1 + 2*cos(t);
y = 4 + 2*sin(t);
5
plot(x, y, 'LineWidth', 2) 4
axis([-3, 1, 1, 7])
axis equal 3
1
-4 -3 -2 -1 0 1 2
grid on
axis off
Las funciones
title('Texto de titulo')
xlabel('Etiqueta del eje OX')
ylabel('Etiqueta del eje OY')
incluyen, respectivamente, un título, una etiqueta asociada al eje OX (horizontal) y una etiqueta
asociada al eje OY (vertical).
La función
inserta en el gráfico una leyenda, usada normalmente para explicar el significado de cada curva. Las
leyendas de la función legend se asignan a las curvas en el orden en que se dibujaron.
El título y las etiquetas de los ejes admiten, entre otras, las propiedades contenidas en la tabla 3.3.
Para obtener información sobre las propiedades de la leyenda, buscar legend properties en el Help de
matlab.
Tabla 3.3: Algunas propiedades de los títulos y las etiquetas de los ejes
Ejemplo 3.30
Dibujar las parábolas y = x2 − 1 e y = −2x2 + x + 5 en el intervalo [−3, 3] usando distintos colores,
marcadores y tipo de líneas y añadiendo título, etiquetas y leyenda.
f = @(x) x.^2 - 1;
g = @(x) -2*x.^2 + x + 5;
x = linspace(-3,3, 40);
plot(x, f(x), 'r<', x, g(x), 'bo-', 'LineWidth', 2)
title('Dos parábolas', 'FontSize', 18, 'FontName', 'Verdana', ...
'FontWeight', 'normal', 'Color', [0.1, 0.3, 0.4])
xlabel('Eje OX', 'FontSize',18, 'Color', [0.1, 0.3, 0.4])
ylabel('Eje OY', 'FontSize',18, 'Color', [0.1, 0.3, 0.4])
legend('Parábola cóncava', 'Parábola convexa', 'Location', 'south')
Dos parábolas
10
0
Eje OY
-5
-10
-15
Parábola cóncava
Parábola convexa
-20
-3 -2 -1 0 1 2 3
Eje OX
hold on
hold off
La primera de ellas anula el borrado previo de la gráfica antes de un nuevo dibujo. La segunda devuelve
el sistema a su estado normal.
La orden hold on debe utilizarse después de haber realizado algún dibujo (el primero) o de haber
utilizado la orden axis([xmin, xmax, ymin, ymax]) para fijar los límites de los ejes en los que se va
a dibujar.
Ejemplo 3.31
x sin(x2 )
Dibujar la curva de ecuación y = para x ∈ [−4, 4], dibujando también los ejes de
2
coordenadas.
f = @(x) 0.5*x.*sin(x.^2);
x = linspace(-4, 4, 300);
axis([-4, 4, -2, 2])
hold on
plot(x, f(x), 'LineWidth', 2)
plot([-4, 4], [0, 0], 'k', 'LineWidth', 1.5)
plot([0,0], [-2, 2], 'k', 'LineWidth', 1.5)
hold off
axis off
Obsérvese que hemos utilizado un número de puntos superior al habitual para realizar la gráfica,
ya que la curva es muy oscilante.
text(x, y, texto_a_añadir)
inserta el texto texto_a_añadir en el punto (x, y) del gráfico. Se puede controlar el tamaño de la
letra con la propiedad ’FontSize’. Para obtener más información sobre otras (muchas) propiedades
aplicables al texto, consultar text properties en el Help de matlab.
Ejemplo 3.32
0.9
axis([0, 1, 0, 1])
0.8
'LineWidth', 2) 0.5
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Ejemplo 3.33
Uso de algunas opciones de text
x = linspace(-10, 10);
y = (x+1)./(x.^2+x+4); x
2
axis([-10, 10, -0.6, 0.6]) x
1
hold on
plot(x, y, 'r', 'LineWidth', 2)
plot([-10,10], [0, 0], 'k')
plot([0,0], [-1, 1], 'k')
plot([1, 1], [0, 1/3], 'k')
plot([-3, -3], [0, -1/5], 'k')
text(0.9, -0.04, 'x_1', 'FontSize', 14)
text(-3.1, 0.04, 'x_2', 'FontSize', 14)
ht = text(3, 0.3, '$f(x) = \frac{x+1}{x^2+x+4}$', 'FontSize', 18);
[Link] = 'latex';
figure
abre una nueva ventana gráfica, cuyo nombre será Figure n, donde n es el primer número entero positivo
que esté libre. También se puede abrir una directamente especificando su número: figure(7).
Cuando hay varias ventanas gráficas abiertas, la manera de designar cuál es la activa (en cuál se
dibujará) es usar la orden figure(n) con el número correspondiente.
La orden
clf
shg
coloca la ventana gráfica activa delante de todas las demás. Muy útil cuando se está escribiendo
programas que dibujan.
scatter(x, y)
dibuja los puntos (xi , yi ) con marcadores (por defecto, con círculos).
La función
stairs(x, y)
stem(x, y)
Ejemplo 3.34 x = e2 t sen(100 t)
Dibujar los puntos definidos por las ecuaciones paramétricas
y = e2 t cos(100 t)
para 300 valores de t comprendidos entre 0 y 1.
t = linspace(0,1,300);
x = exp(2*t).*sin(100*t);
y = exp(2*t).*cos(100*t);
scatter(x,y)
axis equal
axis off
Ejemplo 3.35
Dibujar datos como una función escalonada
x = linspace(-2, 2, 20); 4
stairs(x, g(x)) 3
-1
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
Ejemplo 3.36
Dibujar datos con tallos
x = linspace(-2, 2, 40); 3
stem(x, g(x)) 2
-1
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
subplot(n, m, k)
Esta orden divide el área de la ventana gráfica en una cuadrícula n×m (n filas y m columnas), y
convierte el k-ésimo cuadro en los ejes activos. matlab numera los cuadros de izquierda a derecha y
de arriba hacia abajo.
1 1
0.8 0.8
0.6 0.6
1 2
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
1 1
0.8 0.8
0.6 0.6
3 4
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Cada uno de los ejes así definidos funciona como si fuera una figura diferente: para dibujar en ellos hay
que hacerlos activos con la orden subplot; la orden hold actúa de forma independiente en cada uno
de ellos; cada uno puede tener sus propios título, etiquetas, leyendas, . . .
Ejemplo 3.37
Uso de subplots
x = linspace(0, 2*pi);
subplot(2,2,1)
axis([0, 2*pi, -1.5, 1.5])
hold on
y1 = sin(x);
plot(x, y1, 'Linewidth', 2, 'Color', [0.8, 0.3, 0.1])
subplot(2,2,2)
axis([0, 2*pi, -1.5, 1.5])
hold on
y2 = cos(2*x);
plot(x, y2, 'Linewidth', 2, 'Color', [0.9, 0.7, 0.1])
subplot(2,2,3)
axis([0, 2*pi, -1.5, 1.5])
hold on
y3 = y1.*y2;
scatter(x, y3, 'MarkerEdgeColor', [0, 0.5, 0.7])
subplot(2,2,4)
axis([0, 2*pi, -1.5, 1.5])
hold on
plot(x, y3, 'LineWidth', 2)
stem(x, y1)
stem(x, y2)
1.5 1.5
1 1
0.5 0.5
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
0 2 4 6 0 2 4 6
1.5 1.5
1 1
0.5 0.5
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
0 2 4 6 0 2 4 6
mesh(x, y, z)
siendo x, y y z tres matrices de las mismas dimensiones conteniendo, respectivamente, las coordenadas
(x, y, z) de los nodos de la malla.
Para construir las matrices x, y y z puede ser útil la función
sirve para construir la malla de base, en el plano XY: Si xp es una partición del intervalo [x0 , x1 ] e
yp es una partición del intervalo [y0 , y1 ], entonces x e y son dos matrices, que definen una malla del
rectángulo [x0 , x1 ] × [y0 , y1 ] . La matriz x contiene las coordenadas x de los puntos de la malla y la
matriz y sus coordenadas y.
Ejemplo 3.38
2 2
Dibujar la función z = x e−(x +y ) para (x, y) ∈ [0, π] × [−π/2, π/2] en forma de una malla defor-
mada, con 20 nodos en la dirección x y 30 nodos en la dirección y.
z = x .* exp(-(x.^2+y.^2));
mesh(x, y, z);
axis([0, pi, -pi/2, pi/2, 0, 1.5])
xlabel('X')
ylabel('Y')
zlabel('Z')
Por defecto, las caras de la malla son opacas. Se pueden hacer transparentes usando el parámetro
opcional
mesh(x, y, z, 'FaceAlpha', t)
siendo t un número real positivo en [0, 1]. Para otros parámetros opcionales, véase el Help.
La siguiente variante de mesh hace lo mismo pero, además, dibuja las líneas de nivel de f (x, y) en el
plano XY .
meshc(x, y, z)
Ejemplo 3.39
La misma superficie del Ejemplo 3.38, pero usando las opciones de tansparencia y de líneas de
nivel.
Existe otra forma de proporcionar los datos a la función mesh, en la que x e y son vectores, de longitudes
respectivas n y m, y z es una matriz de dimensión m×n, de forma que z(i,j) es el valor de f en (x(j),
y(i)).
La función surf funciona igual que mesh, pero rellena las caras de la malla de un color sólido, que
varía en función de los valores de f (x, y). También tiene la opción de transparencia y la variante que
dibuja las líneas de nivel.
surf(x, y, z)
surf(x, y, z, 'FaceAlpha', t)
surfc(x, y, z)
Por defecto, surf dibuja las caras de color constante y le superpone la malla. La orden siguiente se
puede usar para cambiar este comportamiento
La función contour dibuja las proyecciones sobre el plano XY de las líneas de nivel (isovalores).
Aquí, las variables x, y y z son iguales que para las funciones mesh y surf. La variable n indica el
número de líneas de nivel a dibujar. Si no se indica, se dibujan 8 líneas.
1.5
0.5
-0.5
-1
-1.5
0 0.5 1 1.5 2 2.5 3
Se puede indicar que se dibujen, solamente, las líneas de nivel correspondientes a unos valores dados:
contour(x, y, z, v)
0.5
-0.5
-1
-1.5
0 0.5 1 1.5 2 2.5 3
Si se desea dibujar solo una línea de nivel, se debe indicar el valor correspondiente por repetido. Por
ejemplo, para dibujar solo la linea f (x, y) = k, se debe escribir v = [k, k].
La opción siguiente escribe, junto a cada línea de nivel, el valor al que corresponde.
1.5
1
2
1.
0.5
1.2
-0.5
1.2
2
1.
-1
-1.5
0 0.5 1 1.5 2 2.5 3
La variante contourf dibuja las líneas de nivel en el plano XY y rellena de color sólido los espacios
entre ellas:
1.5
0.5
-0.5
-1
-1.5
0 0.5 1 1.5 2 2.5 3
Se puede usar la función colorbar para añadir a la gráfica una barra con la graduación de colores y la
correspondencia con los valores de la función representada.
colorbar
Esta función puede ser usada con cualquier otra función gráfica que utilice colores. Dispone de pará-
metros opcionales para indicar su ubicación en la gráfica (ver Help).
En los dibujos tridimensionales, esta función permite elegir el punto de observación del objeto repre-
sentado.
view( v )
siendo v un vector con tres componentes que indica el punto, en el espacio 3D, en el que se sitúa el
observador, siempre suponiendo que el objeto se sitúa en el origen de coordenadas. MATLAB calculará
el azimut y la elevación del vector unitario que apunta al punto v.
colormap(m)
En general, m es una matriz con tres columnas, de modo que la i-ésima fila determina el color, en
coordenadas RGB, del i-ésimo color utilizado. Para más comodidad, MATLAB dispone de una serie
de mapas de colores pre-definidos, que se pueden imponer mediante
colormap('Nombre_del_mapa')
colormap Nombre_del_mapa
3.5 Ejercicios
1. Usando las funciones ez*** adecuadas, representar gráficamente los “objetos” definidos por las
siguientes funciones, expresiones y/o ecuaciones:
1
a) y = cos(2x) + sen(x/2), x ∈ [−π, π]
2
sen2 x
b) y = x ∈ [−2, 2].
x2 − 1
c) x4 + y 3 = 24, x, y ∈ [−6, 6]
d ) y 2 = x4 (1 − x)(1 + x), x, y ∈ [−2, 2]
x = t − 3 sen t,
e)
y = 4 − 3 cos t, t ∈ [0, 10]
x = 4 cos t + cos(4t),
f)
y = 4 sen t − sen(4t), t ∈ [0, 2 π]
x = et/4 sen(2t),
g) y = et/4 cos(2t), t ∈ [−9, 10]
z = t/4
2. Para cada una de las funciones siguientes, definir funciones anónimas que las representen (deben
admitir vectores como argumentos), y dibujar sus gráficas en un intervalo razonable usando la
orden plot. Dibujar también, en cada caso los ejes coordenados.
1
…
1−x e) f (x) =
a) f (x) =
1+x x
x 1
Å ã
b) f (x) = f ) f (x) = ln
ln(x) x−2
√ √ √
c) f (x) = 1 − x + x − 2 g) f (x) = 25 − x2
√
d ) f (x) = 3 h) f (t) = sen2 ( 3 + x)
3. Representar en los mismos ejes los siguientes pares de funciones, en los límites que se indican en
cada caso, de modo que cada curva se distinga perfectamente, tanto en color como en blanco y
negro. Añadir leyendas para identificar las curvas. Añadir también un título y etiquetas en los
ejes.
sen(πx)
h) g(x) = √ y h(x) = 5e−2x ln(x2 + 1) en [0, 2] × [−0.5, 1]
1 + 2x
(x + 10)(x − 10)(18 − x)
i) f (x) = ln(x2 − 27) y g(x) = en [−15, 20] × [−20, 15]
100
4. Reproducir el dibujo siguiente. Las circunferencias grandes tienen radio = 1 y la pequeña tiene
radio = 0.5.
Circunferencia de radio 2
4 B
SE
UNIVER
VILLA
Programación con MATLAB
e an
Ecuaciones Diferenciales
y Análisis Numérico
Este capítulo está dedicado a los conceptos e instrucciones básicas que permiten la escritura de pro-
gramas.
Comenzamos mostrando las instrucciones básicas que permiten que un programa “se comunique” con el
usuario: instrucciones mediante las cuales podemos introducir, de forma interactiva, un dato para uso
del programa; e instrucciones que permiten que el programa pueda escribir información y/o resultados
en la pantalla del ordenador.
Los condicionales y los bucles o repeticiones son la base de la programación estructurada. Sin ellas,
las instrucciones de un programa sólo podrían ejecutarse en el orden en que están escritas (orden
secuencial). Las estructuras de control permiten modificar este orden y, en consecuencia, desarrollar
estrategias y algoritmos para resolver los problemas.
Los condicionales permiten que se ejecuten conjuntos distintos de instrucciones, en función de que
se verifique o no determinada condición.
Los bucles permiten que se ejecute repetidamente un conjunto de instrucciones, ya sea un número
pre-determinado de veces, o bien mientras que se verifique una determinada condición.
var = input('Mensaje')
imprime Mensaje en la pantalla y se queda esperando hasta que el usuario teclea algo en el teclado,
terminado por la tecla Retorno. Lo que se teclea puede ser cualquier expresión que use constantes
y/o variables existentes en el Workspace. Puede ser tambien un vector o matriz. El resultado de esta
expresión se almacenará en la variable var. Si se pulsa la tecla Retorno sin teclear nada se obtendrá
una matriz vacía: [ ].
El caracter de escape \n introducido en el texto del mensaje provoca una salto de línea.
Si lo que se quiere leer son caracteres, hay que encerrarlos entre apóstrofes. O bien se puede añadir el
argumento 's' a la orden input, lo que indica que se van a leer caracteres y, en ese caso, no hay que
encerrar lo que se escriba entre apóstrofes.
113
4. Programación con MATLAB 114
disp(algo)
Si se quieren mezclar, en una línea a imprimir con disp, texto y números, hay que formar con ellos
un vector, pero, para ello, hay que transformar el valor numérico en una cadena de caracteres que
represente ese valor, como en los siguientes ejemplos.
>> x = pi;
>> disp(['El valor de x es: ',num2str(x)])
El valor de x es: 3.1416
disp('')
disp(' HOLA !!! ')
nombre = input('Escribe tu nombre: ', 's');
disp('')
disp([' HOLA ', nombre, ' !!!'])
disp([' Hoy es ', date])
disp(' ')
disp(' >> ¡ADIOS! <<')
Esta orden permite controlar la forma en que se imprimen los datos. Su sintaxis para imprimir en la
pantalla es
donde:
lista_de_datos son los datos a imprimir. Pueden ser constantes y/o variables, separados por comas.
formato es una cadena de caracteres que describe la forma en que se deben imprimir los datos. Puede
contener combibaciones de los siguientes elementos:
Normalmente el formato es una combinación de texto literal y códigos para escribir datos numéricos,
que se van aplicando a los datos de la lista en el orden en que aparecen.
En los ejemplos siguientes se presentan algunos casos simples de utilización de esta instrucción. Para
una comprensión más amplia se debe consultar la ayuda y documentación de MATLAB.
>> x = sin(pi/5);
>> y = cos(pi/5);
>> fprintf('Las coordenadas del punto son x= %10.6f e y=%7.3f \n', x, y)
Las coordenadas del punto son x= 0.587785 e y= 0.809
Observamos que el primer código %10.6f se aplica al primer dato a imprimir, x, y el segundo
código %7.3f al segundo, y.
• el código %3i indica que se escriba un número entero ocupando un total de 3 espacios,
• el código %15.7e indica que se escriba un número en formato exponencial ocupando un total
de 15 espacios, de los cuales 7 son para los dígitos a la derecha del punto decimal.
instruccion-anterior
if expresion
bloque-de-instrucciones
end
instruccion-siguiente
si x < −1,
®
x+1
f (x) =
1 − x2 si x ≥ −1.
instruccion-anterior
if expresion
bloque-de-instrucciones-1
else
bloque-de-instrucciones-2
end
instruccion-siguiente
El siguiente ejercicio es una variante del anterior, en la que se considera una M-función en lugar de
un script. Los datos necesarios (n1 y n2) se obtienen ahora a través de los argumentos de entrada (en
lugar de leerlos del teclado) y el resultado (los datos ordenados) se obtienen a través de la varible de
salida v (en lugar de imprimirlo en la pantalla).
Escribir una M-función que reciba como argumentos de entrada dos números x, y ∈ R y devuelva
un vector cuyas componentes sean los dos números ordenados en orden creciente.
if n1 > n2
v = [n2, n1];
else
v = [n1, n2];
end
Escribir una M-función que, dados dos números x, y ∈ R, devuelva el número del cuadrante del
plano OXY en el que se encuentra el punto (x, y). Si el punto (x, y) está sobre uno de los ejes de
coordenadas se le asignará el número 0.
true
if x*y == 0 x*y == 0 n = 0 fin
n = 0;
false
return
end
false true
n = 3 n = 4
if (x > 0) x > 0
n = 4;
true true
if (y > 0) n = 2 y > 0 y > 0 n = 1
n = 1;
end
else
n = 3;
if (y > 0)
n = 2;
Salida: n
end
end
instruccion-anterior
if expresion-1
bloque-de-instrucciones-1
elseif expresion-2
bloque-de-instrucciones-2
else
bloque-de-instrucciones-3
end
instruccion-siguiente
function Donde(x, a, b)
%
% Donde(x, a, b) escribe en la pantalla la ubicacion
% de x respecto del intervalo [a,b]:
% - si x < a
% - si a <= x < = b
% - si x > b
%
if (a <= x) && (x <= b)
fprintf('El punto %6.2f pertenece al intervalo [%6.2f,%6.2f]\n', x,a,b)
elseif x < a
fprintf('El punto %6.2f esta a la izq. del intervalo [%6.2f,%6.2f]\n', x,a,b)
else
fprintf('El punto %6.2f esta a la dcha. del intervalo [%6.2f,%6.2f]\n', x,a,b)
end
end
siguiente:
instrucciones
Al comienzo, se evalúa expresion. Si el resultado es
false, se termina la repetición. En este caso, no se
ejecuta el bloque de instrucciones.
Si, por el contrario, el resultado de expresion es
true, se ejecuta el bloque de instrucciones. Cuan-
Figura 4.4: Diagrama de flujo de la repetición
do se termina, se vuelve a evalúar la expresion y se
condicional
vuelve a decidir.
instruccion-anterior
while expresion
bloque-de-instrucciones
end
instruccion-siguiente
n = n - 2; true
n = n + 2
end
suma = suma + n
n = n - 2
Salida: n
fin
Ejercicio 4.14 Partiendo de la M-función SumaImpares del ejemplo anterior, escribe una
M-función que reciba como argumento de entrada un número positivo M y devuel-
va el mayor entero n tal que la suma de todos los números naturales impares entre
1 y n sea menor o igual que M .
%---------------------------------------------------------------
% Script Adivina
% Este script genera un numero aleatorio entre 1 y 9
% y pide repetidamente que se teclee un numero hasta acertar
%---------------------------------------------------------------
%
n = randi(9);
fprintf('\n')
fprintf(' Tienes que acertar un numero del 1 al 9 \n') inicio
M = 0; M = 0
while M ~= n
M = input('Teclea un numero del 1 al 9 : '); false
M ~= n
end
beep true
fprintf('\n') Leer M
fprintf(' ¡¡¡¡Acertaste!!!! \n')
Imprimir: ’Acertaste!’
fin
Observación: la orden beep utilizada en el código anterior genera un pitido (si los altavoces están
activados).
Ejercicio 4.16 Partiendo del script Adivina del ejemplo anterior, escribe una M-función que
reciba como argumento de entrada un número natural m, genere de forma aleatoria un número
natural n entre 1 y m, y pida repetidamente al usuario que escriba un número en el teclado hasta
que lo acierte.
Este tipo de estructuras permite implementar bucles en los que se conoce a priori el númeo de veces
que hay que repetir un cierto bloque de instrucciones.
En muchas ocasiones, las repeticiones de un bucle dependen en realidad de una variable entera cuyo
valor se va incrementando hasta llegar a uno dado, momento en que se detienen las repeticiones. Esto
sucede, especialmente, con los algoritmos que manipulan vectores y matrices, que es lo más habitual
cuando se programan métodos numéricos.
En estas ocasiones, para implementar el bucle es preferible utilizar la estructura que se expone en esta
sección: bucle indexado. En MATLAB (igual que en muchos otros lenguajes) este tipo de mecanismos
se implementan mediante una instrucción denominada for.
Por ejemplo, el bucle representado mediante el diagrama siguiente
k = 1
false
k <= n
true
instrucciones
k = k + 1
instruccion-anterior
for k = 1 : n
bloque-de-instrucciones
end
instruccion-siguiente
Se observa que, con esta instrucción, no hay que ocuparse ni de inicializar ni de incrementar dentro del
bucle la variable-índice, k. Basta con indicar, junto a la cláusula for, el conjunto de valores que debe
tomar. Puesto que es la propia instrucción for la que gestiona la variable-índice k, está prohibido
modificar su valor dentro del bloque de instrucciones.
El conjunto de valores que debe tomar la variable-índice k tiene que ser de números enteros, pero
no tiene que ser necesariamente un conjunto de números consecutivos. Son válidos, por ejemplo, los
conjuntos siguientes:
Observación. Siempre que en un bucle sea posible determinar a priori el número de veces que se va a
repetir el bloque de instrucciones, es preferible utilizar la instrucción for, ya que la instrucción while
es más lenta.
k = 1
false
k <= longitud(v)
true
k = k + 1
Salida: v
Observación: Lo que es importante observar y comprender aquí es que en este tipo de bucle, no
hay que ocuparse explícitamente de la inicialización de la variable-índice k ni de su incremento en
cada iteración. La instrucción for se ocupa de todo cuando se proporciona la lista de valores que
debe tomar k.
En el programa Traza se hace uso de la instrucción error, cuya descripción se encuentra en el apar-
tado 4.4.
Los bucles y los condicionales se pueden anidar unos dentro e otros, como en los ejemplos que siguen:
que reciba como argumento una matriz A y devuelva el máximo entre los valores absolutos de todas
sus componentes
Ejercicio 4.21 Modifica la función MayorM del ejemplo anterior, para que devuelva, además, los
números de la fila y columna en que se produce el máximo.
return
que devuelve el control a la ventana de comandos o bien al programa que llamó al actual.
En otras ocasiones es necesario interrumpir la ejecución de un bucle de repetición en algún punto
interno del bloque de instrucciones que se repiten. Lógicamente, ello dependerá de que se verifique o
no alguna condición. Esta interrupción puede hacerse de dos formas:
Las órdenes respectivas en MATLAB, tanto en un bucle while como en un bucle for):
break
continue
true true
instrucciones-1 instrucciones-1
false false
instrucciones-2 instrucciones-2
Figura 4.5: Ruptura de bucle break: Si en algún Figura 4.6: Ruptura de bucle continue: Si en al-
momento de un bucle de repetición (de cualquier gún momento de un bucle de repetición (de cual-
tipo) se llega a una instrucción break, el ciclo quier tipo) se llega a una instrucción continue,
de repetición se interrumpe inmediatamente y la iteración en curso se interrumpe inmediata-
el programa continúa ejecutándose por la ins- mente, y se inicia la siguiente iteración (si pro-
trucción siguiente a la cláusula end del bucle. Si cede).
se trata de un bucle indexado for, la variable-
índice del mismo conserva el último valor que
haya tenido.
warning('Mensaje')
error('Mensaje')
• A los argumentos de salida de una M-función siempre hay que darles algún valor.
• Las instrucciones for, if, while, end, return, ... no necesitan llevar punto y coma al final, ya
que no producen salida por pantalla. Sí necesitanGENERALES
COMENTARIOS llevarlo, en general, las instrucciones incluidas
dentro.
• Las
1. M-funciones
Lo mínimo noquenecesitan
se le pidellevar
a un end al final.
programa De llamarse
para hecho, a así
nivel
es de
queeste
no curso, que lo lleven
produzca
puedeerrores
inducir de
a error. Mejor
ejecución. no ponerlo.
• Todos
2. A los
los programas
argumentos deben incluir
de salida de unas líneas de siempre
una !"#$%&'# comentarioshayque expliquen
que de forma
darles algún sucinta
valor.
su cometido y forma de utilización. En las M-funciones, estas líneas deben ser las siguientes a la
instrucción
3. Hay que comprobar
function, ya que de esta manera
los programas que son utilizadas
se hacen, por MATLAB
dándoles valores acomo el texto de help
los argumentos
de la para verificar
función. En lostodos losdeben
scripts casos serque
las podrían
primeras darse. Esto
líneas del forma parte del proceso de
archivo.
diseño y escritura del programa.
• Es bueno, en la medida de lo posible y razonable, respetar la notación habitual de la formulación
[Link] instrucciones
problema los, &!
al elegir!'( , )*&+,de
nombres , ,#-
las, variables.
(,%"(# ...Ejemplos
no necesitan llevar punto y coma al
final, ya que no producen salida por pantalla. Necesitan llevarlo, en general, las
• Llamar S a una
instrucciones suma y dentro.
incluidas P a un producto,
Ejemplo: mejor que G o N.
• Llamar i, j, k, l, m, n a un número entero, mejor que x ó y.
!"#$%&'##############(#)*+,-#.#/-0$#!++1/12$3!-#
• Llamar T o temp###/4$5&'#
a una temperatura, mejor que P.
1+6'#################(#)*+,-#.#/-0$#!++1/12$3!-#
• Llamar e, err ó error a un error, mejor que vaca.
5. Lashay
• Cuando M-funciones
errores ennoun
necesitan llevar
programa, hay,#- al final.
que De hecho,
leer con a nivel
atención losde este curso,
mensajes que lo
de error que
lleven puede inducir a error. Mejor no ponerlo.
envía matlab. Contienen mucha información, que en la mayoría de los casos es suficiente para
ayudar a detectar el error.
6. Todos los programas deben incluir una líneas de comentarios que expliquen de forma
sucinta
• Hay que su cometido
comprobar y forma de
los programas queutilización. Para insistir
se hacen, dándoles en este
valores a lospunto, en lospara
argumentos exámenes
verificar
se penalizará que no los lleven y/o que no sean adecuados.
que funciona bien en todos los casos que puedan darse. Esto forma parte del proceso de
diseño y escritura
7. Los espacios del programa.
en blanco, en general, hacen el código más legible y por tanto más fácil de
comprender y corregir. Se debe dejar uno o dos espacios en blanco entre el caracter %
• Los espacios en blanco, en general, hacen el código más legible y por lo tanto más fácil de
de cada linea de comentario y el texto del comentario.
comprender y corregir. Se deben dejar uno o dos espacios entre el carácter % de cada línea de
comentario y ellos
Compara texto
dosdel comentario
códigos (comparénse los dos códigos siguientes).
siguientes:
"*+/,!-+78$98:;4816!$2<=># "*+/,!-+#78$98:;4816!$2<=>#
(8$4816!$2<=>##61?*1@?1#@$#016!$#$3!,01,!/$# (HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
(61#@-2#+*013-2#+$,*3$@12#01+-312#-#!:*$@12## (##8$4816!$2<=>##61?*1@?1#@$#016!$#
(A*1#=# (#####$3!,01,!/$#61#@-2#+*013-2#
(78$981;4816!$2<=>#61?*1@?1#,$0&!1+#@$# (#####+$,*3$@12#01+-312#-#!:*$@12#A*1#=#
016!$## (##78$981;4816!$2<=>#61?*1@?1#,$0&!1+#@$#
(:1-01,3!/$# (#####016!$#:1-01,3!/$#
8$4B'# (HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
8:4C'# (#
"-3#!4CD=# ##
8$48$5!'# 8$#4#B'#
8:48:E!'# 8:#4#C'#
1+6# ##
=C4CF='# "-3#!4CD=#
8$48$E=C'# ##8$#4#8$5!'#
8:48:G<=C>'# ##8:#4#8:E!'#
# 1+6#
##
=C##4#CF='#
8$#4#8$E=C'#
8:#4#8:G<=C>'#
#
8. Es bueno respetar, en la medida de lo posible y razonable, la notación habitual de la
formulación de un problema. Ejemplos:
R. Echevarría a. Llamar S a una suma, mejor que G Dpto. EDAN - Universidad de Sevilla
b. Llamar P a un producto, mejor que N
c. Llamar T ó temp a una temperatura, mejor que P
d. Llamar i j k l m n a un número entero, como un índice, mejor que x ó y
4. Programación con MATLAB 134
4.7 Ejercicios
En este listado de ejercicios están incluidos, también, los ejemplos que se han ido mostrando a lo largo
de este capítulo.
En todos estos ejercicios se debe, además, dibujar el diagrama de flujo del algoritmo correspondiente.
1. (HolaMundo.m) Escribir un script que pida el nombre al usuario, escriba un mensaje de saludo
“Hola + nombre”, y escriba la fecha actual.
2. (LeeNumero.m) Escribir un script que lea un numero y escriba un mensaje “El numero leido
es + numero” en una sola linea.
3. (Maximo.m) Usando la instrucción if - end (sin cláusula else), escribir un script que lea dos
números e imprima el máximo entre ellos con un mensaje: “El maximo es + numero”.
4. (Ordena.m) Usando la instrucción if - else - end, escribir un script que lea dos números y
los imprima en pantalla ordenados en orden creciente, con un mensaje “Los numeros ordenados
son + numeros”.
que reciba como argumentos dos numeros y devuelva el mayor de ellos. Se trata del mismo
algoritmo que el ejercicio 3, pero en versión M-función.
que reciba como argumentos dos numeros y devuelva un vector-fila cuyas dos componentes sean
los argumentos de entrada ordenados en orden creciente. Mismo algoritmo del ejercicio 4, en
versión M-función.
7. (Cuadrante.m) Escribe una M-función que reciba como argumentos las coordenadas de un
punto (x, y) del plano y devuelva un número (1,2,3 o 4) que indique en qué cuadrante del plano
se encuentra. El programa debe devolver 0 si el punto se encuentra sobre alguno de los ejes
coordenados.
function Donde(x, a, b)
que reciba como argumentos de entrada un número x ∈ R y los extremos de un intervalo [a, b],
y escriba en la pantalla un mensaje indicando si x < a, x ∈ [a, b] o x > b.
a) Si alguna de las longitudes recibidas toma un valor menor o igual que cero.
b) Si el radicando es negativo, lo cual indica que no existe ningún triángulo con esos lados.
que calcule y devuelva el mayor número impar n para el cual la suma de todos los números
impares entre 1 y n es menor o igual que m.
function Adivina(m)
que genere un número aleatorio entre (por ejemplo) 1 y 25, y le pida repetidamente al usuario
que escriba un número en el teclado, hasta que lo acierte.
que reciba como argumentos de entrada una matriz A y un vector v y devuelva el vector w = A v,
calculado mediente bucles de repetición. El programa debe limitarse a hacer el cálculo y devolverlo
a través de la variable de salida w, sin imprimir nada en la pantalla, a excepción, eventualmente,
de mensajes de error y/o de advertencia.
que reciba como argumento de entrada una matriz A y devuelva su transpuesta, calculada me-
diente bucles de repetición.
16. (Grado1.m) Escribe un script que pida al usuario los coeficientes, a y b, de una ecuación de
primer grado ax + b = 0 e imprima su solución. El script debe asegurarse de que a es distinto de
cero y, en caso contrario, volver a pedirlo.
17. (Grado2.m) Escribe un script que pida al usuario los coeficientes, a, b y c, de una ecuación de
segundo grado ax2 + bx + c = 0 e imprima sus soluciones, indicando cuál de los tres casos posibles
se da:
que reciba como argumento un numero entero positivo n y calcule su factorial por multiplicaciones
sucesivas. La función debe escribir una tabla conteniendo, en cada linea, el valor de la iteración
corriente y el valor calculado hasta dicha iteración.
que reciba como argumento un número N y devuelva dos valores: la media aritmética Ma y la
media geométrica Mg de los números naturales menores o iguales que N .
20. (Fibonacci.m) La sucesión de Fibonacci surgió en el siglo XIII cuando el matemático italiano
Leonardo Pisano Fibonacci estudió un problema relativo a la reproducción de los conejos. Esta
sucesión {fn }n≥1 está definida para n ≥ 1 por fn+2 = fn+1 + fn , siendo f1 = f2 = 1. Escribir
una M-función
que reciba como argumento un número entero N > 2 y devuelva el término N -ésimo de la sucesión
de números de Fibonacci.
Para cada una de las series infinitas de los ejercicios 21 a 24 que siguen, cuyas sumas se indican en
cada caso, se pide escribir una M-función que reciba un valor de x y un valor de N y calcule y devuelva
un valor aproximado, SN , de la suma de la serie, sumando sus N + 1 primeros términos (desde i = 0
hasta i = N ):
XN
SN = ti (ti es el término general de la serie)
i=0
El programa debe, también, imprimir en pantalla una tabla, con una línea por iteración, con los valores
de: (a) número de la iteración,i; (b) valor del término sumado en esa iteración, ti ; (c) valor de la suma
parcial, Si ; (d) error, |Si − f (x)|.
21. (SumaGeom.m)
∞
1 X
= 1 + x + x2 + x3 + · · · + xn + · · · = xi ∀x tal que |x| < 1
1−x
i=0
22. (Sumaxcos.m)
∞
x3 x5 x2n+1 X x2i+1
x cos(x) = x − + + · · · + (−1)n + ··· = (−1)i ∀x ∈ R
2! 4! (2n)! (2i)!
i=0
23. (SumaExp.m)
∞
x2 x3 xn X xi
ex = 1 + x + + + ··· + + ··· = ∀x ∈ R
2! 3! n! i!
i=0
24. (SumaLog.m)
∞
x3 x5 x2n+1 x2i+1
…
1 + x X
ln =x+ + + ··· + + ··· = ∀x tal que |x| < 1
1−x 3 5 2n + 1 2i + 1
i=0
que reciba como argumentos de entrada un vector v de números enteros y un número entero n,
y devuelva el valor lógico true si alguna de las componentes de v es igual a n, y el valor lógico
false en caso contrario. No es necesario que el programa compruebe que los datos de entrada
son, efectivamente, números enteros.
que devuelva, para cada par (x,y), un valor zona de acuerdo con la regla siguiente:
que reciba como argumento de entrada un vector v de números enteros y devuelva dos vectores
w1 y w2 construidos como sigue:
que reciba como argumento de entrada un vector v y proporcione como salida otro vector w
formado por las componentes de v que tengan un valor mayor o igual que la media aritmética de
todas las componentes de v.
que reciba como argumento de entrada un número entero n positivo y devuelva el valor lógico
true si n es primo y false si no lo es. El programa debe comprobar que el valor de n recibido es
realmente un entero positivo y, en caso contrario, emitir un mensaje de “warning” y devolver el
valor lógico false.
Para decidir si un número divide a otro se puede utilizar la función intrínseca de matlab
mod(x,y) que devuelve el resto de la división entera x/y. Para decidir si un número es ente-
ro se puede utilizar la función round(x) que devuelve el número entero más próximo a x.
que reciba como argumento de entrada un vector v y proporcione como salida un entero index
definido como sigue:
• Si todas las componentes del vector v son del mismo signo (todas positivas o todas negativas),
entonces index =0.
• Si no todas las componentes del vector v son del mismo signo, index es el menor valor que
verifica v(index) ∗ v(index + 1) ≤ 0.
Por ejemplo, si v = (2, 5.4, 6.12, 7), entonces debe ser index = 0; si v = (3.2, 1.47, −0.9, 0.75),
entonces debe ser index = 2.
que determine si un año dado es bisiesto o no, teniendo en cuenta que son años bisiestos los
múltiplos de 4, excepto los que son también múltiplos de 100 pero no lo son de 400.
que permita validar una fecha dia-mes-anno. Debe tenerse en cuenta la posibilidad de que sea
un año bisiesto. Para comprobar si lo es, se puede hacer uso de la función EsBisiesto anterior.
function [v]=Divisores(N)
que devuelva, como elementos del vector v, los divisores del numero N.
Escribir un programa que determine si un número natural N dado es perfecto o no. Se dice que
un número es perfecto cuando es igual a la suma de sus divisores, excluido él mismo. Por ejemplo
6 es perfecto, ya que 6 = 1 + 2 + 3.
35. (Goldbach.m) La conjetura de Goldbach es uno de los problemas abiertos más antiguos de la
Teoría de Números. Dice que: “Todo número par mayor que 2 se puede obtener como suma de
dos números primos”.
Escribir un programa
que descomponga un número par N dado en la suma de dos números primos p=[p1,p2] (puede
ser el mismo primo repetido)
36. (Descomponer.m) Escribir un programa que reciba como argumento un número entero y lo
descomponga en unidades, decenas, centenas, etc. (No se pueden usar cadenas de caracteres para
resolver el problema.)
function [d]=Descomponer(N)
El argumento de salida d será un vector fila, de longitud igual al número de dígitos de N conte-
niendo su descomposición. Por ejemplo, si N=36201, debe ser d=[1,0,2,6,3].
37. (Horner.m) Desde el punto de vista de la economía de operaciones, la manera óptima de evaluar
el polinomio
p(x) = cn+1 + cn x + cn−1 x2 + · · · + c2 xn−1 + c1 xn
es utilizando la siguiente expresión:
p(x) = cn+1 + x cn + x cn−1 + . . . x(c2 + c1 x) . . .
que evalúe el polinomio de coeficientes c en el punto x (x puede ser un vector y en ese caso y
debe ser un vector de las mismas dimensiones).
38. (derpol.m) En MATLAB un polinomio se representa mediante un vector que contiene los valores
de sus coeficientes: el polinomio p(x) = a1 xn + a2 xn−1 + · · · + an x + an+1 se representa mediante
el vector-fila c = (a1 , a2 , . . . , an , an+1 ). Escribir una M-función
que recibe como argumento de entrada un vectores, c conteniendo los coeficientes de un polinomio,
y devuelva, en el vector dc, los coeficientes de su derivada.
39. (sumapol.m) En MATLAB un polinomio se representa mediante un vector que contiene los
valores de sus coeficientes: el polinomio p(x) = a1 xn + a2 xn−1 + · · · + an x + an+1 se representa
mediante el vector-fila c = (a1 , a2 , . . . , an , an+1 ). Escribir una M-función
que recibe como argumentos de entrada dos vectores, c1 y c2 (de longitudes desconocidas que
pueden ser distintas) representando dos polinomios, y devuelva, en el vector c, los coeficientes
del polinomio suma de los dos anteriores.
que devuelva bool = true si el punto P pertenece al triangulo XY Z y bool = false si no.
Indicación: La función de MATLAB polyarea sirve para calcular el área de un polígono plano
(ver Help).
41. (Correos.m) Escribir un script que calcule el coste de enviar un paquete de correos, en base a
la siguiente tabla de precios:
Tipo de envío Peso (0-1 kg) Peso (1-5 kg) Peso (5-25 kg)
Tierra 1.50 e 1.50 e + 1 e adicio- 5.50 e + 0.60 e adi-
nal por cada kg o frac- cionales por cada kg o
ción de kg, a partir de fracción de kg, a partir
1 kg de peso. de 5 kg de peso.
Aire 3e 3 e + 1 e adicional 10.20 e + 1.20 e
por cada kg o fracción adicionales por cada kg
de kg, a partir de 1 kg o fracción de kg, a partir
de peso. de 5 kg de peso.
Urgente 18 e 18 e + 6 e adicio- No se realizan envíos ur-
nales por cada kg o frac- gentes de paquetes con
ción de kg, a partir de un peso superior a 5 kg.
1 kg de peso.
El programa pedirá por pantalla los datos peso y tipo de envío. En el caso de un envío por tierra
o aire de más de 25 kg, el programa escribirá en pantalla el mensaje “No se realizan repartos de
más de 25 kg”. En el caso de envíos urgentes de más de 5 kg, el programa escribirá en pantalla
el mensaje “No se realizan repartos urgentes de más de 5 kg”
que reciba como argumentos dos vectores x e y conteniendo respectivamente las abscisas y las
ordenadas de una sucesión de puntos del plano, y devuelva la longitud del camino que se forma
uniendo dichos puntos en el orden dado.
que reciba como argumentos dos vectores x e y conteniendo respectivamente las abscisas y las
ordenadas de una sucesión de puntos del plano, y devuelva la longitud del camino que se forma
uniendo dichos puntos en orden creciente de abscisas.
Esto ejercicio es como el anterior 42, pero es preciso previamente ordenar los puntos de forma
que sus abscisas estén en orden creciente. Para ello hay que usar la orden sort cuyo uso se debe
consultar en el Help de matlab.
5 B
SE
Resolución de sistemas
UNIVER
VILLA
lineales e an
Ecuaciones Diferenciales
y Análisis Numérico
Sea A una matriz cualquiera y sea B otra matriz con el mismo número de filas que A1 . Entonces, la
“solución” del “sistema” (en realidad un sistema lineal por cada columna de B)
AX = B
X = A \ B
X = mldivide(A, B) (equivalente a lo anterior)
Hay que pensar en él como el operador de “división matricial por la izquierda” (A\B ≡ A−1 B).
De forma similar, si C es una matriz con el mismo número de columnas que A, entonces la solución
del sistema
XA=C
se obtiene en MATLAB mediante el operador“/ ” (“división matricial por la derecha”)
X = C / A
X = mrdivide(C, A) (equivalente a lo anterior)
Estos operadores se aplican incluso si A no es una matriz cuadrada. Lo que se obtiene de estas opera-
ciones es, lógicamente, distinto según sea el caso.
De forma resumida, si A es una matriz cuadrada e y es un vector columna, c = A\y es la solución del
sistema lineal Ac=y, obtenida por diferentes algoritmos, en función de las características de la matriz
A (diagonal, triangular, simétrica, etc.). Si A es singular o mal condicionada, se obtendrá un mensaje
de alerta (warning). Si A es una matriz rectangular, A\y devuelve una solución de mínimos cuadrados
del sistema Ac=y.
1
A menos que A sea un escalar, en cuyo caso A\B es la operación de división elemento a elemento, es decir A./B.
141
5. Resolución de sistemas lineales 142
Para una descripción completa del funcionamiento de estos operadores, así como de la elección del
algoritmo aplicado, véase la documentación de MATLAB correspondiente, tecleando en la ventana de
comandos
doc mldivide
Como comprobación calculamos el residuo que, como es de esperar, no es exactamente nulo, debido
a los errores de redondeo:
>> A*x - b
ans =
1.0e-15 *
-0.4441
0
0
det(A)
El uso de abs(det(A)) < epsilon tampoco es recomendable ya que es muy difícil elegir el epsilon
adecuado (véase el Ejemplo 5.4).
Lo aconsejable para testar la singularidad de una matriz es usar su número de condición cond(A).
Ejemplo 5.4 (Matriz no singular, con det. muy pequeño, bien condicionada)
Se considera la matriz 10 × 10 siguiente, que no es singular, ya que es múltiplo de la identidad,
Vemos que su determinante es muy pequeño. De hecho, el test abs(det(A)) < epsilon etique-
taría esta matriz como singular, a menos que se eligiera epsilon extremadamente pequeño. Sin
embargo, esta matriz no está mal condicionada:
>> cond(A)
ans =
1
>> A = diag([24, 46, 64, 78, 88, 94, 96, 94, 88, 78, 64, 46, 24]);
>> S = diag([-13, -24, -33, -40, -45, -48, -49, -48, -45, -40, -33, -24], 1);
>> A = A + S + rot90(S,2);
24 -13 0 0 0 0 0 0 0 0 0 0 0
-24 46 -24 0 0 0 0 0 0 0 0 0 0
0 -33 64 -33 0 0 0 0 0 0 0 0 0
0 0 -40 78 -40 0 0 0 0 0 0 0 0
0 0 0 -45 88 -45 0 0 0 0 0 0 0
0 0 0 0 -48 94 -48 0 0 0 0 0 0
0 0 0 0 0 -49 96 -49 0 0 0 0 0
0 0 0 0 0 0 -48 94 -48 0 0 0 0
0 0 0 0 0 0 0 -45 88 -45 0 0 0
0 0 0 0 0 0 0 0 -40 78 -40 0 0
0 0 0 0 0 0 0 0 0 -33 64 -33 0
0 0 0 0 0 0 0 0 0 0 -24 46 -24
0 0 0 0 0 0 0 0 0 0 0 -13 24
A es singular (la suma de todas sus filas da el vector nulo). Sin embargo, el cálculo de su determi-
nante con la función det da
>> det(A)
ans =
1.0597e+05
cuando debería dar como resultado cero! Esta (enorme) falta de precisión es debida a los errores
de redondeo que se cometen en la implementación del método LU , que es el que MATLAB usa
para calcular el determinante. De hecho, vemos que el número de condición de A es muy grande:
>> cond(A)
ans =
2.5703e+16
5.3 La factorización LU
La factorización LU de una matriz se calcula con MATLAB con la orden
[L, U, P] = lu(A)
Ax = b
utilizando la factorización anterior sólo hay que plantear el sistema, equivalente al anterior,
Lv = P b
P A x = P b ⇐⇒ L U x = P b ⇐⇒
Ux=v
El primero de estos dos sistemas se resuelve fácilmente mediante un algoritmo de bajada, ya que la
matriz del mismo es triangular inferior. Una vez calculada su solución, v, se calcula x como la solución
del sistema U x = v, que se puede calcular mediante un algoritmo de subida, al ser su matriz triangular
superior.
Ejercicio 5.6 Escribir una M-función function [x] = Bajada(A, b) para calcular la solución
x del sistema Ax = b, siendo A una matriz cuadrada triangular inferior.
Algoritmo de bajada
n = dimensión de A
Para cada i = 1, 2, . . . n,
1 i−1
P
xi = bi − Aij xj
Aii j=1
Fin
Para comprobar el funcionamiento del programa, construir una matriz 20 × 20 (por ejemplo) y un
vector columna b de números generados aleatoriamente (con la función rand o bien con randi) y
luego extraer su parte triangular inferior con la función tril(∗) . (Consultar en el help de MATLAB
la utilización de estas funciones).
(∗)Aunque, en realidad esto no es necesario. Obsérvese que en el programa Bajada que sigue, no
se utiliza la parte superior de la matriz A.
%
% Bajada(A, b) es la solucion del sistema lineal de
% matriz triangular inferior Ax = b
%
%------ tolerancia para el test de singularidad
tol = 1.e-10;
%------ inicializaciones
n = length(b);
x = zeros(n,1);
for i = 1:n
Aii = A(i,i);
if abs(Aii) < tol
warning(' La matriz A es singular ')
end
suma = 0; % en version vectorial, sin usar for,
for j = 1:i-1 % se puede calcular la suma con:
suma = suma + A(i,j)*x(j); %
end % suma = A(i,1:i-1)*x(1:i-1)
x(i) = (b(i) - suma)/Aii;
end
Ejercicio 5.7 Escribir una M-función function [x] = Subida(A, b) para calcular la solución
x del sistema Ax = b, siendo A una matriz cuadrada triangular superior.
Algoritmo de subida
n = dimensión de A
Para cada i = n, . . . , 2, 1
1 Pn
xi = bi − Aij xj
Aii j=i+1
Fin
Para comprobar el funcionamiento del programa, construir una matriz A y un vector b de números
generados aleatoriamente (como en el ejercicio anterior) y luego extraer su parte triangular inferior
con la función triu.
Ejercicio 5.8 Escribir una M-función function [x, res] = LU(A, b) que calcule la solución
x del sistema Ax = b y el residuo res=Ax − b siguiendo los pasos siguientes:
— Calcular la factorización LU mediante la función lu de MATLAB
— Calcular la solución del sistema Lv = P b mediante la M-función Bajada
— Calcular la solución del sistema U x = v mediante la M-función Subida
Utilizar la M-función LU para calcular la solución del sistema
Ñ 1 1 0 3 é Ñ x1 é Ñ 1 é
2 1 −1 2 x2 2
3 −1 −1 2 x3 = 3
−1 2 3 −1 x4 4
Observación: Llamamos LU (con mayúsculas) a esta función para que no se confunda con la
función lu (minúsculas) de MATLAB.
1
Ejercicio 5.9 Las matrices de Hilbert definidas por Hij = son un ejemplo notable de
i+j−1
matrices mal condicionadas, incluso para dimensión pequeña.
Utilizar la M-función LU para resolver el sistema Hx = b, donde H es la matriz de Hilbert de
dimensión n = 15 (por ejemplo) y b es un vector de números generados aleatoriamente. Comprobar
que el residuo es grande. Comprobar que la matriz H está mal condicionada calculando su número
de condición.
La función hilb(n) construye la matriz de Hilbert de dimensión n.
Ax = b
utilizando la factorización anterior sólo hay que plantear el sistema, equivalente al anterior,
Lv = b
t
A x = b ⇐⇒ L L x = b ⇐⇒
Lt x = v
Ejercicio 5.10 Escribir una M-función function [x, res] = CHOL(A, b) que calcule la solu-
ción x del sistema con matriz simétrica definida positiva Ax = b y el residuo res=Ax − b siguiendo
los pasos siguientes:
— Calcular la matriz L de la factorización de Cholesky LLt de A mediante la función chol de
MATLAB
— Calcular la solución del sistema Lv = b mediante la M-función Bajada
— Calcular la solución del sistema Lt x = v mediante la M-función Subida
Para comprobar el funcionamiento del programa, se puede generar alguna de las matrices definida
positivas siguientes:
— M = gallery(’moler’, n)
— T = gallery(’toeppd’, n)
— P = pascal(n)
Observación: Por la misma razón que en el ejercicio 5.8, llamamos CHOL (mayúsculas) a esta
función para que no se confunda con la función chol (minúsculas) de MATLAB.
x = linsolve(A, b)
resuelve el sistema lineal Ax = b por factorización LU si A es una matriz cuadrada y regular, y por
factorización QR si no.
La orden
x = linsolve(A, b, opts)
resuelve el sistema lineal Ax = b por el método que resulte más apropiado, dadas las propiedades de
la matriz A, que se pueden especificar mediante el argumento opcional opts (véase la documentación).
0 0 0 0 0
0 1 0 0 0
0 0 0 0 −1
0 0 0 0 0
0 0 3 0 0
Para crear matrices sparse, MATLAB tiene la función sparse, que se puede utilizar de varias formas.
La orden
AS = sparse(A)
A = full(AS)
Ejemplo 5.12
>> A = full(AS)
A =
0 0 1
0 2 0
-1 0 0
La orden
AS = sparse(i, j, s)
crea la matriz sparse AS cuyos elementos no nulos son AS(i(k), j(k)) = s(k)
Ejemplo 5.13
>> i = [2, 3, 5]
>> j = [1, 5, 3]
>> s = [1, -1, 3]
>> AS = sparse(i, j, s)
AS =
(2,1) 1
(5,3) 3
(3,5) -1
Las operaciones con matrices sparse funcionan, con MATLAB, con la misma sintaxis que para las ma-
trices llenas. En general, las operaciones con matrices llenas producen matrices llenas y las operaciones
con matrices sparse producen matrices sparse. Las operaciones mixtas generalmente producen matrices
sparse, a menos que el resultado sea una matriz con alta densidad de elementos distintos de cero. Por
ejemplo, con la misma matriz AS del ejemplo anterior:
Ejemplo 5.14
>> AS(5,3)
ans =
(1,2) 1
(5,3) -1
(3,5) 3
>> AS'
ans =
(1,2) 1
(5,3) -1
(3,5) 3
>> 10*AS
ans =
(2,1) 10
(5,3) 30
(3,5) -10
>> AS(3,5) = 7
AS =
(2,1) 1
(5,3) 3
(3,5) 7
>> AS(1,4) = pi
AS =
(2,1) 1.0000
(5,3) 3.0000
(1,4) 3.1416
(3,5) 7.0000
6 B
SE
Operaciones de lectura y
UNIVER
VILLA
escritura con ficheros e an
Ecuaciones Diferenciales
y Análisis Numérico
En este capítulo veremos operaciones avanzadas de lectura y escritura, concretamente, que operan
sobre el contenido de un fichero.
Una de las cosas que diferencia unos ficheros de otros es el hecho de que su contenido sea o no legible
con un editor de texto plano, también llamado editor ascii 1 . El editor de MATLAB es un editor ascii
y también lo es el Bloc de Notas de Windows y cualquier editor orientado a la programación. Word no
es un editor ascii.
Desde el punto de vista anterior, existen dos tipos de ficheros:
Ficheros formateados o de texto: son aquéllos cuyo contenido está formado por texto plano.
Esto significa que su contenido se reduce a una secuencia de caracteres (los permitidos en el
código de representación de caracteres que se esté utilizando: ASCII, UTF, Unicode, . . . ) y que
no contienen ningún tipo de información sobre el aspecto del texto, como tipos de letra, tamaño,
grosor, color, etc. Los códigos de programa en cualquier lenguaje de programación se escriben
normalmente en ficheros de texto.
1
ASCII = American Standard Code for Information Interchange (Código Estándar Estadounidense para el Intercambio
de Información).
153
6. Operaciones de lectura y escritura con ficheros 154
save nombredefichero
save('nombredefichero') % también se puede escribir así
Esta orden guarda todo el contenido del workspace (nombres de variables y sus valores) en un fichero
de nombre [Link]. Esto puede ser útil si se necesita abandonar una sesión de trabajo
con MATLAB sin haber terminado la tarea: se puede salvar fácilmente a un fichero el workspace en su
estado actual, cerrar MATLAB y luego recuperar todas las variables en una nueva sesión. Esto último
se hace con la orden
load nombredefichero
load('nombredefichero') % equivalente a la anterior
A = rand(5,5);
B = sin(4*A);
x = 1:15;
frase = 'Solo se que no se nada';
who
save memoria
clear
load memoria
6. Comprueba de nuevo, con who, que vuelves a tener todas tus variables.
load memoria
2. Salva, a otro fichero de nombre [Link] el contenido de dos de las variables que tengas, por
ejemplo:
clear A frase
load memo
Los ficheros .mat creados de esta forma son ficheros binarios. Si se intentaran abrir con un editor
de texto plano cualquiera, sólo se verían un montón de signos extraños. Solo matlab “sabe” cómo
interpretar los datos binarios almacenados dentro.
Sin embargo, la orden save también puede ser usada para guardar datos numéricos en un fichero de
texto, usando la opción -ascii:
En este caso, el fichero no tiene la extensión .mat, ya que la misma está reservada para los ficheros
binarios de MATLAB. Se debe especificar la extensión junto con el nombre del fichero. Además, no se
guardan en el fichero los nombres de las variables, solamente sus valores.
load memoria
2. Salva, a un fichero de texto de nombre [Link] el contenido de algunas de las variables que
tengas, por ejemplo:
El cualificador -ascii indica que se salvaguarden los datos en un fichero de texto. Comprueba
que en tu carpeta de trabajo aparece un fichero de nombre [Link].
3. Puedes usar el editor de matlab para «ver» y modificar el contenido de este fichero: con el
cursor sobre el icono del fichero, pulsa el botón derecho del ratón y elige Open as Text.
Es importante hacer notar que los datos guardados de esta manera en un fichero de texto no pueden
ser recuperados en la memoria con la orden load, a menos que se trate de una sola matriz de datos
numéricos, es decir que cada línea del fichero contenga la misma cantidad de datos. En ese caso hay
que usar la versión A=load(’[Link]’) de la orden para guardar los valores de
la matriz en la variable A.
Vamos a escribir un script que llamaremos, por ejemplo, Plot_datos.m, que «cargue» las variables
del fichero y genere la gráfica correspondiente. .
%
% Plot_Datos: lee datos de un fichero binario y hace una grafica
%
fichero = input('Nombre del fichero a cargar : ','s');
load(fichero)
clf
plot(x,y, 'LineWidth',2)
title(titulo)
shg
fichero contiene la temperatura a lo largo del día, medida cada hora, desde las 0 hasta las 23 horas
(en total 24 datos por día). Cada linea contiene las temperaturas de un día. El fichero contiene
cuatro líneas.
Se desea realizar una gráfica que represente, mediante 4 curvas, la evolución de las temperaturas
a lo largo de cada día.
1. Puesto que el fichero [Link] es de texto, se puede editar. Ábrelo para ver su
estructura. Cierra el fichero.
2. Comienza por «cargar» en memoria los datos del fichero [Link], que es un fichero
de texto. Puesto que todas las filas del fichero contienen el mismo número de datos, esto se
puede hacer con la orden load:
Temp = load('[Link]');
que guarda en la variable Temp una matriz con 4 filas y 24 columnas conteniendo los datos
del fichero. Compruébalo. Cada fila de Temp contiene las temperaturas de un día.
hh = 0:23;
4. Se quiere construir una curva con cada fila de la matriz Temp. Esto se puede hacer con la
orden
plot(hh,Temp(1,:),hh,Temp(2,:),hh,Temp(3,:),hh,Temp(4,:),)
plot(hh,Temp)
Consulta en el Help de MATLAB este uso de la función plot en el que el segundo argumento
es una matriz.
5. Añade un título, una etiqueta en el eje OX y una leyenda, para la correcta comprensión de
esta gráfica.
Con frecuencia es necesario utilizar ficheros procedentes de otros usuarios o creados por otros programas
y cuya estructura y formato no se puede elegir o, por el contrario, es preciso generar un fichero con
una estructura determinada.
Para todo ello es necesario recurrir a órdenes de lectura/escritura de «bajo nivel». Estas órdenes son
más complicadas de usar, ya que requieren más conocimiento del usuario pero, a cambio, permiten
elegir todas las caracterísitcas de la operación que se realiza.
Los pasos en el trabajo «de bajo nivel» con ficheros son:
Abrir el fichero e indicar qué tipo de operaciones se van a efectuar sobre él. Esto significa que el
fichero queda listo para poder ser utilizado por nuestro programa.
Leer datos que ya están en el fichero, o bien Escribir datos en el fichero o Añadir datos a los que
ya hay.
Cerrar el fichero.
Sólo hablamos aquí de operaciones con ficheros formateados. No obstante, en algunas ocasiones, como
por ejemplo cuando hay que almacenar grandes cantidades de datos, puede ser conveniente utilizar
ficheros binarios, ya que ocupan menos espacio y las operaciones de lectura y escritura sobre ellos son
más rápidas.
fid=fopen('NombreFichero');
fid=fopen('NombreFichero', 'permiso');
permiso es un código que se utiliza para indicar si se va a utilizar el fichero para leer o para escribir.
Posibles valores para el argumento permiso son:
r indica que el fichero se va a utilizar sólo para lectura. Es el valor por defecto.
(En Windows conviene poner rt).
w indica que el fichero se va a utilizar para escribir en él. En este caso, si existe un fichero con
el nombre indicado, se abre y se sobre-escribe (el contenido previo, si lo había, se pierde).
Si no existe un fichero con el nombre indicado, se crea.
(En Windows conviene poner wt)
a indica que el fichero se va a utilizar para añadir datos a los que ya haya en él.
(En Windows conviene poner at)
Para otros posibles valores del argumento permiso, ver en el Help de MATLAB la descripción
de la función fopen.
Una vez que se ha terminado de utilizar el fichero, hay que cerrarlo, mediante la orden
fclose(fid)
donde fid es el número de identificación que fopen le dió al fichero. Si se tienen varios ficheros abiertos
a la vez se pueden cerrar todos con la orden
fclose('all')
A = fscanf(fid,'formato', dimension);
lee los datos de un fichero de acuerdo con el formato especificado, hasta completar la dimension dada,
o hasta que se encuentre un dato que no se puede leer con dicho formato, o hasta que se termine el
fichero.
formato es una descripción del formato que hay que usar para leer los datos. Se utilizan los mismos
códigos que con la orden fprint.
A es un vector columna en el que se almacenan los datos leídos, en el orden en que están escritos.
dimension es la dimensión de la variable a la que se asignarán los datos (A). Puede ser:
(i) un número entero n; (ii) Inf; (iii) un vector [n, m] o [n, Inf].
Los datos se leen del fichero en el orden natural en que están escritos en el fichero (por filas) y
se almacenan en la matriz A en orden de columnas.
Ejemplo 6.6
El fichero [Link] contine dos columnas de números reales, de longitud desconocida. Se desea
leer dichos datos y guardarlos en una matriz con dos columnas y el número necesario de filas, como
están en el fichero.
fid = fopen('[Link]','r');
En realidad, la opción ’r’ es la opción por defecto, luego no sería necesario indicarlo.
El formato ’ %f’ indica que se van a leer datos numéricos reales, esto es con una parte entera,
un punto decimal y una parte fraccionaria. El número de dígitos de la parte entera y de la
parte fraccionaria vendrá determinado por lo que se lee.
3. Tras la orden anterior, la variable mat contiene un vector columna. Compruébalo. Para que
los datos formen una matriz con 2 columnas como en el fichero hay que reorganizarlos con la
función reshape (ver su descripción más abajo):
mat = reshape(mat,2,[]);
Observa el resultado. Ahora mat es una matriz con dos filas. Lo que deseamos obtener es la
transpuesta de esta matriz:
mat = mat';
4. Cierra el fichero
fclose(fid);
La orden reshape utilizada en el Ejemplo anterior sirve para organizar del datos de una matriz alma-
cenada en memoria como una matriz con otras dimensiones:
asigna los datos de la matriz A a la matriz B de con nf filas y nc columnas. Es necesario que el producto
de nf×nc sea igual al número total de elementos de la matriz A.
También se puede dejar sin especificar una de las dos dimensiones:
Naturalmente, en este caso es necesario que el número total de elementos de A sea un múltiplo de nf
en el primer caso o de nc en el segundo.
Observación: en el Ejemplo 6.6 se podrían haber leído los datos del fichero usando la orden load, ya
que todas las filas están formadas por el mísmo número de datos. Sin embargo la orden fscanf sirve
para otros casos que no se podrían leer con load, como en el Ejemplo siguiente.
Ejemplo 6.7
Se dispone de un fichero [Link] que contiene datos sobre contaminación de las pro-
vincias andaluzas. Cada linea del fichero contiene el nombre de la provincia y tres datos numéricos.
Se quieren recuperar los datos numéricos del fichero en una matriz con tres columnas, es decir,
ignorando el texto de cada línea.
1. Intenta «cargar» el fichero con la orden load. Recibirás un mensaje de error, ya que, como
hemos dicho, load no puede leer ficheros de texto con datos de distintos tipos:
2. Sin embargo, el fichero se puede leer con la orden fscanf. Abre el fichero para lectura:
fid = fopen('[Link]','r');
que indica que se va a leer el contenido del fichero hasta que se termine o bien hasta que los
datos a leer no coincidan con el formato especificado. El formato especifica leer 1 dato real.
Si después de éste siguen quedando datos en el fichero, se «re-visita» el el formato, es decir,
es como si pusiéramos %f %f %f %f %f %f ... hasta que se acabe el fichero.
Verás que el resultado es una matriz vacía. ¿Cual es el error? Que se han encontrado datos
que no se pueden leer con el formato ’ %f’.
4. Escribe ahora la orden:
frewind(fid);
Esta orden «rebobina» el fichero, es decir, coloca el apuntador de lectura al principio del
mismo, ya que, tras la orden anterior puede no estar ahí.
5. Escribe ahora la orden
Esta orden indica que se han de leer del fichero 1 dato tipo string (cadena de caracteres), luego
3 reales y repetir. Ahora obtendrás una matriz con una sola columna. Observa sus valores.
¿Es lo que esperabas? Intenta comprender lo que has obtenido. Quizá te ayude consultar la
Tabla de los caracteres imprimibles del Código ASCII, al final de este Capítulo.
6. Rebobina de nuevo el fichero y escribe la orden:
Esta orden indica que se ignore ( %*s) un dato de tipo string, que luego se lean 3 reales y
luego repetir. Observa el resultado.
7. Para dar a la matriz mat la estructura deseada hay que usar la orden
mat = reshape(mat,3,[])';
que indica que se re–interprete mat como una matriz con 3 filas y el número de columnas que
salgan y luego se transpone.
formato es una descripción del formato que hay que usar para escribir los datos.
lista_de_variables son las variables o expresiones cuyos valores hay que escribir en el fichero.
Ejercicio 6.8 El siglo pasado era costumbre escribir los números de los años con sólo dos dígitos:
76, 77, 78, etc.
Se dispone de un fichero temp_anuales.dat en donde se almacena la temperatura máxima media
de cada mes del año, para varios años del siglo XX. Los años vienen expresados con dos dígitos y
las temperaturas vienen expresadas con una cifra decimal.
Se necesita generar a partir de este fichero otro que contenga los mismos datos, pero con los años
expresados con 4 dígitos, para poder añadir los de este siglo. Las temperaturas deben ser escritas,
también, con un decimal.
Genera también una gráfica que muestre la evolución de las temperaturas a lo largo del año (una
curva para cada año). Dibuja cada curva de un color y añade a cada curva una leyenda con el
número del año al que corresponde, como en la gráfica siguiente:
40
1995
1996
1997
35
1998
1999
30
25
20
15
10
0
1 2 3 4 5 6 7 8 9 10 11 12
1. Puesto que el contenido del fichero temp_anuales.dat es una matriz de números, se puede
leer con la orden load:
datos = load('temp_anuales.dat');
2. Los números de los años forman la primera columna de la matriz datos, así que un vector
conteniendo los años con la numeración de 4 cifras se obtiene con:
anyos = datos(:,1)+1900;
Observación: para escribir los datos en un nuevo fichero no podemos usar la orden save,
ya que los escribiría con formato exponencial (p.e. 9.5000000e+01) y queremos escribir los
años sin decimales y los datos con solo un decimal.
Hay que usar órdenes de «bajo nivel», concretamente la orden fprintf.
3. Abrimos un nuevo fichero para escritura, en el que grabaremos los nuevos datos.
fid = fopen('temp_anuales_new.dat','w');
4. Recuperamos las dimensiones de la matriz datos, para saber cuántas filas tiene:
5. Vamos a realizar un bucle con un índice k variando desde 1 hasta nfils. En cada iteración
del bucle escribiremos una fila del nuevo fichero, mediante las órdenes:
6. Cerramos el fichero
fclose(fid);
plot(1:12,datos(:,2:13),'LineWidth',1.1);
legend(num2str(anyos));
axis([1,12,0,40])
Ejercicio 6.9 Modifica alguno de los scripts de cálculo de series (por ejemplo SumaExp.m) para
que la tabla de las iteraciones se escriba en un fichero en lugar de en la pantalla.
7 B
SE
Interpolación y ajuste de
UNIVER
VILLA
datos e an
Ecuaciones Diferenciales
y Análisis Numérico
7.1 Introducción
En física y otras ciencias con frecuencia es necesario trabajar con conjuntos discretos de valores de
alguna magnitud que depende de otra variable. Pueden proceder de muestreos, de experimentos o
incluso de cálculos numéricos previos.
En ocasiones, para utilizar estos valores en cálculos posteriores es preciso «darles forma» de función,
es decir: es preciso disponer de una función dada por una expresión matemática que «coincida» con
dichos valores.
Existen básicamente dos enfoques para conseguir esto:
165
7. Interpolación y ajuste de datos 166
Interpolación lineal. Por dos puntos dados del plano pasa una sola línea recta.
Más concretamente, dados dos puntos en el plano (x1 , y1 ) y (x2 , y2 ), con x1 6= x2 se trata de
determinar una función polinómica de grado 1
y = ax + b
La solución de este sistema lineal de dos ecuaciones con dos incógnitas proporciona los valores
adecuados de los coeficientes a y b.
(x 2 , y 2 )
(x 1 , y 1 )
x
Interpolación cuadrática. En general, tres puntos del plano determinan una única parábola (poli-
nomio de grado 2).
Dados (x1 , y1 ), (x2 , y2 ) y (x3 , y3 ) con x1 , x2 y x3 distintos dos a dos, se trata de determinar una
función de la forma
y = ax2 + bx + c
(x 1 , y 1 )
(x 3 , y 3 )
(x 2 , y 2 )
p(x) = c1 xN −1 + c2 xN −2 + · · · + cN −1 x + cN
−1 −2
y1 = c1 xN1 + c2 xN
1 + · · · + cN −1 x1 + cN
y2 = c1 xN −1 + c2 xN −2 + · · · + cN −1 x + cN
2 2 2
...
N −1 −2
+ c2 xN
y
N = c1 xN N + · · · + cN −1 xN + cN
−1 −2
xN
1 xN
1 ··· x1 1 c1 y1
N −1 N −2
x2 x2 ··· x2 1 c2 y2
.. .. .. .. .. = ..
. . . .
··· . .
−1 −2
xN
N xN
N · · · xN 1 cN yN
La matriz de este sistema es conocida como matriz de Vandermonde.1 La solución (única) de este
sistema porporciona los coeficientes del polinomio (único) que pasa por los N puntos dados.
1
Alexandre-Théophile Vandermonde (1735-1796) fue un músico, químico y matemático francés. Aunque en realidad
fue Henri-Léon Lebesgue quien en 1940 dió, en su honor, este nombre a las matrices.
2
Joseph Louis Lagrange (1736–1813), fue un matemático, físico y astrónomo italiano nacido en Turín, aunque la mayor
parte de su vida transcurrió entre Prusia y Francia.
(x 2 , y 2 )
(x 3 , y 3 )
(x 1 , y 1 ) (x N , y N )
Interpolación con MATLAB. Dados N puntos (xk , yk ), k = 1, 2, . . . N , con todos los xk distintos,
la instrucción
c=polyfit(x,y,N-1)
calcula los coeficientes del polinomio de grado N-1 que pasa por los N puntos.
x, y son dos vectores de longitud N conteniendo respectivamente las abscisas y las ordenadas
de los puntos.
El vector c devuelto por polyfit contiene los coeficientes del polinomio de interpolación:
p(x) = c1 xN −1 + c2 xN −2 + · · · + cN −1 x + cN .
Ejercicio 7.1 Calcular el polinomio de interpolación de grado 2 que pasa por los puntos
Representar su gráfica en el intervalo [0, 7], señalando con marcadores los puntos interpolados y
dibujando también los ejes coordenados.
1. Crea dos vectores x e y conteniendo respectivamente las abscisas y las ordenadas de los
puntos:
x = [1,2,3];
y = [-3,1,3];
2. Calcula con la orden polyfit los coeficientes de la parábola que pasa por estos puntos:
c = polyfit(x,y,2);
3. Para dibujar la gráfica de esta función comenzamos por crear un vector z de puntos en el
intervalo [0, 9]:
z = linspace(0,7);
4. Necesitamos ahora calcular los valores del polinomio p(x) en todos estos puntos. Para ello
usamos la función polyval:
p = polyval(c,z);
5. Dibuja la parábola:
plot(z,p)
6. Añade ahora los ejes de coordenadas y los marcadores de los puntos del soporte de interpo-
lación:
hold on
plot([-1,7],[0,0], 'k','LinWidth',1.1)
plot([0,0],[-10,6],'k','LineWidth',1.1)
plot(x,y,'r.','MarkerSize',15)
axis([-2,8,-11,7])
hold off
−2
−4
−6
−8
−10
−2 −1 0 1 2 3 4 5 6 7 8
Ejercicio 7.2 La temperatura del aire cerca de la tierra depende de la concentración K del ácido
carbónico (H2 CO3 ) en él. En la tabla siguiente se recoge, para diferentes latitudes L sobre la tierra y
para el valor de K = 0.67, la variación δK de la temperatura con respecto a una cierta temperatura
de referencia:
L 65 35 5 -25 -55
δK -3.1 -3.32 -3.02 -3.2 -3.25
x = [-55,-25,5,35,65];
y = [-3.25,-3.2,-3.02,-3.32,-3.1];
c = polyfit(x,y,4);
z = linspace(-60,70);
p = polyval(c,z);
plot(z,p,'Color',[0,0,1],'LineWidth',1.2)
hold on
plot(x,y,'r.','MarkerSize',15)
4. Usa el polinomio que acabas de obtener para calcular δK para L = 10 e imprime su valor:
L = 10;
delta=polyval(c,10);
fprintf('Para L =%3i delta = %6.3f \n ',L, delta)
−2.8
−2.9
−3
−3.1
−3.2
−3.3
−3.4
−3.5
−60 −40 −20 0 20 40 60 80
Ejercicio 7.3 (Propuesta) Calcula el polinomio de grado 10 que interpola los valores:
x 0 2 3 5 6 8 9 11 12 14 15
y 10 20 30 -10 10 10 10.5 15 50 60 85
150
100
50
−50
−100
−150
−200
−250
−300
−350
−400
0 2 4 6 8 10 12 14 16
Esta práctica pretende mostrar que el procedimiento de interpolación global es, en general
inestable, ya que los polinomios tienden a hacerse oscilantes al aumentar su grado y eso puede
producir grandes desviaciones sobre los datos.
Ejercicio 7.4 (Propuesta) Calcula el polinomio de grado 5 que interpola los valores:
s1=interp1(x,y,z)
calcula el valor en el/los punto(s) z de la interpolante lineal a trozos que pasa por los puntos (x,y).
x y son dos vectores de la misma dimensión que contienen las abscisas y las ordenadas de los
puntos dados.
z son las abscisas de los puntos a interpolar, es decir, los puntos en los cuales queremos evaluar la
interpolante. Puede ser una matriz.
100
80
60
40
20
−20
0 2 4 6 8 10 12 14 16
s1 son los valores calculados, es decir, los valores de la función interpolante en z. s1 tendrá las
mismas dimensiones que z
Calcula, a partir del polinomio de interpolación lineal a trozos, el valor interpolado para x = 1
y compáralo con el obtenido mediante un polinomio de interpolación global de grado 10 (el de la
Práctica 7.3)
Dibuja juntas las gráficas del polinomio de interpolación lineal a trozos y del polinomio de inter-
polación de grado 10.
1. Crea dos variables x e y con las abscisas y las ordenadas que vas a interpolar:
2. Calcula a partir del polinomio de interpolación lineal a trozos el valor interpolado para x = 1:
s1 = interp1(x,y,1)
c = polyfit(x,y,10);
s = polyval(c,1) % = -245.7853
3. Dibuja el polinomio de interpolación lineal a trozos y el polinomio de grado 10 junto con los
puntos a interpolar:
t = linspace(0,15);
p = polyval(c,t);
plot(x,y,x,y,'o',t,p)
axis([-1,16,-400,150])
Observa la gráfica y compara los valores los valores s1 y c1 ¿Qué puedes decir? ¿Cuál de los
dos valores da una aproximación a priori más acertada?
150
100
50
−50
−100
−150
−200
−250
−300
−350
−400
0 2 4 6 8 10 12 14 16
s = spline(x,y,z)
Los argumentos de entrada x, y y z tienen el mismo significado que en el comando interp1 (ver la
Sección 7.3), el argumento de salida almacena en el vector s los valores de la función spline en los
puntos arbitrarios guardados en el vector z.
Ejercicio 7.7 Calcula el spline cúbico que interpola los datos de la Práctica 7.3. Dibuja en la
misma ventana el spline calculado junto con el polinomo de interpolación lineal a trozos.
1. Crea dos variables x e y con las abscisas y las ordenadas que vas a interpolar:
z=linspace(0,15);
s=spline(x,y,z);
3. Dibuja el spline cúbico y la interpolante lineal a trozos junto con los datos:
plot(x,y,'or',z,s,x,y,'LineWidth',2)
90
80
70
60
50
40
30
20
10
−10
−2 0 2 4 6 8 10 12 14 16 18
Ejercicio 7.8 (Propuesta) El fichero [Link] contiene una matriz con dos columnas, que
corresponden a las abscisas y las ordenadas de una serie de datos.
Hay que leer los datos del fichero, y calcular y dibujar juntos el polinomio de interpolación global
y el spline cúbico que interpolan dichos valores, en un intervalo que contenga todos los puntos del
soporte.
−2.5
Interpolante global
Spline cubico
−3
−3.5
−4
−60 −40 −20 0 20 40 60 80
Existen distintos tipos de splines que se diferencian en la forma que toman en los extremos. Para más
información consulta en el Help de MATLAB.
Ejercicio 7.9 (Para ampliar conocimientos) Cuando se calcula un spline cúbico con la función
spline es posible cambiar la forma en que éste se comporta en los extremos. Para ello hay que
añadir al vector y dos valores extra, uno al principio y otro al final. Estos valores sirven para imponer
el valor de la pendiente del spline en el primer punto y en el último. El spline así construido se
denomina sujeto.
Naturalmente, todos los procedimientos de interpolación antes explicados permiten aproximar fun-
ciones dadas: basta con interpolar un soporte de puntos construido con los valores exactos de una
función.
En esta práctica se trata de calcular y dibujar una aproximación de la función sen(x) en el intervalo
[0, 10] mediante la interpolación con dos tipos distintos de spline cúbico y comparar estos resultados
con la propia función. Hay por lo tanto que dibujar tres curvas en [0, 10]:
1. La curva y = sen(x).
0.5
−0.5
−1
0 2 4 6 8 10
Cuando lo que se minimiza es la suma de las distancias de los puntos a la curva (medidas como se
muestra en la figura) hablamos de ajuste por mínimos cuadrados. La descripción detallada de este
método se escapa de los objetivos de esta práctica. Veremos solamente cómo se puede hacer esto con
MATLAB en algunos casos sencillos.
c=polyfit(x,y,m)
x y son dos vectores de la misma dimensión que contienen respectivamente las abscisas y las orde-
nadas de los N puntos.
Ejercicio 7.10 Calcula y dibuja los polinomios de ajuste de grado 1, 2, 3 y 6 para los siguientes
datos:
(0.9, 0.9) (1.5, 1.5) (3, 2.5) (4, 5.1) (6, 4.5) (8, 4.9) (9.5, 6.3)
Una vez calculados, escribe las expresiones analíticas de los polinomios que has obtenido.
p1=polyfit(x,y,1);
subplot(2,2,1)
plot(x,y,'o',[0,10],[polyval(p1,0),polyval(p1,10)])
title('Ajuste por una recta')
p2=polyfit(x,y,2);
xp=linspace(0,10);
subplot(2,2,2)
plot(x,y,'o',xp,polyval(p2,xp));
title('Ajuste por una parabola')
6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 2 4 6 8 10 0 2 4 6 8 10
6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 2 4 6 8 10 0 2 4 6 8 10
Al realizar esta práctica dibujando cada curva en un cuadro distinto, como se hace aquí, debes
tener cuidado de fijar en cada subplot los mismos ejes con el comando axis, para poder comparar
bien. Si no lo haces los resultados te pueden resultar confusos.
x 0 2 3 5 6 8 9 11 12 14 15
y 10 20 30 -10 10 10 10.5 15 50 60 85
100
80
60
40
20
−20
−5 0 5 10 15 20
Ejercicio 7.12 (Propuesto) En el fichero [Link] está recogido el precio de una determi-
nada acción de la Bolsa española a lo largo de 88 días, tomado al cierre de cada día.
Queremos ajustar estos datos por una función que nos permita predecir el precio de la acción para
un corto intervalo de tiempo más allá de la última cotización.
Lee los datos del fichero y represéntalos. Calcula los polinomios de ajuste de grados 1, 2 y 4 y
represéntalos también.
4
x 10
2.9
Interp. lineal a trozos
Recta de regresion
2.85 Ajuste parabola
Ajuste polin. grado 4
2.8
2.75
2.7
2.65
2.6
2.55
2.5
2.45
2.4
0 10 20 30 40 50 60 70 80 90 100
se ve que se pueden encontrar ln(b) y m ajustando los datos ln(x), ln(y) mediante un polinomio
de grado 1.
de donde se pueden encontrar ln(b) y m ajustando los datos x y ln(y) mediante una recta.
Función logarítmica: y = m ln(x) + b Está claro que hay ajustar los datos ln(x) e y mediante
una recta.
1
Función y = La anterior relación se puede escribir también:
mx + b
1 1
y= ⇔ mx + b =
mx + b y
1
lo que muestra que x y se relacionan linealmente. Por lo tanto se pueden calcular m y b
y
1
ajustando x y mediante una recta de regresión.
y
mediante una función potencial y = bxm y dibujar la función obtenida así como los datos con
marcadores.
podemos calcular la recta de regresión y = αx+β para los datos (ln(x), ln(y)) y luego tomar b = eβ
y m = α.
c = polyfit(log(x),log(y),1);
Obtendrás c = [0.6377, -1.4310], lo que significa, por lo dicho antes, que la función que
buscamos es:
y = e−1.4310 · x0.6377 = 0.2391 · x0.6377
plot(x,y,'rx','LineWidth',1.5)
hold on
xs = linspace(1.5,5.5);
ys = exp(c(2)) * xs.^(c(1));
plot(xs,ys,'LineWidth',1.5)
hold off
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
Ejercicio 7.14 (Propuesta) Ajustar los datos contenidos en el fichero [Link] mediante
una función exponencial y = b emx y dibujar la función obtenida así como los datos.
0.06
0.05
0.04
0.03
0.02
0.01
0
2 2.5 3 3.5 4 4.5 5
Ejercicio 7.15 Determinar una función, de las indicadas antes, que se ajuste lo mejor posible (a
simple vista) a los datos de la siguiente tabla:
x=0:0.5:5;
y=[6, 4.83, 3.7, 3.15, 2.41, 1.83, 1.49, 1.21, 0.96, 0.73, 0.64];
plot(x,y,'rx','LineWidth',1.2)
Observa, mirando la gráfica, que una función lineal no proporcionaría el mejor ajuste porque
los puntos claramente no siguen una línea recta. De las restantes funciones, la logarítmica
también se excluye, ya que el primer punto es x = 0. Lo mismo sucede con la función potencia
ya que ésta se anula en x = 0 y nuestros datos no.
Vamos a realizar el ajuste con las funciones exponencial y inversa. A la vista de la gráfica
que vamos a obtener veremos cuál de ellas se ajusta mejor a los datos.
2. Comienza por calcular la función exponencial y = bemx . Utiliza la función polyfit con x y
log(y) para calcular los coeficientes b y m:
c=polyfit(x,log(y),1);
m1=c(1);
b1=exp(c(2));
3. Determina ahora la función inversa y = 1/(mx + b), utilizando de nuevo polyfit con x y
1./y:
p=polyfit(x,1./y,1);
m2=p(1);
b2=p(2);
t1=linspace(0,5);
t2=linspace(0.1,5)
y1=b1*exp(m1*t1);
y2=1./(m2*t2+b2);
axis([-1,6,0,7])
plot(t1,y1,t2,y2)
legend('Datos','Ajuste por exponencial','Ajuste por hipérbola')
7
Datos
Ajuste exponencial
Ajuste hiperbola
6
0
−1 0 1 2 3 4 5 6
A la vista de la gráfica obtenida, ¿cuál de las dos funciones se ajusta mejor a los datos?
8 B
SE
Resolución de ecuaciones
UNIVER
VILLA
no lineales e an
Ecuaciones Diferenciales
y Análisis Numérico
8.1 Introducción
f (x) = 0. (8.1)
A las soluciones de esta ecuación, también se les suele llamar ceros de la función f , por ser los puntos
en los que f = 0.
Desde el punto de vista geométrico, las soluciones de la ecuación 8.1 son los puntos en los que la gráfica
de la función f «toca» al eje OX, aunque no necesariamente lo atraviesa (véanse las figuras)
y y
α x α x
Figura 8.1: La gráfica «toca» al eje de abscisas en el punto α, lo que significa que α es un cero
de la función f , aunque no en ambos casos la función cambia de signo en α.
En las aplicaciones surgen frecuentemente problemas que conducen a ecuaciones del tipo (8.1) cuyas
soluciones no se pueden calcular explícitamente.
Por ejemplo, la ecuación
2
ex + e−x =
cos x
que aparece, por ejemplo, cuando se quieren determinar las frecuencias de las oscilaciones transversales
de una viga con extremos empotrados y sometida a un golpe. Las soluciones de esta ecuación no se
pueden calcular por métodos analíticos.
Incluso para ecuaciones tan aparentemente sencillas como las polinómicas
es bien conocido que para n ≥ 5, no existe una fórmula explícita de sus soluciones.
183
8. Resolución de ecuaciones no lineales 184
s = roots(p)
que calcula todas las raíces (reales y complejas) de un polinomio y puede ser utilizado para calcular
las soluciones reales de la ecuación 8.2.
Ejemplo 8.1
Calcular las soluciones de la ecuación polinómica:
x3 − 9x2 − x + 5 = 0.
Recuerda que hay que incluir los coeficientes nulos, si los hay.
roots(p)
obtendras las raíces: 9.0494, -0.7685 y 0.7190 que, puesto que son todas reales, son las
soluciones de la ecuación.
Recuerda también, siempre, que estamos realizando cálculos numéricos con el ordenador y que, por
consiguientes, todo lo que calculemos es aproximado.
Ejemplo 8.2
Calcular las soluciones de la ecuación polinómica:
2x2 (x + 2) = −1.
1. Comenzamos por escribir la ecuación en forma homogénea (con segundo miembro cero) y
desarrollar el polinomio, para disponer de todos sus coeficientes:
p = [2, 4, 0, 1];
roots(p)
obtendras las raíces: -2.1121, 0.0560 + 0.4833i y 0.0560 - 0.4833i. Por consiguiente,
en el campo real, la única solución de la ecuación es x = −2.1121.
Ejercicio 8.3 Se quiere determinar el volumen V ocupado por un gas a temperatura T y presión
p a partir de la ecuación de estado
ñ Å ã2 ô
N
p+a (V − N b) = kN T,
V
f (V ) = pV 3 − (pbN + kN T )V 2 + aN 2 V − abN 3 = 0
10 5
1
0.8
0.6
0.4
0.2
funcion f(V)
0
-0.2
-0.4
-0.6
-0.8
-1
-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1
V: volumen
cuando la ecuación f (x) = 0 es no lineal (en el caso lineal su resolución es inmediata) son en general
algoritmos iterados.
Esto significa que, a partir de un punto o intervalo iniciales (dependiendo del algoritmo), se construye
una sucesión de aproximaciones (calculadas una a partir de la anterior) cuyo límite es la solución
buscada. Estos algoritmos, cuando convergen, lo hacen a una de las soluciones de la ecuación. Si ésta
tuviera más de una solución, habría que utilizar el algoritmo una vez para cada una, cambiando el
punto inicial.
En la práctica, lógicamente, sólo se efectúan un número finito de iteraciones y el algoritmo se detiene
cuando se verifica algún criterio, previamente establecido.
Para resolver el Problema 8.3, MATLAB dispone de la función
donde
funcion es un manejador de la función que define la ecuación, f . Puede ser el nombre de una función
anónima dependiente de una sola variable, o también un manejador de una M-función, en cuyo
caso se escribiría @funcion. Ver los ejemplos a continuación.
Ejemplo 8.4
La ecuación:
x
x + ln =0
3
tiene una solución cerca de x = 1. Calcularla.
1. Comienza por definir una función anónima que evalúe la expresión del primer miembro:
fzero(fun,1)
Tambien podríamos haber usado una M-función para definir la función, en lugar de una función
anónima. Mostramos a continuación cómo se haría.
Salvamos el fichero y lo cerramos. Recuerda que el fichero tiene que llamarse igual que la
M-función.
fzero(@mifuncion,1)
Como se ha dicho antes, los algoritmos de aproximación de raíces de ecuaciones no lineales necesitan
que el punto inicial que se tome esté «cerca» de la solución. ¿Cómo de cerca? Esta pregunta no tiene
una respuesta fácil.
Hay condiciones matemáticas que garantizan la convergencia del algoritmo si el punto inicial se toma
en un entorno adecuado de la solución. Estas condiciones no son inmediatas de comprobar y requieren
un análisis de cada caso.
El ejemplo siguiente pone de manifiesto la necesidad de elegir un punto cercano a la solución.
Ejemplo 8.5
Analizar la influencia del punto inicial en la convergencia de fzero al aproximar la solución de
x
x + ln =0
3
2. Prueba ahora eligiendo un punto inicial más alejado de la solución, por ejemplo:
fzero(fun,15)
El algoritmo utilizado por fzero comienza por «localizar» un intervalo en el que la función cambie de
signo. No funcionará, pues, con ceros en los que no suceda esto, como pasa en el Ejercicio E 8.1.
Ejemplo 8.6
Calcular, si existe, una solución positiva de la ecuación
sen(x) − 2 cos(2x) = 2 − x2
sen(x) − 2 cos(2x) + x2 − 2 = 0
La gráfica te mostrará que esta función tiene dos ceros: uno positivo cerca de x = 1 y otro
negativo cerca de x = −1.
fzero(fun,1)
Como se ha dicho antes, en los casos en que la ecuación tenga varias soluciones, habrá que utilizar
fzero una vez para cada solución que interese calcular.
Ejemplo 8.7
Calcular las soluciones de la ecuación
x √ 1
sen cos x = en [0, π].
2 5
1. El primer paso es escribir la ecuación en forma homogénea (f (x) = 0). Por ejemplo
x √ 1
f (x) = sen cos x − =0
2 5
Comprobarás que hay una solución en el intervalo [0.4, 0.7] y otra en el intervalo [1.5, 2].
Cuando se utiliza la gráfica por ordenador para «localizar» las soluciones de una ecuación es preciso
prestar especial atención a los factores de escala de los ejes en el dibujo y hacer sucesivos dibujos
«acercándose» a determinadas zonas de la gráfica para comprender la situación. Asimismo, siempre
es recomendable hacer previo análisis teórico de la función. Como ejemplo de esta situación, véase
el Ejercicio E 8.3 y el Ejercicio E 8.6 al final del capítulo.
Para calcular (de forma aproximada) en qué punto(s) se anula una tal función se pueden interpolar sus
valores por alguno de los procedimientos estudiados en el Tema anterior, definir una función anónima
con la interpolante y calcular con fzero el cero de esta función.
Ejercicio 8.8 El fichero [Link] contiene, en dos columnas, las abscisas y las ordenadas de un
conjunto de puntos correspondientes a los valores de una función.
Se desea interpolar dichos valores mediante una interpolante lineal a trozos y calcular el valor para
el que dicha interpolante toma el valor 0.56.
mat = load('[Link]');
xs = mat(:,1);
ys = mat(:,2);
2. Si denotamos por g(x) a la interpolante, lo que tenemos que hacer es resolver la ecuación
Necesitamos definir una función anónima que evalúe el primer miembro de la anterior ecuación
3. Ahora utilizamos fzero para calcular un cero de esta función partiendo, por ejemplo, del
punto medio del intervalo ocupado por los xs (se podría elegir cualquier otro).
fzero(fun,mean(x))
Ejercicio 8.9 Repite la práctica anterior, pero interpolando mediante un spline cúbico.
mat = load('[Link]');
xs = mat(:,1);
ys = mat(:,2);
Esta segunda opción es mejor porque con ella los coeficientes se calculan una sola vez.
fzero(fun,mean(x))
cuando la ecuación f (x) = 0 es no lineal (en el caso lineal su resolución es inmediata) son en general
algoritmos iterados. Esto significa que, a partir de un punto o intervalo iniciales (dependiendo del
algoritmo), se construye una sucesión de aproximaciones (calculadas una a partir de la anterior) cuyo
límite es la solución buscada. Estos algoritmos, cuando convergen, lo hacen a una de las soluciones de
la ecuación. Si ésta tuviera más de una solución, habría que utilizar el algoritmo una vez para cada
una, cambiando el punto inicial. En la práctica, lógicamente, sólo se efectúan un número finito de
iteraciones y el algoritmo se detiene cuando se verifica algún criterio, previamente establecido.
1) Subdividir en dos partes el intervalo en que se sabe que la función cambia de signo y tiene una
sola raíz.
2) Averiguar, utilizando el Teorema de Bolzano, en cual de las dos mitades se encuentra la raíz y
descartar la otra mitad del intervalo.
4) Continuar este proceso hasta que el subintervalo elegido tenga una longitud lo suficientemente
pequeña como para que cualquiera de sus puntos sea una aproximación aceptable de la solución.
La elección óptima como aproximación es, entonces, el punto medio del intervalo.
α α x2
a0 x0 b0 a1 x1 b1 a2 b2
En el algoritmo anterior, cabe plantearse la posibilidad de que, en alguna de las iteraciones, sea f (x) =
0. En ese caso, lo lógico sería detener el algoritmo y dar x como solución.
Habría, pues, que introducir, después del paso b1), un test para detectar si f (x) = 0 y, en caso
afirmativo, parar.
Observación: cuando este algoritmo se implementa en un programa, para ser ejecutado en un ordena-
dor, el test mencionado antes no debe consistir en preguntar si f (x) es igual a 0 ya que, debido
a los errores de redondeo, es prácticamente imposible que un número real que resulte de hacer cálculos
sea exactamente igual a cero: será muy pequeño en valor absoluto.
Estas consideraciones dan lugar a la siguiente versión (mejorada) del algoritmo:
Ejercicio 8.10 (Método de bisección) Escribir una M-función que aproxime la solución de
f (x) = 0 en el intervalo [a, b] utilizando el método de bisección.
fa = fun(a);
fb = fun(b);
if ( fa*fb > 0 )
error('La funcion debe cambiar de signo en el intervalo');
end
e = (b-a)*0.5;
x = (a+b)*0.5;
A las soluciones de x = g(x) se les llama puntos fijos de la función g. Geométricamente hallar un
punto fijo de g es determinar la abscisa del punto de corte de las gráficas de y = g(x) e y = x en [a, b].
y =x
y =g(x)
a α x2 x1 x0 b
a α b a x1 x3 α x2 x0 b
El método de aproximaciones sucesivas es un método iterativo que consiste en tomar una apro-
ximación inicial x0 ∈ [a, b] y calcular los demás términos de la sucesión {xn }n≥0 mediante la relación
xn+1 = g(xn ).
Es obvio que cualquier ecuación f (x) = 0 puede escribirse en la forma x = g(x) (por ejemplo x =
x + f (x), pero también de muchas otras formas). Bajo ciertas condiciones sobre la función g el método
de aproximaciones sucesivas es convergente con orden de convergencia lineal (orden 1).
Cuando se utiliza el método de aproximaciones sucesivas para calcular una aproximación del punto fijo,
se suelen terminar las iteraciones cuando el valor absoluto de la diferencia entre dos puntos sucesivos
sea menor que una tolerancia pre-establecida, ε:
|xn − xn−1 | ≤ ε.
Por otra parte, para prevenir el caso en que el algoritmo no converja por alguna razón, es preciso
imponer que el programa se detenga tras realizar un número máximo prefijado de iteraciones sin que
se haya obtenido la convergencia.
b) Dados n ≥ 0 y xn .
Ejercicio 8.11 (Método de aproximaciones sucesivas) Escribir una M-función que aproxi-
me una solución de x = g(x) utilizando el método de aproximaciones sucesivas.
x = xcero;
for k = 1:Nmax
x0 = x;
x = fun(x);
if ( abs(x-x0) < epsi )
return
end
end
%
f (x) = 0
consiste en generar una sucesión {xn }n≥0 construida a partir de un valor inicial x0 mediante el método
iterado:
f (xn )
xn+1 = xn −
f 0 (xn )
Obviamente, el método de Newton necesita del conocimiento de la derivada f 0 (x) y que esta no se anule
en ningún término de la sucesión. La interpretación geométrica del método de Newton es la siguiente:
xn+1 es la abscisa del punto de intersección con el eje OX de la tangente a la curva y = f (x) en el
punto (xn , f (xn ))
a b
α x2 x1 x0
Tomando el punto inicial x0 cercano a la solución, bajo ciertas condiciones sobre la función f , el método
de Newton es convergente con orden de convergencia cuadrático (de orden 2).
La forma de detener las iteraciones de este método para obtener una aproximación de la raíz es similar
al método de aproximaciones sucesivas:
|xn−1 − xn | < ε.
Por otra parte, para prevenir el caso en que el algoritmo no converja por alguna razón, es preciso
imponer que el programa se detenga tras realizar un número máximo prefijado de iteraciones sin que
se haya obtenido la convergencia.
b) Dados n ≥ 0 y xn .
f (xn )
b.1) Hacer xn+1 = xn − .
f 0 (xn )
b.2) Si |xn+1 − xn | < tol o bien |f (xn+1 )| < ε, parar y devolver xn+1 como aproximación.
Ejercicio 8.12 Escribir una M-función que aproxime la solución de f (x) = 0 utilizando el método
de Newton.
8.7 Ejercicios
E 8.1 Utiliza el comando fzero para aproximar la solución de la ecuación
x2 = 0
E 8.5 Calcular todos los puntos de corte de la curva y = cos3 (3x) con la circunferencia de centro
(π/2, 0) y radio 1.
E 8.6 Las frecuencias naturales de la vibración de una viga homogénea sujeta por un extremo son las
soluciones de :
f (x) := cos(x) cosh(x) + 1 = 0
Se desea saber qué raíces tiene f en el intervalo [0, 15]. Calcular dichas raíces utilizando la función
MATLAB fzero. Realizar las gráficas necesarias para «localizar» los ceros de la función y luego
calcularlos, uno a uno, partiendo en cada ocasión de un punto suficientemente próximo.
(Solución: x1 = 1.8751, x2 = 4.6941, x3 = 7.8548, x4 = 10.9955, x5 = 14.1372)
E 8.7 (Caida.m) Se desea saber si un tanque con desperdicios tóxicos que se deja caer en el océano
alcanza una determinada velocidad crítica al tocar fondo, ya que, de ser así, podría romperse,
liberando su carga y contaminando el océano. Aplicando las leyes de Newton, se obtiene la
siguiente ecuación para la velocidad de caída, v:
gd v p − b p − b − cv
= − − 2 ln
p c c p−b
que, a partir de los valores de p (peso) y d (altura), calcule la velocidad v con la que el tanque
llega al fondo, considerando los siguientes valores del resto de los parámetros:
g = 32.174 ft/s2 , b = 470.327 lb, c = 0.08.
Determinar si tiene o no peligro de rotura un tanque de 527.436 lb de peso que se deja caer en
un punto con 300 ft de profundidad, teniendo en cuenta que la velocidad crítica es vc = 40 ft/s.
`2 `1
L= + l2
sen(π − γ − α) sen α
E 8.9 (PesoMaximo.m) Determinar el peso máximo que puede tener el tanque del Ejercicio E 8.7
para que, con los mismos datos, no se rompa al tocar fondo, es decir, que llegue al fondo con
una velocidad v ≤ vc . (Indicación: hallar el cero de la función g : p → v(p) − vc , donde v(p) es la
velocidad con que toca fondo un tanque con peso p.
E 8.13 Modificar el programa Biseccion de forma que imprima una tabla con la evolución de las itera-
ciones, es decir, con una linea por iteración indicando (a) el número k de la iteración; (b) el valor
xk de la aproximación; (c) el valor f (xk ) de la función en ese punto.
E 8.14 Modificar el programa AASS de forma que devuelva, además de la aproximación de la solución, el
número de iteraciones realizadas:
E 8.15 Modificar el programa AASS de forma que imprima una tabla con la evolución de las iteraciones,
es decir, con una linea por iteración indicando (a) el número k de la iteración; (b) el valor xk de
la aproximación; (c) el valor f (xk ) de la función en ese punto.
E 8.16 Modificar el programa AASS de forma que para el test de parada se utilice el error relativo en
lugar del error absoluto.
E 8.17 Modificar el programa AASS de forma que se produzca un mensaje de warning cuando se alcance
el número máximo de iteraciones.
E 8.18 (Newton.m) Escribir una M-función Newton.m que implemente el método de Newton para
aproximar un cero de una función.
E 8.19 Modificar el programa Newton de forma que devuelva, además de la aproximación de la solución,
el número de iteraciones realizadas:
E 8.20 Modificar el programa Newton de forma que imprima una tabla con la evolución de las iteraciones,
es decir, con una linea por iteración indicando (a) el número k de la iteración; (b) el valor xk de
la aproximación; (c) el valor f (xk ) de la función en ese punto.
E 8.21 Modificar el programa Newton de forma que para el test de parada se utilice el error relativo en
lugar del error absoluto:
|xn+1 − xn |
<ε
|xn |
E 8.22 Modificar el programa Newton de forma que se produzca un warning cuando se alcance el número
máximo de iteraciones y otro warning cuando se llegue a un xk en que la derivada f 0 (xk ) sea
muy pequeña en valor absoluto .
9 B
SE
UNIVER
VILLA
Integración numérica
e an
Ecuaciones Diferenciales
y Análisis Numérico
Si se conoce una primitiva, F , de la función f , es bien sabido que el valor de la integral definida se
puede calcular mediante la Regla de Barrow:
Z b
f (x) dx = F (b) − F (a).
a
En la mayoría de los casos, sin embargo, no se puede utilizar esta fórmula, ya que no se conoce dicha
primitiva. Es posible, por ejemplo, que no se conozca la expresión matemática de la función f , sino
sólo sus valores en determinados puntos. Pero también hay funciones (de apariencia sencilla) para las
que se puede demostrar que no tienen ninguna primitiva que pueda escribirse en términos de funciones
2
elementales, como por ejemplo f (x) = e−x .
La integración numérica es una herramienta de las matemáticas que proporciona fórmulas y téc-
nicas para calcular aproximaciones de integrales definidas. Gracias a ella se pueden calcular, aunque
sea de forma aproximada, valores de integrales definidas que no pueden calcularse analíticamente y,
sobre todo, se puede realizar ese cálculo en un ordenador.
integral(fun,a,b);
1
En antiguas versiones de MATLAB, la función básica para el cálculo numérico de integrales definidas era quad, que
ha sido sustituida por integral, que, aunque tiene más funcionalidades, funciona igual que quad en los casos en que esta
última funcionaba.
203
9. Integración numérica 204
f = @(x) x.*sin(4*log(x));
q = integral(f,0.2,3)
Es preciso que la función del integrando esté vectorizada. Debes obtener el valor -0.2837.
También se podría escribir, directamente,
Ejercicio 9.2 Calcular la siguiente integral definida utilizando una M-función para representar el
integrando Z 8
0.8
(x e−x + 0.2) dx
0
Usando la función area, representar gráficamente el área por debajo de la curva de la función a
integrar, las rectas x = 0 y x = 8 y el eje OX. Consultar para ello el help de MATLAB.
1. Escribe, en un fichero de nombre mifun.m, una M-función que evalúe la función a integrar:
q = integral(@mifun, 0,8);
Recuerda que, puesto que el integrando viene definido mediante una M-función, es necesario
poner el símbolo @ delante del nombre al transmitirlo a otra función.
3. La gráfica pedida se puede hacer con las órdenes siguientes. Incluímos un título que utilizamos
para mostrar el valor calculado de la integral.
0.8
La función f (x) = x e−x + 0.2 es positiva en el intervalo [0, 8], por lo tanto el valor de la
integral calculado es el área de la región mostrada en la figura.
x = linspace(0,8);
y = mifun(x);
area(x, y, 'FaceColor', [0.1,0.8,0])
title(['Integral=', num2str(q,'%12.5f')], 'FontSize', 14);
Integral = 3.16043
0.6
0.5
0.4
0.3
0.2
0.1
0
0 1 2 3 4 5 6 7 8
Ejercicio 9.3 (Propuesta) Calcular la integral definida del Ejercicio 9.1, usando una M-función
Ejercicio 9.5 Calcular el área de la región plana delimitada por la curva de ecuación
y = sen(4 ln(x)), el eje OX y las rectas verticales x = 1 y x = 2.
f = @(x) sin(4*log(x));
x1 = linspace(0.95,2.3);
y1=f(x1);
plot(x1,y1)
grid on
A = integral(f,1,2)
Si la función f es negativa en [a, b], entonces el área de la región delimitada por la curva y = f (x), el
eje OX y las rectas verticales x = 1 y x = b es
Z b
A=− f (x) dx
a
Ejercicio 9.6 Calcular el área de la región plana delimitada por la curva de ecuación
1
y= , el eje OX y las rectas verticales x = 3 y x = 4.
x(1 − ln(x))
f = @(x) 1./(x.*(1-log(x)));
x1 = linspace(3, 4);
y1 = f(x1);
plot(x1, y1)
grid on
A = - integral(f,3,4)
Ejercicio 9.7 Calcular el área total de la región plana delimitada por la curva de ecuación
y = sen(4 ln(x + 1)) y el eje OX entre los puntos de abscisas x = 0 y x = 9.
f = @(x) sin(4*log(x+1)); y
x1 = linspace(0,10);
y1 = f(x1);
plot(x1,y1)
grid on
x
2. La observación de la gráfica nos muestra que la función f cambia de signo en dos puntos del
intervalo [0, 9], x1 y x2 , y que es positiva en [0, x1 ], negativa en [x1 , x2 ] y de nuevo positiva
en [x2 , 9].
Tenemos que calcular las tres áreas por separado y para ello lo primero es calcular los valores
x1 y x2 :
A1 = integral(f,0,xcero1); % A1 = 0.7514
A2 = - integral(f,xcero1,xcero2); % A2 = 1.6479
A3 = integral(f,xcero2,9); % A3 = 3.5561
A = A1 + A2 + A3; % A = 5.9554
Ejercicio
9.8 Calcular el área de la región del primer cuadrante delimitada por la curva y =
x
x sen 6 log , el eje OX y las rectas verticales x = 1 y x = 10.
2
Ejercicio 9.10 Calcular los puntos de corte de las curvas siguientes, así como el área de la región
plana encerrada entre ellas
1. Tenemos que comenzar por identificar la zona en cuestión y hacernos una idea de la locali-
zación de los puntos de corte. Para ello vamos a dibujar ambas curvas, por ejemplo, entre
x = −5 y x = 5 (podrían ser otros valores):
y=f(x)
y=g(x)
30
y
f = @(x) x.^2-4; 20
yf = f(x); 0
x
yg = g(x);
plot(x,yf,x,yg) -10
grid on -20
-30
Las ordenadas de los puntos de corte son los valores en estos puntos de la función f (o g, ya
que toman el mismo valor en ellos):
Así pues, las dos curvas se cortan en los puntos (−1.4932, −1.7703) y (2.6043, 2.7826).
3. El área de la región que queda encerrada entre las dos curvas es, ya que en [x1 , x2 ] se tiene
g(x) ≥ f (x): Z x2 Z x2
A= (g(x) − f (x)) dx = − (f (x) − g(x)) dx
x1 x1
A = - integral(fg,x1,x2); % A = 20.6396
y=f(x)
y=g(x)
y
Ejercicio 9.11 Calcular los puntos de corte de las curvas siguientes, así como el área de la región
plana encerrada entre ellas
cos x 1
y= = f (x) y y = 0.3 x − = g(x)
x2 + 1 2
Ejercicio 9.12 Calcular el área de la región plana delimitada por las rectas verticales x = 2 y
x = 3 y las gráficas de las funciones
x
f (x) = tg(x) y g(x) = −
3
se cortan en el punto (0, 0) y en otro punto (x̄, ȳ) del primer cuadrante. Se pide calcular las
coordenadas de dicho punto y el área de la región encerrada entre ambas curvas.
se cortan en dos puntos del intervalo [π, 2π]. Se pide calcular las coordenadas de dichos puntos así
como el área de la región encerrada entre ambas curvas.
Ejercicio 9.15 Calcula el área de la o las regiones planas delimitadas por las curvas
La longitud de este arco de curva desde (a, f (a)) hasta (b, f (b)) viene dada por
Z bq
2
L= 1 + f 0 (x) dx (9.1)
a
entonces la longitud del arco de curva desde (x(t0 ), y(t0 )) hasta (x(t1 ), y(t1 )) viene dada por
Z t1 q
2 2
L= f 0 (t) + g 0 (t) dt (9.2)
t0
Ejercicio 9.16 Calcular la longitud del arco de la curva y = sen 4 ln(x + 1) entre x = 0 y x = 9.
4
f 0 (x) =
cos 4 ln(x + 1)
x+1
La longitud del arco indicado se puede calcular con las órdenes:
df = @(x) 4*cos(4*log(x+1))./(x+1);
fint = @(x) sqrt(1+df(x).^2);
L = integral(fint, 0, 9)
Ejercicio 9.17 Calcular la longitud del arco de la curva definida por las ecuaciones paramétricas
®
x = cos(80t) − cos3 (t)
t ∈ [0, π]
y = sen(80t) − sen3 (t)
Denotando f (t) = cos(80t) − cos3 (t) y g(t) = sen(80t) − sen3 (t), se tiene
El volumen del sólido de revolución generado por rotación de la región determinada por la curva
y = f (x), las rectas verticales x = a y x = b y el eje OX en torno a dicho eje viene dado por
Z b
V =π f (x)2 dx
a
El área de la superficie de revolución generada al girar el arco de curva y = f (x), x ∈ [a, b] alrededor
del eje OX viene dada por
Z b q
2
S = 2π f (x) 1 + f 0 (x) dx
a
Ejercicio 9.18 Calcular el volumen del sólido de revolución generado por rotación alrededor del
eje OX de la región plana delimitada por la curva
4
y= , t ∈ [0, 7],
1 + 8e−t
las rectas verticales x = 0 y x = 7 y el eje OX.
f = @(x) 4./(1+8*exp(-x));
V = pi*integral(@(x) f(x).^2, 0, 7); % V = 197.4628
Ejercicio 9.19 Calcular el volumen del sólido de revolución (elipsoide) generado por rotación
alrededor del eje OX de la región plana encerrada dentro de la elipse de ecuación
(x − 3)2 + 3y 2 = 4.
4 − (x − 3).2
y=±
3
El sólido mencionado se obtiene por rotación alrededor del eje OX de la parte positiva
4 − (x − 3).2
y=+
3
Se puede calcular su volumen con las órdenes
Ejercicio 9.20 Calcular el volumen del toroide generado por rotación alrededor del eje OX de la
región plana encerrada dentro de la elipse de ecuación
x2 + 3(y − 3)2 = 4.
Como hemos mencionado antes, muchas veces es preciso trabajar con funciones de las que sólo se cono-
cen sus valores en un conjunto de puntos. Para calcular (de forma aproximada) la integral definida de
una tal función se pueden interpolar sus valores por alguno de los procedimientos vistos anteriormente
y luego integrar dicha interpolante.
La posibilidad más sencilla es hacer una interpolación lineal a trozos. Entonces la integral en todo el
intervalo es la suma de las integrales en cada uno de los subintervalos, que es fácil de hacer, ya que en
cada uno de ellos la función a integrar es un polinomio de grado uno.
q = trapz(x,y);
x,y son dos vectores de la misma dimensión y representan las coordenadas de los puntos.
q es el valor de la integral
Ejercicio 9.21 Calcula la integral definida de una función que viene dada a través del siguiente
conjunto de puntos:
x 0 2 3 5 6 8 9 11 12 14 15
y 10 20 30 -10 10 10 10.5 15 50 60 85
1. Define dos vectores que contienen las abscisas y las ordenadas de los puntos:
q = trapz(x,y)
Si se lleva a cabo otro tipo de interpolación, el procedimiento a seguir es definir una función (por
ejemplo anónima) y luego utilizar integral.
Ejercicio 9.22 Se consideran los siguientes datos correspondientes a los valores de una función:
x 0 2 3 5 6 8 9 11 12 14 15
y 10 20 30 -10 10 10 10.5 15 50 60 85
Se desea interpolar dichos valores mediante un spline cúbico s = s(x), representarlo gráficamente
y calcular la integral definida: Z 14
s(x) dx
10
1. Crea dos variables x e y con las abscisas y las ordenadas que vas a interpolar:
2. Calcula con la orden spline los coeficientes del spline que pasa por estos puntos:
c = spline(x,y);
3. Define una función anónima s = s(x) que evalúa el spline en los puntos x:
s = @(x) ppval(c,x);
a = integral(s,10,14)
plot(x,y,'rx','LineWidth',1.5)
hold on
xs = linspace(min(x),max(x));
ys = s(xs);
plot(xs,ys,'LineWidth',1.5)
hold off
100
80
60
40
20
-20
0 5 10 15
Ejercicio 9.23 Usando los datos del Ejercicio 9.22, calcular el area de la región del cuarto cua-
drante (x ≥ 0, y ≤ 0) delimitada por el spline cúbico que interpola estos datos y el eje OX.
1. Crea dos variables x e y con las abscisas y las ordenadas de los puntos a interpolar:
c = spline(x,y);
s = @(x) ppval(c,x);
xs = linspace(0,15);
ys = s(xs);
plot(xs,ys)
grid on
100
80
60
40
20
-20
0 5 10 15
3. La observación de la gráfica muestra que el spline cúbico cambia de signo en tres puntos del
intervalo [0, 15]: x1 , x2 y x3 y que las regiones que nos interesan son las que quedan entre
x = 0 y x = x1 por un lado, y entre x = x2 y x = x3 por otro. Lo primero es calcular los
valores x1 , x2 y x3 :
x1 = fzero(s,1) % x1 = 1.1240
x2 = fzero(s,4) % x2 = 4.2379
x3 = fzero(s,6) % x3 = 5.5685
A1 = -integral(s,0,x1) % A1 = 4.7719
A2 = -integral(s,x2,x3); % A2 = 8.6921
A = A1 + A2 % A = 13.4640
Ejercicio 9.24 Se considera el siguiente conjunto de puntos correspondientes a los valores de una
función:
x 12 13 14 15 16 17 18 19 20
y -11.43 -7.30 8.86 13.82 -1.76 -16.73 -8.06 13.87 17.24
Se pide calcular el area de la región delimitada por el spline cúbico que intepola estos datos, el eje
OX y las rectas verticales x = 12 y x = 20.
Ejercicio 9.25 El fichero [Link] contiene una matriz con dos columnas, que corresponden
a las abscisas y las ordenadas de una serie de datos.
Se pide leer los datos del fichero, y calcular y dibujar el spline cúbico que interpola dichos valores,
en un intervalo que contenga todos los puntos del soporte. También se desea calcular el área de la
región del primer cuadrante delimitada por el spline cúbico.
1. Puesto que el fichero [Link] es de texto, se puede editar. Ábrelo (como texto) para ver
su estructura. Cierra el fichero.
2. Dado que todas las líneas del fichero tienen el mismo número de datos, se puede leer con la
orden load:
datos = load('[Link]');
que guarda en la variable datos una matriz 50 × 2 cuya primera columna contiene los valores
de las abscisas y la segunda los de las ordenadas. Compruébalo.
x = datos(:,1);
y = datos(:,2);
c = spline(x,y);
f = @(x) ppval(c,x);
xs = linspace(min(x),max(x));
ys = f(xs);
plot(xs,ys)
grid on
0.8
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
7. En la gráfica se observa que el spline cúbico cambia de signo en tres puntos del intervalo
[1, 3]: x1 y x2 . La región del primer cuadrante delimitada por el spline es la que queda entre
x = x1 y x = x2 . Comienza por calcular los valores de x1 y x2 :
x1 = fzero(f,1.5) % x1 = 1.5713
x2 = fzero(f,2.5) % x2 = 2.3566
A = integral(f,x1,x2) % A = 0.5002
x 0 1 2 3 4 5 6 7 8 9 10
y 0 0.28 -1.67 -2.27 1.63 4.95 0.92 -6.32 -5.33 4.72 9.64
Se desea interpolar dichos valores mediante una interpolante lineal a trozos s(x), calcular la integral
definida Z 2π
s(x) dx
π/2
y realizar una representación gráfica de los datos junto con la interpolante lineal a trozos.
Ejercicio 9.27 El fichero [Link] contiene una matriz con dos columnas, que corresponden
a las abscisas y las ordenadas de una serie de datos.
Se pide leer los datos del fichero, y calcular la función logarítmica f (x) = m ln(x) + b que mejor se
ajusta a estos valores, y dibujarla en un intervalo que contenga todos los puntos del soporte.
También se pide calcular el área de la región del plano delimitada por la gráfica de la función f (x),
el eje OX y las rectas verticales x = 1.5 y x = 3.5.
Ejercicio 9.28 Escribir una M-función para calcular la siguiente función definida a trozos:
si x ≥ 0
®
0.5x2
f (x) =
x3 si no
y calcular el área de la región delimitada por la curva y = f (x), el eje OX y las rectas verticales
x = −2 y x = 2.5.
1. En un primer intento, escribiríamos probablemente algo como lo que sigue para calcular la
función f :
2. Para que una M-función como la anterior sea «vectorial» hay que calcular los valores de y
componente a componente, como se muestra a continuación:
Escribe el código anterior en un fichero y sálvalo con el nombre mifun.m. Comprueba, desde
la ventana de comandos, que el programa funciona correctamente, dando distintos valores a
la variable de entrada.
x = linspace(-3, 3);
y = mifun(x);
plot(x, y, 'LineWidth', 1.5)
-2
2.5
A1 = -integral(@mifun,-2,0) % A1 = 4
A2 = integral(@mifun,0,2.5) % A2 = 2.6042
A = A1 + A2; % A = 6.6042
Observación 1
Las órdenes para realizar esta práctica se pueden escribir todas en un solo fichero de nombre,
por ejemplo, practica23.m, de la siguiente manera (para ejecutarlo hay que invocar el programa
principal, es decir, el que lleva el nombre del fichero):
function practica23
%-------------------------------------------------
% Practica numero 23 : funcion "principal"
%
x = linspace(-2.5, 3);
y = mifun(x);
plot(x, y, 'LineWidth', 1.5)
axis([-2.5, 3, -15, 10])
grid on
%
A1 = - integral(@mifun, -2, 0);
A2 = integral(@mifun, 0, 2.5);
A = A1 + A2;
fprintf(' El area calculada es: %12.6f \n',A)
%-------------------------------------------------
% mifun : funcion auxiliar utilizada por practica23
% No es "visible" fuera de este programa
%
Observación 2
En realidad, utilizando adecuadamente las operaciones vectoriales de Matlab y los operadores de
comparación, es posible escribir la función f (x) mediante una sóla expresión y, por lo tanto, como
una función anónima, de la siguiente manera:
En esta expresión se hace uso del hecho de que, cuando se combinan en una expresión aritmética
datos numéricos y datos lógicos, Matlab interpreta estos últimos como ceros y unos (números).
10 B
SE
UNIVER
VILLA
iniciales para ecuaciones
e an
Ecuaciones Diferenciales
diferenciales ordinarias y Análisis Numérico
Una ecuación diferencial es una ecuación en que la incógnita es una función y que, además, involucra
también las derivadas de la función hasta un cierto orden. La incógnita no es el valor de la función en
uno o varios puntos, sino la función en sí misma.
Cuando la incógnita es una función de una sola variable se dice que la ecuación es ordinaria, debido a
que la o las derivadas que aparecen en ella son derivadas ordinarias (por contraposición a las derivadas
parciales de las funciones de varias variables).
Comenzamos con la resolución de ecuaciones diferenciales ordinarias (edo) de primer orden, llamadas
así porque en ellas sólo aparecen las derivadas de primer orden.
Más adelante trataremos la resolución de sistemas de ecuaciones diferenciales ordinarias (sdo) de primer
orden. Con ello podremos, en particular, abordar la resolución de ecuaciones diferenciales de orden
superior a uno, ya que éstas siempre se pueden reducir a la resolución de un sistema de primer orden.
al problema de determinar, de entre todas las soluciones de la ecuación diferencial y 0 = f (t, y), aquélla
que pasa por el punto (t0 , y0 ), es decir, que verifica y(t0 ) = y0 .
Sólo para algunos (pocos) tipos muy especiales de ecuaciones diferenciales es posible encontrar la
expresión de sus soluciones en términos de funciones elementales. En la inmensa mayoría de los casos
prácticos sólo es posible encontrar aproximaciones numéricas de los valores de una solución en algunos
puntos.
Así, una aproximación numérica de la solución del problema de Cauchy
y 0 = f (t, y), t ∈ [t0 , tf ]
(10.2)
y(t0 ) = y0
1
En realidad basta con que f sea continua y localmente Lipschitziana con respecto de y en Ω.
223
10. Problemas de valores iniciales para edo’s 224
yk ≈ y(tk ), k = 0, 1, . . . , n,
yk ≈ y(tk ), k = 0, 1, . . .
Cada una de ellas implementa un método numérico diferente, siendo adecuado usar unas u otras en
función de las dificultades de cada problema en concreto. Para problemas no demasiado «difíciles» será
suficiente con la función ode45 ó bien ode23.
Exponemos aquí la utilización de la función ode45, aunque la utilización de todas ellas es similar, al
menos en lo más básico. Para más detalles, consúltese la documentación.
donde:
odefun es un manejador de la función que evalúa el segundo miembro de la ecuación, f (t, y). Puede
ser el nombre de una función anónima dependiente de dos variables, siendo t la primera de ellas
e y la segunda, o también una M-función, en cuyo caso se escribiría @odefun. Ver los ejemplos a
continuación.
[t0,tf] es el intervalo en el que se quiere resolver la ecuación, i.e. t0= t0 , tf= tf .
1. Comenzamos por definir una función anónima que represente el segundo miembro de la
ecuación f (t, y):
f = @(t,y) 5*y;
ode45(f, [0,1], 1)
100
50
3. Para comparar con la gráfica de la solución exacta, dibujamos en la misma ventana, sin borrar
la anterior:
hold on
t = linspace(0,1);
plot(t, exp(5*t), 'r')
shg
4. Para hacer lo mismo utilizando una M-función en vez de una función anónima, tendríamos
que escribir, en un fichero (el nombre odefun es sólo un ejemplo):
y guardar este texto en un fichero de nombre odefun.m (el mismo nombre que la función).
Este fichero debe estar en la carpeta de trabajo de MATLAB.
Después, para resolver la ecuación, usaríamos ode45 en la forma:
ode45(@odefun, [0,1], 1)
Si se desean conseguir los valores tk e yk de la aproximación numérica para usos posteriores, se debe
usar la función ode45 en la forma:
obteniéndose así los vectores t e y de la misma longitud que proporcionan la aproximación de la solución
del problema:
y(k) ≈ y(t(k)), k = 1, 2, . . . , length(y).
Calculamos la solución en el intervalo [0, 1] y luego calculamos el valor en t = 0.632 utilizando una
interpolación lineal a trozos:
f = @(t,y) t.*exp(t/y);
[t,y] = ode45(f, [0,1], 1);
v = interp1(t, y, 0.632)
Obsérvese que t e y, las variables de salida de ode45, son dos vectores columna de la misma
longitud.
En ocasiones puede ser necesario calcular las aproximaciones yk en unos puntos tk pre-determinados.
Para ello, en lugar de proporcionar a ode45 el intervalo [t0,tf], se le puede proporcionar un vector
con dichos valores, como en el ejemplo siguiente:
Ejercicio 10.3 Calcular (aproximaciones de) los valores de la solución del problema
y 0 = −2y
y(0) = 10
en los puntos: 0, 0.1, 0.2, . . . , 1.9, 2. Comparar (gráficamente) con la solución exacta y = 10 e−2t .
f = @(t,y) -2*y;
t = 0:0.1:2;
[t,y] = ode45(f, t, 10);
plot(t, y, 'LineWidth', 1.1)
hold on
plot(t, 10*exp(-2*t), 'r.')
legend('Sol. aproximada', 'Sol. exacta')
shg
10
Sol. aproximada
Sol. exacta
9
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
1. Comenzamos por visualizar la gráfica de la solución, para obtener, por inspección, una pri-
mera aproximación:
f = @(t,y) 0.5*(10*t-log(y+1));
ode45(f, [0,1], 1)
grid on
shg
3.5
2.5
1.5
0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Vemos así que el valor y = 1.5 se produce para un valor de t en el intervalo [0.5, 0.6].
2. Usaremos ahora la función fzero para calcular la solución de la ecuación y(t) = 1.5 escrita
en forma homogénea, es decir, y(t) − 1.5 = 0, dando como aproximación inicial (por ejemplo)
t = 0.5. Para ello, necesitamos construir una función que interpole los valores que nos devuelva
ode45. Usamos interpolación lineal a trozos:
Gracias a las características vectoriales del lenguaje de MATLAB, estos problemas se resuelven con
las mismas funciones (ode45, ode23, etc.) que los problemas escalares. Basta escribir el problema en
forma vectorial. Para ello, denotamos:
y1 f1 (t, y1 , y2 , . . . , yn ) y10
0
y2 f2 (t, y1 , y2 , . . . , yn ) y2
Y =
, F (t, Y ) =
, Y0 =
,
... ... ...
yn fn (t, y1 , y2 , . . . , yn ) yn0
Con esta notación el problema planteado se escribe:
Y 0 = F (t, Y )
(10.5)
Y (t ) = Y
0 0
y se puede resolver (numéricamente) con la función ode45, de forma similar al caso escalar, con las
adaptaciones oportunas, como en el ejemplo siguiente.
y1 (0) = 0
y2 (0) = 1
y (0) = 1
3
2. Definimos una función anónima para F (t, Y ) y el dato inicial (ambos son vectores columna
con tres componentes):
3. Utilizamos ode45 para calcular y dibujar la solución. Añadimos una leyenda para identificar
correctamente la curva correspondiente a cada componente: y1 (t), y2 (t), y3 (t).
0.5
−0.5
−1
−1.5
0 2 4 6 8 10 12 14 16
4. Vamos ahora a utilizar ode45 recuperando las variables de salida de la aproximación numérica:
Obsérvese que t es un vector columna y que y es una matriz con tantas filas como t y tres
columnas: cada columna es una de las componentes de la solución Y . Es decir:
y(k,i) ≈ yi (tk )
5. Si, por ejemplo, sólo quisiéramos dibujar la gráfica de la segunda componente y2 (t) usaríamos
la orden plot con la segunda columna de y:
plot(t, y(:,2))
legend('y_2(t)')
shg
1
y2(t)
0.8
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
0 2 4 6 8 10 12 14 16
Ejercicio 10.6 Calcular (una aproximación de) la solución del sistema siguiente para t ∈ [0, 15]
y dibujar la solución en el plano de fases:
y10 = 0.5y1 − 0.2y1 y2
y 0 = −0.5y + 0.1y y
2 2 1 2
y1 (0) = 4
y (0) = 4
2
2. Para ver cómo se haría, vamos a escribir la función del segundo miembro como una M-función.
También se podría hacer como una función anónima, como en el ejercicio anterior.
Observaciones:
a) Esta función debe ser guardada en un fichero de nombre odefun.m, que debe estar en la
carpeta de trabajo de MATLAB para poder ser encontrada.
b) El nombre odefun es sólo un ejemplo.
c) La orden dy = zeros(2,1); en la función odefun tiene el objetivo de crear la variable dy
desde un principio con las dimensiones adecuadas: un vector columna con dos componentes.
3. Calculamos la solución:
plot(y(:,1), y(:,2))
xlabel('y_1')
ylabel('y_2')
4.5
3.5
3
y2
2.5
1.5
1
2 3 4 5 6 7 8 9
y1
En esta sección nos planteamos la resolución numérica de problemas de Cauchy para ecuaciones dife-
renciales ordinarias de orden superior a uno, es decir, problemas como:
y (n) = f (t, y, y 0 , . . . , y (n−1) )
y(t0 ) = y0
y 0 (t0 ) = y1 (10.7)
...
y (n−1) (t ) = y
0 n−1
En efecto, se tiene:
z10 = y 0 = z2
z1 (t0 ) = y(t0 ) = y0
z 0 = y 00 = z3
z2 (t0 ) = y 0 (t0 ) = y1
2
...
...
z 0 = y (n) = f (t, y, y 0 , . . . , y (n−1) ) = f (t, z , z , . . . , z
z (t ) = y (n−1) (t ) = y
n 1 2 n−1 ) n 0 0 n−1
Siguiendo los pasos antes indicados, este problema se puede reducir al siguiente:
Z 0 = F (t, Z)
Z(0) = Z0
con Ñ é Ñ é
z2 0
F (t, Z) = , Z0 = .
−7 sen z1 − 0.1 cos t 1
1. Comenzamos por usar ode45 para dibujar la solución del sistema (2 componentes)
0.5
−0.5
−1
−1.5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Ñ é
z1 (t)
En esta gráfica están representadas las dos componentes de la solución Z(t) =
z2 (t)
2. A nosotros sólo nos interesa la primera de las componentes de Z(t), ya que es y(t) = z1 (t). Por
lo tanto usaremos ode45 para recuperar los valores de la solución y dibujaremos únicamente
la primera componente:
0.3
0.2
0.1
−0.1
−0.2
−0.3
−0.4
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Ejercicio 10.8 Calcular (una aproximación de) la solución del sistema siguiente para t ∈ [0, 25]
y dibujar la solución en el plano de fases:
0
y1 = ay1 − by1 y2
y 0 = −cy + dy y
2 2 1 2
y1 (0) = 4
y (0) = 4
2
siendo a, b, c, d > 0.
3. Para la llamada a ode45, damos valores a los parámetros y definimos la siguiente función
anónima:
a = 0.5;
b = 0.2;
c = 0.5;
d = 0.1;
odefun = @(t,y) odeparam(t, y, a, b, c, d)
[t, y] = ode45(odefun, [0, 25], [4; 4]);
params.a = 0.5;
params.b = 0.2;
params.c = 0.5;
params.d = 0.1;
odefun = @(t,y) odeparam(t, y, params)
[t, y] = ode45(odefun, [0, 25], [4; 4]);
11 B
SE
UNIVER
VILLA
sistemas diferenciales
e an
Ecuaciones Diferenciales
ordinarios y Análisis Numérico
Exponemos aquí la forma de utilizar las herramientas de MATLAB para resolver problemas de
contorno relativos a sistemas diferenciales ordinarios de primer orden. En estos problemas se bus-
ca encontrar una solución del sistema, definida en un intervalo [a, b], que verifique unas condiciones
adicionales que involucran los valores de la solución en ambos extremos del intervalo.
Concretamente, queremos calcular una solución del problema
Como es conocido, el problema (11.1) puede no tener solución, tener un número finito de soluciones o
incluso tener infinitas soluciones.
Las mismas herramientas que permiten calcular (cuando la haya) una solución de (11.1) permiten tam-
bién resolver problemas de contorno relativos a ecuaciones diferenciales ordinarias de orden superior, ya
que, mediante un cambio de variable adecuado, estos problemas se pueden reducir al problema (11.1).
y, también, la ecuación diferencial en los extremos y el punto medio de cada subintervalo [xi , xi+1 ].
Esto da lugar a un sistema (no lineal) de ecuaciones que es resuelto por métodos iterados.
Por esta razón, y también porque el problema (11.1) puede tener más de una solución, es necesario
proporcionar ciertas aproximaciones iniciales.
El uso básico de la función bvp4c es el siguiente:
donde:
odefun es un manejador de la función que evalúa el segundo miembro de la ecuación, f (x, y). Como
en otros casos, puede ser el nombre de una función anónima o un manejador de una M-función.
235
11. Problemas de contorno para sistemas diferenciales ordinarios 236
initsol es una estructura que contiene las aproximaciones iniciales. Esta estructura se crea mediante
la función bvpinit, cuya descripción está más adelante.
sol es una estructura1 que contiene la aproximación numérica calculada, así como algunos datos
informativos sobre la resolución. Los campos más importantes de la estructura sol son:
xinit es un vector que describe una malla inicial, i.e., una partición inicial a = x0 < x1 < . . . <
xN = b del intervalo [a, b]. Normalmente bastará con una partición con pocos puntos creada, por
ejemplo, con xinit = linspace(a,b,5). La función bvp4c refinará esta malla inicial, en caso
necesario. Sin embargo, en casos difíciles, puede ser necesario proporcionar una malla que sea
más fina cerca de los puntos en los que se prevea que la solución cambia rápidamente.
yinit es una aproximación inicial de la solución. Puede ser, o bien un vector, o bien una función.
Se usará la primera opción cuando se desee proporcionar una aproximación inicial de la solución
que sea constante. Si se desea proporcionar, como aproximación inicial, una función y = y0 (x),
se suministrará el nombre de una función que evalúe y0 (x).
1. Comenzamos por definir una función anónima que represente el segundo miembro de la
ecuación f (t, y):
2. Definimos también una función anónima que evalúe la función que define las condiciones de
contorno
1
En MATLAB, una estructura es una colección de datos, organizados en campos, que se identifican mediante un
nombre de campo. Ello permite manipular de forma compacta conjuntos de datos de distintos tipos, a diferencia de las
matrices, en que todas las componentes deben ser del mismo tipo.
xinit = linspace(0,1,5);
yinit = [0; pi];
solinit = bvpinit(xinit, yinit);
Ahora, para dibujar las curvas de la solución (y1 (x), y2 (x)), utilizamos la estructura sol:
— sol.x es un vector fila que contiene los puntos a = x0 < x1 < . . . < xN = b del soporte
de la solución numérica
— sol.y es una matriz con dos filas; la primera fila sol.y(1,:) son los valores de y1 (x) en
los puntos sol.x; la segunda fila sol.y(2,:) son los valores de y2 (x) en los puntos sol.x.
En consecuencia, para obtener la gráfica de la solución ponemos:
plot(sol.x, sol.y)
legend('y_1(x)', 'y_2(x)', 'Location', 'Best')
4
−2
−4
−6
−8
y1
y2
−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Observación: la forma en que se ha construido la gráfica anterior sólo utiliza los valores de
la solución numérica en los puntos del soporte sol.x. No hace uso del hecho de que la
aproximación obtenida es un polinomio cúbico a trozos. Para tener esto en cuenta,
habría que usar la función deval, que se explica a continuación.
Para evaluar la solución obtenida con bvp4c en un punto o vector de puntos xp, se debe utilizar la
función deval:
yp = deval(sol, xp);
Ejercicio 11.1 (sigue. . . ) Repetimos el ejercicio, utilizando una partición más fina del intervalo
[0, 1] para realizar la gráfica, de forma que las curvas se vean más «suaves». Utilizamos la función
deval para calcular los valores de la solución aproximada en los puntos de la partición.
function ejercicio1
%
dydx = @(x,y) [ 3*y(1) - 2*y(2); -y(1) + 0.5*y(2)];
condcon = @(ya, yb) [ ya(1); yb(2)-pi];
xinit = linspace(0,1,5);
yinit = [0; pi];
solinit = bvpinit(xinit, yinit);
xx = linspace(0,1);
yy = deval(sol, xx);
plot(xx,yy)
legend('y_1', 'y_2', 'Location', 'Best')
shg
4
−2
−4
−6
−8
y1
y2
−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Comenzamos por escribir el problema en forma de un problema de contorno para un SDO de primer
orden. Hacemos el cambio de variable z1 = y, z2 = y 0 y el problema (11.2) se transforma en:
0
z1 = z2
z20 = 100z1
z1 (0) = 1, z1 (1) = 2
function ejercicio2
%
dydx = @(t,y) [ y(2); 100*y(1)];
condcon = @(ya, yb) [ ya(1)-1; yb(1)-2];
xinit = linspace(0,1,5);
yinit = [1.5; 0];
solinit = bvpinit(xinit, yinit);
xx = linspace(0,1);
yy = deval(sol, xx);
plot(xx, yy(1,:), 'LineWidth', 1.1)
shg
2
1.8
1.6
1.4
1.2
0.8
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
La función bvp4c permite, mediante un argumento opcional, modificar algunas de los parámetros que
utiliza el algoritmo por defecto. Para ello hay que usar la función en la forma
donde
options es una estructura que contiene los valores que se quieren dar a los parámetros opcionales.
Esta estructura se crea mediante la función bvpset que se explica a continuación.
nombre1, nombre2 son los nombres de las propiedades cuyos valores por defecto se quieren modificar.
Las propiedades no mencionadas conservan su valor por defecto.
valor1, valor2 son los valores que se quieren dar a las propiedades anteriores.
Para una lista de las propiedades disponibles, véase el help de MATLAB. Los valores actuales se pueden
obtener usando la función bvpset sin argumentos. Mencionamos aquí la propiedad que vamos a usar en
el ejercicio siguiente, que nos permite determinar el nivel de precisión que bvp4c utilizará para detener
las iteraciones y dar por buena una aproximación:
RelTol es el nivel de tolerancia de error relativo que se aplica a todas las componentes del vector
residuo ri = S 0 (xi ) − f (xi , S(xi )). Esencialmente, la aproximación se dará por buena cuando
S 0 (xi ) − f (xi , S(xi ))
≤ RelT ol
máx{|f (xi , S(xi ))|}
El valor por defecto de RelTol es 1.e-3.
Stats si tiene el valor ’on’ se muestran, tras resolver el problema, estadísticas e información sobre
la resolución. El valor por defecto es ’off’.
para el valor por defecto de RelTol y para RelTol=1.e-5 y comparar con la solución exacta.
Comenzamos por escribir el problema en forma de un problema de contorno para un SDO de primer
orden. Hacemos el cambio de variable z1 = y, z2 = y 0 y el problema (11.2) se transforma en:
0
z1 = z2
0
3pz1
z2 = − 2
p + t2
−0.1 0.1
z1 (−0.1) = √
, z1 (0.1) = √
p + 0.01 p + 0.01
function ejercicio3
% Parametros y constantes
%
p = 1.e-5;
yenb = 0.1/sqrt(p+0.01);
% Funciones anonimas
%
dydt = @(t,y) [ y(2); -3*p*y(1)/(p+t^2)^2];
% Aproximaciones iniciales
% La aproximacion inicial se ha obtenido promediando
% las condiciones de contorno
%
tguess = linspace(-0.1, 0.1, 10);
yguess = [0; 10];
solinit = bvpinit(tguess, yguess);
% Graficas
% Dibujamos solo la primera componente de la solucion del sistema
%
xx = linspace(-0.1, 0.1);
yy1 = deval(sol1, xx, 1);
yy2 = deval(sol2, xx, 1);
yye = exacta(xx);
plot(xx, yy1,'b', xx, yy2,'g', xx, yye,'r*')
shg
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
−0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1
Ejercicio 3: zoom
−0.6
Sol. con RelTol estandar
Sol. con RelTol=1.e−5
Sol. exacta
−0.65
−0.7
−0.75
−0.8
−0.85
−0.9
−0.95
−1
−0.04 −0.035 −0.03 −0.025 −0.02 −0.015 −0.01 −0.005 0
12 B
SE
UNIVER
VILLA
Problemas de minimización e an
Ecuaciones Diferenciales
y Análisis Numérico
Exponemos aquí brevemente la forma de utilizar las herramientas de MATLAB para resolver proble-
mas de minimización con restricciones.
Concretamente, queremos calcular una solución del problema
Minimizar J(x)
sujeto a x ∈ Rn
(12.1)
Gi (x) = 0 i = 1, . . . , me
G (x) ≤ 0 i = m + 1, . . . , m
i e
• Cotas: ci ≤ x ≤ cs
• Igualdades lineales: Ax = b, donde A es una matriz
• Desigualdades lineales: Cx ≤ d, donde C es una matriz
• Igualdades nolineales: f (x) = 0
• Desigualdades no lineales: g(x) ≤ 0
Cada tipo de restrición habrá que incorporarla de manera distinta a la función fmincon.
Los usos más sencillos de la función fmincon son los siguientes:
x = fmincon(funJ, x0, C, d)
x = fmincon(funJ, x0, C, d, A, b)
x = fmincon(funJ, x0, C, d, A, b, ci, cs)
x = fmincon(funJ, x0, C, d, A, b, ci, cs, resnolin)
x = fmincon(funJ, x0, C, d, A, b, ci, cs, resnolin, options)
donde:
243
12. Problemas de minimización 244
funJ es un manejador de la función a minimizar J(x). Como en otros casos, puede ser el nombre de
una función anónima o un manejador de una M-función.
C, d son, respectivamente, la matriz y el segundo miembro que definen las restricciones lineales de
desigualdad Cx ≤ d.
A, b son, respectivamente, la matriz y el segundo miembro que definen las restricciones lineales
de igualdad Ax = b. Si hay que utilizar A y b y nuestro problema no tiene restricciones con
desigualdades lineales, se pone C=[] y b=[].
ci, cs son, respectivamente, las cotas inferior y superior ci ≤ x ≤ cs . Si no hay cota inferior para
alguna de las variables se especificará -Inf. Análogamente, si no hay cota superior, se pondrá
Inf. Como antes, si alguno de los argumentos anteriores (C, d, A, b) no se usa, se utilizarán
matrices vacías ([]) en su lugar.
resnolin es un manejador de la función que define las restricciones no lineales, tanto de igualdad
f (x) = 0 como de desigualdad g(x) ≤ 0. Puede ser el nombre de una función anónima o un
manejador de una M-función. Esta función debe aceptar como argumento de entrada un vector
x y devolver dos vectores: c = g(x), ceq = f (x):
options es una estructura que permite modificar algunos de los parámetros internos de la función
fmincon. Esta estructura se crea mediante la función:
Aunque son numerosas las opciones que se pueden modificar, mencionamos aquí sólo dos:
Algorithm permite elegir entre cuatro algoritmos, en función de las caracterísitcas y/o di-
ficultad del problema a resolver. Los posibles valores son: 'interior-point', 'sqp',
'active-set' y 'trust-region-reflective'.
Display permite elegir el nivel de verbosity, es decir, de información suministrada durante la
ejecución. Posibles valores son: 'off', 'iter', 'final' y 'notify'.
Se trata en ambos casos de restricciones lineales de desigualdad, que se pueden expresar en la forma
Ñ é Ñ é
−1 −2 −3 0
Cx ≤ d con C = , d=
1 2 3 72
0 ≤ x1 + 2x2 + 3x3 ≤ 72
x2 + x2 + x2 ≤ 500
1 2 3
Las restricciones lineales son las mismas del ejercicio anterior. La restricción no lineal hay que
escribirla en forma homogénea
x21 + x22 + x23 − 500 ≤ 0
y escribir una M-función para definir dichas restricciones:
13 B
SE
UNIVER
VILLA
iniciales y de contorno para
e an
edp’s unidimensionales Ecuaciones Diferenciales
y Análisis Numérico
Dedicamos este capítulo a describir brevemente cómo resolver con MATLAB problemas de valores
iniciales y de contorno para sistemas de EDP parabólicas/elípticas unidimensionales.
donde:
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
donde:
247
13. Problemas de valores iniciales y de contorno para edp’s unidimensionales 248
pdefun es un manejador de la función que proporciona los coeficientes del sistema (funciones c, f y
s). Puede ser el nombre de una función anónima o un manejador de una M-función. Esta función
responde a la siguiente declración:
∂u
Recibe los valores x, t y u=u(x,t), dudx= (x,t) y devuelve los valores de c, f y s para estos
∂x
argumentos. Más concretamente, devuelve tres vectores c, f y s, de la misma dimensión que u
y dudx, de forma que f = f(x, t, u, dudx), s = s(x, t, u, dudx) y c es la diagonal de la matriz
c(x, t, u, dudx).
function u0 = icfun(x)
xmesh es un vector que define una partición del intervalo espacial [a, b]: xmesh = [x1, x2, ...xN],
con x1 = a < x2 ...< xN = b y N ≥ 3. Son los puntos del intervalo [a, b] en los que se va a
calcular la aproximación de la solución.
tspan es un vector que define una discretización del intervalo temporal [tO , tf ]: [t1, t2, ...tM],
con t0 = t1 < t2 ...< tM = tf . Son los tiempos en los que se van a calcular las aproximaciones
de la solución u del problema.
El producto cartesiano xmesh × tspan define una malla en cuyos nodos se va a calcular una
aproximación de la solución u(x, t):
Los vectores xmesh y tspan juegan diferentes papeles en pdepe. En particular, el coste y la
precisión de la solución dependen en gran medida de la longitud del vector xmesh. Sin embargo,
el cálculo es mucho menos sensible a los valores del vector tspan.
sol es una matriz tridimensional (esto es, con tres subíndices) que contiene la aproximación de
la solución: sol(:,:,k) son los valores de la k-ésima componente de u en los nodos de la
Ejemplo 13.1
Calcular una aproximación de la solución de:
2
2 ∂u = ∂ u
π (x, t) ∈ [0, 1] × [0, 2],
∂t ∂x2
∂u (13.2)
u(0, t) = 0, (1, t) = −πe−t , t ∈ [0, 2]
∂x
u(x, 0) = sen(πx), x ∈ [0, 1],
Así pues, la EDP en este ejemplo puede ser representada por la M-función
function u0 = icfun(x)
%
% funcion que proporciona el dato inicial para t = 0
%
u0 = sin(pi*x);
end
Creamos ahora los vectores que definen la malla. Para este problema usaremos una malla con 20
puntos equi-espaciados en el intervalo espacial [0, 1] y 10 puntos en el intervalo temporal [0, 2]:
xmesh = linspace(0,1,20);
tspan = linspace(0,2,10);
m = 0;
u = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
En este caso la solución, u, solo tiene una componente, de manera que u es una matriz con dos
subíndices. Podemos ahora representar la solución usando la función surf (ver su descripción en
la sección 3.4.2:
surf(xmesh,tspan,u)
title('Solucion numerica calculada con una malla de 20 puntos')
xlabel('Distacia x')
ylabel('Tiempo t')
Ñ é
u1
Ejercicio 13.2 Calcular la solución u = del siguiente sistema de reacción-difusión:
u2
∂u ∂2u
= D + Au, (x, t) ∈ [0, 1] × [0, 2],
∂x2
∂t
∂u1
)
(0, t) = 0, u2 (0, t) = 0,
∂x t ∈ [0, 2] (13.3)
∂u2
u (1, t) = 1, (1, t) = 0,
1
∂x
u (x, 0) = x2 , u (x, 0) = 0, x ∈ [0, 1],
1 2
En este ejemplo tenemos un problema de valores iniciales y de contorno para un sistema de dos
EDP parabólicas. La EDP en (13.3) se escribe, en forma desarrollada:
Ü ê Ü 2 ê Ü ê
∂u1 Ñ é ∂ u1 Ñ éÑ é ∂u1 Ñ é
∂t 1 0 ∂x 2 0 1 u1 ∂ ∂x u2
= + = +
∂u2 0 2
2
∂ u2 1 0 u2 ∂x ∂u2 u1
2
∂t ∂x 2 ∂x
Por lo tanto, se puede escribir en la forma general de la EDP de 13.1 para m = 0 y:
Ü ê
Ñ é ∂u1 Ñ é
∂u 1 0 ∂u ∂x ∂u u2
c x, t, u, = , f x, t, u, = , s x, t, u, =
∂x 0 1 ∂x ∂u2 ∂x u1
2
∂x
El código de la función pdefun se escribe (observemos que ninguno de los coeficientes c, f y s
depende explícitamente de x ni de t):
Las condiciones iniciales en este caso se escriben (recordemos que u es un vector columna de dos
componentes):
function u0 = icfun(x)
%
% funcion que proporciona el dato inicial
%
u0 = [x^2; 0];
end
y la M-función sería (recordamos que se devuelven las diagonales de q(0, t) y de q(1, t)):
De nuevo, usaremos una malla con 200 puntos equi-espaciados en el intervalo espacial [0, 1] y 50
puntos en el intervalo temporal [0, 2]:
xmesh = linspace(0,1,200);
tspan = linspace(0,2,50);
m = 0;
sol = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
La solución es devuelta por pdepe como una matriz 3D: sol(:,:,1) son los valores de u1 y
sol(:,:,2) son los valores de u2 :
subplot(1,2,1)
contour(xmesh,tspan,u1,50)
title('Primera componente de la solucion: u_1')
xlabel('Distacia x')
ylabel('Tiempo t')
colorbar
%
subplot(1,2,2)
contour(xmesh,tspan,u2,50)
title('Segunda componente de la solucion: u_2')
xlabel('Distacia x')
ylabel('Tiempo t')
colorbar
1.8 1.8
0.9
1.4 1.4
0.7
1.2 1.2
0.15
0.6
Tiempo t
Tiempo t
1 1
0.5
0.8 0.8
0.4 0.1
0.6 0.6
0.3
0.4 0.4
0.2 0.05
0.2 0.2
0.1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Distacia x Distacia x