Apuntes Matlab Metodos Numericos Sevilla 2020
Apuntes Matlab Metodos Numericos Sevilla 2020
Apuntes Matlab Metodos Numericos Sevilla 2020
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
6 de febrero de 2020
IDAD DE
S
SE
UNIVER
VILLA
Índice e an
Ecuaciones Diferenciales
y Análisis Numérico
3
2.4. Forma en que se muestran los resultados . . . . . . . . . . . . . . . . . . . . . . 51
2.5. Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.5.1. Instrucciones de asignación . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.5.2. Variables pre-definidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.6. Funciones matemáticas elementales . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.7. Operaciones de comparación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.8. Operadores lógicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.9. Orden general de evaluación de expresiones . . . . . . . . . . . . . . . . . . . . . 57
2.10. Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.10.1. Construcción de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.10.2. Acceso a los elementos de vectores y matrices . . . . . . . . . . . . . . . 62
2.10.3. Operaciones con vectores y matrices . . . . . . . . . . . . . . . . . . . . . 64
2.10.4. Operaciones de comparación y lógicas con matrices . . . . . . . . . . . . 66
2.11. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3. Gráficos 73
3.1. Generalidades sobre la representación gráfica de funciones . . . . . . . . . . . . . 73
3.1.1. Representación gráfica de funciones de una variable real . . . . . . . . . . 74
3.1.2. Curvas planas definidas mediante ecuaciones paramétricas . . . . . . . . 75
3.1.3. Curvas planas en coordenadas polares . . . . . . . . . . . . . . . . . . . . 75
3.1.4. Curvas en tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.1.5. Gráficas de funciones de dos variables: superficies . . . . . . . . . . . . . 79
3.1.6. Superficies definidas mediante ecuaciones paramétricas . . . . . . . . . . 80
3.1.7. Representación mediante curvas de nivel de una función de dos variables 81
3.2. Funciones gráficas de MATLAB fáciles de usar . . . . . . . . . . . . . . . . . . . 83
3.2.1. Gráficas de funciones de una variable real . . . . . . . . . . . . . . . . . 83
3.2.2. Curvas planas definidas implícitamente . . . . . . . . . . . . . . . . . . . 85
3.2.3. Curvas planas definidas por ecuaciones paramétricas . . . . . . . . . . . 86
3.2.4. Curvas planas en coordenadas polares . . . . . . . . . . . . . . . . . . . . 87
3.2.5. Sobre los nombres de las variables independientes . . . . . . . . . . . . . 88
3.2.6. Curvas en tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.2.7. Superficies definidas mediante ecuaciones explícitas . . . . . . . . . . . . 90
3.2.8. Superficies definidas mediante ecuaciones paramétricas . . . . . . . . . . 91
3.2.9. Representación mediante curvas de nivel de una función de dos variables 92
3.2.10. Observación importante . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.3. Funciones básicas para el dibujo de curvas . . . . . . . . . . . . . . . . . . . . . 94
3.3.1. La orden plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.3.2. Color y tipo de línea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.3.3. Personalizar las propiedades de las líneas . . . . . . . . . . . . . . . . . . 98
3.3.4. Personalizar los ejes: orden axis . . . . . . . . . . . . . . . . . . . . . . . 100
3.3.5. Anotaciones: título, etiquetas y leyendas . . . . . . . . . . . . . . . . . . 101
3.3.6. Dibujar cosas encima de otras . . . . . . . . . . . . . . . . . . . . . . . . 103
3.3.7. Añadir texto a la gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
3.3.8. Gestión de la ventana gráfica . . . . . . . . . . . . . . . . . . . . . . . . 105
3.3.9. Otros tipos de gráficos de datos planos . . . . . . . . . . . . . . . . . . . 105
3.3.10. Varios ejes en la misma gráfica . . . . . . . . . . . . . . . . . . . . . . . . 107
3.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
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 . . . . 97
3.2. Algunas propiedades de las líneas, controlables con la orden plot . . . . . . . . 99
3.3. Algunas propiedades de los títulos y las etiquetas de los ejes . . . . . . . . . . . 102
9
IDAD DE
S
SE
UNIVER
VILLA
Índice de figuras e an
Ecuaciones Diferenciales
y Análisis Numérico
11
4.1. Condicional simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.2. Condicional doble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.3. Condicional doble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
4.4. Repetición condicional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.5. Ruptura de bucle break . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.6. Ruptura de bucle continue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.1. Región delimitada por la gráfica de la función y = f (x), el eje de abscisas y las
rectas x = a y x = b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
IDAD DE
S
1 B
SE
Introducción a los
UNIVER
VILLA
ordenadores e an
Ecuaciones Diferenciales
y Análisis Numérico
El modelo básico de los ordenadores actuales se atribuye a von Neumann1 . Consta principalmen-
te 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.
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
MEMORIA CENTRAL
DATOS PROGRAMA
UNIDAD DE CONTROL
ENTRADA/SALIDA MEMORIA
REGISTRO DE
SECUNDARIA
INSTRUCCIONES (disco duro)
UNIDAD ARITMÉTICO-LÓGICA
REGISTRO DE
OPERACIONES
La unidad aritmético-lógica o ALU (Arithmetic Logic Unit), que realiza las opera-
ciones elementales que constituyen el programa (sumar, multiplicar, comparar, etc).
1.1.2 Memoria
Son los componentes de hardware en los que se almacena la información procesada por el
ordenador. 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).
Con el objetivo de que el procesador pueda obtener los datos de la memoria central más rápi-
damente, 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.
La memoria central es un circuito que se puede imaginar como una enorme tabla que almacena
informació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 correla-
tivamente, comenzando 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.
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
Memoria secundaria es una memoria más barata donde, a cambio de mayores tiempos de
acceso, se puede conservar gran cantidad de información (disco duro, CD (Compact Disk), DVD
(Digital Versatile Disk), memoria flash, etc). En ella se almacena todo lo que la unidad central
de procesos (véase la Sección 1.1.1) no necesita urgentemente. La principal característica de
este tipo de memorias es que su contenido no se pierde al apagar el ordenador.
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 ordena-
dor. Son el sistema operativo, los compiladores e intérpretes y los utilitarios.
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.
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
instrucciones 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
ensamblador 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:
2
Fue usado ampliamente en el pasado para el desarrollo de software, pero actualmente sólo se utiliza en
contadas ocasiones, especialmente cuando se requiere la manipulación directa del hardware o se pretenden
rendimientos inusuales del los equipos.
A = B
Ejemplo 1.1
5 22 -1
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)
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 valo-
res 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.
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.
3
Variaciones con repetición de 2 elementos (0, 1) tomados de n en n = 2n .
Ejemplo 1.8
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
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)
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 fraccio-
nario, 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
1 1 k
términos 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
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
53 3
− 1 0 1 = 10 +
0 0 1 1 5 5
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
| {z } 0010
1010001000001110 = 1010 | {z } 0000
| {z } 1110
| {z } = A20E
A 2 0 E
Notación octal La notación octal es similar a la hexadecimal, pero con el sistema de numera-
ción en base 8 (8 = 22 ). 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
Observación Distintos lenguajes de programación tienen diversas formas de escribir los núme-
ros binarios, 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 representa-
bles, 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
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 comple-
mento 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
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 nume-
ració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 − .
4
Las operaciones bitwise o bit a bit son las que se realizan sobre los números binarios actuando directamente
sobre sus bits, uno a uno. Son operaciones primitivas, que llevan a cabo dispositivos electrónicos específicos, y
son muy rápidas de realizar, incluso para procesadores sencillos de bajo costo. La operación de complemento
a 1 es la operación bitwise NOT.
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
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 comple-
mento 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 nume-
ració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
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, denominado 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úme-
ro 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
1 25 32 1 25 32 0 25 0 25 1 25 32
0 26 0 26 1 26 64 1 26 64 0 26
1 27 128 1 27 128 0 27 0 27 1 27 128
1 28 256 1 28 256 0 28 0 28 1 28 256
1 29 512 1 29 512 0 29 0 29 1 29 512
0 210 0 210 1 210 1024 1 210 1024 0 210
0 211 0 211 1 2 11
2048 1 2 11
2048 0 211
1 212 4096 1 212 4096 0 212 0 212 1 212 4096
0 213 0 213 1 213 8192 1 213 8192 0 213
0 214 0 214 1 214 16384 1 214 16384 0 214
1 215 32768 1 − 1 − 1 − 1 215 32768
37805 −5037 −27730 −27731 37805-32767=5038
Como flotante La coma flotante (en inglés, floating point) es un modo “normalizado” de
representació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 signi-
ficativa (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.
5
IEEE son las siglas del Institute of Electrical and Electronics Engineers.
Ejemplo 1.33
Sea, pues, un número R dado por la siguiente representación binaria en coma flotante:
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 existe el bit implícito.
Ejemplo 1.35
Representar en simple precisión el número R = −5.625(10) = −1.01101 · 210
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. Introducción
1 1 a0 los
0 0ordenadores
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 40
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 .
El mayor exponente positivo que se puede almacenar en doble precisión es, como hemos
visto, 21 0 = 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 1 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
exponente 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.
R1 = |1.00000000000000000000000
{z }
24 dígitos
R2 = |1.00000000000000000000001
{z }
24 dígitos
Ejemplo 1.36
R = 0.1(10) = 0.00011
˘
(2)
= 0.00011001100110011 . . .
=4
z}|{
− 100
| {z } 1100 · · · × 2
= = 1.10011001100110011001100
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
{z } ×2
R = |1.10011001100110011001101 .
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 )
2 B
SE
Introducción a MATLAB.
UNIVER
VILLA
Operaciones elementales e an
Ecuaciones Diferenciales
y Análisis Numérico
2.1 Introducción
Al iniciar matlab nos aparecerá una ventana más o menos como la de la Figura 2.1 (depen-
diendo 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 salva-
guardar 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
comunicaremos 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.
43
2. Introducción a MATLAB. Operaciones elementales 44
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, consti-
tuyen 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 expre-
sió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:
Las funciones anónimas pueden tener varios argumentos y hacer uso de variables previamente
definidas:
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.)
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,
utilizando el punto decimal (no la coma):
Ejemplo 2.7
23. -10.1 22.0765
2.e-2 = 2 × 10-2 = 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 /
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.
variable = expresión
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
>> (2+3i)*(1-i) % Producto de números complejos
ans =
5.0000 + 1.0000i
>> 2^2000 % Número no almacenable en el ordenador
ans =
Inf
>> 0/0 % Número imposible de calcular
ans =
NaN
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 corres-
pondientes 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 nu-
mé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 |
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.
Ejemplos 2.21
>> (1 > 5) & (5 < 10)
ans =
0 (logical)
>> (0 < 5) | (0 > 5)
ans =
A B A& B A| B
true true true true
true false false true
false true false true
false false false false
1 (logical)
Cuando aparecen valores numéricos en una operación lógica, matlab transforma el valor nu-
mé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)
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.
Observación.
El hecho de que, al introducirlas, se escriban las matrices por filas no significa que interna-
mente, 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
espaciados (ú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
Se puede acceder a una fila o columna completa sustituyendo el índice correspondiente por :
(dos puntos)
>> A(2, :)
ans =
6 3 1
>> A(3, :) = [5, 5, 5]
A =
1 2 3
6 3 1
5 5 5
2 0 4
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 ma-
trices 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 .
5
8
9
11
14
15
2.11 Ejercicios
1. Comprobar el valor de las siguientes expresiones:
π
r
π 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
! !
4 5
2 3 1 3
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
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π . . . ) f ) w6 = (1, 3, . . . , n − 1, n, n − 2, . . . , 4, 2)
c) w3 = (n, n − 1, n − 2, . . . , 2, 1) (suponiendo, aquí, que n es par)
d) w4 = (2n, 2n − 2, . . . , 4, 2)
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 si-
guientes 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 = ... .. .. . 0
. . , . . , 0
C= 1 0
.. .. .. ..
0 2 0 0 1 π . .
. .
0 0 ··· 0 4 0 0 ··· 0 1 0 ··· ··· 1 0
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
adecuadas).
c) Resolver el sistema lineal de ecuaciones M x = b, siendo b = (1, 2 . . . n)t . (Obser-
vació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
73
3. Gráficos 74
mucho más grande de dificultad. Por un lado, es preciso utilizar técnicas de geometría proyectiva
para determinar la perspectiva y conseguir impresión de tridimensionalidad. Por otro, aparece
la necesidad de utilizar algoritmos y técnicas complejas para determinar partes ocultas. Y, aún
más, iluminación, transparencias, aplicación de texturas, etc. Todo ello queda fuera del ámbito
de este curso.
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:
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).
0
−6 −4 −2 0 2 4 6 8
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).
θ, que es el ángulo que forma el segmento OP con una semirrecta fija de origen O deno-
minada eje polar.
t x y
10
0 0 1
Y
8
1 -1.5 2.4
6 t=10 2 -0.7 5.2
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 para-
8 5.0 4.4
mé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
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
En tal sistema de coordenadas, se dice que el par (r, θ) son las coordenadas polares del punto
P (ver Figura 3.7).
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
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 pola- Figura 3.8: Coordenadas cartesianas y po-
res. lares.
x = r cos(θ), y = r sen(θ);
p y
r = x2 + y 2 , θ = arctan
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.
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 coorde- Figura 3.10: Curva de ecuación, en coor-
nadas polares, r = θ, θ ∈ [0, 9π/2] denadas polares, r = 2 sen(6θ), θ ∈
[0, 2π].
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 las coordenadas de los puntos que definen la curva:
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
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 Ω.
Un mallado en rectángulos de un dominio rectangular es fácil de construir a partir de sendas
particiones de sus lados. Un mallado en triángulos es más complicado y precisa de algoritmos
y programas especializados.
Y
Y
X
X
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 dis-
poner 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 represen-
tación de la superficie como una red deformada, como en las Figuras 3.16 y 3.17.
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.
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 co- Figura 3.20: Mallado deformado con caras
lor plano. de color constante.
Figura 3.21: Mallado deformado con co- Figura 3.22: Mallado deformado con caras
lor plano en cada cara, dependiente de la de color interpolado a partir de los valores
altura. en los vértices.
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.
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
ezplot(f)
ezplot(f,[a,b])
Ejemplo 3.1
Representar, usando ezplot, la gráfica de la función
x sen(x2 )
f (x) = en [−2π, 2π]
2
Por defecto, la función ezplot hace variar la variable independiente en el intervalo [−2π, 2π].
Por lo tanto, basta escribir en la línea de comandos
ezplot('sin(x^2)*x/2')
3
−1
−2
−3
−6 −4 −2 0 2 4 6
x
Ejemplo 3.2
Representar, usando ezplot, la gráfica de la función
f (x) = ln(x)
ezplot('log(x)')
2
1.5
0.5
−0.5
−1
−1.5
−2
0 1 2 3 4 5 6
x
Ejemplo 3.4
Representar, usando ezplot, la gráfica de la función
x sen(x2 )
f (x) = para x ∈ [0, π]
2
f=@(x) x*sin(x^2)/2;
ezplot(f, [0,pi])
1.5
0.5
−0.5
−1
ezplot(f)
ezplot(f,[a,b])
ezplot(f,[a,b,c,d])
f es la función (de dos variables) que describe la ecuación, y puede ser especificada como
sigue:
[a,b] (opcional) la curva se dibuja en el cuadrado del plano [a, b] × [a, b]. Si no se especifica,
ezplot usa por defecto [a, b] = [−2π, 2π].
[a,b,c,d] (opcional) la curva se dibuja en el rectángulo del plano [a, b] × [c, d].
Ejemplo 3.5
Dibujar, usando ezplot, la circunferencia x2 + y 2 = 7
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
ezplot('x^2+y^2-7')
y
−2
−4
−6
−6 −4 −2 0 2 4 6
x
Utilizando una función anónima, hubiéramos tenido que definir la función de dos variables
f (x, y) = x2 + y 2 − 7, y luego usar ezplot
ezplot(@(x,y) x^2+y^2-7)
ezplot(f, g)
ezplot(f, g, [a,b])
f, g son las dos funciones que definen la curva. Pueden ser especificadas, como antes, me-
diante sus expresiones o mediantes manejadores.
Ejemplo 3.7
Dibujar, usando ezplot, la curva definida por las ecuaciones (paramétricas)
ezplot('sin(3*t)', 'cos(t)')
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
y
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
x x
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)
ezpolar(f)
ezpolar(f, [a,b])
[a, b] (opcional) son los límites del intervalo en el que varía θ, es decir donde varía el
ángulo. Si no se especifica, ezpolar toma por defecto [a, b] = [0, 2π].
Ejemplo 3.9
Representar, usando ezpolar, la gráfica de la siguiente función definida en coordenadas
polares
r = sen(2θ) cos(3θ), para θ ∈ [0, π]
ezpolar('sin(2*t)*cos(3*t)', [0,pi])
90
1
120 60
0.8
0.6
150 30
0.4
0.2
180 0
210 330
240 300
270
r = sin(2 t) cos(3 t)
ezpolar('nombre*sin(nombre)')
Por otra parte, cuando se trata de una función de dos variables, y si no se especifica el orden,
MATLAB tomará como primer argumento la primera en orden alfabético y como segundo,
la segunda. Como ejemplo, compruébese cómo las dos órdenes siguientes producen resultados
distintos:
ezplot('1+v/(sin(a)+2)')
ezplot('1+v/(sin(w)+2)')
Esto no sucede si se utilizan funciones anónimas, ya que se indica expresamente el orden de las
variables.
ezplot3(f, g, h)
ezplot3(f, g, h, [a,b])
f,g,h son las expresiones de las tres funciones dependientes de t. Se proporcionan en cual-
quiera de las formas habituales
Ejemplo 3.11
Dibujar, usando ezplot3, la siguiente curva tridimensional:
√
x = 3 cos(t), y = t sen(t2 ), z = t, t ∈ [0, 2π]
2.5
1.5
z
0.5
0
10
5 3
2
0 1
0
−5 −1
−2
−10 −3
y
x
ezmesh(f)
ezmesh(f, [a,b])
ezmesh(f, [a,b,c,d])
f es la expresión de la función (de dos variables) f (x, y). Se puede proporcionar en cualquiera
de las formas habituales.
[a,b] (opcional) establece el rectángulo para dibujar la superficie: [a, b] × [a, b]. Si no se
especifica, ezmesh toma [a, b] = [−π, π].
Ejemplo 3.13
Dibujar la superficie z = x e−(x +y )
2 2
ezmesh('x*exp(-x^2-y^2)')
0.5
−0.5
3
2
1 3
2
0 1
−1 0
−1
−2 −2
−3 −3
y
x
ezmesh(f, g, h)
ezmesh(f, g, h, [a,b])
ezmesh(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, ezmesh
toma [a, b] = [−2π, 2π].
[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.14
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
ezmesh('cos(s)*(2+cos(t))','sin(s)*(2+cos(t))','t',[0,2*pi])
4
z
0
3
2
3
1 2
0 1
−1 0
−1
−2 −2
−3 −3
y
x
ezcontour(f)
La función ezcontourf actúa como el anterior, pero rellena con colores los espacios entre curvas
ezcontourf(f)
La función ezmeshc actúa igual que ezmesh , pero dibuja también las curvas de nivel
ezmeshc(f)
Ejemplo 3.16
Dibujar usando ezcontourf la curvas de nivel de la función
2 +y 2 )
f (x, y) = x e−(x
ezcontourf('x*exp(-x^2 - y^2)')
3
0
y
−1
−2
−3
−3 −2 −1 0 1 2 3
x
Ejercicio 3.17 Dibujar usando ezmeshc 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).
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.18
x2 + 2
Usando la función plot, dibujar la curva y = para x ∈ [−2, 3]
x+5
y = (x.^2+2)./(x+5); 1.6
plot(x,y) 1.4
1.2
resultar equívoco.
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.19
1
Usando la función plot, dibujar las curvas y = x + x2 − x e y = ex + sen(x2 ) en el
e
intervalo [−π, π].
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.20
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
-2
yg = g(xg); -4
-8
-6 -4 -2 0 2 4 6
cos(x/4)
f (x) = ,
ln(1 + x2 )
Ejemplo 3.22
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)
t = linspace(0, 2*pi);
5.5
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.25
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.
Para información sobre el resto de propiedades de las líneas, se puede consultar, en el Help de
matlab, Line Properties.
Tabla 3.2: Algunas propiedades de las líneas, controlables con la orden plot
Ejemplo 3.26
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
'MarkerEdgeColor', 'blue',... 1
'MarkerFaceColor', 'g',... 0
'LineWidth', 2) -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.27
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.28
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); 5
y = 4 + 2*sin(t);
plot(x, y, 'LineWidth', 2) 4
axis([-3, 1, 1, 7]) 3
axis equal
2
1
-4 -3 -2 -1 0 1 2
grid on
axis off
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.29
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.30
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.31
1
0.9
'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.32
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);
ht.Interpreter = '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.33 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.34
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.35
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.36
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
3.4 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 represen-
ten (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
r
1−x e) f (x) =
a) f (x) = x
1+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.
a) f (x) = 2 sen3 (x) cos2 (x) y g(x) = ex − 2x − 3, x ∈ [−1.5, 1.5]
b) f (x) = log(x + 1) − x y g(x) = 2 − 5x, x ∈ [0, 5]
c) f (x) = 6 sen(x) y g(x) = 6x − x3 , x ∈ [−π/2, π/2]
−x2 +2
d ) f (x) = e sen(x/2) y g(x) = −x3 + 2x + 2, x ∈ [−1, 2]
√
e) f (x) = r2 − x2 , para r = 1 y r = 4
√ !
x 2 + 4x
f ) g(x) = 3x + 2x − 6 y h(x) = sen en [1, 2] × [0, 1.5]
x3
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
programas.
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.
111
4. Programación con MATLAB 112
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! <<')
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
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.
if x*y == 0
n = 0; true
n = 0
x*y == 0 fin
return
false
end
if (y > 0)
true true
n = 1; n = 2 y > 0 y > 0 n = 1
end
else
n = 3;
if (y > 0)
n = 2;
end Salida: n
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
Naturalmente, este mecanismo precisa que, dentro del bloque de instrucciones se modifique,
en alguna de las repeticiones, el resultado de evaluar expresion. En caso contrario, el progra-
ma entraria en un bucle infinito. Llegado este caso, se puede detener el proceso pulsando la
combinación de teclas CTRL + C.
instruccion-anterior
while expresion
bloque-de-instrucciones
end
instruccion-siguiente
suma = suma + n;
end suma <= 5000
false
n = n - 2;
true
n = n + 2
end
suma = suma + n
n = n - 2
Salida: n
fin
%---------------------------------------------------------------
% 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')
fprintf(' Pulsa CTRL + C si te rindes ... \n') inicio
M = 0;
while M ~= n M = 0
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.
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
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, ade-
má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 Figura 4.6: Ruptura de bucle continue: Si
algún momento de un bucle de repetición en algún momento de un bucle de repetición
(de cualquier tipo) se llega a una instrucción (de cualquier tipo) se llega a una instrucción
break, el ciclo de repetición se interrumpe continue, la iteración en curso se interrum-
inmediatamente y el programa continúa eje- pe inmediatamente, y se inicia la siguiente
cutándose por la instrucción siguiente a la iteración (si procede).
cláusula end del bucle. Si se trata de un bu-
cle indexado for, la variable-índice del mis-
mo 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í necesitan llevarlo, en general, las
instrucciones incluidas dentro.
• Las M-funciones no necesitan llevar end al final. De hecho, a nivel de este curso, que lo
lleven puede inducir a error. Mejor no ponerlo.
• Todos los programas deben incluir unas líneas de comentarios que expliquen de forma
sucinta su cometido y forma de utilización. En las M-funciones, estas líneas deben ser las
siguientes a la instrucción function, ya que de esta manera son utilizadas por MATLAB
como el texto de help de la función. En los scripts deben ser las primeras líneas del
archivo.
• Cuando hay errores en un programa, hay que leer con atención los mensajes de error
que envía matlab. Contienen mucha información, que en la mayoría de los casos es
suficiente para ayudar a detectar el error.
• Hay que comprobar los programas que se hacen, dándoles valores a los argumentos para
verificar que funciona bien en todos los casos que puedan darse. Esto forma parte
del proceso de diseño y escritura del programa.
• Los espacios en blanco, en general, hacen el código más legible y por lo tanto más fácil
de comprender y corregir. Se deben dejar uno o dos espacios entre el carácter % de cada
línea de comentario y el texto del comentario (comparénse los dos códigos 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:
a. Llamar S a una suma, mejor que G
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
e. Llamar e ó err ó error a un error, en lugar de vaca
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 correspon-
diente.
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 componen-
tes 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)
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
mediente 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:
• Dos soluciones reales y distintas
• Una solución real doble
• Dos soluciones complejas conjugadas
18. (Factorial.m) Escribe una M-función
que reciba como argumento un numero entero positivo n y calcule su factorial por mul-
tiplicaciones 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.
19. (Medias.m) Escribe una M-funció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 ):
N
X
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)! i=0
(2i)!
23. (SumaExp.m)
∞
x x2 x3 xn X xi
e =1+x+ + + ··· + + ··· = ∀x ∈ R
2! 3! n! i=0
i!
24. (SumaLog.m)
r ∞
1 + x x3 x5 x2n+1 X x2i+1
ln =x+ + + ··· + + ··· = ∀x tal que |x| < 1
1−x 3 5 2n + 1 i=0
2i + 1
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
entero 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.
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)
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.
40. (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 adi- 5.50 e + 0.60 e
cional por cada kg o adicionales por cada
fracción de kg, a par- kg o fracción de kg, a
tir de 1 kg de peso. partir de 5 kg de peso.
Aire 3e 3 e + 1 e adicio- 10.20 e + 1.20 e
nal por cada kg o frac- adicionales por cada
ción de kg, a partir de kg o fracción de kg, a
1 kg de peso. partir de 5 kg de peso.
Urgente 18 e 18 e + 6 e adi- No se realizan envíos
cionales por cada kg o urgentes de paquetes
fracción de kg, a partir con un peso superior a
de 1 kg de peso. 5 kg.
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”
41. (Camino.m) Escribir una M-función
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 41, 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
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
operaciones 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
>> A*x - b
ans =
1.0e-15 *
-0.4441
0
0
det(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
etiquetarí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
determinante 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)
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.
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 n
P
xi = bi − Aij xj
Aii j=i+1
Fin
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
i+j−1
de 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.
Una vez calculada la matriz L, para calcular la solución del sistema lineal de ecuaciones
Ax = b
utilizando la factorización anterior sólo hay que plantear el sistema, equivalente al anterior,
Lv = b
A x = b ⇐⇒ L Lt x = b ⇐⇒
Lt x = v
Ejercicio 5.10 Escribir una M-función function [x, res] = CHOL(A, b) que calcule
la solució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)
x = linsolve(A, b, opts)
resuelve el sistema lineal Ax = b por el método que resulte más apropiado, dadas las propie-
dades de la matriz A, que se pueden especificar mediante el argumento opcional opts (véase la
documentación).
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 matrices llenas. En general, las operaciones con matrices llenas producen matrices llenas y
las operaciones con matrices sparse producen matrices sparse. Las operaciones mixtas general-
mente 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
UNIVER
VILLA
escritura con ficheros en
e an
Ecuaciones Diferenciales
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 per-
mitidos 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.
155
6. Operaciones de lectura y escritura con ficheros en MATLAB 156
de Word), ficheros que contienen libros de hojas de cálculo (Excel), ficheros de imágenes,
ficheros que contienen programas ejecutables (es decir, escritos en código-máquina), etc.
Sólo se pueden utilizar con programas adecuados para cada fichero.
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 nombredefichero.mat. 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 memo.mat 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
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(’nombredefichero.extension’) 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
1. Puesto que el fichero temperaturas.dat es de texto, se puede editar. Ábrelo para ver
su estructura. Cierra el fichero.
2. Comienza por «cargar» en memoria los datos del fichero temperaturas.dat, 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('temperaturas.dat');
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)
5. Añade un título, una etiqueta en el eje OX y una leyenda, para la correcta comprensión
de esta gráfica.
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 datos3.dat 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('datos3.dat','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.
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
almacenada 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 contaminacion.dat que contiene datos sobre contaminación de
las provincias 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('contaminacion.dat','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’.
frewind(fid);
Esta orden indica que se han de leer del fichero 1 dato tipo string (cadena de carac-
teres), 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.
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
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 en MATLAB 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 va-
lores 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:
169
7. Interpolación y ajuste de datos en MATLAB 170
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
(polinomio 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
que pase por dichos puntos, es decir tal que
y1 = ax21 + bx1 + c
y2 = ax22 + bx2 + c
y = ax2 + bx + c
3 3 3
(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 xN 1 + c2 x N
1 + · · · + cN −1 x1 + cN
−1 −2
= c1 x N + c2 x N + · · · + cN −1 x2 + cN
y
2 2 2
...
N −1 −2
+ c2 x N + · · · + cN −1 xN + cN
y
N = c1 x N N
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 )
c=polyfit(x,y,N-1)
calcula los coeficientes del polinomio de grado N-1 que pasa por los N 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 interpola-
dos y dibujando también los ejes coordenados.
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
interpolació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
6
−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)
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).
100
80
60
40
20
−20
0 2 4 6 8 10 12 14 16
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.
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
interpolación de grado 10.
1. Crea dos variables x e y con las abscisas y las ordenadas que vas a interpolar:
s1 = interp1(x,y,1)
c = polyfit(x,y,10);
s = polyval(c,1) % = -245.7853
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)
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 Datos7.dat 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 aproxi-
mar funciones 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).
2. El spline que calcula MATLAB por defecto (denominado not-a-knot).
3. El spline sujeto con pendiente = −1 en x = 0 y pendiente = 2 en x = 10.
Observación: Este ejercicio no es imprescindible.
2
sen(x)
spline not−a−knot
spline sujeto
1.5
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
ordenadas 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);
3. Dibuja la recta de regresión (recuerda que para dibujar una recta bastan dos puntos):
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);
5. Dibuja la parábola
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
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 me-
diante 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
y
1
b 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 datos13.dat me-
diante 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.
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
Dada f : [a, b] ⊂ R 7→ R, continua, se plantea el problema de encontrar soluciones de la
ecuació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
187
8. Resolución de ecuaciones no lineales 188
es bien conocido que para n ≥ 5, no existe una fórmula explícita de sus soluciones.
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.
p = [2, 4, 0, 1];
roots(p)
obtendras las raíces: -2.1121, 0.0560 + 0.4833i y 0.0560 - 0.4833i. Por consi-
guiente, en el campo real, la única solución de la ecuación es x = −2.1121.
f (V ) = pV 3 − (pbN + kN T )V 2 + aN 2 V − abN 3 = 0
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.4, MATLAB dispone de la función
donde
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)
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 ??.
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 suce-
sivos 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 ?? y el ?? al final del capítulo.
Ejercicio 8.8 El fichero Datos.dat 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('Datos.dat');
xs = mat(:,1);
ys = mat(:,2);
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('Datos.dat');
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.
4) Continuar este proceso hasta que el subintervalo elegido tenga una longitud lo suficiente-
mente 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
ordenador, 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
|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 .
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 aproxima-
ción.
Ejercicio 8.12 Escribir una M-función que aproxime la solución de f (x) = 0 utilizando el
método de Newton.
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 elementales, como por ejemplo f (x) = e−x .
2
integral(fun,a,b);
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 repre-
sentar 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.
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.
La función f (x) = x e−x + 0.2 es positiva en el intervalo [0, 8], por lo tanto el valor
0.8
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
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
Ejercicio9.8 Calcular
x el área de la región del primer cuadrante delimitada por la curva
y = 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 localizació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
plot(x,yf,x,yg)
grid on -20
-30
2. La observación de la gráfica nos muestra que ambas curvas se cortan en dos puntos,
-6 -4 -2 0 2 4 6
(x1 , y1 ) y (x2 , y2 ).
Para calcularlos comenzamos por calcular sus abscisas x1 y x2 , que son las soluciones
de la ecuación
f (x) = g(x) equivalente a f (x) − g(x) = 0
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
y= = f (x) 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 para-
mé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
(
f 0 (t) = −80 sen(80t) + 3 cos2 (t) sen(t)
g 0 (t) = 80 cos(80t) − 3 sen2 (t) cos(t)
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
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 alrede-
dor 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 rota-
ción alrededor del eje OX de la región plana encerrada dentro de la elipse de ecuación
(x − 3)2 + 3y 2 = 4.
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
conocen 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)
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áfi-
camente 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:
x = [0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15];
y = [0, 20, 30, -10, 10, 10, 10.5, 15, 50 , 60 , 85];
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)
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
cuadrante (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:
x = [0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15];
y = [0, 20, 30, -10, 10, 10, 10.5, 15, 50 , 60 , 85];
2. Comienza por representar gráficamente el spline para analizar su signo:
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
A1 = -integral(s,0,x1) % A1 = 4.7719
A2 = -integral(s,x2,x3); % A2 = 8.6921
5. El área total buscada es finalmente:
A = A1 + A2 % A = 13.4640
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 datos21.dat contiene una matriz con dos columnas, que corres-
ponden 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 datos21.dat 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('datos21.dat');
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
1
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 datos22.dat contiene una matriz con dos columnas, que corres-
ponden 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:
(
0.5x2 si x ≥ 0
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 cal-
cular la función f :
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
%
function [y] = mifun(x)
y = zeros(size(x));
for k = 1:length(x)
xk = x(k);
if xk>=0
y(k) = 0.5*xk^2;
else
y(k) = xk^3;
end
end
Observación 2
En realidad, utilizando adecuadamente las operaciones vectoriales de Matlab y los opera-
dores 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
sistemas diferenciales e an
Ecuaciones Diferenciales
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.
y 0 = f (t, y), f : Ω ⊂ R2 → R
admite, en general, infinitas soluciones. Si, por ejemplo, f ∈ C 1 (Ω; R)1 , por cada punto (t0 , y0 ) ∈
Ω pasa una única solución ϕ : I → R, definida en un cierto intervalo I ⊂ R, que contiene a t0 .
Se denomina Problema de Cauchy (PC) o Problema de Valor Inicial (PVI)
y 0 = f (t, y),
(10.1)
y(t ) = y
0 0
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.
1
En realidad basta con que f sea continua y localmente Lipschitziana con respecto de y en Ω.
225
10. Resolución de ecuaciones y sistemas diferenciales ordinarios 226
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.
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
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, ten-
drí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)
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
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 f (t, y1 , y2 , . . . , yn ) y
Y = , F (t, Y ) = 2 , Y0 = 2 , (10.5)
... ... ...
0
yn fn (t, y1 , y2 , . . . , yn ) yn
Con esta notación el problema planteado se escribe:
Y 0 = F (t, Y )
(10.6)
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.
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
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 efecto, se tiene:
z10 = y 0 = z2
z1 (t0 ) = y(t0 ) = y0
z 0 = y 00 = z
z (t ) = y 0 (t ) = y
2 3 2 0 0 1
...
...
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) = Z
0
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)
0.4
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
11 B
SE
UNIVER
VILLA
contorno para sistemas e an
Ecuaciones Diferenciales
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 busca 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
(
y 0 = f (x, y), x ∈ [a, b]
(11.1)
g y(a), y(b) = 0.
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) per-
miten tambié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).
237
11. Resolución de problemas de contorno para sistemas diferenciales ordinarios 238
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.
ccfun es un manejador de la función que evalúa las condiciones de contorno g(ya , yb )
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:
sol.x contiene la malla final del intervalo [a, b].
sol.y contiene la aproximación numérica de la solución y(x) en los puntos sol.x.
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
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.
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 condi-
ciones de contorno
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
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». Uti-
lizamos 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
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áme-
tros 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 opcio-
nales. 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
0
S (xi ) − f (xi , S(xi ))
máx{|f (xi , S(xi ))|}
≤ RelT ol
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.
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];
condcon = @(ya, yb) [ya(1)+yenb; yb(1)-yenb ];
exacta = @(t) t./sqrt(p+t.^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
Resolución de problemas de
UNIVER
VILLA
minimización e an
Ecuaciones Diferenciales
y Análisis Numérico
Exponemos aquí brevemente la forma de utilizar las herramientas de MATLAB para resolver
problemas 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
Gi (x) ≤ 0 i = me + 1, . . . , m
• 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)
245
12. Resolución de problemas de minimización 246
donde:
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.
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.
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:
0 ≤ x1 + 2x2 + 3x3 ≤ 72
2
x1 + x22 + x23 ≤ 500
Las restricciones lineales son las mismas del ejercicio anterior. La restricción no lineal hay
que escribirla en forma homogénea
Cálculo Numérico II
Curso 2015/16
Índice
3
0. 4
4. Integración numérica 57
4.1. La función integral (quad en versiones antiguas de MATLAB) . . . . . . . . . . . . . . 57
4.2. La función trapz (fórmula de los trapecios) . . . . . . . . . . . . . . . . . . . . . . . . 61
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: http://www.mathworks.es
1.1 Comenzando
Al iniciar MATLAB nos aparecerá una ventana más o menos como la de la Figura 1.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
5
1. Introducción y elementos básicos 6
interés. La mayoría de las funciones tienen mas y/o distintas funcionalidades que las que se exponen
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 Tabla 1.1)
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 recu-
perar 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:
Tabla 1.1: Algunas constantes pre-definidas en MATLAB. Sus nombres son reservados y no deberían
ser usados para variables ordinarias.
Se puede utilizar MATLAB como simple calculadora, escribiendo expresiones aritméticas y terminando
por RETURN (<R>). Se obtiene el resultado inmediatamente a través de la variable del sistema ans
(de answer ). Si no se desea que MATLAB escriba el resultado en el terminal, debe terminarse la orden
por punto y coma (útil, sobre todo en programación).
Ejemplos 1.1
>> 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
1.2.4 Variables
En MATLAB las variables no son nunca declaradas: su tipo y su tamaño cambian de forma dinámica
de acuerdo con los valores que le son asignados. Así, una misma variable puede ser utilizada, por
ejemplo, para almacenar un número complejo, a continuación una matriz 25 × 40 de números enteros
y luego para almacenar un texto. Las variables se crean automáticamente al asignarles un contenido.
Asimismo, es posible eliminar una variable de la memoria si ya no se utiliza.
Ejemplos 1.2
>> a=10
a =
10.
>> pepito=exp(2.4/3)
pepito =
2.2255
>> pepito=a+pepito*(4-0.5i)
pepito =
18.9022 - 1.1128i
>> clear pepito
Para conocer en cualquier instante el valor almacenado en una variable basta con teclear su nombre
(Atención: recuérdese que las variables AB, ab, Ab y aB son distintas, ya que MATLAB distingue entre
mayúsculas y minúsculas).
Otra posibilidad es hojear el Workspace ó espacio de trabajo, abriendo la ventana correspondiente. Ello
nos permite ver el contenido de todas las variables existentes en cada momento e, incluso, modificar
su valor.
Algunos comandos relacionados con la inspección y eliminación de variables se describen en la tabla
siguiente:
1.2.5 Formatos
Por defecto, MATLAB muestra los números en formato de punto fijo con 5 dígitos. Se puede modificar
este comportamiento mediante el comando format
>> 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.
1.4.1 Scripts
En términos generales, en informática, un script (guión o archivo por lotes) es un conjunto de instruc-
ciones (programa), usualmente simple, guardadas en un fichero (usualmente de texto plano) que son
ejecutadas normalmente mediante un intérprete. Son útiles para automatizar pequeñas tareas. También
puede hacer las veces de un "programa principal" para ejecutar una aplicación.
Así, para llevar a cabo una tarea, en vez de escribir las instrucciones una por una en la línea de
comandos de MATLAB, se pueden escribir las órdenes una detrás de otra en un fichero. Para ello se
puede utilizar el Editor integrado de MATLAB. Para iniciarlo, basta pulsar el icono hoja en blanco (New
script) de la barra de MATLAB, o bien
File → New → Script
UN script de MATLAB debe guardarse en un fichero con sufijo .m para ser reconocido. El nombre
del fichero puede ser cualquiera razonable, es decir, sin acentos, sin espacios en blanco y sin caracteres
«extraños».
Para editar un script ya existente, basta hacer doble-click sobre su icono en la ventana Current Folder.
Para ejecutar un script que esté en el directorio de trabajo, basta escribir su nombre (sin el sufijo) en
la linea de comandos.
1.4.2 M-Funciones
Una función (habitualmente denominadas M-funciones en MATLAB), es un programa con una «inter-
faz» de comunicación con el exterior mediante argumentos de entrada y de salida.
Las funciones MATLAB responden al siguiente formato de escritura:
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 1.3
El siguiente código debe guardarse en un fichero de nombre areaequi.m.
La primera linea 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 1.4
>> areaequi(1.5)
ans =
0.9743
>> rho = 4 * areaequi(2) +1;
>> sqrt(areaequi(rho))
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 1.5
>> help areaequi
areaequi(long) devuelve el area del triangulo
equilatero de lado = long
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:
Las funciones anónimas pueden tener varios argumentos y hacer uso de variables previamente definidas:
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.)
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 linea 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 linea 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 de una funcion pertenezca al Base Workspace, hay que declararla
global: la orden
en una función hace que las variables a, suma y error pertenezcan al Base Workspace y por lo tanto,
seguirán estando disponibles cuando finalize la ejecución de la función.
1.6 Matrices
Como ya se ha dicho, 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
exhaustiva (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.
Se pueden también utilizar los vectores/matrices como objetos para construir otras matrices (bloques):
>> v1 = 1:4;
Las siguientes funciones generan vectores de elementos regularmente espaciados, útiles en muchas
circunstancias, especialmente para creación de gráficas.
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 → Mathematics → Arrays and Matrices
Se puede acceder a una fila o columna completa sustituyendo el índice adecuado por : (dos puntos)
También se puede acceder a una “submatriz” especificando los valores de las filas y columnas que la
definen:
Los operadores aritméticos representan las correspondientes operaciones matriciales siempre que tengan
sentido.
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.
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 .
Sea A una matriz cualquiera y sea b un vector columna con el mismo número de filas que A. Entonces,
la solución del sistema
Ax = b
se calcula en MATLAB mediante el operador no estándar \ denominado backward slash, (barra inversa
en español)
x = A \ b
x = mldivide(A, b) (equivalente a lo anterior)
Como comprobación calculamos el resíduo 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
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 cons-
truirla, ya que los segmentos serán imperceptibles.
Para dibujar una curva plana en MATLAB se usa el comando
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.
Ejemplo 1.18
x2 + 2
Dibujar la curva y = para x ∈ [−2, 3]
x+5
Se pueden dibujar dos o más curvas de una sóla vez, proporcionando al comando plot varios pares de
vectores abscisas-ordenadas, como en el ejemplo siguiente.
Ejemplo 1.19
Dibujar las curvas y = 2 sen3 (x) cos2 (x) e y = ex − 2x − 3 para x ∈ [−1.5, 1.5]
A cada par de vectores abscisas-ordenadas en el comando plot se puede añadir un argumento opcional
de tipo cadena de caracteres, que modifica el aspecto con que se dibuja la curva. Para más información,
hojear el help (help plot).
Ejemplo 1.20
Las siguientes órdenes dibujan la curva y = sen3 (x) en color rojo (r, de red) y con marcadores ∗,
en lugar de una línea contínua.
>> x = linspace(0,pi,30);
>> plot(x,sin(x).^3,'r*')
La orden siguiente dibuja la curva y = sen3 x en color negro (k, de black), con marcadores ∗ y con
línea punteada, y la curva y = cos2 x en color azul (b, de blue) y con marcadores +
Ademas, mediante argumentos opcionales, es posible modificar muchas otras propiedades de las curvas.
Esto se hace siempre mediante un par de argumentos en la orden de dibujo que indican
Esta propiedad afectará a todas las curvas que se dibujen con la misma orden.
Ejemplo 1.21
Para establecer un grosor determinado de las líneas se usa la propiedad LineWidth (atención a las
mayúsculas):
Se pueden añadir elementos a la gráfica, para ayudar a su comprensión. Para añadir una leyenda que
identifique cada curva se usa el comando siguiente, que asigna las leyendas a las curvas en el orden en
que han sido dibujadas.
legend('Leyenda1', 'Leyenda2')
Para añadir etiquetas a los ejes que aclaren el significado de las variables se usan los comandos
grid on
También es muy útil la orden, siguiente, que define las coordenadas mínimas y máxima del rectángulo
del plano OXY que se visualiza en la gráfica.
Cada nueva orden de dibujo borra el contenido previo de la ventana gráfica, si existe. Para evitar esto
existen la órdenes
hold on
hold off
La orden hold on permanece activa hasta que se cierre la ventana gráfica o bien se dé la orden hold
off
Ejemplos 1.22
Las siguientes órdenes darán como resultado la gráfica de la figura 1.3
x = linspace(0,pi,30);
axis([-0.5,pi+0.5,-0.5,1.5])
hold on
plot(x,sin(x).^3,'g', x, cos(x).^2,'b+', 'LineWidth', 1.1)
x = linspace(-0.95, pi);
plot(x, log(x+1)/2, 'r', 'LineWidth', 1.1)
plot([-5,5], [0, 0], 'k', 'LineWidth', 1)
plot([0, 0], [-5,5], 'k', 'LineWidth', 1)
legend('Leyenda1', 'Leyenda2', 'Leyenda3')
xlabel('Etiqueta del eje OX')
ylabel('Etiqueta del eje OY')
hold off
shg
1.5
Leyenda1
Leyenda2
Leyenda3
1
Etiqueta del eje OY
0.5
−0.5
−0.5 0 0.5 1 1.5 2 2.5 3 3.5
Etiqueta del eje OX
1.8 Ejercicios
3. Almacenar, en una variable de nombre n, un número entero par n ≥ 4. Generar los vectores de
longitud n variable que se indican:
cos(x/4) x
(4.a) f (x) = (4.c) f (x) =
ln(2 + x2 ) ln(2 + x)
1
(4.b) f (x) = e(−x
2 +2)
sin(x/2) (4.d) f (x) = x + x2 − x
e
5. Definir funciones anónimas para calcular las funciones siguientes. Debe hacerse de forma que x
pueda ser un vector.
1
r
1−x (5.e) f (x) =
(5.a) f (x) = x
1+x
x 1
(5.b) f (x) = (5.f) f (x) = ln
ln(x) x−2
√
√ √ (5.g) f (x, r) = r2 − x2
(5.c) f (x) = 1 − x + x − 2 2
sen (4t)
(5.d) f (x) = 3 (5.h) f (t) =
cos(5t)
6. Dibujar las curvas definidas por las funciones del Ejercicio 5, utilizando para ello las funciones
anónimas allí definidas y la orden plot. (Dibujar la función del apartado (5.g) para r = 5: f (5, x).
Dibujar las dos componentes de la función del apartado (5.h)).
7. Representar en la misma ventana gráfica las siguientes funciones, de modo que cada curva tenga
un color diferente:
8. Representar en una misma ventana las funciones siguientes, incluyendo una leyenda y etiquetas
en los ejes.
√ !
x2 + 4x
(8.a) g(x) = 3 + 2x − 6 y h(x) = sen
x en el cuadro [1, 2] × [0, 1.5]
x3
√
x+2
(8.b) g(x) = cos(x + 2x + 1) y h(x) = log
2 en el cuadro [1, 2] × [−1.5, 2]
x3
sen(πx)
(8.c) g(x) = √ y h(x) = 5e−2x ln(x2 + 1) en el cuadro [0, 2] × [−0.5, 1]
1 + 2x
(x + 10)(x − 10)(18 − x)
(8.d) f (x) = ln(x2 − 27) y g(x) = en el cuadro [−15, 20] × [−20, 15]
100
9. Resolver, si es posible, los siguientes sistemas lineales, comprobando que la solución es correcta:
2x1 + x2 − x3 = −1 x + 2y + 3z = 2
(9.a) 2x1 − x2 + 3x3 = −2 (S.C.D.) (9.d) 2x + y + 3z = 1 (S.C.I.)
3x1 − 2x2 = 1 x + y + 2z = 1
2x + 2y + t = 1
2x + 3y = 8
2x − 2y + z = −2
(9.b) 3x − y = −2 (S.C.D.) (9.e) (S.I.)
x − z + t = 0
−3x + y + z = 0
−4x + 4y − 2z = 1
x1 + x2 + x3 = 1 2x + 3y = 5
(9.c) 2x1 − x2 + x3 = 2 (S.C.I.) (9.f) x−y = 2 (S.I.)
x1 − 2x2 = 1 3x + y = 6
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.
./0%!)
Figura 2.1: Estructura condicional simple: Si Figura 2.2: Estructura condicional doble: Si
expresión toma el valor lógico true, se ejecu- expresión toma el valor lógico true, se ejecutan
tan las instrucciones del bloque instrucciones las instrucciones del bloque instrucciones 1.
y, después, el programa se continúa ejecutando Si, por el contrario, la expresión toma el
por la instrucción siguiente. Si, por el contrario, valor lógico false, se ejecuta el bloque de
la expresión toma el valor lógico false, no se instrucciones 2. En ambos casos, el programa
ejecutan las instrucciones, y el programa con- continúa ejecutándose por la instrucción siguien-
tinúa directamente por la instrucción siguiente. te.
En casi todos los lenguajes de programación, este tipo de estructuras se implementan mediante una ins-
trucción (o super-instrucción) denominada if, cuya sintaxis puede variar ligeramente de unos lenguajes
a otros. En MATLAB, concretamente, su forma es la siguiente:
27
2. Programación con MATLAB 28
instruccion-anterior
if expresion
instruccion-anterior bloque-de-instrucciones-1
if expresion else
bloque-de-instrucciones bloque-de-instrucciones-2
end end
instruccion-siguiente instruccion-siguiente
Observación. La orden de MATLAB sort([x,y]) tiene el mismo efecto que esta M-función.
Estas estructuras se pueden “anidar”, es decir, se puede incluir una estructura condicional dentro de
uno de los bloques de instrucciones de uno de los casos de otra.
En ocasiones interesa finalizar la ejecución de un programa en algún punto “intermedio” del mismo,
es decir, no necesariamente despues de la última instrucción que aparece en el código fuente. Esto se
hace mediante la orden
return
if (x > 0)
if (y > 0)
n=1; Entrada: x, y
else Salida: n
n=4;
end x*y == 0
true
n = 0 fin
else
false
if (y > 0)
n=2;
else true false true true
fin
Se pueden construir estructuras condicionales más complejas (con más casos). En MATLAB, estas
estructuras se implementan mediante la versión más completa de la instrucción if:
instruccion-anterior
true if expresion-1
expresión-1 instrucciones-1
bloque-de-instrucciones-1
false
elseif expresion-2
expresión-2
true
instrucciones-2 bloque-de-instrucciones-2
else
false
bloque-de-instrucciones-3
instrucciones-3 end
instruccion-siguiente
Este mecanismo de programación permite repetir un grupo de instrucciones mientras que se verifique
una cierta condición. Su diagrama de flujo y su sintaxis en MATLAB son los siguientes:
false instruccion-anterior
expresión while expresion
true
bloque-de-instrucciones
end
instrucciones instruccion-siguiente
en un bucle infinito. Llegado este caso, se puede detener el proceso pulsando la combinación de teclas
CTRL + C.
false
k <= n
true
suma = suma + k
k = k + 1
fin
Ejercicio.
Partiendo de la M-función SumaNat del ejemplo anterior, escribe una M-función que calcule y
devuelva la suma de todos los números naturales impares entre 1 y n.
%---------------------------------------------------------------
% 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);
disp(' ')
disp(' Tienes que acertar un numero del 1 al 9')
disp(' Pulsa CTRL + C si te rindes ...')
disp(' ')
M = 0;
while M~=n
M=input('Teclea un numero del 1 al 9 : ');
end
beep
disp(' ')
disp(' ¡¡¡¡Acertaste!!!! ')
Ejercicio.
Partiendo del script Adivina del ejemplo anterior, escribe una M-función que reciba como argu-
mento 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.
k = 1
false
instruccion-anterior
k <= n for k = 1 : n
bloque-de-instrucciones
true
end
instrucciones instruccion-siguiente
k = k + 1
Estructura de repetición for: Para cada valor de la variable k desde 1 hasta n, se ejecuta una
vez el bloque-de-instrucciones. Cuando se termina, el programa se continúa ejecutando por la
instrucción-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.
Las órdenes respectivas en MATLAB son (tanto en un bucle while como en un bucle for):
break
continue
instrucciones-1
true
instrucciones
condición !"#$%
para break
romper true condición
continue para
false romper
false
instrucciones-2
instrucciones
warning('Mensaje')
error('Mensaje')
imprime ??? Mensaje en la pantalla y detiene en ese punto la ejecución del programa.
Además de estas funciones y sus versiones más completas, MATLAB dispone de otras órdenes para
gestionar la detección de errores. Véase Error Handling en la documentación.
La instrucción input permite almacenar en una variable un dato que se introduce a través del teclado.
La orden
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 return. Lo que se teclea puede ser cualquier expresión que use constantes y/o
variables existentes en el Workspace. El resultado de esta expresión se almacenará en la variable var.
Si se pulsa la tecla return sin teclear nada se obtendrá una matriz vacía: [ ].
La instrucción disp permite imprimir en la pantalla el valor de una (matriz) constante o variable, sin
imprimir el nombre de la variable o ans y sin dejar líneas en blanco. Su utilidad es muy limitada.
disp(algo)
>> disp(pi)
3.1416
>> x = pi;
>> disp(['El valor de x es: ',num2str(x)])
El valor de x es: 3.1416
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 combinaciones de los siguientes elementos:
Códigos de conversión: formados por el símbolo %, una letra (como f, e, i, s) y eventual-
mente unos números para indicar el número de espacios que ocupará el dato a imprimir.
Texto literal a imprimir
Caracteres de escape, como \n.
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.
• Las instrucciones for, if, while, end, return, ... no necesitan llevar punto y coma al final, ya
que
1. no
Loproducen
mínimo salida
que sepor pantalla.
le pide Sí necesitan
a un programa llevarlo,
para en general,
llamarse las no
así es que instrucciones
produzcaincluidas
errores de ejecución.
dentro.
• Las
2. M-funciones no necesitan
A los argumentos llevar
de salida end !"#$%&'#
de una al final. De hecho, ahay
siempre nivel
quededarles
este curso,
algúnque lo lleven
valor.
puede inducir a error. Mejor no ponerlo.
3. Hay que comprobar los programas que se hacen, dándoles valores a los argumentos
• Todospara verificar todos
los programas losincluir
deben casosunas
que líneas
podrían darse. Esto forma
de comentarios parte deldeproceso
que expliquen forma de
sucinta
diseño y escritura del programa.
su cometido y forma de utilización. En las M-funciones, estas líneas deben ser las siguientes a la
instrucción function, ya que de esta manera son utilizadas por MATLAB como el texto de help
4. Las instrucciones !'(, &!, )*&+,, ,#-, (,%"(# ... no necesitan llevar punto y coma al
de la final,
función. En los scripts deben ser las primeras líneas del archivo.
ya que no producen salida por pantalla. Necesitan llevarlo, en general, las
instrucciones incluidas dentro. Ejemplo:
• Es bueno, en la medida de lo posible y razonable, respetar la notación habitual de la formulación
de un problema al elegir los nombres de las variables. Ejemplos
!"#$%&'##############(#)*+,-#.#/-0$#!++1/12$3!-#
###/4$5&'#
• Llamar S a una1+6'#################(#)*+,-#.#/-0$#!++1/12$3!-#
suma y P a un producto, mejor que G o N.
• Llamar i, j, k, l, m, n a un número entero, mejor que x ó y.
5. Las M-funciones no necesitan llevar ,#- al final. De hecho, a nivel de este curso, que lo
• Llamar T o temp
lleven puede a una
inducir temperatura,
a error. mejor
Mejor no que P.
ponerlo.
• Llamar e, err ó error a un erro, mejor que vaca.
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 funciona quetodos
bien en no loslos
lleven y/o que
casos que no sean adecuados.
puedan darse. Esto forma parte del proceso de
diseño y escritura del programa.
7. Los espacios 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
Comparay ellos
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:
a. Llamar S a una suma, mejor que G
b. Llamar P a un producto, mejor que N
c. Llamar T ó temp a una temperatura, mejor que P
Cálculo Numérico d. Llamar
II - Grado en iMatemáticas
j k l m n a un número entero, como un índice, mejor
Dpto. EDAN que x ó dey Sevilla
- Universidad
e. Llamar e ó err ó error a un error, en lugar de vaca
2. Programación con MATLAB 39
2.8 Ejercicios
1. Escribir una M-función function [N] = Multiplo(k, h) que reciba como argumentos de en-
trada dos números enteros k y h y devuelva:
• N = 1 si k + h es múltiplo de 2,
• N = 2 si, además, es múltiplo de 3,
• N = 0 en otro caso.
2. Escribir una M-función de nombre function [flag] = Divide(m, k, h) que reciba como ar-
gumento de entrada tres números m, k y h, y proporcione como salida un entero flag definido
como sigue:
• N = 1 si el punto (x, y) está (estrictamente) dentro del círculo de centro (1, 1) y radio r = 1,
• N = 2 si, además, está (estrictamente) dentro del círculo de centro (0, 0) y radio r = 1,
• N = 0 en otro caso.
4. Escribir una M-función de nombre function [A] = Areatri(a, b, c) que calcule el área de
un triángulo a partir de las longitudes de sus lados (fórmula de Herón):
a+b+c
p(p − a)(p − b)(p − c), con p =
p
A= .
2
El programa debe emitir un error con un mensaje explicativo en los casos eventuales en que no
se pueda calcular el área: (a) Si alguna de las longitudes es menor que cero; (b) Si el radicando
el negativo, lo cual indica que no existe ningún triángulo que tenga esos lados.
si x < 0,
−x
f (x) =
ex − 1 si x ≥ 0.
−x + 2 si x < 2,
f (x) =
x − 2 si x ≥ 2.
7. Calcula la suma de los 100 primeros números pares (empezando en el 2), utilizando (a) un bucle
for y (b) un bucle while.
log(k + 1) + 2k
8. Calcula e imprime en pantalla los 10 primeros términos de la sucesión xk =
2k 2 − 1
9. Construye un vector con los diez primeros términos de la sucesión
x1 = 0.25
2x − 1
xn+1 = n , n≥1
x2n + 1
11. Escribe una M-función function [u] = ProdAb(A, b, n) que reciba como datos de entrada
una matriz cuadrada A de dimensión n, y un vector b de longitud n, y devuelva el producto
A b calculado elemento a elemento mediante bucles, es decir, sin usar el producto matricial de
MATLAB.
12. Escribe una M-función de nombre function [H] = fhilbert(n) que reciba como argumento
de entrada un entero positivo n y devuelva la matriz de Hilbert de dimensión n, es decir, la
matriz definida por
n 1
H = hij i,j=1 , hij = .
i+j−1
Estas matrices son un ejemplo notable de mal condicionamiento, incluso para dimensión pequeña.
13. Escribe una M-función function [w] = SuperaMedia(v) que reciba como argumento de en-
trada 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. Por
ejemplo, si v=(10,1,7), entonces w=(10,7); si v=(6,5,6,7), entonces w=(6,6,7).
14. Escribir una M-función function [v]=Inter(x,y) que reciba como argumento de entrada dos
vectores fila x e y, y devuelva como salida otro vector v que contenga los elementos comunes a x
e y, es decir, la intersección de ambos conjuntos. Se supondrá que ninguno de los vectores x e y
tiene elementos repetidos, cosa que no es necesario verificar.
15. Escribir una M-función function [x] = Bajada(A, b) para calcular la solución x del sistema
Ax = b en el caso en que A es una matriz cuadrada triangular inferior utilizando el algoritmo
de bajada:.
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).
16. Escribir una M-función function [x] = Subida(A, b) para calcular la solución x del sistema
Ax = b en el caso en que A es una matriz cuadrada triangular superior utilizando el algoritmo
de subida:.
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.
17. La fórmula de los rectángulos es una fórmula de integración numérica para aproximar el valor
de una integral definida. Dada una partición a = x1 < x2 . . . < xn = b del intervalo [a, b] tal
que todos los subintervalos [xi , xi+1 ] tienen la misma amplitud, h = (b − a)/(n − 1), la integral
definida de f entre a y b se puede aproximar por:
Z b n−1
X
f (x) ≈ h f (xi )
a i=1
Z 1
2
Utilizando esta fórmula con 15 subintervalos, aproximar el valor de la integral e−x dx.
0
18. La fórmula de los puntos medios es una fórmula de integración numérica para aproximar el
valor de una integral definida. Dada una partición a = x1 < x2 . . . < xn+1 = b del intervalo
[a, b] tal que todos los subintervalos [xi , xi+1 ] tienen la misma amplitud, h = (b − a)/n, la integral
definida de f entre a y b se puede aproximar por:
b n
xi + xi+1
Z X
f (x) ≈ h f
a 2
i=1
19. La fórmula de los trapecios es una fórmula de integración numérica para aproximar el valor
de una integral definida. Dada una partición a = x1 < x2 . . . < xn+1 = b del intervalo [a, b] tal
que todos los subintervalos [xi , xi+1 ] tienen la misma amplitud, h = (b − a)/h, la integral definida
de f entre a y b se puede aproximar por:
b n
hh
Z X i
f (x) ≈ f (x1 ) + 2 f (xi ) + f (xn+1 )
a 2
i=2
20. La fórmula de Simpson es una fórmula de integración numérica para aproximar el valor de
una integral definida. Dada una partición a = x1 < x2 . . . < xn+1 = b del intervalo [a, b] tal que
todos los subintervalos [xi , xi+1 ] tienen la misma amplitud, h = (b − a)/n, la integral definida de
f entre a y b se puede aproximar por:
b n n−1
hh xi + xi+1
Z X X i
f (x) ≈ f (x1 ) + 2 f (xi ) + 4 f + f (xn+1 )
a 6 2
i=2 i=2
Escribir una M-función function [v] = FSimpson(fcn,a,b,n) que calcule el valor de la inte-
gral de fcn en [a, b] utilizando la fórmula de Simpson con n subintervalos.
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.
α α x2
a0 x0 b0 a1 x1 b1 a2 b2
43
3. Resolución de ecuaciones no lineales 44
Algoritmo de dicotomía
a+b
c) Hacer x = , e = e/2 y n = n + 1 y volver al paso b).
2
Ejercicio 3.1 (Método de dicotomía) Escribir una M-función que aproxime la solución de
f (x) = 0 en el intervalo [a, b] utilizando el método de dicotomía.
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;
fa = fx;
else
b = x;
fb = fx;
end
x = (a+b)*0.5;
e = e*0.5;
end
Ejercicio 3.2 Modificar el programa Dicotomia de forma que devuelva, además de la aproxima-
ción de la solución, el número de iteraciones realizadas:
Ejercicio 3.3 Modificar el programa Dicotomia 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.
Este método tiene más en cuenta las propiedades de la función f , es capaz de detectar la cercanía del
cero a uno de los extremos del intervalo y, en definitiva, cuando el intervalo es pequeño, la secante
por los puntos (a, f (a)) y (b, f (b)) aproxima a la derivada de f . No es extraño, por lo tanto, que se
obtengan mejores resultados que con el método de dicotomía.
b.1) Si |b − a| ≤ ε o bien |f (x)| < tol, parar y devolver x como aproximación de la solución.
b.2) Si f (a)f (x) < 0, hacer b = x.
b.3) Si no, hacer a = x
f (b)a − f (a)b
c) Hacer x = y volver al paso b).
f (b) − f (a)
Ejercicio 3.4 (Método de Regula Falsi) Escribir una M-función que aproxime la solución de
f (x) = 0 en el intervalo [a, b] utilizando el método de la falsa posición.
fa = fun(a);
fb = fun(b);
if ( fa*fb > 0 )
error('La funcion debe cambiar de signo en el intervalo');
end
x = (fb*a-fa*b)/(fb-fa);
Ejercicio 3.5 Modificar el programa RegulaFalsi de forma que devuelva, además de la aproxi-
mación de la solución, el número de iteraciones realizadas:
Ejercicio 3.6 Modificar el programa RegulaFalsi de forma que imprima una tabla con la evo-
lució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.
Ejercicio 3.7 Modificar el programa RegulaFalsi de forma que para el test de parada se utilice
el error relativo en lugar del error relativo.
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 α b a x1 x3 α x2 x0 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 3.8 (Método de aproximaciones sucesivas) Escribir una M-función que aproxime
la solución de x = g(x) en el intervalo [a, b] 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
%
Ejercicio 3.9 Modificar el programa AASS de forma que devuelva, además de la aproximación de
la solución, el número de iteraciones realizadas:
Ejercicio 3.10 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.
Ejercicio 3.11 Modificar el programa AASS de forma que para el test de parada se utilice el error
relativo en lugar del error absoluto.
Ejercicio 3.12 Modificar el programa AASS de forma que se produzca un mensaje de warning
cuando se alcance el número máximo de iteraciones.
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 − 0
f (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 .
|xn − xn−1 |
b.1) Si ≤ ε, n ≥ 1, parar y devolver xn como aproximación.
|xn |
f (xn )
b.2) Hacer xn+1 = xn − 0 y repetir el paso b).
f (xn )
Ejercicio 3.13 Escribir una M-función que aproxime la solución de f (x) = 0 utilizando el método
de Newton.
Ejercicio 3.14 Modificar el programa Newton de forma que devuelva, además de la aproximación
de la solución, el número de iteraciones realizadas:
Ejercicio 3.15 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.
Ejercicio 3.16 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 |
Ejercicio 3.17 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.
solucion=fzero(funcion,xcero)
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.
xcero es un valor «cercano» a la solución, a partir del cual el algoritmo iterado de búsqueda de la
solución comenzará a trabajar.
solucion es el valor (aproximado) de la solución encontrado por el algoritmo.
Ejercicio 3.18 La ecuación siguiente tiene una solución cerca de x = 1. Calcular una aproxima-
ción. x
x + ln =0
3
1. Definimos una función anónima que evalúe la expresión del primer miembro y usamos el
comando fzero tomando x=1 como valor inicial:
2. Probamos ahora eligiendo un punto inicial más alejado de la solución, por ejemplo:
fzero(fun,15)
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. En la ventana de comandos, para calcular el cero de la función escribimos:
fzero(@mifuncion,1)
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 ejemplo siguiente:
x2 = 0
La determinación teórica de un punto inicial cercano a la solución puede ser una tarea difícil. Si ello
es posible, un estudio gráfico previo puede resultar de gran ayuda.
Ejercicio 3.20 Calcular, si existe, una solución positiva de la ecuación sen(x)−2 cos(2x) = 2−x2
determinando un punto inicial a partir de la gráfica de la función.
La gráfica mostrará que esta función tiene dos ceros: uno positivo cerca de x = 1 y otro
negativo cerca de x = −1.
fzero(fun,1)
ans =
0.8924
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.
Comprobaremos que hay una solución en el intervalo [−1, 1], otra en [2, 4] y otra en [6, 7].
fzero(fun,1) % sol. x = 0
fzero(fun,2) % sol. x = 2.4674
fzero(fun,6) % sol. x = 6.2832
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.
También será necesario un análisis teórico para no caer en equivocaciones.
2. Define una función anónima para f y dibuja su gráfica en el intervalo [−1, 20] (por ejemplo)
Resulta obvio, a la vista de la gráfica y del análisis del apartado anterior, que no hay ningún
cero de f a la derecha de x = 10. Sin embargo, entre x = −1 y x = 20 la situación resulta
confusa, ya que la escala del eje OY es demasiado pequeña: observa que la parte visible es
[−1, 20] × [−10000, 2000]. Vuelve a dibujar la gráfica, pero ahora en el intervalo (por ejemplo)
[0, 10].
3. Con esta nueva gráfica parece claro que hay un cero en el intervalo [0, 2] y otro en el intervalo
[6, 8]. Tratamos de aproximarlos con fzero.
x1 = fzero(fun,0); % x1 = 0.7815
x2 = fzero(fun,6); % x2 = 7.1686
Ejercicio 3.24 Las frecuencias naturales de la vibración de una viga homogénea sujeta por un
extremo son las soluciones de :
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 de un punto suficientemente próximo.
Ejercicio 3.25 Comparar los resultados que se obtienen al calcular (numéricamente) la solución
de la ecuación
f (x) = cos(x) cosh(x) + 1 = 0
en el intervalo [1, 2] utilizando la función AASS del Ejercicio 3.6 y la función fzero.
integral(fun,a,b);
La función integral utiliza una fórmula de integración adaptativa (i.e. una fórmula en la que los puntos
del soporte de integración son elegidos de acuerdo con ciertos criterios con el objetivo de minimizar el
error) y utiliza parámetros de tolerancia estándar para determinar la precisión (aunque sus valores se
pueden elegir utilizando ciertos argumentos opcionales).
En el caso escalar, la función a integrar (anónima o M-función) debe ser escrita de forma vectorizada,
esto es, que admita como argumento de entrada un vector y devuelva un vector de la misma dimensión
con los valores de la función.
Z 3
Ejercicio 4.1 Calcular la integral definida x sen(4 ln(x)) dx
0.2
1
En realidad, la función integral permite también calcular (aproximar) integrales impropias y ciertas integrales de
línea de funciones de variable compleja, aunque en estas prácticas nos limitamos a ejemplos con funciones reales de una
variable real.
57
4. Integración numérica 58
Z 8
0.8
Ejercicio 4.2 Calcular la integral definida (x e−x + 0.2) dx describiendo el integrando
0
mediante una M-función.
Escribe una M-función de nombre (por ejemplo) mifun, en un fichero de nombre mifun.m:
Después, para calcular la integral, habría que escribir (en la ventana de comandos o en otro pro-
grama)
integral(@mifun, 0, 8)
Comparar con el valor exacto de la integral (imprimir los números con todos sus decimales).
1
(Una primitiva de f (x) = arc tg(x + 4) es F (x) = − ln(x2 + 8x + 17) + (x + 4) arc tg(x + 4)).
2
Ejercicio 4.4 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
1 2
A = integral(f,1,2)
A =
0.7166
f = @(x) sin(4*log(x+1)); y
x1 = linspace(0,10);
y1=f(x1);
plot(x1,y1) A1
A3
grid on x1
x2 9 x
A2
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 4.6 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):
x = linspace(-5,5);
yf = f(x);
yg = g(x);
plot(x,yf,x,yg)
(x 2 , y 2 )
grid on
(x 1 , y 1 )
2. La observación de la gráfica nos muestra que ambas curvas se cortan en dos puntos, (x1 , y1 )
y (x2 , y2 ). Para calcularlos comenzamos por calcular sus abscisas x1 y x2 :
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
Ejercicio 4.7 Calcular el área de la región que queda encerrada entre las curvas
√
y = f (t) = t sen( 5t) e y = g(t) = 0.5 − t.
trapz(x,y);
x,y son dos vectores de la misma dimensión y representan las coordenadas de los puntos. Obviamente,
los puntos deben estar ordenados en orden creciente según sus abscisas.
Ejercicio 4.8 Calcular la integral definida entre 0 y 15 de la función discreta dada por el siguiente
conjunto de datos:
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
Ejercicio 4.9 Aproximar la integral del Ejercicio 4.3 mediante la fórmula de los trapecios (función
trapz), utilizando varios soportes, con distinto número de puntos ( 5, 15, 30, 100, 1000). Comparar
con el resultado obtenido con integral. Imprimir los números con todos sus decimales.
Ejercicio 4.10 Se define el promedio integral de una función f en un intervalo [a, b] como:
b
1
Z
f¯ = f (x) dx
b−a a
Calcular el promedio integral de la función discreta definida por los datos almacenados en el fichero
discreta1.txt.
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. En Matemáticas también es necesario en ocasiones,
sustituir o aproximar una función dada por una más simple o con alguna estructura especial, de forma
que ambas tomen los mismos valores en un cierto conjunto de puntos.
En algunas ocasiones, lo anterior no es posible o el resultado no posee las propiedades deseadas. En
estos casos se puede recurrir al ajuste
que verifique
p(xi ) = yi ∀i = 1, 2, . . . n
Para calcular los coeficientes de este polinomio con MATLAB se usa la orden
coef = polyfit(x,y,n-1)
La orden polyfit calcula los coeficientes del polinomio de interpolación resolviendo el sistema de
Vandermonde.
Una vez calculados los coeficientes, se puede evaluar el polinomio en un punto z mediante la orden
pz = polyval(coef,z)
63
5. Interpolación y ajuste polinómico 64
coef es un vector que contiene los coeficientes de un polinomio en orden de mayor a menor potencia.
pz la salida de esta función es el valor del polinomio en el punto z (si z es un vector, pz será un
vector de la misma dimensión que z).
Ejercicio 5.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. (Escribir un script con las órdenes siguientes).
%
% Calcular y representar graficamente el polinomio de interpolacion
% global que pasa por tres puntos dados
%
x = [1, 2, 3]
y = [-3, 1, 3];
coef = polyfit(x, y, 2);
z = linspace(0,7);
pz = polyval(coef, z);
plot(z, pz, 'LineWidth',1.1)
hold on
plot([-1,7],[0,0], 'k','LineWidth',1.1)
plot([0,0],[-10,6],'k','LineWidth',1.1)
plot(x,y,'r.','MarkerSize',20)
axis([-2, 8, -11, 7])
hold off
Ejercicio 5.2 Este ejercicio 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 (fenómeno de Runge).
Dibuja su gráfica, así como los puntos con marcadores, y observa las inestabilidades cerca de los
extremos. Escribe un script o M-función que lleve a cabo todo lo que se pide.
y3
y2
y8
pz
y1 y5 y6 y7
y4
x1 x2 x3 z x4 x5 x6 x7 x8
pz = interp1(x,y,z)
pz la salida de esta función es el valor del polinomio en el punto z (si z es un vector, pz será un
vector de la misma dimensión que z).
Ejercicio 5.3 El fichero de datos censo.txt contiene dos columnas que corresponden al censo de
EEUU entre los años 1900 y 1990 (en millones de personas).
Representar gráficamente la evolución del censo en esos años y estimar la población que había en
el año 1956. (Escribir un script con las órdenes siguientes).
%
% Interpolacion lineal a trozos de un conjunto de datos
% del censo de los EEUU durante el siglo XX
%
Para ello es necesario, naturalmente, aumentar los grados de libertad de que disponemos, es decir, el
número de coeficientes que podemos elegir, i.e. el grado de los polinomios.
Concretamente, mediante splines cúbicos, se pueden obtener funciones interpolantes de clase C 2 .
El problema ahora es:
s(xi ) = yi , i = 1, 2, . . . n
y tal que, restringida a cada subintervalo [xi , xi+1 ], sea un polinomio de grado tres.
Para calcular splines cúbicos, MATLAB dispone de la función spline, que se puede usar de dos formas
distintas.
pz = spline(x, y, z)
pz la salida de esta función es el valor del spline en el punto z (si z es un vector, pz será un vector
de la misma dimensión que z).
coef = spline(x, y)
coef es una estructura de datos que contiene los coeficientes que definen el spline (no es necesario,
a los efectos de este curso, saber más sobre esta estructura)
Una vez calculados los coeficientes del spline con esta función, para evaluarlo en un punto z hay que
usar la orden
pz = ppval(coef, z)
pz la salida de esta función es el valor del spline en el punto z (si z es un vector, pz será un vector
de la misma dimensión que z).
Ejercicio 5.6 El fichero DatosSpline.dat contiene una matriz con dos columnas, que correspon-
den a las abscisas y las ordenadas de una serie de datos.
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.
(Escribir un script con las órdenes siguientes).
% ------------- graficas
plot(z,p,'b--','LineWidth',1.1);
hold on
plot(z,s,'g', 'LineWidth',1.1);
plot(x,y,'r.', 'MarkerSize',15)
legend('Interpolante global', 'Spline cubico', 'Puntos del soporte')
hold off
shg
Ejercicio 5.7 En este ejemplo se muestra el uso de la función spline en la segunda de las formas
que se ha explicado antes: en una primera etapa se calculan los coeficientes del spline y se almacenan
en una variable, y en una segunda etapa se evalúa el spline en los puntos deseados, utilizando la
función ppval. Esto permite no repetir el cálculo de los coeficientes (que siempre son los mismos)
cada vez que se desea evaluar el spline.
Con los datos del mismo fichero DatosSpline.dat se pide:
a) Definir una función anónima que represente el spline cúbico s(x) que interpola dichos valores,
es decir, que calcule s(x) para cualquier x.
b) Calcular
Z 40
V = s(x) dx
0
(Escribir un script con las órdenes siguientes).
Ejercicio 5.8 Se consideran los mismos valores que en el ejercicio 5.2 y el ejercicio 5.4.
Representa gráficamente (juntos) el polinomio de interpolación global, el interpolante lineal a
trozos y el spline cúbico. Representa también los puntos del soporte de interpolación, mediante
marcadores. Añade las leyendas adecuadas para que se pueda identificar cada curva adecuadamente.
Calcula el valor interpolado para z = 1 por cada uno de los procedimientos.
Escribe un script o M-función que lleve a cabo todo lo que se pide.
Ejercicio 5.9 (Prescindible. 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.
En este ejercicio 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, utilizando para ello un soporte regular con 8 puntos. Hay por lo tanto que
dibujar tres curvas en [0, 10]:
1. La curva y = sen(x).
La técnica de interpolación que hemos explicado antes requiere que la función que interpola los datos
pase exactamente por los mismos. En ocasiones esto no da resultados muy satisfactorios, por ejemplo
si se trata de muchos datos. También sucede con frecuencia que los datos vienen afectados de algún
error, por ejemplo porque provienen de mediciones. No tiene mucho sentido, pues, obligar a la función
que se quiere construir a «pasar» por unos puntos que ya de por sí no son exactos.
Otro enfoque diferente es construir una función que no toma exactamente los valores dados, sino que
«se les parece» lo más posible, por ejemplo minimizando el error, medido éste de alguna manera.
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. Veremos solamente cómo se
puede hacer esto con MATLAB en algunos casos sencillos.
La función polyfit usada ya para calcular el polinomio de interpolación global sirve también para
ajustar unos datos por un polinomio de grado dado:
coef = 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 5.10 Calcula y dibuja (en el intervalo [0.5, 10]) los polinomios de ajuste de grado 1, 2,
3 y 6 para los siguientes datos:
% ------------ graficas
z = linspace(0.5, 10);
plot(z, p1(z), z, p2(z), z, p3(z), z, p6(z))
hold on
plot(x, y, 'm.', 'MarkerSize', 15)
legend('Recta de regresion', 'Parábola de regresion', ...
'Cubica de regresion','Pol. interpolacion global',...
'Location', 'SouthEast')
hold off
shg
Ejercicio 5.11 En el fichero Finanzas.txt está recogido el precio de una determinada 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 recta 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 mediante una linea poligonal, para observar las oscilaciones
de las cotizaciones. Calcula y representa gráficamente la recta de regresión para estos datos.
% ------- graficas
plot(t,c,'r.')
hold on
plot(t, fc(t), 'LineWidth', 1.1)
hold off
shg
function FunInterp(f, a, b, n)
y 0 = f (t, y), f : Ω ⊂ R2 → R
admite, en general, infinitas soluciones. Si, por ejemplo, f ∈ C 1 (Ω; R)1 , por cada punto (t0 , y0 ) ∈ Ω
pasa una única solución ϕ : I → R, definida en un cierto intervalo I ⊂ R, que contiene a t0 .
Se denomina Problema de Cauchy (PC) o Problema de Valor Inicial (PVI)
(
y 0 = f (t, y),
(6.1)
y(t0 ) = y0
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 ]
(6.2)
y(t0 ) = y0
yk ≈ y(tk ), k = 0, 1, . . . , n,
yk ≈ y(tk ), k = 0, 1, . . .
1
En realidad basta con que f sea continua y localmente Lipschitziana con respecto de y en Ω.
75
6. Resolución de ecuaciones y sistemas diferenciales ordinarios 76
MATLAB dispone de toda una familia de funciones para resolver (numéricamente) ecuaciones diferen-
ciales:
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.
y = 5y en [0, 1]
0
y(0) = 1
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)
150
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).
Si se utiliza la función ode45 en esta forma, no se obtiene ninguna gráfica.
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 6.3 Calcular (aproximaciones de) los valores de la solución del problema
0
y = −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: