1.1: Nueva página

Última actualización
Guardar como PDF

Page ID: 118775

$ \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } $ $ \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} $$\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$ $\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$$\newcommand{\AA}{\unicode[.8,0]{x212B}}$

Modelado y Simulación

Este libro trata sobre modelar y simular sistemas físicos. Antes de que podamos construir cualquier modelo, ayudará tener una comprensión de alto nivel de lo que es un modelo. También necesitaremos familiarizarnos con las herramientas que utilizamos para construirlas. En este capítulo, veremos el proceso de modelado e introduciremos MATLAB, el lenguaje de programación que usaremos para representar modelos y ejecutar. Al final del capítulo encontrarás ejercicios que puedes utilizar para poner a prueba tus conocimientos.

Modelado

Cuando digo “modelar”, estoy hablando de algo como Figura [fig:modeling]. En la esquina inferior izquierda de la figura está el sistema, algo en el mundo real que nos interesa. A menudo, es algo complicado, así que tenemos que decidir qué detalles se pueden dejar fuera; eliminar detalles se llama abstracción.

El proceso de modelado

El resultado de la abstracción es un modelo, mostrado en la parte superior izquierda; un modelo es una descripción del sistema que incluye solo las características que creemos que son esenciales. Un modelo se puede representar en forma de diagramas y ecuaciones, que pueden ser utilizados para el análisis matemático. También se puede implementar en forma de un programa de computadora, que puede ejecutar simulaciones.

El resultado del análisis y la simulación podría ser una predicción sobre lo que hará el sistema, una explicación de por qué se comporta de la manera en que lo hace, o un diseño específico diseñado para satisfacer un requisito u optimizar el rendimiento.

Podemos validar predicciones y diseños de pruebas tomando medidas del mundo real y comparando los datos que obtenemos con los resultados de análisis y simulación.

Para cualquier sistema físico hay muchos modelos posibles, cada uno incluyendo y excluyendo diferentes características o incluyendo diferentes niveles de detalle. El objetivo del proceso de modelado es encontrar el modelo que mejor se adapte a su propósito (predicción, explicación o diseño).

A veces el mejor modelo es el más detallado. Si incluimos más características, el modelo es más realista y esperamos que sus predicciones sean más precisas. Pero muchas veces un modelo más simple es mejor. Si incluimos solo las características esenciales y dejamos fuera el resto, obtenemos modelos con los que es más fácil trabajar, y las explicaciones que brindan pueden ser más claras y convincentes.

Como ejemplo, supongamos que alguien te preguntó por qué la órbita de la Tierra es casi elíptica. Si modelas la Tierra y el Sol como masas puntuales (ignorando su tamaño real), calculas la fuerza gravitacional entre ellos usando la ley de gravitación universal de Newton, y calculas la órbita resultante usando las leyes de movimiento de Newton, puedes demostrar que el resultado es una elipse.

Por supuesto, la órbita real de la Tierra no es una elipse perfecta, debido a las fuerzas gravitacionales de la Luna, Júpiter y otros objetos en el sistema solar, y porque las leyes del movimiento de Newton son solo aproximadamente ciertas (no toman en cuenta los efectos relativistas).

Pero agregar estas características al modelo no mejoraría la explicación; más detalle sólo sería una distracción de la causa fundamental. No obstante, si el objetivo es predecir la posición de la Tierra con gran precisión, podría ser necesario incluir más detalles.

Elegir el mejor modelo depende de para qué sirve el modelo. Por lo general, es una buena idea comenzar con un modelo simple, aunque sea probable que sea demasiado simple, y probar si es lo suficientemente bueno para su propósito. Entonces puedes agregar características poco a poco, comenzando por las que esperas que sean más esenciales. Este proceso se llama modelado iterativo.

La comparación de los resultados de modelos sucesivos proporciona una forma de validación interna para que pueda detectar errores conceptuales, matemáticos y de software. Y al agregar y eliminar características, puede saber cuáles tienen el mayor efecto en los resultados, y cuáles pueden ignorarse. La comparación de resultados con datos del mundo real proporciona validación externa, que generalmente es la prueba más fuerte.

La figura [fig:modeling] muestra que los modelos pueden ser utilizados tanto para el análisis como para la simulación; en este libro haremos algunos análisis, pero el foco está en la simulación. Y la herramienta que usaremos para construir simulaciones es MATLAB. Así que comencemos.

Una calculadora glorificada

MATLAB es un lenguaje de programación con características que soportan modelado y simulación. Tiene muchas características, por lo que puede parecer abrumador, pero en el fondo MATLAB es una calculadora glorificada. Cuando inicie MATLAB debería ver una ventana titulada MATLAB que contiene ventanas más pequeñas tituladas Carpeta actual, Ventana de comandos y Espacio de trabajo. En Octave, Carpeta Actual se llama Explorador de Archivos.

El Intérprete

La Ventana de Comandos ejecuta el intérprete, lo que le permite ingresar comandos; una vez ingresados, el intérprete ejecuta el comando e imprime el resultado. Inicialmente, la Ventana de Comandos contiene un mensaje de bienvenida con información sobre la versión del software que estás ejecutando, seguido de un prompt, que se ve así:

>>

El símbolo >> le pide que introduzca un comando. El tipo de comando más simple es una expresión matemática, como 2 + 1. Si escribe una expresión y luego presiona (o), MATLAB evalúa la expresión e imprime el resultado.

>> 2 + 1
ans = 3

Solo para que quede claro: en este ejemplo, MATLAB mostró >>; escribí 2 + 1 y luego presioné, y MATLAB mostró ans = 3.

En esta expresión, el signo más es un operador y los números 2 y 1 son operandos. Una expresión puede contener cualquier número de operadores y operandos. No hay que poner espacios entre ellos; algunas personas sí y otras no. He aquí un ejemplo sin espacios:

>> 1+2+3+4+5+6+7+8+9
ans = 45

Hablando de espacios, es posible que haya notado que MATLAB pone una línea en blanco entre ans = y el resultado. En mis ejemplos lo dejaré afuera para ahorrar espacio.

Los otros operadores aritméticos son más o menos lo que cabría esperar. La resta se denota con un signo menos (-), la multiplicación se designa con un asterisco (*), la división se denota con una barra hacia adelante (/).

>> 2*3 - 4/5
ans = 5.2000

Otro operador común es la exponenciación, que utiliza el símbolo ^, a veces llamado “intercalación” o “sombrero”. Entonces, 2 elevado al 16 poder es

>> 2^16
ans = 65536

El orden de las operaciones es lo que esperarías del álgebra básica: la exponenciación ocurre antes de la multiplicación y división, y la multiplicación y división ocurren antes de la suma y resta. Si desea anular el orden de las operaciones, puede usar paréntesis.

>> 2 * (3-4) / 5
ans = -0.4000

Cuando agregué los paréntesis, eliminé algunos espacios para hacer más clara la agrupación de operandos a un lector humano. Esta es la primera de muchas pautas de estilo que recomendaré para que tus programas sean más fáciles de leer. El estilo no cambia lo que hace el programa; el intérprete de MATLAB no comprueba el estilo. Pero los lectores humanos sí, y el humano más importante que leerá tu código eres tú.

Y eso nos lleva al Primer Teorema de Depuración:

El código legible es código depurable.

Vale la pena dedicar tiempo a hacer que tu código sea bonito; ¡te ahorrará tiempo!

Funciones matemáticas

MATLAB sabe cómo calcular casi todas las funciones matemáticas de las que ha oído hablar. Por ejemplo, conoce todas las funciones trigonométricas, así es como las usa:

>> sin(1)
ans = 0.8415

Este comando es un ejemplo de una llamada a una función. El nombre de la función es pecado, que es la abreviatura habitual para el seno trigonométrico. El valor entre paréntesis se llama argumento.

Las funciones trigonométricas pecado, cos y tan —entre muchas otras— funcionan en radianes, por lo que el argumento del ejemplo se interpreta como 1 radián. MATLAB también proporciona funciones trig que funcionan en grados: sind, cosd y tand.

Algunas funciones toman más de un argumento, en cuyo caso los argumentos están separados por comas. Por ejemplo, atan2 calcula la tangente inversa, que es el ángulo en radianes entre el eje x positivo y el punto con las coordenadas x e y dadas.

>> atan2(1,1)
ans = 0.7854

Si ese poco de trigonometría no te resulta familiar, no te preocupes por ello. Es sólo un ejemplo de una función con múltiples argumentos.

MATLAB también proporciona funciones exponenciales, como exp, que calcula$e$ elevado a la potencia dada. Entonces exp (1) es solo$e$:

>> exp(1)
ans = 2.7183

La inversa de exp es log, que calcula la base logarítmica$e$:

>> log(exp(3))
ans = 3

Este ejemplo también demuestra que las llamadas a funciones se pueden anidar; es decir, se puede usar el resultado de una función como argumento para otra.

De manera más general, puede usar una llamada a función como operando en un.

>> sqrt(sin(0.5)^2 + cos(0.5)^2)
ans = 1

Como probablemente adivinaste, sqrt calcula la raíz cuadrada.

Hay muchas otras funciones matemáticas, pero esto no pretende ser un manual de referencia. Para conocer otras funciones, debes leer la documentación.

Variables

Por supuesto, MATLAB es bueno para algo más que evaluar expresiones. Una de las características que hace que MATLAB sea más potente que una calculadora es la capacidad de darle un nombre a un valor. Un valor con nombre se llama variable.

MATLAB viene con algunas variables predefinidas. Por ejemplo, el nombre pi se refiere a la cantidad matemática$\pi$, que es aproximadamente esta:

>> pi
ans = 3.1416

Y si haces algo con números complejos, te puede resultar conveniente que tanto i como j estén predefinidos como la raíz cuadrada de$-1$.

Puede usar un nombre de variable en cualquier lugar donde pueda usar un número, por ejemplo, como operando en una expresión,

>> pi * 3^2
ans = 28.2743

o como argumento a una función:

>> sin(pi/2)
ans = 1

Siempre que evalúe una expresión, MATLAB asigna el resultado a una variable llamada ans. Puede usar ans en un cálculo posterior como taquigrafía para “el valor de la expresión anterior”.

>> 3^2 + 4^2
ans = 25

>> sqrt(ans)
ans = 5

Pero ten en cuenta que el valor de ans cambia cada vez que evalúas una expresión.

Declaraciones de Asignación

Puedes crear tus propias variables, y darles valores, con una declaración de asignación. El operador de asignación es el signo igual (=), usado así:

>> x = 6 * 7
x = 42

Este ejemplo crea una nueva variable llamada x y le asigna el valor de la expresión 6 * 7. MATLAB responde con el nombre de la variable y el valor calculado.

Existen algunas reglas a la hora de asignar un valor a las variables. En cada declaración de asignación, el lado izquierdo tiene que ser un nombre de variable legal. El lado derecho puede ser cualquier expresión, incluyendo llamadas a funciones. Casi cualquier secuencia de letras minúsculas y mayúsculas es un nombre de variable legal. Algunos signos de puntuación también son legales, pero el guión bajo (_) es la única no letra de uso común. Los números están bien, pero no al principio. No se permiten espacios. Los nombres de variables son sensibles a mayúsculas, por lo que x y X son variables diferentes.

Veamos algunos ejemplos de declaraciones de asignación.

>> fibonacci0 = 1;

>> LENGTH = 10;

>> first_name = 'bob'
first_name = 'bob'

Los dos primeros ejemplos demuestran el uso del punto y coma, que suprime la salida de un comando. En este caso MATLAB crea las variables y les asigna valores pero no muestra nada.

El tercer ejemplo demuestra que no todo en MATLAB es un número. Una secuencia de caracteres entre comillas simples es una cadena.

Aunque i, j y pi están predefinidos, usted es libre de reasignarlos. Es común usar i y j para otros fines, pero es raro asignar un valor diferente a pi.

Variables en el espacio de trabajo

Al crear una nueva variable, ésta aparece en la ventana Espacio de trabajo y se agrega al espacio de trabajo, que es un conjunto de variables y sus valores.

El comando who imprime los nombres de las variables en el espacio de trabajo:

>> x = 5;
>> y = 7;
>> z = 9;
>> who

Your variables are:

x  y  z

El comando clear elimina las variables especificadas del espacio de trabajo:

>> clear x
>> who

Your variables are:

y z

Pero ten cuidado: si no especificas ninguna variable, clear las elimina todas.

Para mostrar el valor de una variable, puede usar la función disp:

>> disp(z)
     9

pero es más fácil simplemente escribir el nombre de la variable:

>> z
z = 9

Ahora que ya has visto cómo usarlos, demos un paso atrás y pensemos por qué usaríamos variables.

¿Por qué Variables?

Hay una serie de razones para usar variables. Uno grande es evitar recalcular un valor que usas repetidamente. Por ejemplo, si su cálculo usa$e$ con frecuencia, es posible que desee calcularlo una vez y guardar el resultado.

>> e = exp(1)
e = 2.7183

Las variables también hacen más evidente la conexión entre el código y las matemáticas subyacentes. Si está calculando el área de un círculo, es posible que desee usar una variable llamada r:

>> r = 3
r = 3

>> area = pi * r^2
area = 28.2743

De esa manera, tu código se asemeja a la fórmula familiar$a = \pi r^2$.

También puede usar una variable para romper un cálculo largo en una secuencia de pasos. Supongamos que estás evaluando una expresión grande y peluda como esta:

p = ((x - theta) * sqrt(2 * pi) * sigma)^-1 * ...
exp(-1/2 * (log(x - theta) - zeta)^2 / sigma^2)

Puede usar puntos suspensivos para romper la expresión en varias líneas. Sólo tienes que entrar... al final de la primera línea y continuar a la siguiente. Pero a menudo es mejor dividir el cálculo en una secuencia de pasos y asignar resultados intermedios a las variables:

shiftx = x - theta
denom = shiftx * sqrt(2 * pi) * sigma
temp = (log(shiftx) - zeta) / sigma
exponent = -1/2 * temp^2
p = exp(exponent) / denom

Los nombres de las variables intermedias explican su papel en el cálculo: shiftx es el valor de x desplazado por theta, no debería sorprender que exponente sea el argumento de exp, y denom termine en el denominador. Elegir nombres informativos hace que el código sea más fácil de leer y entender, lo que facilita la depuración.

Errores

Todo error es una oportunidad de aprendizaje. Siempre que aprendas una nueva característica, debes intentar cometer tantos errores como sea posible, lo antes posible. Cuando cometes errores deliberados, ves cuáles son los mensajes de error. Posteriormente, cuando cometas errores accidentales, sabrás lo que significan los mensajes.

Veamos algunos errores comunes. Uno grande para principiantes es dejar fuera el * para multiplicar, como en este ejemplo:

area = pi r^2

Este código debería producir el siguiente mensaje de error:

area = pi r^2 | Error: Expresión no válida. Verifique si falta operador de multiplicación, delimitadores faltantes o desequilibrados, u otro error de sintaxis. Para construir matrices, use corchetes en lugar de paréntesis.

El mensaje indica que la expresión no es válida y sugiere varias cosas que podrían estar mal. En este caso, una de sus conjeturas es acertada: nos falta un operador de multiplicación.

Otro error común es dejar fuera los paréntesis alrededor de los argumentos de una función. Por ejemplo, en notación matemática es común escribir algo así como$\sin \pi$, pero en MATLAB si escribes

sin pi

debería recibir el siguiente mensaje de error:

Función indefinida 'sin' para argumentos de entrada de tipo 'char'.

El problema es que cuando dejas fuera los paréntesis, MATLAB trata el argumento como una cadena de caracteres (que tienen tipo 'char'). En este caso el mensaje de error es útil, pero en otros casos los resultados pueden ser desconcertantes. Por ejemplo, si llamas abs, que calcula valores absolutos, y olvidas los paréntesis, obtienes un resultado sorprendente:

>> abs pi
ans =  112   105

No voy a explicar este resultado; por ahora, sólo voy a sugerir que se ponga paréntesis alrededor de los argumentos.

Aquí hay otro error común. Si estuvieras traduciendo la expresión matemática\[\frac{1}{2 \sqrt \pi}\] a MATLAB, podrías tener la tentación de escribir esto:

1 / 2 * sqrt(pi)

Pero eso estaría mal por el orden de las operaciones. La división y la multiplicación se evalúan de izquierda a derecha, por lo que esta expresión multiplicaría 1/2 por sqrt (pi).

Para mantener sqrt (pi) en el denominador, podrías usar paréntesis,

1 / (2 * sqrt(pi))

o hacer explícita la división,

1 / 2 / sqrt(pi)

Los dos últimos ejemplos nos llevan al Segundo Teorema de Depuración:

Lo único peor que recibir un mensaje de error es no recibir un mensaje de error.

Los programadores principiantes a menudo odian los mensajes de error y hacen todo lo posible para que los mensajes desaparezcan. Los programadores experimentados saben que los mensajes de error son tu amigo. Pueden ser difíciles de entender, e incluso engañosas, pero merece la pena el esfuerzo para entenderlas.

Documentación

MATLAB viene con dos formas de documentación, help y doc. El comando help funciona en la Ventana de Comandos; simplemente ingrese help seguido del nombre de un comando.

>> help sin

Deberías ver una salida como esta:

sin Seno de argumento en radianes. sin (X) es el seno de los elementos de X.

Ver también asin, sind, sinpi.

Alguna documentación usa vocabulario que aún no hemos cubierto. Por ejemplo, los elementos de X podrían no tener sentido hasta que lleguemos a vectores y matrices unos capítulos a partir de ahora.

Las páginas de doc suelen ser mejores. Si ingresa doc sin, aparece una ventana del navegador con información más detallada sobre la función, incluyendo ejemplos de cómo usarla. Los ejemplos suelen utilizar vectores y matrices, por lo que puede que aún no tengan sentido, pero puedes obtener una vista previa de lo que viene.

Revisión del Capítulo

Este capítulo proporciona una visión general del proceso de modelado, incluyendo abstracción, análisis y simulación, medición y validación.

También introdujo MATLAB, el lenguaje de programación que usaremos para escribir simulaciones. Hasta ahora, hemos visto variables y valores, operaciones aritméticas y funciones matemáticas.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

El intérprete es el programa que lee y ejecuta MATLAB o código. Imprime un prompt para indicar que está esperando que escribas un comando, que es una línea de código ejecutada por el intérprete.

Un operador es un símbolo, como * o +, que representa una operación matemática. Un operando es un número o variable que aparece en una expresión junto con operadores. Una expresión es una secuencia de operandos y operadores que especifica un cálculo matemático y produce un valor.

Una función es un cálculo con nombre; por ejemplo, log10 es el nombre de una función que calcula logaritmos en base 10. Una llamada a función es un comando que hace que una función se ejecute y calme un resultado. Un argumento es una expresión que aparece en una llamada a una función para especificar el valor sobre el que opera la función.

Una variable es un valor con nombre. Una instrucción de asignación es un comando que crea una nueva variable (si es necesario) y le da un valor. Un espacio de trabajo es un conjunto de variables y sus valores.

Finalmente, una cadena es un valor que consiste en una secuencia de caracteres (a diferencia de un número).

En el siguiente capítulo, comenzarás a escribir programas más largos y aprenderás sobre los números de coma flotante.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

[centavo] Quizás hayas escuchado que un centavo caído desde lo alto del Empire State Building iría tan rápido cuando golpeara el pavimento que quedaría incrustado en el concreto o que si golpeaba a una persona le rompería el cráneo.

Podemos probar este mito haciendo y analizando un modelo. Para comenzar, asumiremos que el efecto de la resistencia al aire es pequeño. Esto va a resultar una mala suposición, pero tengan cuidado conmigo.

Si la resistencia del aire es insignificante, la fuerza primaria que actúa sobre el centavo es la gravedad, lo que hace que el centavo se acelere hacia abajo.

Si la velocidad inicial es 0, la velocidad después de$t$ segundos es$a t$, y la distancia que ha caído el centavo es\[h = a t^2 / 2\] Usando álgebra, podemos resolver para$t$:\[t = \sqrt{ 2 h / a}\] Taponar la aceleración de la gravedad, $a = \SI{9.8}{\meter\per\second\squared}$, y la altura del Empire State Building$h = \SI{381}{\meter}$,, obtenemos$t = \SI{8.8}{\second}$. Entonces, la computación$v = a t$ obtenemos una velocidad de impacto de$\SI{86}{\meter\per\second}$, que es de aproximadamente 190 millas por hora. Eso suena como que podría doler.

Utilice MATLAB para realizar estos cálculos y compruebe que obtiene el mismo resultado.

El resultado en el ejercicio anterior no es exacto porque ignora la resistencia del aire. En realidad, una vez que el centavo llega a unos 18 m, la fuerza ascendente de la resistencia aérea equivale a la fuerza descendente de la gravedad, por lo que el centavo deja de acelerar. Después de eso, no importa hasta dónde caiga el centavo; golpea la acera a unos 18 m, mucho menos de 86 m.

Como ejercicio, computar el tiempo que tarda el centavo en llegar a la acera si asumimos que acelera con aceleración constante$a = \SI{9.8}{\meter\per\second\squared}$ hasta alcanzar la velocidad terminal, luego cae con velocidad constante hasta llegar a la acera.

El resultado que obtienes no es exacto, pero es una aproximación bastante buena.

Guiones

Hasta ahora hemos escrito todos nuestros programas “al instante”. Si solo estás escribiendo unas pocas líneas, esto no es tan malo. Pero, ¿y si estás escribiendo cien? Volver a escribir cada línea de código cada vez que quiera cambiar o probar su programa llevará mucho tiempo y será tedioso. Por suerte, no tienes que hacerlo. En este capítulo, veremos una manera de ejecutar muchas líneas a la vez: scripts.

Tu Primer Guión

Un script es un archivo que contiene código MATLAB. Cuando ejecuta un script, MATLAB ejecuta los comandos en él, uno tras otro, exactamente como si los hubiera escrito en el prompt. Los scripts también se denominan a veces M-files porque usan la extensión .m, abreviatura de MATLAB.

Puede crear scripts con cualquier editor de texto o procesador de textos, pero la forma más sencilla es hacer clic en el botón Nuevo script en la esquina superior izquierda de la interfaz de MATLAB, que abre un editor de texto diseñado para MATLAB.

Para probarlo, cree un nuevo script e ingrese el siguiente código:

x = 5

A continuación, pulsa el botón Guardar. Debería aparecer una ventana de diálogo donde pueda elegir el nombre del archivo y la carpeta donde irá su script. Cambia el nombre a myscript.m y guárdalo en cualquier carpeta que quieras.

Ahora haz clic en el botón verde Ejecutar. Es posible que reciba un mensaje que diga que el script no se encuentra en la carpeta actual. Si es así, haga clic en el botón que dice Cambiar carpeta y debería ejecutarse.

También puede ejecutar un script escribiendo su nombre en la Ventana de Comandos y presionando. Por ejemplo, si ingresa myscript, MATLAB debe ejecutar su script y mostrar el resultado:

>> myscript
x = 5

Hay algunas cosas a tener en cuenta al usar scripts. Primero, no debes incluir la extensión .m cuando ejecutas un script. Si lo haces, recibirás un mensaje de error como este:

>> myscript.m
Undefined variable "myscript" or class "myscript.m".

Segundo, cuando nombras un nuevo guión, trata de elegir algo significativo y memorable. No elija un nombre que ya esté en uso; si lo hace, reemplazará una de las funciones de MATLAB por la suya propia (al menos temporalmente). Puede que no te des cuenta de inmediato, pero podrías tener algún comportamiento confuso más tarde.

Además, el nombre del script no puede contener espacios. Si crea un archivo llamado my script.m, MATLAB se quejará cuando intente ejecutarlo:

>> my script
Undefined function or variable 'my'.

Puede ser difícil recordar en qué carpeta se encuentra un script. Para que las cosas sean simples, por ahora, te sugiero poner todos tus scripts en una carpeta.

¿Por qué Scripts?

Hay algunas buenas razones para usar un guión. Cuando escribes más de un par de líneas de código, podría tomar algunos intentos para que todo esté bien. Poner tu código en un script hace que sea más fácil de editar que escribirlo en el prompt. Asimismo, si estás ejecutando un script repetidamente, ¡es mucho más rápido escribir el nombre del script que reescribir el código! Y es posible que pueda reutilizar un guión de un proyecto a otro, ahorrándole un tiempo considerable en todos los proyectos.

Pero el gran poder de los scripts viene con una gran responsabilidad: hay que asegurarse de que el código que está ejecutando sea el código que cree que está ejecutando. Siempre que inicies un nuevo script, comienza con algo simple, como x = 5, que produzca un efecto visible. Después ejecuta tu guión y confirma que obtienes lo que esperas. Cuando escribe el nombre de un script, MATLAB busca el script en una ruta de búsqueda, que es una secuencia de carpetas. Si no encuentra el script en la primera carpeta, busca en la segunda, y así sucesivamente. Si tienes scripts con el mismo nombre en diferentes carpetas, podrías estar mirando una versión y ejecutando otra.

Si el código que estás ejecutando no es el código que estás viendo, ¡encontrarás que depurar es un ejercicio frustrante! Entonces no es de extrañar que este sea el Tercer Teorema de Depuración:

Asegúrate de que el código que estás ejecutando es el código que crees que estás ejecutando.

Ahora que ya has visto cómo escribir un guión, usemos uno para hacer algo un poco más complicado.

La secuencia de Fibonacci

La secuencia de Fibonacci, denotada$F$, es una secuencia de números donde cada número es la suma de los dos anteriores. Se define por las ecuaciones$F_1 = 1$,$F_2 = 1$, y, para$i > 2$,$F_{i} = F_{i-1} + F_{i-2}$. La siguiente expresión calcula el número de Fibonacci$n$ th:\[F_n = \frac{1}{\sqrt{5}} \left[ \left( \frac{1 + \sqrt{5}}{2} \right)^{n} - \left( \frac{1 - \sqrt{5}}{2} \right)^{n} \right]\]

Podemos traducir esta expresión en MATLAB, así:

s5 = sqrt(5);
t1 = (1 + s5) / 2;
t2 = (1 - s5) / 2;
diff = t1^n - t2^n;
ans = diff / s5

Utilizo variables temporales como t1 y t2 para hacer que el código sea legible y explícito el orden de las operaciones. Las primeras cuatro líneas tienen punto y coma al final, por lo que no muestran nada. La última línea asigna el resultado a ans.

Si guardamos este script en un archivo llamado fibonacci1.m, podemos ejecutarlo así:

>> n = 10
>> fibonacci1
ans = 55.0000

Antes de llamar a este script, hay que asignar un valor a n. Si no se define n, se obtiene un error:

>> clear n
>> fibonacci1
Undefined function or variable 'n'.

Error in fibonacci1 (line 9)
diff = t1^n - t2^n;

Este script solo funciona si hay una variable llamada n en el espacio de trabajo; de lo contrario, debería obtener un error. MATLAB le dirá en qué línea del script se encuentra el error y mostrará la línea.

Los mensajes de error pueden ser útiles, ¡pero ten cuidado! En este ejemplo, el mensaje dice que el error está en fibonacci, pero el problema real es que no hemos asignado un valor a n. Y eso nos lleva al Cuarto Teorema de Depuración:

Los mensajes de error indican dónde se descubrió el problema, no dónde fue causado.

A menudo hay que trabajar hacia atrás para encontrar la línea de código (o línea faltante) que causó el problema.

Números de punto flotante

En el apartado anterior, el resultado que calculamos fue 55.0000. Dado que los números de Fibonacci son enteros, es posible que se haya sorprendido al ver los ceros después del punto decimal.

Están ahí porque MATLAB realiza cálculos utilizando números de coma flotante. Con números de coma flotante, los enteros se pueden representar exactamente, pero la mayoría de las fracciones no pueden.

Por ejemplo, si calculas la fracción 2/3, el resultado es solo aproximado: la respuesta correcta tiene un número infinito de 6s:

>> 2/3
ans = 0.6666

No es tan malo como parece este ejemplo: MATLAB usa más dígitos de los que muestra por defecto. Puede usar el comando format para cambiar el formato de salida:

>> format long
>> 2/3
ans = 0.666666666666667

En este ejemplo, los primeros 14 dígitos son correctos; el último ha sido redondeado.

Los números grandes y pequeños se muestran en notación científica. Por ejemplo, si usamos la función factorial incorporada para calcular$100!$, obtenemos el siguiente resultado:

>> factorial(100)
ans = 9.332621544394410e+157

La e en esta notación no es el número trascendental conocido como$e$; es solo una abreviatura de “exponente”. Entonces esto significa que$100!$ es aproximadamente$9.33 \times 10^{157}$. La solución exacta es un entero de 158 dígitos, pero con números de punto flotante de doble precisión, solo conocemos los primeros 16 dígitos.

Se pueden introducir números usando la misma notación.

>> speed_of_light = 3.0e8
speed_of_light = 300000000

Aunque el formato de coma flotante puede representar números muy grandes y pequeños, hay límites. Las variables predefinidas realmax y realmin contienen los números más grandes y pequeños que MATLAB puede manejar.

>> realmax
ans = 1.797693134862316e+308

>> realmin
ans = 2.225073858507201e-308

Si el resultado de un cálculo es demasiado grande, MATLAB “redondea” a.

>> factorial(170)
ans = 7.257415615307994e+306

>> factorial(171)
ans = Inf

División por cero también devuelve Inf.

>> 1/0
ans = Inf

Para operaciones que no están definidas, MATLAB devuelve NaN, que significa “no un número”.

>> 0/0
ans = NaN

Comentarios

Los programas cortos y simples son fáciles de leer, pero a medida que se hacen más grandes y complejos, se vuelve más difícil averiguar qué hacen y cómo. Para eso están los comentarios.

Un textual es una línea de texto que se agrega a un programa para explicar cómo funciona. No tiene efecto en la ejecución del programa; está ahí para los lectores humanos. Los buenos comentarios hacen más legibles los programas; los malos comentarios son inútiles en el mejor de los casos y engañosos en

Para escribir un comentario, se utiliza el símbolo de porcentaje (%) seguido del texto del comentario.

>> speed_of_light = 3.0e8     % meters per second
speed_of_light = 300000000

El comentario va desde el símbolo de porcentaje hasta el final de la línea. En este caso especifica las unidades del valor. En un mundo ideal, MATLAB haría un seguimiento de las unidades y las propagaría a través del cálculo, pero por ahora esa carga recae sobre el programador.

Evite comentarios que sean redundantes con el código:

>> x = 5        % assign the value 5 to x

Los buenos comentarios proporcionan información adicional que no está en el código, como unidades en el ejemplo anterior, o el significado de una variable:

>> p = 0         % position from the origin in meters
>> v = 100       % velocity in meters / second
>> a = -9.8      % acceleration of gravity in meters / second^2

Si usas nombres de variables más largos, es posible que no necesites comentarios explicativos, pero hay una compensación: los nombres de variables más largos son más claros, pero el código más largo puede volverse más difícil de leer. Además, si estás traduciendo de matemáticas que usan nombres cortos de variables, puede ser útil para que tu programa sea consistente con tus matemáticas.

Documentación

Todo script debe proporcionar documentación, que es un comentario que explique qué hace el script y cuáles son sus requisitos.

Por ejemplo, podría poner algo como esto al inicio de fibonacci1.m:

% Computes a numerical approximation of the nth Fibonacci number.
% Precondition: you must assign a value to n before running this script.
% Postcondition: the result is stored in ans.

Una condición previa es algo que debe ser cierto cuando se inicia el guión para que funcione correctamente. Una condición posterior es algo que será cierto cuando se complete el guión.

Si hay un comentario al principio de un script, MATLAB asume que es la documentación del script. Entonces, si escribes ayuda fibonacci1, obtienes el contenido del comentario (sin los signos porcentuales).

>> help fibonacci1
  Computes a numerical approximation of the nth Fibonacci number.
  Precondition: you must assign a value to n before running this script.
  Postcondition: the result is stored in ans.

De esa manera, los scripts que escribes se comportan igual que los scripts predefinidos. Incluso puedes usar el comando doc para ver tu comentario en la Ventana de Ayuda.

Asignación e Igualdad

Para los programadores principiantes, una fuente común de confusión es la asignación y el uso del signo igual.

En matemáticas, el signo igual significa que los dos lados de la ecuación tienen el mismo valor. En MATLAB, una sentencia de asignación parece una igualdad matemática, pero no lo es.

Una diferencia es que los lados de una declaración de asignación no son intercambiables. El lado derecho puede ser cualquier expresión legal, pero el lado izquierdo tiene que ser una variable, que se llama el objetivo de la asignación. Entonces esto es legal:

>> y = 1;
>> x = y + 1
x = 2

Pero esto no es:

>> y + 1 = x
 y + 1 = x
       |
Error: Incorrect use of '=' operator.
To assign a value to a variable, use '='.
To compare values for equality, use '=='.

En este caso el mensaje de error no es muy útil. El problema aquí es que la expresión del lado izquierdo no es un objetivo válido para una asignación.

Otra diferencia entre asignación e igualdad es que una igualdad matemática es verdadera (o falsa) por toda la eternidad; una declaración de asignación es temporal. Cuando asignas x = y + 1, obtienes el valor actual de y. Si y cambia más tarde, x no se actualiza.

Una tercera diferencia es que una igualdad matemática es una afirmación que puede o no ser cierta. En matemáticas,$y = y+1$ es una afirmación que pasa a ser falsa para todos los valores de$y$. En MATLAB, y = y + 1 es una sentencia de asignación sensible y útil. Lee el valor actual de y, agrega 1, y reemplaza el valor antiguo por el nuevo valor.

>> y = 1;
>> y = y + 1
y = 2

Cuando lee el código de MATLAB, puede resultarle útil pronunciar el signo igual como “obtiene” en lugar de “igual”. Entonces x = y + 1 se pronuncia “x obtiene el valor de y más uno”.

Revisión del Capítulo

En este capítulo se presentaron guiones y se sugirieron razones para utilizarlos. Calculamos elementos de una secuencia de Fibonacci, pero debido a que usamos números de coma flotante, los resultados a veces solo fueron aproximados. Y vimos cómo agregar comentarios a un programa para documentar lo que hace y explicar cómo funciona.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

Un archivo M es un archivo que contiene un script, que es una secuencia de comandos MATLAB/Octave. La ruta de búsqueda es la lista de carpetas donde el intérprete busca M-files.

Una condición previa es algo que debe ser cierto cuando se inicia el script para que funcione correctamente; una condición posterior es algo que será cierto cuando se complete el script.

El objetivo de una sentencia de asignación es la variable en el lado izquierdo.

El punto flotante es una forma de representar y almacenar números en una computadora. es un formato para escribir y mostrar números grandes y pequeños; por ejemplo, 3.0e8 representa$3.0 \times 10^8$ o 300,000,000.

Un textual es parte de un programa que proporciona información adicional sobre el programa, pero no afecta su ejecución.

En el siguiente capítulo, aprenderás a escribir programas que realizan tareas repetitivas usando bucles.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Para poner a prueba su comprensión de las declaraciones de asignación, escriba algunas líneas de código que intercambien los valores de x e y. Pon tu código en un script llamado swap.m y pruébelo.

Si funciona correctamente, deberías poder ejecutarlo así:

>> x = 1, y = 2
x = 1
y = 2

>> swap

>> x, y
x = 2
y = 1

[Bikegame]

Imagina que eres el operador de un sistema de bicicletas compartidas con dos ubicaciones: Boston y Cambridge.

Observa que todos los días el 5 por ciento de las bicicletas en Boston se dejan en Cambridge, y el 3 por ciento de las bicicletas en Cambridge se dejan en Boston. A principios de mes, hay 100 bicicletas en cada ubicación.

Escribe un guión llamado bikeupdate.m que actualice el número de bicicletas en cada ubicación de un día a otro. La condición previa es que las variables b y c contengan el número de bicicletas en cada ubicación al inicio del día. La condición posterior es que b y c se han modificado para reflejar el movimiento neto de las bicicletas.

Para probar su programa, inicialice b y c en el prompt y luego ejecute el script. El script debe mostrar los valores actualizados de b y c, pero no ninguna variable intermedia.

Recuerda que las bicicletas son cosas contables, por lo que b y c siempre deben ser valores enteros. Es posible que desee usar la función round para calcular el número de bicicletas que se mueven cada día.

Si ejecutas tu script repetidamente, puedes simular el paso del tiempo de un día a otro (puedes repetir un comando presionando la flecha y luego).

¿Qué pasa con las bicis? ¿Todos terminan en un solo lugar? ¿El sistema alcanza un equilibrio, oscila o hace otra cosa?

En el siguiente capítulo, veremos cómo ejecutar tu script automáticamente y cómo trazar los valores de b y c a lo largo del tiempo.

Loops

Los programas que hemos visto hasta ahora son de código de línea recta; es decir, ejecutan una instrucción tras otra de arriba a abajo. Este capítulo introduce una de las características más importantes del lenguaje de programación, el bucle for, que permite a programas simples realizar tareas complejas y repetitivas. Este capítulo también introduce los conceptos matemáticos de secuencia y serie, y un proceso para escribir programas, desarrollo incremental.

Empezaremos revisando el ejercicio del capítulo anterior; si no lo hiciste, es posible que quieras echarle un vistazo antes de continuar.

Actualización de Variables

En la Sección [juego de bicicletas], te pedí que escribieras un programa que modele un sistema de bicicletas compartidas con bicicletas moviéndose entre dos estaciones. Cada día, el 5 por ciento de las bicicletas en Boston se dejan en Cambridge, y el 3 por ciento de las bicicletas en Cambridge se dejan en Boston.

Para actualizar el estado del sistema, es posible que hayas tenido la tentación de escribir algo como

b = b - 0.05*b + 0.03*c
c = c + 0.05*b - 0.03*c

Pero eso estaría mal, así que muy mal. ¿Por qué? El problema es que la primera línea cambia el valor de b, así que cuando se ejecuta la segunda línea, obtiene el valor antiguo de c y el nuevo valor de b. En consecuencia, el cambio en b no siempre es lo mismo que el cambio en c, ¡que viola el Principio de Conservación de Bicicletas!

Una solución es usar variables temporales como b_new y c_new:

b_new = b - 0.05*b + 0.03*c
c_new = c + 0.05*b - 0.03*c
b = b_new
c = c_new

Esto tiene el efecto de actualizar las variables simultáneamente; es decir, lee ambos valores antiguos antes de escribir cualquiera de los valores nuevos.

La siguiente es una solución alternativa que tiene la ventaja añadida de simplificar el cálculo:

b_to_c = 0.05*b - 0.03*c
b = b - b_to_c
c = c + b_to_c

Es fácil mirar este código y confirmar que obedece a Conservación de Bicicletas. Incluso si el valor de b_to_c es incorrecto, al menos el número total de bicicletas es correcto. Y eso nos lleva al Quinto Teorema de Depuración:

La mejor manera de evitar un error es hacerlo imposible.

En este caso, eliminar la redundancia también elimina la oportunidad de un error.

Taxonomía de errores

Cuanto más entiendas los errores, mejor estarás en la depuración. Hay cuatro tipos de errores:

Error de sintaxis: Ha escrito un comando que no puede ejecutarse porque viola una de las reglas de sintaxis del lenguaje. Por ejemplo, en MATLAB, no puede tener dos operandos seguidos sin un operador, por lo que pi r^2 contiene un error de sintaxis. Cuando el intérprete encuentra un error de sintaxis, imprime un mensaje de error y deja de ejecutar su programa.
Error de tiempo de ejecución: Tu programa empieza a funcionar, pero algo sale mal en el camino. Por ejemplo, si intentas acceder a una variable que no existe, eso es un error de tiempo de ejecución. Cuando el intérprete detecta el problema, imprime un mensaje de error y se detiene.
Error lógico: Tu programa se ejecuta sin generar ningún mensaje de error, pero no hace lo correcto. El problema en la sección anterior, donde cambiamos el valor de b antes de leer el valor antiguo, es un error lógico.
Error numérico: La mayoría de los cálculos en MATLAB son aproximadamente correctos. La mayoría de las veces los errores son lo suficientemente pequeños como para que no nos importe, pero en algunos casos los errores de redondeo son un problema.

Los errores de sintaxis suelen ser los más fáciles de tratar. A veces los mensajes de error son confusos, pero MATLAB suele indicarle dónde está el error, al menos aproximadamente.

Los errores de tiempo de ejecución son más difíciles porque, como mencioné antes, MATLAB puede indicarle dónde detectó el problema, pero no qué lo causó.

Los errores lógicos son difíciles porque MATLAB no puede ayudar en absoluto. Desde el punto de vista de MATLAB no hay nada malo en el programa; solo tú sabes lo que se supone que debe hacer el programa, así que solo tú puedes verificarlo.

Los errores numéricos pueden ser complicados porque no está claro si el problema es tu culpa. Para la mayoría de los cálculos simples, MATLAB produce el valor de coma flotante más cercano a la solución exacta, lo que significa que los primeros 15 dígitos significativos deben ser correctos.

Pero algunos cálculos están mal condicionados, lo que significa que incluso si tu programa es correcto, los errores de redondeo se acumulan y el número de dígitos correctos puede ser menor. A veces MATLAB te puede advertir que esto está sucediendo, ¡pero no siempre! La precisión (el número de dígitos en la respuesta) no implica exactitud (el número de dígitos que son correctos).

Error absoluto y relativo

Hay dos formas de pensar sobre los errores numéricos. El primero es el error absoluto, o la diferencia entre el valor correcto y la aproximación. A menudo escribimos la magnitud del error, ignorando su signo, cuando no importa si la aproximación es demasiado alta o demasiado baja.

La segunda forma de pensar sobre los errores numéricos es el error relativo, donde el error se expresa como una fracción (o porcentaje) del valor exacto.

Por ejemplo, podríamos querer estimar$9!$ usando la fórmula$\sqrt {18 \pi} ( 9 / e)^9$. La respuesta exacta es$9 \cdot 8 \cdot 7 \cdot 6 \cdot 5 \cdot 4 \cdot 3 \cdot 2 \cdot 1 = 362,880$. La aproximación es$359,536.87$. Entonces el error absoluto es$3,343.13$.

A primera vista, eso suena como mucho —estamos fuera por tres mil— pero debemos considerar el tamaño de lo que estamos estimando. Por ejemplo, los $3,000 importan mucho si estamos hablando de un salario anual, pero para nada si estamos hablando de la deuda nacional.

Una forma natural de manejar este problema es usar el error relativo. En este caso, dividiríamos el error por$362,880$, rindiendo$0.00921$, que es apenas menos del 1 por ciento. Para muchos propósitos, estar apagado en un 1 por ciento es bastante bueno.

para Loops

Volvamos al ejemplo de bicicletas compartidas. En el capítulo anterior, escribimos un guión bikeupdate.m para simular un día en la vida de un sistema de bicicletas compartidas. Para simular un mes entero, tendríamos que ejecutar el script 30 veces. Podríamos ingresar el mismo comando 30 veces, pero es más sencillo usar un loop, que es un conjunto de sentencias que se ejecuta repetidamente.

Para crear un bucle, podemos usar la sentencia for, así:

for i=1:30
    bike_update
end

La primera línea incluye lo que parece una declaración de asignación, y es como una instrucción de asignación, excepto que se ejecuta más de una vez. La primera vez, crea la variable i y le asigna el valor 1. La segunda vez, i obtiene el valor 2, y así sucesivamente, hasta e incluyendo 30.

El operador de dos puntos (:) especifica un rango de números enteros. Puede crear un rango en el símbolo del sistema:

>> 1:5
ans =  1     2     3     4     5

La variable que se utiliza en la sentencia for se llama la variable de bucle. Es común usar los nombres i, j y k como variables de bucle.

Las declaraciones dentro del bucle se llaman el cuerpo. Por convención, se les aplica sangría para mostrar que están dentro del bucle, pero la sangría no afecta la ejecución del programa. La sentencia end marca el final del bucle.

Para ver un bucle en acción, puede ejecutar uno que muestre la variable de bucle:

>> for i=1:5
    i
end

i = 1
i = 2
i = 3
i = 4
i = 5

Como muestra este ejemplo, se puede ejecutar un bucle for desde la línea de comandos, pero es más común ponerlo en un script.

Cree un script llamado bikeloop.m que utilice un bucle for para ejecutar bikeupdate.m 30 veces. Antes de ejecutarlo, hay que asignar valores a b y c. Para este ejercicio, comience con los valores b = 100 y c = 100.

Si todo va bien, su guión mostrará una larga secuencia de números en la pantalla. Probablemente sea demasiado largo para encajar, e incluso si encajara, sería difícil de interpretar. ¡Una gráfica sería mucho mejor!

Trazar

Si la salida de tu programa es un largo flujo de números, puede ser difícil ver lo que está sucediendo. Trazar los resultados puede hacer las cosas más claras.

La función plot es una herramienta versátil para trazar gráficos bidimensionales. Desafortunadamente, es tan versátil que puede ser difícil de usar (y difícil de leer la documentación). Empezaremos simples y trabajaremos nuestro camino hacia arriba.

Para trazar un solo punto, escriba

>> plot(1, 2, 'o')

Una Ventana de Figura debe aparecer con una gráfica y un solo círculo azul en la$x$ posición 1 y$y$ la posición 2.

La letra entre comillas simples es una cadena de estilo que especifica cómo se debe trazar el punto; o indica un círculo. Otras formas incluyen +, *, x, s (para un cuadrado), d (para un diamante) y ^ (para un triángulo).

También puede especificar el color iniciando la cadena de estilo con un código de color:

>> plot(1, 2, 'ro')

Aquí, r significa rojo; los otros colores incluyen g para verde, b para azul, c para cian, m para magenta, y para amarillo y k para negro.

Cuando usas la trama de esta manera, solo puede trazar un punto a la vez. Si vuelves a ejecutar la trama, borra la cifra antes de hacer la nueva trama. El comando hold le permite anular ese comportamiento: hold on le dice a MATLAB que no borre la figura cuando hace una nueva gráfica; hold off devuelve al comportamiento predeterminado.

Prueba esto:

>> clf
>> hold on
>> plot(1, 1, 'ro')
>> plot(2, 2, 'go')
>> plot(3, 3, 'bo')
>> hold off

El comando clf borra la figura antes de comenzar a trazar.

Si ejecutas el código anterior, deberías ver una figura con tres círculos. MATLAB escala la gráfica automáticamente para que los ejes se ejecuten desde los valores más bajos de la gráfica hasta los más altos.

Modifique bikeloop.m para que borre la figura antes de ejecutar el bucle. Entonces, cada vez que pasa por el bucle, se debe trazar el valor de b versus el valor de i con un círculo rojo.

Una vez que lo consigas funcionando, modifícalo para que grafique los valores de c con diamantes azules.

Secuencias

Ahora que tenemos la capacidad de escribir bucles, podemos utilizarlos para explorar secuencias y series, que son útiles para describir y analizar sistemas que cambian con el tiempo.

En matemáticas, una secuencia es un conjunto de números que corresponde a los enteros positivos. Los números en la secuencia se llaman elementos. En notación matemática, los elementos se denotan con subíndices, por lo que el primer elemento de la serie$A$ es$A_1$, seguido de$A_2$, y así sucesivamente.

Un bucle for es una forma natural de calcular los elementos de una secuencia. Como ejemplo, en una secuencia geométrica, cada elemento es un múltiplo constante del elemento anterior. Como ejemplo más específico, veamos la secuencia con$A_1 = 1$ y la relación$A_{i+1} = A_i/2$, para todos$i$. En otras palabras, cada elemento es la mitad de grande que el anterior a él.

El siguiente bucle calcula los primeros 10 elementos de$A$:

a = 1
for i=2:10
    a = a/2
end

La primera línea inicializa la variable a con el primer elemento de la secuencia,$A_1$. Cada vez que pasa por el bucle, encontramos el siguiente valor de a dividiendo el valor anterior por 2, y asignamos el resultado de nuevo a a. Al final, a contiene el décimo elemento.

Los demás elementos se muestran en la pantalla, pero no se guardan en una variable. Posteriormente, veremos cómo guardar los elementos de una secuencia en un vector.

Este bucle calcula la secuencia de manera recurrente, lo que significa que cada elemento depende del anterior. Para esta secuencia, también es posible calcular el elemento$i$ th directamente, en función de$i$, sin usar el elemento anterior.

En notación matemática,$A_i = A_1 r^{i-1}$, donde$r$ está la relación de elementos sucesivos. En el ejemplo anterior,$A_{i+1} = A_i/2$, entonces$r = 1/2$.

Escriba un script llamado sequence.m que utilice un bucle para calcular elementos de$A$.

Serie

En matemáticas, una serie es la suma de los elementos de una secuencia. Es un nombre terrible, porque en el inglés común, “sequence” y “series” significan prácticamente lo mismo, pero en matemáticas, una secuencia es un conjunto de números, y una serie es una expresión (una suma) que tiene un solo valor. En notación matemática, una serie a menudo se escribe usando el símbolo de suma$\sum$.

Por ejemplo, la suma de los primeros 10 elementos de$A$ es\[\sum_{i=1}^{10} A_i\]

Un bucle for es una forma natural de calcular el valor de esta serie:

A1 = 1;
total = 0;
for i=1:10
    a = A1 * (1/2)^(i-1);
    total = total + a;
end
ans = total

Vamos a recorrer lo que está pasando aquí. A1 es el primer elemento de la secuencia, por lo que le asignamos para que sea 1; también creamos total, que almacenará la suma acumulativa. Cada vez que pasa por el bucle, establecemos a al elemento$i$ th y agregamos a al total. Al final, fuera del bucle, almacenamos total como ans.

La forma en que estamos usando total se llama acumulador, es decir, una variable que acumula el resultado poco a poco.

Este ejemplo calcula los términos de la serie directamente. Como ejercicio, escribir un guión llamado series.m que compute la misma suma calculando los elementos de manera recurrente. Tendrás que tener cuidado con por dónde comienzas y detengas el bucle.

Generalización

Tal y como está escrito, el ejemplo anterior siempre suma los primeros 10 elementos de la secuencia, pero podríamos tener curiosidad por saber qué pasa al total a medida que aumentamos el número de términos en la serie. Si has estudiado series geométricas, podrías saber que esta serie converge en 2; es decir, a medida que el número de términos va al infinito, la suma se acerca a 2 asintóticamente.

Para ver si eso es cierto para nuestro programa, podemos reemplazar la constante 10 en Listing [lst:series_10] con una variable llamada n:

A1 = 1;
total = 0;
for i=1:n
    a = A1 * (1/2)^(i-1);
    total = total + a;
end
ans = total

El código en Listing [lst:series_n] ahora puede calcular cualquier número de términos, con la condición previa de que se tiene que establecer n antes de ejecutar el código. Puse este código en un archivo llamado series.m, luego lo ejecuté con diferentes valores de n:

>> n=10; series
total = 1.99804687500000

>> n=20; series
total = 1.99999809265137

>> n=30; series
total = 1.99999999813735

>> n=40; series
total = 1.99999999999818

Seguro que parece que está convergiendo en 2.

Reemplazar una constante por una variable se llama generalización. En lugar de computar un número fijo y específico de términos, el nuevo script es más general; puede calcular cualquier número de términos. Esta es una idea importante a la que volveremos cuando hablemos de funciones.

Desarrollo Incremental

A medida que comienzas a escribir programas más largos, es posible que pases más tiempo depurando. Cuanto más código escribas antes de comenzar a depurar, más difícil es encontrar el problema.

El desarrollo incremental es una forma de programación que trata de minimizar el dolor de la depuración mediante el desarrollo y la prueba en pequeños pasos. Los pasos fundamentales son:

Siempre comience con un programa de trabajo. Si tienes un ejemplo de un libro, o un programa que escribiste que sea similar a lo que estás trabajando, empieza con eso. De lo contrario, empieza con algo que sepas que es correcto, como x = 5. Ejecute el programa y confirme que estás ejecutando el programa que crees que estás ejecutando. Este paso es importante, porque en la mayoría de los ambientes hay pequeñas cosas que pueden tropezar cuando comienzas un nuevo proyecto. Quítelos del camino para que puedas concentrarte en la programación.
Haga un cambio pequeño y comprobable a la vez. Un cambio comprobable es aquel que muestra algo en la pantalla (o tiene algún otro efecto) que puedes verificar. Idealmente, debes saber cuál es la respuesta correcta o poder verificarla realizando otro cálculo.
Ejecute el programa y vea si el cambio funcionó. Si es así, vuelve al paso 2. Si no, tendrás que hacer alguna depuración, pero si el cambio que hiciste fue pequeño, no debería tardar mucho en encontrar el problema.

Con el desarrollo incremental, es más probable que tu código funcione la primera vez, y si no lo hace, es más probable que el problema sea obvio. Y eso nos lleva al Sexto Teorema de Depuración:

El mejor tipo de depuración es el tipo que no tienes que hacer.

En la práctica, hay dos problemas con el desarrollo incremental. Primero, a veces hay que escribir código extra para generar salida visible que pueda verificar. Este código extra se llama andamiaje porque lo usas para construir el programa y luego eliminarlo cuando termines. Pero el tiempo que ahorras en la depuración casi siempre vale la pena el tiempo que inviertes en andamios.

El segundo problema es que cuando estás empezando, es posible que no sepas cómo elegir pasos que van desde x = 5 hasta el programa que intentas escribir. Veremos un ejemplo extendido en el Capítulo 6.

Si te encuentras escribiendo más de unas pocas líneas de código antes de comenzar a probar y estás gastando mucho tiempo depurando, deberías probar el desarrollo incremental.

Revisión del Capítulo

En este capítulo, utilizamos un bucle para realizar una computación repetitiva, actualizar un modelo 30 veces, y para calcular secuencias y series. Además, se utilizó la función plot para visualizar los resultados.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

El error absoluto es la diferencia entre una aproximación y una respuesta exacta. El error relativo es la misma diferencia expresada como una fracción o porcentaje de la respuesta exacta.

Un bucle es parte de un programa que se ejecuta repetidamente. Una variable de bucle es una variable a la que se le asigna un valor diferente cada vez que pasa por el bucle. Un rango es una secuencia de valores asignados a la variable de bucle, a menudo especificada con el operador de dos puntos, por ejemplo, 1:5. El cuerpo de un bucle es el conjunto de declaraciones dentro del bucle que se ejecuta repetidamente. Un acumulador es una variable que se utiliza para acumular un resultado poco a poco.

En matemáticas, una secuencia es un conjunto de números que corresponden a los enteros positivos. Se llama a los números que componen la secuencia. Una serie es la suma de una secuencia de elementos. A veces calculamos los elementos de una secuencia de manera recurrente, lo que significa que cada nuevo elemento depende de elementos anteriores. En ocasiones podemos calcular los elementos directamente, sin utilizar elementos anteriores.

La generalización es una forma de hacer que un programa sea más versátil, por ejemplo, reemplazando un valor específico por una variable que pueda tener cualquier valor. El desarrollo incremental es una forma de programación mediante la realización de una serie de pequeños cambios comprobables. El andamio es un código que escribes para ayudarte a programar o depurar pero que no forma parte del programa terminado.

En el siguiente capítulo, usaremos un vector para almacenar los resultados de un bucle.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Hace años estaba en una tienda de dulces y vi un letrero que decía “Compra una libra de caramelo, obtén otro cuarto de libra gratis”. Eso es bastante simple.

Pero si dirigía la tienda de caramelos, ofrecería un trato especial a cualquiera que pudiera resolver el siguiente problema:

Si compras una libra de dulce de azúcar, te daremos otro cuarto de libra gratis. Y luego te daremos un cuarto de cuarto de libra, o un decimosexto. Y luego te daremos una cuarta parte de eso, y así sucesivamente. ¿Cuánto fudge obtendrías en total?

Escribe un script llamado fudge.m que resuelva este problema. Pista: comenzar con series.m y generalizarlo reemplazando la relación 1/2 por una variable, r.

[fib2]

Ya hemos visto la secuencia de Fibonacci$F$, que se define de manera recurrente como

\[\mathrm{for}~i \ge 3, \quad F_{i} = F_{i-1} + F_{i-2}\]

Para comenzar, hay que especificar los dos primeros elementos, pero una vez que los tenga, puede calcular el resto. La secuencia de Fibonacci más común comienza con$F_1 = 1$ y$F_2 = 1$.

Escribe un script llamado fibonacci2.m que utilice un bucle for para calcular los primeros 10 elementos de esta secuencia de Fibonacci. Como postcondición, tu script debe asignar el décimo elemento a ans.

Ahora generaliza tu script para que calcule el elemento$n$ th para cualquier valor de n, con la condición previa de que tengas que establecer n antes de ejecutar el script. Para mantener las cosas simples por ahora, se puede suponer que n es mayor que 0.

Sugerencia: tendrás que usar dos variables para realizar un seguimiento de los dos elementos anteriores de la secuencia. Es posible que desee llamarlos prev1 y prev2. Inicialmente, prev1 =$F_1$ y prev2 =$F_2$. Al final del ciclo, tendrás que actualizar prev1 y prev2; ¡piensa detenidamente en el orden de las actualizaciones!

Vectores

En el capítulo anterior se utilizó un bucle para calcular los elementos de una secuencia, pero sólo pudimos almacenar el último elemento. En este capítulo, usaremos un vector para almacenar todos los elementos. También aprenderemos cómo seleccionar elementos de un vector, y cómo realizar operaciones aritméticas vectoriales y vectores comunes como reducir y aplicar.

Creación de vectores

Un vector en MATLAB es una secuencia de números. Hay varias formas de crear vectores; una de las más comunes es poner una secuencia de números entre corchetes:

>> [1 2 3]
ans = 1     2     3

Otra forma de crear un vector es el operador colon, que ya hemos utilizado para crear un rango de valores en un bucle for.

>> 1:3
ans = 1     2     3

En general, cualquier cosa que puedas hacer con un número, también puedes hacerlo con un vector. Por ejemplo, puede asignar un valor vectorial a una variable:

>> X = [1 2 3]
X = 1     2     3

Tenga en cuenta que las variables que contienen vectores suelen ser letras mayúsculas. Eso es solo una convención; MATLAB no lo requiere, pero es una forma útil de recordar qué variables son vectores.

Aritmética vectorial

Al igual que con los números, se puede hacer aritmética con vectores. Si agrega un número a un vector, MATLAB incrementa cada elemento del vector:

>> Y = X + 5
Y = 6     7     8

El resultado es un nuevo vector; el valor original de X no ha cambiado.

Si agrega dos vectores, MATLAB agrega los elementos correspondientes de cada vector y crea un nuevo vector que contiene las sumas:

>> Z = X + Y
Z = 7     9    11

Pero agregar vectores sólo funciona si los operandos son del mismo tamaño. De lo contrario se obtiene un error:

>> W = [1 2]
W = 1     2

>> X + W
Matrix dimensions must agree.

Este mensaje de error puede resultar confuso, porque pensamos en X y W como vectores, no matrices. Pero en MATLAB, un vector es una especie de matriz. Volveremos a las matrices más tarde, pero mientras tanto, es posible que veas el término en mensajes de error y documentación.

Si divide dos vectores, podría sorprenderse con el resultado:

>> X / Y
ans = 0.2953

MATLAB está realizando una operación llamada división derecha, que no es lo que esperábamos. Para dividir los elementos de X por los elementos de Y, hay que usar. /, que es la división por elementos:

>> X ./ Y
ans = 0.1667    0.2857    0.3750

La multiplicación tiene el mismo problema. Si usa *, MATLAB hace multiplicación matricial. Con estos dos vectores, la multiplicación matricial no está definida, por lo que se obtiene un error:

>> X * Y
Error using  *
Incorrect dimensions for matrix multiplication.
Check that the number of columns in the first matrix
matches the number of rows in the second matrix.
To perform element-wise multiplication, use '.*'.

En este caso, el mensaje de error es bastante útil. Como sugiere, puedes usar. * para realizar multiplicación por elementos:

>> X .* Y
ans = 6    14    24

Como ejercicio, vea qué sucede si usa el operador de exponenciación (^) con un vector.

Selección de elementos

Puede seleccionar un elemento de un vector con paréntesis:

>> Y = [6 7 8 9]
Y = 6    7     8     9

>> Y(1)
ans = 6

>> Y(4)
ans = 9

Esto significa que el primer elemento de Y es 6 y el cuarto elemento es 9. El número entre paréntesis se llama índice porque indica qué elemento del vector desea.

El índice puede ser un nombre de variable o una expresión matemática:

>> i = 1;
>> Y(i)
ans = 6
>> Y(i+1)
ans = 7

Podemos usar un bucle para mostrar los elementos de Y:

for i=1:4
     Y(i)
end

Cada vez que a través del bucle usamos un valor diferente de i como índice en Y.

En el ejemplo anterior tuvimos que conocer el número de elementos en Y. Podemos hacerlo más general usando la función length, que devuelve el número de elementos en un vector:

for i=1:length(Y)
     Y(i)
end

Esta versión funciona para un vector de cualquier longitud.

Errores de indexación

Un índice puede ser cualquier tipo de expresión, pero el valor de la expresión tiene que ser un entero positivo, y tiene que ser menor o igual que la longitud del vector. Si es cero o negativo, obtendrás un error:

>> Y(0)
Array indices must be positive integers or logical values.

Si no es un número entero, se obtiene un error:

>> Y(1.5)
Array indices must be positive integers or logical values.

Si el índice es demasiado grande, también obtienes un error:

>> Y(5)
Index exceeds the number of array elements (4).

Los mensajes de error usan la matriz de palabras en lugar de matriz, pero significan lo mismo, al menos por ahora.

Vectores y Secuencias

Los vectores y las secuencias van muy bien unidos. Por ejemplo, otra forma de evaluar la secuencia de Fibonacci del Capítulo 2 es almacenar valores sucesivos en un vector. Recuerde que la definición de la secuencia de Fibonacci es$F_1 = 1$,$F_2 = 1$, y$F_{i} = F_{i-1} + F_{i-2}$ para$i > 2$.

Listing [lst:fib_vec] muestra cómo podemos calcular los elementos de esta secuencia y almacenarlos en un vector, usando una letra mayúscula para el vector F y letras minúsculas para los enteros i y n.

F(1) = 1
F(2) = 1
for i=3:n
    F(i) = F(i-1) + F(i-2)
end

Si tuviste algún problema con, tienes que apreciar la sencillez de este guión. La sintaxis de MATLAB es similar a la notación matemática, lo que facilita la comprobación de la corrección.

Si solo quieres el número$n$ th Fibonacci, almacenar toda la secuencia desperdicia algo de espacio. Pero si desperdiciar espacio hace que tu código sea más fácil de escribir y depurar, probablemente esté bien.

Trazado de vectores

Si llama a plot con un vector como argumento, MATLAB traza los índices en el eje x y los elementos en el eje y. Para trazar los números de Fibonacci que calculamos en Listing [lst:fib_vec], usaríamos

plot(F)

La figura [fig:fibonacci] muestra el resultado.

Los primeros 10 elementos de la secuencia de Fibonacci

Esta forma de ver un vector suele ser útil para depurar, especialmente si es lo suficientemente grande como para que mostrar los elementos en la pantalla sea difícil de manejar.

Operaciones de vectores comunes

Hemos cubierto algunas de las características básicas de los vectores. Veamos ahora algunos patrones comunes que usamos para trabajar con datos almacenados en vectores.

Reducir

Frecuentemente usamos bucles para correr a través de los elementos de un vector y sumarlos, multiplicarlos juntos, calcular la suma de sus cuadrados, y así sucesivamente. Este tipo de operación se llama reducir, porque reduce un vector con múltiples elementos hasta un solo número.

Por ejemplo, el bucle en Listing [lst:vec_reduce] suma los elementos de un vector llamado X (que suponemos que se ha definido).

total = 0
for i=1:length(X)
    total = total + X(i)
end
ans = total

El uso del total como acumulador es similar a lo que vimos en el Capítulo 3. Nuevamente, usamos la función length para encontrar el límite superior del rango, por lo que este bucle funcionará independientemente de la longitud de X. Cada vez que pasa por el bucle, agregamos en el i-ésimo elemento de X, así al final del bucle total contiene la suma de los elementos.

MATLAB proporciona funciones que realizan algunas operaciones de reducción. Por ejemplo, la función sum calcula la suma de los elementos en un vector y prod calcula el producto.

Aplicar

Otro uso común de un bucle es correr a través de los elementos de un vector, realizar alguna operación sobre los elementos y crear un nuevo vector con los resultados. Esta operación se llama apply, porque se aplica la operación a cada elemento en el vector.

Por ejemplo, el bucle en Listing [lst:vec_apply] crea un vector Y que contiene los cuadrados de los elementos de X (suponiendo, nuevamente, que X ya está definido).

for i=1:length(X)
    Y(i) = X(i)^2
end

Muchas operaciones de aplicación se pueden hacer con operadores por elementos. La siguiente declaración es más concisa que el bucle en Listing [lst:vec_apply].

Y = X .^ 2

¡También corre más rápido!

Revisión del Capítulo

En este capítulo, se utilizó un vector para almacenar los elementos de una secuencia. Aprendimos a seleccionar elementos de un vector y realizar aritmética vectorial. Realizamos operaciones de reducción y aplicación usando bucles for, funciones de MATLAB y operaciones por elementos.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

Un vector es una secuencia de valores, que es una especie de matriz, también llamada matriz en alguna documentación de MATLAB.

Un índice es un valor entero utilizado para indicar uno de los elementos en un vector o matriz (también llamado subíndice en alguna documentación de MATLAB).

Una operación es por elementos si actúa sobre los elementos individuales de un vector o matriz (a diferencia de algunas operaciones de álgebra lineal).

Puede aplicar una operación a todos los elementos de un vector, y puede un vector a un solo valor, por ejemplo, calculando la suma de los elementos.

En el siguiente capítulo, conoceremos la idea más importante en programación informática: ¡funciones!

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Escribe un bucle que compute los primeros$n$ elementos de la secuencia geométrica$A_{i+1} = A_i/2$ con$A_1 = 1$. Observe que la notación matemática pone$A_{i+1}$ en el lado izquierdo de la igualdad. Cuando traduzca a MATLAB, es posible que desee reescribirlo con$A_{i}$ en el lado izquierdo.

Escribir una expresión que compute la raíz cuadrada de la suma de los cuadrados de los elementos de un vector, sin usar un bucle.

[fibratio]

La relación de números consecutivos de Fibonacci,$F_{n+1}/F_{n}$, converge a un valor constante a medida que$n$ aumenta. Escribir un script que compute un vector con los primeros$n$ elementos de una secuencia de Fibonacci (asumiendo que la variable n está definida) y luego calcula un nuevo vector que contiene las proporciones de números consecutivos de Fibonacci. Trace este vector para ver si parece converger. ¿En qué valor converge?

El siguiente conjunto de ecuaciones se basa en un famoso ejemplo de un sistema caótico, el atractor Lorenz (ver https://greenteapress.com/matlab/lorenz):\[\begin{aligned} x_{i+1} &=& x_i + \sigma \left( y_i - x_i \right) dt \\ y_{i+1} &=& y_i + \left[ x_i (r - z_i) - y_i \right] dt \\ z_{i+1} &=& z_i + \left( x_i y_i - b z_i \right) dt\end{aligned}\]

Escribe un script que compute los primeros 10 elementos de las secuencias$X$$Y$,$Z$ y los almacene en vectores llamados X, Y y Z.

Utilice los valores iniciales$X_1 = 1$,$Y_1 = 2$, y$Z_1 = 3$ con los valores$\sigma = 10$,$b = 8/3$,$r = 28$, y$dt = 0.01$.
Lee la documentación para plot3 y comet3, y traza los resultados en tres dimensiones.
Una vez que el código esté funcionando, use punto y coma para suprimir la salida y luego ejecute el programa con longitudes de secuencia de 100, 1,000 y 10,000.
Ejecute el programa nuevamente con diferentes condiciones de inicio. ¿Qué efecto tiene en el resultado?
Ejecute el programa con diferentes valores para$\sigma$$b$, y$r$, y vea si puedes tener una idea de cómo afecta cada variable al sistema.

El mapa logístico (ver https://greenteapress.com/matlab/logistic) se describe mediante la siguiente ecuación:

\[X_{i+1} = r X_i (1-X_i)\]

donde$X_i$ es un número entre 0 y 1, y$r$ es un número positivo.

Escribe un script llamado logmap.m que compute los primeros 50 elementos de$X$ con r = 3.9 y X1 = 0.5, donde r es el parámetro del mapa logístico y X1 es el valor inicial.
Trazar los resultados para un rango de valores$r$ de 2.4 a 4.0. ¿Cómo cambia el comportamiento del sistema a medida que varías$r$?

Trabajar con Condiciones

Algunos de los bucles del capítulo anterior no funcionan si el valor de n no se establece correctamente antes de que se ejecute el bucle. Por ejemplo, este bucle calcula la suma de los primeros n elementos de una secuencia geométrica:

A1 = 1;
total = 0;
for i=1:n
    a = A1 * 0.5^(i-1);
    total = total + a;
end
ans = total

Funciona para cualquier valor positivo de n, pero ¿y si n es negativo? En ese caso, obtienes:

total = 0

¿Por qué? Porque la expresión 1: -1 significa “todos los números del 1 al -1, contando hacia arriba por 1”. No es inmediatamente obvio lo que eso debería significar, pero la interpretación de MATLAB es que no hay ningún número que se ajuste a esa descripción, por lo que el resultado es

>> 1:-1
ans = 1x0 empty double row vector

El resultado es un vector vacío con una fila y 0 columnas. Si haces un bucle sobre un vector vacío, el bucle nunca se ejecuta en absoluto, razón por la cual en este ejemplo el valor de total es cero para cualquier valor negativo de n.

Si estás seguro de que nunca cometerás un error, y que las condiciones previas de tus funciones siempre estarán satisfechas, entonces no tienes que comprobarlo. Pero para el resto de nosotros, es peligroso escribir un guión, como este, que silenciosamente produzca la respuesta incorrecta (o al menos una respuesta sin sentido) si el valor de entrada es negativo. Una mejor alternativa es usar una sentencia if.

if Declaraciones

La sentencia if le permite verificar ciertas condiciones y ejecutar sentencias si se cumplen las condiciones. En el ejemplo anterior, podríamos escribir:

if n<0
    ans = NaN
end

La sintaxis es similar a un bucle for. La primera línea especifica la condición que nos interesa; en este caso estamos preguntando si n es negativo. Si lo es, MATLAB ejecuta el cuerpo de la sentencia, que es la secuencia sangrada de sentencias entre el if y el final.

MATLAB no requiere sangrar el cuerpo de una sentencia if, pero hace que su código sea más legible, por lo que debe hacerlo. No me hagas decírtelo otra vez.

En este ejemplo, lo “correcto” a hacer si n es negativo es establecer ans = NaN, que es una forma estándar de indicar que el resultado es indefinido (no un número).

Si no se cumple la condición, no se ejecutan las declaraciones en el cuerpo. A veces hay sentencias alternativas para ejecutar cuando la condición es falsa. En ese caso se puede extender la declaración if con una cláusula else.

La versión completa del ejemplo anterior podría verse así:

if n<0
    ans = NaN
else
    A1 = 1;
    total = 0;
    for i=1:n
        a = A1 * 0.5^(i-1);
        total = total + a;
    end
    ans = total
end

Las declaraciones como si y para eso contienen otras declaraciones se denominan declaraciones compuestas. Todas las declaraciones compuestas terminan con... final.

En este ejemplo, una de las declaraciones de la cláusula else es un bucle for. Poner una declaración compuesta dentro de otra es legal y común, y a veces se llama anidación.

Operadores Relacionales

Los operadores que comparan valores, como se llaman operadores relacionales porque prueban la relación entre dos valores. El resultado de un operador relacional es uno de los valores lógicos: ya sea 1, que representa “verdadero”, o 0, que representa “falso”.

Los operadores relacionales suelen aparecer en declaraciones if, pero también puede evaluarlos en el prompt:

>> x = 5;
>> x < 10
ans = 1

Puede asignar un valor lógico a una variable:

>> flag = x > 10
flag = 0

Una variable que contiene un valor lógico a menudo se llama bandera porque marca el estado de alguna condición.

Los otros operadores relacionales son =, que se explican por sí mismos, ==, para “igual” y ~=, para “no igual”.

No olvides que == es el operador que prueba la igualdad, y = es el operador de asignación. Si intentas usar = en una instrucción if, obtienes un error:

>> if x=5
 if x=5
     |
Error: Incorrect use of '=' operator. 
To assign a value to a variable, use '='.
To compare values for equality, use '=='.
 
Did you mean:
>> x = 5

En este caso, el mensaje de error es bastante útil.

Operadores lógicos

Para probar si un número cae en un intervalo, es posible que tengas la tentación de escribir algo como 0 < x < 10, pero eso estaría mal, así que muy mal. Desafortunadamente, en muchos casos, obtendrás la respuesta correcta por la razón equivocada. Por ejemplo:

>> x = 5;
>> 0 < x < 10            % right for the wrong reason
ans = 1

¡Pero no te dejes engañar!

>> x = 17
>> 0 < x < 10            % just plain wrong
ans = 1

El problema es que MATLAB está evaluando los operadores de izquierda a derecha, así que primero comprueba si 0<x. Lo es, entonces el resultado es 1. Entonces compara el valor lógico 1 (no el valor de x) con 10. Desde 1<10, el resultado es verdadero, aunque x no esté en el intervalo.

Para los programadores principiantes, ¡esto es un bicho malvado, malvado!

Una forma de evitar este problema es usar una instrucción if anidada para verificar las dos condiciones por separado:

ans = 0
if x>0
    if x<10
        ans = 1
    end
end

Pero es más conciso usar el operador AND, &&, para combinar las condiciones.

>> x = 5;
>> x>0 && x<10
ans = 1

>> x = 17;
>> x>0 && x<10
ans = 0

El resultado de AND es cierto si ambos operandos son ciertos. El operador OR, ||, es verdadero si uno o ambos operandos son verdaderos.

...

No obstante, hay que tener cuidado con el alcance del bucle. En la versión anterior, el bucle va de 3 a n, y cada vez asignamos un valor al i-ésimo elemento.

También funcionaría para “desplazar” el índice por dos, ejecutando el bucle de 1 a n-2:

F(1) = 1
F(2) = 1
for i=1:n-2
    F(i+2) = F(i+1) + F(i)
end

Cualquiera de las dos versiones está bien, pero hay que elegir un enfoque y ser consistente. Si combinas elementos de ambos, te confundirás. Prefiero la versión en Listing [lst:fib_vec] que tiene F (i) en el lado izquierdo de la asignación, de manera que cada vez que pasa por el bucle asigna el i-ésimo elemento.

Todo es una Matriz

En matemáticas (específicamente en álgebra lineal) un vector es unidimensional y una matriz es bidimensional.

Pero en MATLAB, todo es una matriz (excepto cadenas). Esto se puede ver si usa el comando whos para mostrar las variables en el espacio de trabajo. quién es similar a quién, pero también muestra el tamaño de cada valor y otra información.

Para demostrarlo, voy a hacer uno de cada tipo de valor:

>> number = 5
number = 5

>> vector = [1 2 3 4 5]
vector = 1     2     3     4     5

>> matrix = ones(2,3)
matrix =
     1     1     1
     1     1     1

Las funciones integradas construyen una nueva matriz con el número dado de filas y columnas, y establece todos los elementos en 1. Ahora veamos qué tenemos.

>> whos
  Name        Size            Bytes  Class     Attributes

  number      1x1                 8  double              
  vector      1x5                40  double               
  matrix      2x3                48  double

Según MATLAB, las tres son matrices. Tienen la misma clase, doble, lo que significa que contienen números de coma flotante de doble precisión.

Pero tienen tamaños de diferencia: el tamaño del número es 1x1, lo que significa que tiene 1 fila y 1 columna; vector tiene 1 fila y 5 columnas; y matriz tiene 2 filas y 3 columnas.

Entonces, puedes pensar en tus valores como números, vectores y matrices, pero en MATLAB son todas matrices.

Buscar

Otro uso más de los bucles es buscar los elementos de un vector y devolver el índice del valor que estás buscando (o el primer valor que tiene una propiedad en particular).

Por ejemplo, el bucle en Listing [lst:vec_search] encuentra el índice del elemento 0 en X:

for i=1:length(X)
    if X(i) == 0
        ans = i
    end
end

Algo curioso de este bucle es que sigue persiguiendo que encuentra lo que está buscando. Eso podría ser lo que quieres; si el valor objetivo aparece más de uno, este bucle proporciona el índice del último.

Pero si quieres el índice del primero (o sabes que solo hay uno), puedes guardar algún bucle innecesario usando la sentencia break.

for i=1:length(X)
    if X(i) == 0
        ans = i
        break
    end
end

break hace más o menos lo que parece. Termina el bucle y pasa inmediatamente a la siguiente declaración después del bucle (en este caso, no hay una, así que el código termina).

Funciones

Este capítulo introduce la idea más importante en la programación informática:! Para explicar por qué las funciones son tan importantes, comenzaré explicando uno de los problemas que resuelven: las colisiones de nombres.

Colisiones de nombres

Todos los scripts se ejecutan en el mismo espacio de trabajo, por lo que si un script cambia el valor de una variable, todos los demás scripts ven el cambio. Con una pequeña cantidad de scripts simples, eso no es un problema, pero eventualmente las interacciones entre scripts se vuelven inmanejables.

Por ejemplo, el siguiente script calcula la suma de los primeros n términos en una secuencia geométrica, pero también tiene el efecto secundario de asignar valores a A1, total, i y a.

A1 = 1;
total = 0;
for i=1:10
    a = A1 * 0.5^(i-1);
    total = total + a;
end
ans = total

Si estaba usando alguno de esos nombres de variables antes de llamar a este script, podría sorprenderse al descubrir, después de ejecutar el script, que sus valores habían cambiado. Si tienes dos scripts que usan los mismos nombres de variables, es posible que encuentres que funcionan por separado y luego se rompen cuando intentas combinarlos. Este tipo de interacción se llama colisión de nombres.

A medida que aumenta el número de scripts que escribes, y se hacen más largos y complejos, las colisiones de nombres se vuelven más problemáticas. Evitar este problema es una de varias motivaciones para las funciones.

Definición de Funciones

Una función es como un script, excepto que cada función tiene su propio espacio de trabajo, por lo que cualquier variable definida dentro de una función solo existe mientras la función se está ejecutando y no interfiere con variables en otros espacios de trabajo, aunque tengan el mismo nombre. Las entradas y salidas de funciones se definen cuidadosamente para evitar interacciones inesperadas.

Para definir una nueva función, se crea un archivo M con el nombre que desea y le pone una definición de función. Por ejemplo, para crear una función llamada myfunc, cree un archivo M llamado myfunc.m y ponga la siguiente definición en él (Listing [lst:function_def]):

function res = myfunc(x)
    s = sin(x)
    c = cos(x)
    res = abs(s) + abs(c)
end

La primera palabra sin comentario del archivo tiene que ser función, porque así es como MATLAB dice la diferencia entre un script y un archivo de función.

Una definición de función es una declaración compuesta. La primera línea se llama la firma de la función; define las entradas y salidas de la función. En Listing [lst:function_def] la variable de entrada se llama x. Cuando se llama a esta función, el argumento proporcionado por el usuario se asignará a x.

La variable de salida se denomina res, que es la abreviatura de result. Puedes llamar a la variable de salida como quieras, pero como convención, me gusta llamarla res. Por lo general, lo último que hace una función es asignar un valor a la salida.

Una vez que haya definido una nueva función, la llama de la misma manera que llama a las funciones integradas de MATLAB. Si llama a la función como una declaración, MATLAB pone el resultado en ans:

>> myfunc(1)

s = 0.84147098480790

c = 0.54030230586814

res = 1.38177329067604

ans = 1.38177329067604

Pero es más común (y mejor estilo) asignar el resultado a una variable:

>> y = myfunc(1)

s = 0.84147098480790

c = 0.54030230586814

res = 1.38177329067604

y = 1.38177329067604

Mientras está depurando una nueva función, es posible que desee mostrar resultados intermedios como este, pero una vez que esté funcionando, querrá agregar punto y coma para que sea una función silenciosa. Una función silenciosa calcula un resultado pero no muestra nada (excepto a veces mensajes de advertencia). La mayoría de las funciones integradas son silenciosas.

Cada función tiene su propio espacio de trabajo, que se crea cuando se inicia la función y se destruye cuando termina la función. Si intentas acceder (leer o escribir) las variables definidas dentro de una función, encontrarás que no existen.

>> clear
>> y = myfunc(1);
>> who
Your variables are: y

>> s
Undefined function or variable 's'.

El único valor de la función a la que se puede acceder es el resultado, que en este caso se le asigna a y.

Si tienes variables llamadas s o c en tu espacio de trabajo antes de llamar a myfunc, seguirán ahí cuando se complete la función.

>> s = 1;
>> c = 1;
>> y = myfunc(1);
>> s, c

s = 1
c = 1

Así que dentro de una función puedes usar cualquier nombre de variable que quieras sin preocuparte por colisiones.

Documentación de funciones

Al inicio de cada archivo de función, debe incluir un comentario que explique lo que hace la función:

% res = myfunc(x)
% Compute the Manhattan distance from the origin to the
% point on the unit circle with angle (x) in radians.

function res = myfunc(x)
% this is not part of documentation given by help function

    s = sin(x);
    c = cos(x);
    res = abs(s) + abs(c);
end

Cuando pides ayuda, MATLAB imprime el comentario que proporcionas.

>> help myfunc
  res = myfunc(x)
  Compute the Manhattan distance from the origin to the
  point on the unit circle with angle (x) in radians.

Hay muchas convenciones sobre lo que se debe incluir en estos comentarios. Entre otras cosas, es una buena idea incluir lo siguiente:

Firma: La firma de la función, que incluye el nombre de la función, las variables de entrada y las variables de salida.
Descripción: Una descripción clara, concisa, abstracta de lo que hace la función. Una descripción abstracta es aquella que deja de lado los detalles de cómo funciona la función e incluye solo información que alguien que usa la función necesita conocer. Puedes poner comentarios adicionales dentro de la función que expliquen los detalles.
Variables: Una explicación de lo que significan las variables de entrada; por ejemplo, en este caso es importante señalar que x se considera un ángulo en radianes.
Condiciones: Cualquier condición previa y postcondición.

Funciones de Nombramiento

Hay algunas “gotchas” que surgen cuando comienzas a definir funciones. La primera es que el nombre “real” de tu función está determinado por el nombre del archivo, no por el nombre que pones en la firma de la función. Por cuestión de estilo, debes asegurarte de que siempre sean los mismos, pero si cometes un error, o si cambias el nombre de una función, es fácil confundirte.

En el espíritu de cometer errores a propósito, edita myfunc.m y cambia el nombre de la función de myfunc a something_else, así:

function res = something_else (x)
    s = sin(x);
    c = cos(x);
    res = abs(s) + abs(c);
end

Ahora llama a la función desde la Ventana de Comandos, así:

>> y = myfunc(1)
y = 1.3818

La función todavía se llama myfunc, porque ese es el nombre del archivo. Si intentas llamarlo así:

>> y = something_else(1)
Undefined function or variable 'something_else'.

No funciona. El nombre del archivo es lo que importa; se ignora el nombre de la función.

El segundo “gotcha” es que el nombre del archivo no puede tener espacios. Por ejemplo, si cambias el nombre del archivo a mi func.m e intentas ejecutarlo, obtienes:

>> y = my func(1)
 y = my func(1)
        |
Error: Unexpected MATLAB expression.

Esto falla porque MATLAB piensa que my y func son dos nombres de variables diferentes.

El tercer “gotcha” es que los nombres de sus funciones pueden colisionar con las funciones integradas de MATLAB. Por ejemplo, si crea un archivo M llamado sum.m y luego llama a sum, MATLAB podría llamar a su nueva función, ¡no a la versión incorporada! A cuál se llama realmente depende del orden de los directorios en la ruta de búsqueda y (en algunos casos) de los argumentos. Como ejemplo, pon el siguiente código en un archivo llamado sum.m:

function res = sum(x)
   res = 7;
end

Y luego prueba esto:

>> sum(1:3)

ans = 6

>> sum

ans = 7

En el primer caso MATLAB utilizó la función incorporada; en el segundo caso ejecutó tu función! Este tipo de interacción puede ser muy confusa. Antes de crear una nueva función, verifique si ya hay una función MATLAB con el mismo nombre. Si lo hay, ¡elige otro nombre!

Variables de entrada múltiple

Las funciones pueden tomar más de una variable de entrada. Por ejemplo, la siguiente función en Listing [lst:hyp_function] toma dos variables de entrada, a y b:

function res = sum_squares(a, b)
    res = a^2 + b^2;
end

Esta función calcula la suma de cuadrados de dos números, a y b.

Si lo llamamos desde la Ventana de Comandos con los argumentos 3 y 4, podemos confirmar que la suma de sus cuadrados es 25.

>> ss = sum_squares(3, 4)
ss = 25

Los argumentos que proporcione se asignan a las variables de entrada en orden, por lo que en este caso 3 se asigna a a y 4 se asigna a b. MATLAB comprueba que usted proporciona el número correcto de argumentos; si proporciona muy pocos, obtiene

>> ss = sum_squares(3)
Not enough input arguments.

Error in sum_squares (line 4)
    res = a^2 + b^2;

Este mensaje de error puede ser confuso, porque sugiere que el problema está en sum_squares en lugar de en la llamada a la función. Téngalo en mente cuando esté depurando.

Si proporcionas demasiados argumentos, obtienes

ss = sum_squares(3, 4, 5)
Error using sum_squares
Too many input arguments.

Ese es un mejor mensaje de error, porque está claro que el problema no está en la función, está en la forma en que estamos usando la función.

Revisión del Capítulo

Ahora que conocemos las funciones, y todas las formas en que pueden salir mal, vamos a darles un buen uso. En el siguiente capítulo desarrollaremos un programa que utiliza varias funciones para buscar triples pitagóricos (y voy a explicar cuáles son esos).

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

Una función es una secuencia con nombre de instrucciones almacenadas en un archivo M. Una función puede tener una o más variables de entrada, que obtienen sus valores cuando se llama a la función, y variables de salida, que devuelven un valor de la función a la persona que llama.

La primera línea de una definición de función es su firma, que especifica el nombre de la función, las variables de entrada y las variables de salida.

Una función silenciosa no muestra nada, genera una figura ni tiene ningún otro efecto que no sea devolver valores de salida.

Ejercicio

Antes de continuar, es posible que desee trabajar en el siguiente ejercicio.

[Hypotenuse_exercise] Escribe una función llamada hipotenusa que tome dos parámetros, a y b, que representen las longitudes de dos lados de un triángulo rectángulo. Se debe asignar a res la longitud del tercer lado del triángulo, dada por la fórmula

\[c = \sqrt{a^2 + b^2}\]

Condicionales

En este capítulo, usaremos funciones y una nueva característica, declaraciones condicionales, para buscar triples pitagóricos. Un triple pitagórico es un conjunto de enteros, como 3, 4 y 5, que son las longitudes de los lados de un triángulo rectángulo. Matemáticamente, es un conjunto de enteros$a$,$b$, y$c$ tal que$a^2 + b^2 = c^2$. Este ejemplo también demostrará el proceso de desarrollo incremental del que hablamos en Capítulo.

Operadores Relacionales

Supongamos que tenemos tres variables, a, b y c, y queremos comprobar si forman un triple pitagórico. Podemos usar el operador de igualdad (==) para comparar dos valores:

>> a = 3;
>> b = 4;
>> c = 5;
>> a^2 + b^2 == c^2

ans = logical 1

El resultado es un valor lógico, lo que significa que es 1, que significa “verdadero”, o 0, que significa “falso”. Aquí hay un ejemplo donde el resultado es falso:

>> c = 6;
>> a^2 + b^2 == c^2
ans = logical 0

Es un error común usar el operador de asignación (=) en lugar del operador de igualdad (==). Si lo haces, obtienes un error:

>> a^2 + b^2 = c^2
 a^2 + b^2 = c^2
           |
Error: Incorrect use of '=' operator.
To assign a value to a variable, use '='.
To compare values for equality, use '=='.

El operador de igualdad es uno de varios operadores relacionales, llamados así porque prueban las relaciones entre valores. Por ejemplo, x < 10 es true (1) si el valor de x es menor que 10 o false (0) si es contrario. Y x > 0 es verdadero si x es mayor que 0.

Los otros operadores relacionales son <= para “menor o igual”, >= para “mayor o igual” y ~= para “no igual”.

if Declaración

Ahora supongamos que cuando encontremos un triple pitagórico queremos mostrar un mensaje. La sentencia if le permite verificar ciertas condiciones y ejecutar sentencias si se cumplen las condiciones. Por ejemplo:

if a^2 + b^2 == c^2
    disp("Yes, that is a Pythagorean triple.")
end

La sintaxis es similar a un bucle for. La primera línea especifica la condición que nos interesa. Si la condición es verdadera, MATLAB ejecuta el cuerpo de la sentencia, que es la secuencia sangrada de sentencias entre el if y el final.

MATLAB no requiere sangrar el cuerpo de una sentencia if, pero hace que su código sea más legible, por lo que debe hacerlo.

Si no se cumple la condición, las declaraciones en el cuerpo no lo son.

A veces hay sentencias alternativas para ejecutar cuando la condición es falsa. En ese caso, puede extender la declaración if con una cláusula else.

La versión completa del ejemplo anterior podría verse así:

if a^2 + b^2 == c^2
    disp("Yes, that is a Pythagorean triple.")
else
    disp("No, that is not a Pythagorean triple.")
end

Las declaraciones como si y para eso contienen otras declaraciones se denominan declaraciones compuestas. Todas las declaraciones compuestas terminan con el final.

Desarrollo Incremental

Ahora que tenemos operadores relacionales y si declaraciones, comencemos a escribir el programa.

Estos son los pasos que seguiremos para desarrollar el programa de manera incremental:

Escribe un script llamado DelmarFindTriples.m y comienza con un bucle que enumera los valores de a y los muestra.
Escribe un segundo bucle que enumera los valores de b y un tercer bucle que enumera los valores de c.
Utilice la declaración if de la sección anterior para verificar si a, b y c forman un triple pitagórico.
Mostrar los valores que pasan la prueba.
Transforme el script en una función y haga que tome una variable de entrada que especifique el rango a buscar.

En el camino, optimizaremos el programa para eliminar el trabajo innecesario.

Funciones lógicas

El primer paso es crear una función lógica, que es una función que devuelve un valor lógico. La siguiente función toma tres variables de entrada, a, b y c, y devuelve true (1) si forman un triple pitagórico y false (0) en caso contrario.

function res = is_pythagorean(a, b, c)
    if a^2 + b^2 == c^2
        res = 1;
    else
        res = 0;
    end
end

Podemos usar esta función así:

>> is_pythagorean(3, 4, 5)
ans = 1

Pero podemos escribir la misma función de manera más concisa, así:

function res = is_pythagorean(a, b, c)
    res = a^2 + b^2 == c^2;
end

El resultado del operador de igualdad es un valor lógico, que podemos asignar directamente a res.

Poner esta función en un archivo llamado is_pythagorean.m, para que podamos usarla como parte de nuestro programa.

Bucles anidados

El siguiente paso es escribir bucles que enumeren diferentes valores de a, b y c. Crea un nuevo archivo llamado find_triples.m donde desarrollaremos el resto del programa.

Empezaremos con un bucle para un:

for a=1:3
    a
end

Puede parecer una tontería comenzar con un programa tan simple, pero este es un elemento esencial del desarrollo incremental: empezar simple y probar a medida que avanza.

La salida es como se esperaba.

1
2
3

Ahora agregaremos un segundo bucle para b. Podría ser tentador escribir algo como esto:

for a=1:3
    disp(a)
end
for b=1:4
    disp(b)
end

Pero eso recorre los valores de a y luego recorre los valores de b, y eso no es lo que queremos.

En cambio, queremos considerar cada par de valores posibles, así:

for a=1:3
    for b=1:4
        disp([a,b])
    end
end

Ahora un bucle está dentro del otro. El bucle interno se ejecuta tres veces, una por cada valor de a, así que así es como se ve la salida (he ajustado el espaciado para aclarar la estructura):

>> find_triples
     1     1
     1     2
     1     3
     1     4
     2     1
     2     2
     2     3
     2     4
     3     1
     3     2
     3     3
     3     4

La columna izquierda muestra los valores de a y la columna derecha muestra los valores de b.

El siguiente paso es buscar valores de c que puedan hacer un triple pitagórico. El mayor valor posible para c es a + b, porque de lo contrario no podríamos formar un triángulo (ver https://greenteapress.com/matlab/triangle).

for a=1:3
    for b=1:4
        for c=1:a+b
            disp([a,b,c])
        end
    end
end

Después de cada pequeño cambio, vuelva a ejecutar el programa y verifique la salida.

Poniéndolo Juntos

Ahora, en lugar de mostrar todos los triples, agregaremos una sentencia if y mostraremos solo triples pitagóricos:

for a=1:3
    for b=1:4
        for c=1:a+b
            if is_pythagorean(a, b, c)
                disp([a,b,c])
            end
        end
    end
end

El resultado es solo un triple:

>> find_triples
     3     4     5

Podrías notar que estamos desperdiciando un poco de esfuerzo aquí. Después de verificar el caso cuando a es 1 y b es 2, no tiene sentido verificar el caso cuando a es 2 y b es 1. Podemos ahorrar el trabajo extra ajustando el rango de b:

for b=a:4

Podemos ahorrar aún más trabajo ajustando el rango de c:

for c=b:a+b

Aquí está la versión final:

for a=1:3
    for b=a:4
        for c=b:a+b
            if is_pythagorean(a, b, c)
                disp([a,b,c])
            end
        end
    end
end

Encapsulación y Generalización

Como script, este programa tiene el efecto secundario de asignar valores a a, b y c, lo que sería malo si alguno de esos nombres estuviera en uso. Al envolver el código en una función, podemos evitar colisiones de nombres; este proceso se llama encapsulación porque aísla este programa del espacio de trabajo.

El primer borrador de la función no toma variables de entrada:

function res = find_triples()
    for a=1:3
        for b=a:4
            for c=b:a+b
                if is_pythagorean(a, b, c)
                    disp([a,b,c])
                end
            end
        end
    end
end

Los paréntesis vacíos en la firma no son necesarios, pero hacen evidente que no hay variables de entrada. Del mismo modo, es una buena idea al llamar a la nueva función para usar paréntesis como recordatorio de que es una función, no un script:

>> find_triples()

Tampoco es necesaria la variable de salida; nunca se le asigna un valor. Pero lo puse ahí como cuestión de hábito y así mis firmas de funciones todas tienen la misma estructura.

El siguiente paso es generalizar esta función agregando variables de entrada. La generalización natural consiste en reemplazar los valores constantes 3 y 4 por una variable para que podamos buscar un rango arbitrariamente grande de valores.

function res = find_triples(n)
    for a=1:n
        for b=a:n
            for c=b:a+b
                if is_pythagorean(a, b, c)
                    disp([a,b,c])
                end
            end
        end
    end
end

Aquí están los resultados para el rango de 1 a 15:

>> find_triples(15)
     3     4     5
     5    12    13
     6     8    10
     8    15    17
     9    12    15

Los triples$5,12,13$ y$8,15,17$ son nuevos, pero los otros son solo múltiplos del$3,4,5$ triángulo.

Agregar una declaración continue

Como mejora final, modifiquemos la función para que solo muestre el “más bajo” de cada triple pitagórico, y no los múltiplos.

La forma más sencilla de eliminar los múltiplos es verificar si a y b comparten un factor común. Si lo hacen, al dividir ambos por el factor común se obtiene un triángulo más pequeño y similar que ya se ha comprobado.

MATLAB proporciona una función gcd que calcula el mayor divisor común de dos números. Si gcd (a, b) es mayor que 1, a y b comparten un factor común y podemos usar la sentencia continue para saltar al siguiente par. Listing [lst:triples_function] contiene la versión final de esta función:

function res = find_triples(n)
    for a=1:n
        for b=a:n
            for c=b:a+b
                if gcd(a,b) > 1
                    continue
                end
                if is_pythagorean(a, b, c)
                    disp([a,b,c])
                end
            end
        end
    end
end

La instrucción continue hace que el programa termine la iteración actual inmediatamente, salte a la parte superior del bucle y “continue” con la siguiente iteración.

En este caso, como hay tres bucles, puede que no sea obvio a qué bucle saltar, pero la regla es saltar al bucle más interno (que es lo que queremos).

Aquí están los resultados con n = 40:

>> find_triples(40)
     3     4     5
     5    12    13
     7    24    25
     8    15    17
     9    40    41
    12    35    37
    20    21    29

Cómo funcionan las funciones

Revisemos la secuencia de pasos que ocurren cuando llamas a una función:

Antes de que la función comience a ejecutarse, MATLAB crea un nuevo espacio de trabajo para ella.
MATLAB evalúa cada uno de los argumentos y asigna los valores resultantes, en orden, a las variables de entrada (que viven en el nuevo espacio de trabajo).
Se ejecuta el cuerpo del código. En algún lugar del cuerpo se asigna un valor a la variable de salida.
Se destruye el espacio de trabajo de la función; lo único que queda es el valor de la variable de salida y cualquier efecto secundario que tenga la función (como mostrar valores).
El programa se reanuda desde donde lo dejó. El valor de la llamada a la función es el valor de la variable de salida.

Cuando estás leyendo un programa y llegas a una llamada de función, hay dos formas de interpretarlo. Puedes pensar en el mecanismo que acabo de describir, y seguir la ejecución del programa en la función y volver, o puedes asumir que la función funciona correctamente, y pasar a la siguiente declaración después de la llamada a la función.

Cuando usas una función incorporada, es natural asumir que funciona, en parte porque normalmente no tienes acceso al código en el cuerpo de la función. Pero cuando empiezas a escribir tus propias funciones, podrías encontrarte siguiendo el “flujo de ejecución”. Esto puede ser útil mientras estás aprendiendo, pero a medida que adquieras experiencia, deberías sentirte más cómodo con la idea de escribir una función, probarla para asegurarte de que funciona y luego olvidarte de los detalles de cómo funciona.

Olvidar los detalles se llama abstracción; en el contexto de las funciones, abstracción significa olvidarse de cómo funciona una función y simplemente asumir (después de pruebas apropiadas) que funciona. Para muchas personas, lleva algún tiempo sentirse cómodo con las funciones. Si eres uno de ellos, podrías sentirte tentado a evitar funciones, y en ocasiones puedes arreglártelas sin ellas.

Pero los programadores experimentados utilizan funciones extensamente, por varias buenas razones. Primero, cada función tiene su propio espacio de trabajo, por lo que el uso de funciones ayuda a evitar colisiones de nombres. Las funciones también se prestan al desarrollo incremental: puedes depurar primero el cuerpo de la función (como script), luego encapsularlo como una función y luego generalizarlo agregando variables de entrada.

Además, las funciones permiten dividir un problema grande en piezas pequeñas, trabajar en las piezas una a la vez y luego armar una solución completa.

Una vez que tengas una función funcionando, podrás olvidarte de los detalles de cómo funciona y concentrarte en lo que hace. Este proceso de abstracción es una herramienta importante para gestionar la complejidad de los programas grandes.

Revisión del Capítulo

En este capítulo, encontramos operadores relacionales y declaraciones if, y los usamos para desarrollar un programa que busca triples pitagóricos. Escribimos una función lógica, que es una función que devuelve un valor lógico (1 para “true” o 0 para “false”).

También vimos un ejemplo de desarrollo incremental, o desarrollando programas gradualmente, agregando solo unas pocas líneas de código a la vez y probando a medida que avanza. Si desarrollas programas de esta manera, tendrás menos errores y los encontrarás más rápidamente.

Este capítulo definió dos nuevos términos: encapsulación es el proceso de envolver parte de un programa en una función con el fin de limitar las interacciones (incluidas las colisiones de nombres) entre la función y el resto del programa; la abstracción es el proceso de ignorar los detalles de cómo una función funciona con el fin de enfocarse en un modelo más simple de lo que hace la función.

El siguiente capítulo introduce una nueva herramienta, llamada fzero, que utilizaremos para resolver ecuaciones no lineales.

Ejercicio

Antes de continuar, es posible que desee trabajar en el siguiente ejercicio.

Existe una interesante conexión entre los números de Fibonacci y los triples pitagóricos. Si$F$ es una secuencia de Fibonacci, entonces

\[\big(F_i F_{i+3}, \, 2 F_{i+1} F_{i+2}, \, F_{i+1}^2 + F_{i+2}^2 \big)\]

es un triple pitagórico, para todos$i \ge 1$.

Escribe una función llamada fib_triple que tome n como variable de entrada, calme los primeros n números de Fibonacci, los almacene en un vector y verifique si esta fórmula produce triples pitagóricos para los números en el.

Búsqueda cero

En este capítulo usaremos la función fzero de MATLAB para resolver ecuaciones no lineales. Las ecuaciones no lineales son útiles para modelar sistemas físicos; por ejemplo, en uno de los ejercicios al final de este capítulo, se puede utilizar fzero para encontrar el punto de equilibrio de un objeto flotando sobre el agua. Usar fzero también es una oportunidad para aprender sobre los manejadores de funciones, que necesitaremos para los siguientes capítulos.

Resolver ecuaciones no lineales

¿Qué significa “resolver” una ecuación? Esa puede parecer una pregunta obvia, pero tomemos un minuto para pensarlo, empezando por una simple.

Supongamos que queremos saber el valor de una variable,$x$, pero todo lo que sabemos de ella es la relación$x^2 = a$. Si has tomado álgebra, probablemente sepas cómo resolver esta ecuación: tomas la raíz cuadrada de ambos lados y obtienes$x = \pm \sqrt{a}$. Entonces, con la satisfacción de un trabajo bien hecho, pasas al siguiente problema.

Pero, ¿qué has hecho realmente? La relación que derivaste es equivalente a la relación con la que iniciaste —contienen la misma información sobre$x$ — entonces, ¿por qué es preferible la segunda a la primera?

Hay dos razones. Una es que la relación ahora es explícita en$x$: porque$x$ está completamente sola en el lado izquierdo, podemos tratar al lado derecho como una receta para la computación$x$, asumiendo que conocemos el valor de$a$.

La otra razón es que la receta está escrita en términos de operaciones que sabemos realizar. Suponiendo que sabemos calcular raíces cuadradas, podemos calcular el valor de$x$ para cualquier valor de$a$.

Cuando la gente habla de resolver una ecuación, lo que suele significar es algo así como “encontrar una relación equivalente que sea explícita en una de las variables”. En el contexto de este libro, eso es lo que llamaremos un analítico, para distinguirlo de una solución numérica, que es lo que vamos a hacer a continuación.

= = Para demostrar una solución numérica, considere la ecuación$x^2 - 2x = 3$. You could solve this analytically, either by factoring it or by using the quadratic formula, and you would discover that there are two solutions, $x=3$ and $x=-1$. Alternatively, you could solve it numerically by rewriting it as $x = \pm \sqrt{2x+3}$.

Esta ecuación no es explícita, ya que$x$ aparece en ambos lados, por lo que no está claro que esta jugada haya hecho algún bien en absoluto. Pero supongamos que teníamos razones para esperar una solución cerca del 4. Podríamos comenzar con$x=4$ como un valor inicial y luego usar la ecuación$x = \sqrt{2x+3}$ para calcular aproximaciones sucesivas de la solución. (Para entender por qué funciona esto, consulte https://greenteapress.com/matlab/fixed.)

Esto es lo que sucede:

>> x = 4;
>> x = sqrt(2*x+3)
x = 3.3166

>> x = sqrt(2*x+3)
x = 3.1037

>> x = sqrt(2*x+3)
x = 3.0344

>> x = sqrt(2*x+3)
x = 3.0114

>> x = sqrt(2*x+3)
x = 3.0038

Después de cada iteración, x está más cerca de la respuesta correcta, y después de cinco iteraciones el error relativo es de aproximadamente 0.1 por ciento, lo que es lo suficientemente bueno para la mayoría de los propósitos.

Las técnicas que generan soluciones numéricas se denominan numéricas. Lo bueno del método que acabamos de usar es que es sencillo. Pero no siempre funciona, y no se usa a menudo en la práctica. Veremos una mejor alternativa en la siguiente sección.

Búsqueda cero

La función fzero de MATLAB que utiliza métodos numéricos para buscar soluciones a ecuaciones no lineales. Para poder usarlo, tenemos que reescribir la ecuación como una función de error, así:

\[f(x) = x^2 - 2x -3\]

El valor de la función error es 0 si$x$ es una solución y distinto de cero si no lo es. Esta función es útil porque podemos utilizar valores de$f(x)$, evaluados a diversos valores de$x$, para inferir la ubicación de las soluciones. Y eso es lo que hace fzero. Los valores de$x$ donde$f(x) = 0$ se llaman ceros de la función o raíces.

Para usar fzero hay que definir una función MATLAB que calcule la función error, así:

function res = error_func(x)
    res = x^2 - 2*x -3;
end

Puedes llamar a error_func desde la Ventana de Comandos y confirmar que$3$ y$-1$ son ceros:

>> error_func(3)
ans = 0

>> error_func(-1)
ans = 0

Pero pretendamos que no sabemos dónde están las raíces; sólo sabemos que una de ellas está cerca de 4. Entonces podríamos llamar a fzero así:

>> fzero(@error_func, 4)
ans = 3.0000

¡Éxito! Encontramos uno de los ceros.

El primer argumento es un manejador de función que especifica la función de error. El símbolo @ nos permite nombrar la función sin llamarla. Lo interesante aquí es que en realidad no estás llamando directamente a error_func; solo le estás diciendo a fzero dónde está. A su vez, fzero llama a tu función de error, más de una vez, de hecho.

El segundo argumento es el valor inicial. Si proporcionamos un valor diferente, obtenemos una raíz diferente (al menos a veces).

>> fzero(@error_func, -2)
ans = -1

Alternativamente, si conoces dos valores que paréntesis la raíz, puedes proporcionar ambos.

>> fzero(@error_func, [2,4])
ans = 3

El segundo argumento es un vector que contiene dos elementos.

Quizás tengas curiosidad por saber cuántas veces fzero llama a tu función y dónde. Si modificas error_func para que muestre el valor de x cuando se llama y luego vuelvas a ejecutar fzero, obtienes

>> fzero(@error_func, [2,4])
x = 2
x = 4
x = 2.75000000000000
x = 3.03708133971292
x = 2.99755211623500
x = 2.99997750209270
x = 3.00000000025200
x = 3.00000000000000
x = 3
x = 3
ans = 3

No en vano, comienza por la computación$f(2)$ y$f(4)$. Después calcula un punto en el intervalo,$2.75$, y evalúa$f$ ahí. Después de cada iteración, el intervalo se hace más pequeño y la conjetura se acerca a la raíz verdadera. La función fzero se detiene cuando el intervalo es tan pequeño que el cero estimado es correcto a aproximadamente 15 dígitos.

Si quieres saber más sobre cómo funciona fzero, consulta el Capítulo 2.

¿Qué podría salir mal?

El problema más común que tiene la gente con fzero es dejar fuera el @. En ese caso, obtienes algo así:

>> fzero(error_func, [2,4])
Not enough input arguments.

Error in error_func (line 2)
    res = x^2 - 2*x -3;

El error se produce porque MATLAB trata el primer argumento como una llamada a una función, por lo que llama a error_func sin argumentos.

Otro problema común es escribir una función de error que nunca asigne un valor a la variable de salida. En general, las funciones siempre deben asignar un valor a la variable de salida, pero MATLAB no impone esta regla, por lo que es fácil de olvidar.

Por ejemplo, si escribes

function res = error_func(x)
    y = x^2 - 2*x -3
end

y luego llamarlo desde la Ventana de Comando,

>> error_func(4)
y = 5

parece que funcionó, pero no se deje engañar. Esta función asigna un valor a y, y muestra el resultado, pero cuando la función termina, y desaparece junto con el espacio de trabajo de la función. Si intentas usarlo con fzero, obtienes

>> fzero(@error_func, [2,4])
y = -3

Error using fzero (line 231)
FZERO cannot continue because user-supplied function_handle ==>
error_func failed with the error below.

Output argument "res" (and maybe others) not assigned during call
to "error_func".

Si lo lees detenidamente, este es un mensaje de error bastante bueno, siempre que entiendas que “argumento de salida” y “variable de salida” son lo mismo.

Habrías visto el mismo mensaje de error al llamar a error_func desde el intérprete, si hubieras asignado el resultado a una variable:

>> x = error_func(4)
y = 5

Output argument "res" (and maybe others) not assigned during
call to "error_func".

Otra cosa puede salir mal: si proporcionas un intervalo para la suposición inicial y en realidad no contiene una raíz, obtienes

>> fzero(@error_func, [0,1])
Error using fzero (line 272)
The function values at the interval endpoints must differ in sign.

Hay otra cosa que puede salir mal cuando usas fzero, pero esta es menos probable que sea tu culpa. Es posible que fzero no pueda encontrar una raíz.

Generalmente, fzero es robusto, por lo que es posible que nunca tengas un problema, pero debes recordar que no hay garantía de que fzero funcione, especialmente si proporcionas un solo valor como suposición inicial. Incluso si proporciona un intervalo que entre paréntesis a una raíz, las cosas aún pueden salir mal si la función de error es discontinua.

Elegir un Valor Inicial

Cuanto mejor sea tu valor inicial, más probable es que fzero funcione, y menos iteraciones necesitará.

Cuando estás resolviendo problemas en el mundo real, normalmente tendrás cierta intuición sobre la respuesta. Esta intuición suele ser suficiente para proporcionar una buena suposición inicial.

Si no, otra forma de elegir una suposición inicial es trazar la función de error y aproximar visualmente los ceros. Si tienes una función como error_func que toma una variable de entrada escalar y devuelve una variable de salida escalar, puedes trazarla con ezplot:

>> ezplot(@error_func, [-2,5])

El primer argumento es un manejador de función; el segundo es el intervalo en el que desea trazar la función. Al examinar la parcela, se pueden estimar las ubicaciones de las dos raíces.

Funciones de Vectorización

Cuando llamas a ezplot, podrías recibir la siguiente advertencia (o error, si estás usando Octave):

Warning: Function failed to evaluate on array inputs;
vectorizing the function may speed up its evaluation and
avoid the need to loop over array elements.

Esto significa que MATLAB intentó llamar a error_func con un vector, y falló. El problema es que utiliza los operadores * y ^. Con vectores, esos operadores no hacen lo que queremos, que es multiplicación y exponenciación por elementos (ver “” en la página).

Si reescribes error_func así:

function res = error_func(x)
    res = x.^2 - 2.*x -3;
end

el mensaje de advertencia desaparece y ezplot corre más rápido.

Depuración

Cuando empiezas a escribir programas más largos, es posible que dediques más tiempo a depurar. Entonces terminaremos este capítulo con algunos consejos de depuración.

Más colisiones de nombres

Las funciones y variables ocupan el mismo espacio de trabajo, lo que significa que cada vez que aparece un nombre en una expresión, MATLAB comienza buscando una variable con ese nombre; si no hay uno, busca una función.

Como resultado, si tienes una variable con el mismo nombre que una función, la variable sombrea la función; metafóricamente, no puedes ver la función porque la variable está en el camino.

Por ejemplo, si asignas un valor a sin y luego intentas usar la función sin, podrías obtener un error:

>> sin = 3;
>> sin(5)
Index exceeds the number of array elements (1).

'sin' appears to be both a function and a variable.
If this is unintentional, use 'clear sin' to remove
the variable 'sin' from the workspace.

Dado que el valor que asignamos a sin es un número, y un número se considera una$1 \times 1$ matriz, MATLAB intenta acceder al quinto elemento de la matriz y encuentra que no hay uno.

En este caso, MATLAB es capaz de detectar el error, y el mensaje de error es bastante útil. Pero si el valor del pecado fuera un vector, o si el argumento fuera menor, estarías en problemas. Por ejemplo,

>> sin = 3;
>> sin(1)
ans = 3

Solo para revisar, ¡el seno de 1 no es 3!

Puede evitar estos problemas eligiendo cuidadosamente los nombres de las funciones. Use nombres largos y descriptivos para las funciones, no letras simples como f. Para ser aún más claros, usa nombres de funciones que terminen en func. Y antes de definir una función, compruebe si MATLAB ya tiene una función con el mismo nombre.

Depurando tu cabeza

Cuando estés trabajando con una nueva función o una nueva función de idioma por primera vez, debes probarla de forma aislada antes de meterla en tu programa.

Por ejemplo, supongamos que sabes que x es el seno de algún ángulo y quieres encontrar el ángulo. Encuentra la función ASIN de MATLAB, y estás bastante seguro de que calcula la función sinusoidal inversa. Bastante seguro que no es lo suficientemente bueno; quieres estar muy seguro.

Ya que sabemos$\sin 0 = 0$, podríamos intentar

>> asin(0)
ans = 0

lo cual es correcto. Ahora, también sabemos que el seno de$90$ es 1, así que si intentamos asin (1), esperamos que la respuesta sea 90, ¿verdad?

>> asin(1)
ans = 1.5708

¡Uy! Olvidamos que las funciones trigonométricas en MATLAB funcionan en radianes, no grados. La respuesta que obtuvimos es$\pi/2$, que es$90$, en radianes.

Con este tipo de pruebas, realmente no estás comprobando errores en tu programa; estás comprobando tu comprensión. Si cometes un error porque estás confundido acerca de cómo funciona MATLAB, puede tomar mucho tiempo encontrarlo, porque cuando miras el código, se ve bien.

Lo que nos lleva al Séptimo Teorema de Depuración:

Los peores bichos no están en tu código; están en tu cabeza.

Revisión del Capítulo

Este capítulo introdujo fzero, una función que podemos utilizar para resolver ecuaciones no lineales. Para usar fzero, hay que escribir una función de error y pasar un manejador de función como argumento. Usar funciones de esta manera puede ser complicado al principio, pero ponernos cómodos con ella, porque la vamos a usar mucho.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

Si podemos resolver una ecuación realizando operaciones algebraicas y derivando una forma explícita de calcular un valor, el resultado es una solución analítica. De lo contrario, podemos usar un método numérico, que encuentra una solución numérica a la ecuación, que suele ser una aproximación.

Para resolver ecuaciones no lineales, a menudo las reescribimos como funciones y luego encontramos uno o más ceros de la función, es decir, argumentos que hacen el valor de la función$0$.

Un manejador de función es una forma de referirse a una función por su nombre (y pasarla como argumento) sin llamarla.

Finalmente, el sombreado es una especie de colisión de nombres en la que una nueva definición hace que una definición existente se vuelva invisible. En MATLAB, los nombres de variables pueden sombrear funciones integradas (con resultados divertidos).

En el siguiente capítulo, vamos a escribir funciones que toman vectores como entradas y vectores de retorno como salidas.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Escribe una función llamada cheby6 que evalúe el sexto polinomio de Chebyshev. Se debe tomar una variable de entrada,$x$, y devolver

\[32 x^6 - 48 x^4 + 18 x^2 - 1\]
Utilice ezplot para mostrar una gráfica de esta función en el intervalo de$-1$ a$1$. Estimar la ubicación de cualquier ceros en este rango.
Usa fzero para encontrar tantas raíces diferentes como puedas. ¿Fzero siempre encuentra la raíz que está más cerca del valor inicial?

[pato]

Cuando un pato está flotando en el agua, ¿cuánto de su cuerpo está sumergido? ¹

Para estimar una solución a este problema, asumiremos que la parte sumergida de un pato está bien aproximada por una sección de una esfera. Si una esfera con radio$r$ se sumerge en el agua a una profundidad$d$, el volumen de la esfera debajo de la línea de flotación es

\[V = \frac{\pi}{3} (3r d^2 - d^3) \quad \mbox{as long as} \quad d < 2 r\]

También asumiremos que la densidad de un pato es$\rho = 0.3$ g/cm$^3$ (0.3 veces la densidad del agua) y que su masa es$\frac{4}{3} \pi r^3 \rho$ g.

Por último, según la ley de flotabilidad, un objeto flota al nivel donde el peso del agua desplazada es igual al peso total del objeto.

Aquí hay algunas sugerencias sobre cómo proceder:

Escribir una ecuación que$\rho$ relacione$d$,, y$r$.
Reorganice la ecuación para que el lado derecho sea cero. Nuestro objetivo es encontrar valores de$d$ que sean raíces de esta ecuación.
Escribe una función MATLAB que evalúe esta función. Pruebalo, luego conviértalo en una función silenciosa.
Hacer una conjetura sobre el valor de$d_0$ usar como valor inicial.
Usa fzero para encontrar una raíz cerca$d_0$.
Verifique para asegurarse de que el resultado tenga sentido. En particular, ¡comprueba eso$d < 2 r$, porque de lo contrario la ecuación de volumen no funciona!
Prueba diferentes valores de$\rho$ y$r$ y mira si obtienes el efecto que esperas. ¿Qué pasa a medida que$\rho$ aumenta? ¿Va al infinito? ¿Va a cero? ¿Qué pasa a medida que$r$ aumenta? ¿Va al infinito? ¿Va a cero?

Funciones de los vectores

Ahora que tenemos funciones y vectores, los pondremos juntos para escribir funciones que toman vectores como variables de entrada y vectores de retorno como variables de salida. También verás dos patrones para la computación con vectores: la cuantificación existencial y la universal.

Funciones y Vectores

En esta sección veremos patrones comunes que involucran funciones y vectores, y aprenderás a escribir una sola función con la que pueda funcionar así como escalares.

Vectores como variables de entrada

Dado que muchas de las funciones incorporadas toman vectores como argumentos, no debería sorprendernos que puedas escribir funciones que tomen vectores como entrada. Aquí hay un ejemplo simple (pero no muy útil):

function res = display_vector(X)
    for i=1:length(X)
        display(X(i))
    end
end

No hay nada especial en esta función. La única diferencia con las funciones escalares que hemos visto es que ésta usa una letra mayúscula para recordarnos que X es un vector.

El uso de display_vector en realidad no devuelve un valor; solo muestra los elementos del vector que obtiene como una variable de entrada:

>> display_vector(1:3)
    1
    2
    3

Aquí hay un ejemplo más interesante que encapsula el código de Listing [lst:vec_reduce] en la página para sumar los elementos de un vector:

function res = mysum(X)
    total = 0;
    for i=1:length(X)
        total = total + X(i);
    end
    res = total;
end

Llamé a esta función mysum para evitar una colisión con la función incorporada sum, que hace más o menos lo mismo.

Así es como lo llamas desde la Ventana de Comandos:

>> total = mysum(1:3)
total = 6

Debido a que esta función tiene una variable de salida, me propuse asignarla a una variable.

Vectores como variables de salida

Tampoco hay nada de malo en asignar un vector a una variable de salida. Aquí hay un ejemplo que encapsula el código de Listing [lst:vec_apply] en la página:

function res = mysquare(X)
    for i=1:length(X)
        Y(i) = X(i)^2;
    end
    res = Y;
end

Esta función cuadra cada elemento de X y lo almacena como un elemento de Y. Después asigna Y a la variable de salida, res. Así es como usamos esta función:

>> V = mysquare(1:3)
V = 1     4     9

La variable de entrada es un vector con los elementos 1,2,3. La variable de salida es un vector con los elementos 1,4,9.

Funciones de Vectorización

Las funciones que trabajan en vectores casi siempre funcionarán también en escalares, porque MATLAB considera que un escalar es un vector con longitud 1.

>> mysum(17)
ans = 17

>> mysquare(9)
ans = 81

Desafortunadamente, lo contrario no siempre es cierto. Si escribes una función con entradas escalares en mente, es posible que no funcione en vectores.

¡Pero podría! Si los operadores y funciones que usas en el cuerpo de tu función funcionan en vectores, entonces tu función probablemente funcionará en vectores. Por ejemplo, aquí está la primera función que escribimos:

function res = myfunc(x)
    s = sin(x);
    c = cos(x);
    res = abs(s) + abs(c);
end

¡Y lo! Resulta trabajar en vectores:

>> Y = myfunc(1:3)
Y = 1.3818    1.3254    1.1311

Algunas de las otras funciones que escribimos no funcionan en vectores, pero se pueden parchear con solo un poco de esfuerzo. Por ejemplo, aquí está la hipotenusa de la Sección [hipotenuse_exercise]:

function res = hypotenuse(a, b)
    res = sqrt(a^2 + b^2);
end

Esto no funciona en vectores porque el operador ^ intenta hacer exponenciación matricial, que solo funciona en matrices cuadradas.

>> hypotenuse(1:3, 1:3)
Error using  ^  (line 51)
Incorrect dimensions for raising a matrix to a power.
Check that the matrix is square and the power is a scalar.
To perform element-wise matrix powers, use '.^'.

Pero si reemplaza ^ con el operador element-wise (. ^), ¡funciona!

>> A = [3,5,8];
>> B = [4,12,15];
>> C = hypotenuse(A, B)

C = 5    13    17

La función coincide con los elementos correspondientes de los dos vectores de entrada, por lo que los elementos de C son las hipotenusas de los pares$(3,4)$$(5,12)$,, y$(8,15)$, respectivamente.

En general, si escribes una función usando solo operadores elementales y funciones que trabajan en vectores, la nueva función también funcionará en vectores.

Sumas y diferencias

Otra operación de vector común es la suma acumulativa, que toma un vector como entrada y calcula un nuevo vector que contiene todas las sumas parciales del original. En notación matemática, si$V$ es el vector original, los elementos de la suma acumulativa,$C$, son

\[C_i = \sum_{j=1}^i V_j\]

En otras palabras, el elemento$i$ th de$C$ es la suma de los primeros$i$ elementos de$V$. MATLAB proporciona una función llamada cumsum que calcula sumas acumuladas:

>> X = 1:5

X = 1     2     3     4     5

>> C = cumsum(X)

C = 1     3     6    10    15

La operación inversa de cumsum es diff, que calcula la diferencia entre elementos sucesivos del vector de entrada.

>> D = diff(C)

D = 2     3     4     5

Observe que el vector de salida es más corto en uno que el vector de entrada. Como resultado, la versión de diff de MATLAB no es exactamente la inversa de cumsum. Si lo fuera, esperaríamos que cumsum (diff (X)) fuera X:

>> cumsum(diff(X))

ans = 1     2     3     4

Pero no lo es.

Escribe una función llamada mydiff que calcule la inversa de cumsum para que cumsum (mydiff (X)) y mydiff (cumsum (X)) ambos devuelvan X.

Productos y Ratios

La versión multiplicativa de cumsum es cumprod, que calcula el producto acumulado. En notación matemática, eso es

\[P_i = \prod_{j=1}^i V_j\]

En MATLAB, eso parece

>> V = 1:5

V = 1     2     3     4     5

>> P = cumprod(V)

P = 1     2     6    24   120

MATLAB no proporciona la versión multiplicativa de diff, que se llamaría ratio, y que calcularía la relación de elementos sucesivos del vector de entrada.

Escribe una función llamada myratio que calcule la inversa de cumprod, de modo que cumprod (myratio (X)) y myratio (cumprod (X)) ambos devuelvan X.

Puedes usar un bucle, o si quieres ser inteligente, puedes aprovechar el hecho de que$e^{\ln a + \ln b} = a b$.

Computación con vectores

En esta sección, veremos dos patrones comunes para trabajar con vectores y los conectaremos a las ideas correspondientes desde la matemática, la cuantificación existencial y universal. Y aprenderás sobre los vectores lógicos, que contienen los valores booleanos 0 y 1.

Cuantificación Existencial

A menudo es útil verificar los elementos de un vector para ver si hay alguno que satisfaga una condición. Por ejemplo, tal vez quieras saber si hay algún elemento positivo. En términos matemáticos, verificar si algo existe se llama cuantificación existencial, y se denota con el símbolo$\exists$, que se pronuncia “existe”. Por ejemplo,\[\exists x \mbox{~in~} S: x>0\] significa, “existe algún elemento$x$ en el conjunto$S$ tal que”$x>0$. En MATLAB, es natural expresar esta idea con una función lógica, como existe, que devuelve 1 si hay tal elemento y 0 si no lo hay.

function res = exists(X)
    for i=1:length(X)
        if X(i) > 0
            res = 1;
            return
        end
    end
    res = 0;
end

No hemos visto antes la declaración return; es similar a break excepto que se rompe de toda la función, no solo del loop. Ese comportamiento es lo que queremos aquí porque en cuanto encontremos un elemento positivo, conocemos la respuesta (¡existe!) y podemos terminar la función inmediatamente sin mirar el resto de los elementos.

Si llegamos al final del bucle, eso significa que no encontramos lo que buscábamos, por lo que el resultado es 0.

Cuantificación Universal

Otra operación común sobre vectores es verificar si todos los elementos satisfacen una condición, que se denomina cuantificación universal, denotada con el símbolo$\forall$ y pronunciada “para todos”. Entonces la expresión\[\forall x \mbox{~in~} S: x>0\] significa “para todos$x$ los elementos del conjunto”$S$,$x>0$.

Una forma de evaluar esta expresión en MATLAB es reducir el problema a la cuantificación existencial, es decir, reescribir

\[\forall x \mbox{~in~} S: x>0\]

a lo siguiente:\[{\sim} \exists x \mbox{~in~} S: x \le 0\] donde${\sim} \exists$ significa “no existe”. Es decir, verificar que todos los elementos sean positivos es lo mismo que verificar que no hay elementos que sean no positivos.

Escribe una función llamada forall que tome un vector y devuelva 1 si todos los elementos son positivos y 0 si hay elementos no positivos.

Vectores Lógicos

Cuando se aplica un operador lógico a un vector, el resultado es un vector lógico: un vector cuyos elementos son los valores lógicos 1 y 0. Veamos un ejemplo:

>> V = -3:3

V = -3    -2    -1     0     1     2     3

>> L = V>0

L =  0     0     0     0     1     1     1

En este ejemplo, L es un vector lógico cuyos elementos corresponden a los elementos de V. Para cada elemento positivo de V, el elemento correspondiente de L es 1.

Los vectores lógicos se pueden usar como banderas para almacenar el estado de una condición. Y a menudo se usan con la función find, que toma un vector lógico y devuelve un vector que contiene los índices de los elementos que son “verdaderos”.

Aplicando find a L a partir del ejemplo anterior rendimientos

>> find(L)

ans = 5     6     7

lo que indica que los elementos 5, 6 y 7 tienen el valor 1.

Si no hay elementos “verdaderos”, el resultado es un vector vacío.

>> find(V>10)

ans = Empty matrix: 1x0

Este ejemplo calcula el vector lógico y lo pasa como argumento para encontrar sin asignarlo a una variable intermedia. Puedes leer esta versión de manera abstracta como “encuentra los índices de elementos de V que son mayores a 10”.

También podemos usar find para escribir existe de manera más concisa:

function res = exists(X)
    L = find(X>0)
    res = length(L) > 0
end

Escribe una versión de forall usando find.

Depuración en Cuatro Actos

Cuando estás depurando un programa, y especialmente si estás trabajando en un error duro, hay cuatro cosas para probar:

Lectura: Examina tu código, léelo para ti mismo y comprueba que significa lo que querías decir.
Correr: Experimenta haciendo cambios y ejecutando diferentes versiones. A menudo, si muestra lo correcto en el lugar correcto del programa, el problema se vuelve obvio, pero es posible que tenga que invertir tiempo construyendo andamios.
Rumiando: ¡Tómate un tiempo para pensar! ¿Qué tipo de error es: sintaxis, tiempo de ejecución o lógico? ¿Qué información se puede obtener de los mensajes de error o de la salida del programa? ¿Qué tipo de error podría causar el problema que estás viendo? ¿Qué cambiaste último, antes de que apareciera el problema?
Retirándose: En algún momento, lo mejor que puedes hacer es retroceder, deshacer los cambios recientes, hasta que vuelvas a un programa que funcione y que entiendas. Entonces puedes comenzar a reconstruir.

Los programadores principiantes a veces se atascan en una de estas actividades y se olvidan de las demás. Cada actividad viene con su propio modo de falla. Por ejemplo, leer tu código podría ayudar si el problema es un error tipográfico, pero no si el problema es un malentendido conceptual. Si no entiendes lo que hace tu programa, puedes leerlo 100 veces y nunca ver el error, porque el error está en tu cabeza.

Ejecutar experimentos puede ayudar, especialmente si realizas pruebas pequeñas y simples. Pero si ejecutas experimentos sin pensar o leer tu código, podrías caer en un patrón al que llamo “programación de caminata aleatoria”, que es el proceso de hacer cambios aleatorios hasta que el programa haga lo correcto. No hace falta decir que la programación de caminata aleatoria puede llevar mucho tiempo.

La salida es tomar más tiempo para pensar. La depuración es como una ciencia experimental. Deberías tener al menos una hipótesis sobre cuál es el problema. Si hay dos o más posibilidades, trata de pensar en una prueba que elimine a una de ellas.

Tomar un descanso a veces ayuda con el pensamiento. También lo hace hablar. Si le explicas el problema a otra persona (o incluso a ti mismo), a veces encontrarás la respuesta antes de terminar de hacer la pregunta.

Pero incluso las mejores técnicas de depuración fallarán si hay demasiados errores o si el código que intentas corregir es demasiado grande y complicado. A veces la mejor opción es retirarse, simplificando el programa hasta llegar a algo que funcione, y luego reconstruir.

Los programadores principiantes suelen ser reacios a retirarse, porque no pueden soportar eliminar una línea de código (aunque esté mal). Si te hace sentir mejor, copia tu programa en otro archivo antes de comenzar a desnudarlo. Entonces puedes volver a pegar las piezas, poco a poco.

Para resumir, aquí está el Octavo Teorema de Depuración:

Encontrar un error duro requiere leer, correr, rumiar y, a veces, retroceder. Si te quedas atascado en alguna de estas actividades, prueba las otras.

Revisión del Capítulo

Este capítulo presenta patrones para trabajar con vectores, incluyendo la cuantificación existencial y universal. Aprendimos a escribir funciones que toman vectores como variables de entrada y vectores de retorno como variables de salida. Y aprendimos sobre los vectores lógicos, que contienen los valores 1 y 0 para representar verdadero y falso.

Algunas de las funciones de este capítulo no son idiomáticas de MATLAB; muchas de ellas se pueden hacer más simplemente usando operadores y funciones integradas de MATLAB, en lugar de escribirlas usted mismo. Pero estos ejemplos demuestran conceptos que necesitarás saber cuando trabajes en problemas más complicados.

En el siguiente capítulo, aplicaremos las herramientas que hemos aprendido hasta ahora al objetivo central de este libro, modelar sistemas físicos.

Ecuaciones diferenciales ordinarias

En el capítulo anterior, encontramos el punto de equilibrio donde un pato flotaría sobre el agua. Este tipo de problema se llama estático porque no se mueve. Este capítulo introduce problemas dinámicos, que involucran cosas que cambian con el tiempo.

Además, aprenderás sobre una herramienta matemática para describir sistemas físicos, ecuaciones diferenciales y dos herramientas computacionales para resolverlos, el método de Euler y ode45.

Pero primero tengo una sugerencia rápida sobre la organización del código en archivos.

Funciones y archivos

Hasta ahora sólo hemos puesto una función en cada archivo. También es posible poner más de una función en un archivo, pero solo la primera, la función de nivel superior, se puede llamar desde la Ventana de Comandos. Las otras funciones de ayuda se pueden llamar desde cualquier lugar dentro del archivo, pero no desde la Ventana de Comandos o cualquier otro archivo.

Mantener múltiples funciones en un archivo es conveniente, pero dificulta la depuración porque no se puede llamar a las funciones de ayuda desde la ventana de comandos.

Para ayudar con este problema, a menudo uso la función de nivel superior para desarrollar y probar mis funciones auxiliares. Por ejemplo, para escribir un programa para Section [duck], crearía un archivo llamado duck.m y comenzaría con una función de nivel superior llamada duck que no toma variables de entrada y no devuelve ningún valor de salida.

Entonces escribiría una función llamada error_func para evaluar la función de error para fzero. Para probar error_func, lo llamaría de pato y luego llamaría pato desde la Ventana de Comando.

Así es como podría ser mi primer borrador:

function res = duck()
    error = error_func(10)
end

function res = error_func(d)
    rho = 0.3;      % density in g / cm^3
    r = 10;         % radius in cm
    res = d;
end

Este programa no está completo, pero es suficiente código para probar. Una vez que este programa esté funcionando, terminaría de escribir error_func. Y una vez que hubiera terminado y probado error_func, modificaría pato para usar fzero.

Este problema podría requerir solo dos funciones, pero si hubiera más, podría escribirlas y probarlas una a la vez y luego combinarlas en un programa de trabajo.

Ahora, volvamos a las ecuaciones diferenciales.

Ecuaciones diferenciales

Una ecuación diferencial (DE) es una ecuación que describe las derivadas de una función desconocida. “Resolver un DE” significa encontrar una función cuyas derivadas satisfagan la ecuación.

Por ejemplo, supongamos que nos gustaría predecir la población de levadura que crece en una solución nutritiva. Supongamos que sabemos que la población inicial es de 5 mil millones de células de levadura. Cuando la levadura crece en condiciones particularmente amigables con las levaduras, la tasa de crecimiento en cualquier momento es proporcional a la población actual. Si$y(t)$ definimos que es la población a la vez$t$, podemos escribir la siguiente ecuación para la tasa de crecimiento:\[\frac{dy}{dt}(t) = a y(t)\] donde$\frac{dy}{dt}(t)$ es la derivada de$y(t)$ y$a$ es una constante que caracteriza la rapidez con que crece la población. Esta ecuación es diferencial porque relaciona una función con una de sus derivadas.

Se trata de una ecuación diferencial ordinaria (ODE) porque todas las derivadas involucradas se toman con respecto a una misma variable. Si relacionara derivadas con respecto a diferentes variables (derivadas parciales), sería una ecuación diferencial parcial (PDE).

Esta ecuación es de primer orden porque involucra solo primeras derivadas. Si se tratara de segundas derivadas, sería de segundo orden, y así sucesivamente.

Por último, es lineal porque cada término involucra$t$$y$, o$dy/dt$ elevado a la primera potencia; si alguno de los términos involucra productos o poderes de$t$$y$, o$dy/dt$ sería no lineal.

Ahora supongamos que queremos predecir la población de levaduras en el futuro. Podemos hacerlo usando el método de Euler.

Método de Euler

Aquí tienes una prueba para ver si eres tan listo como Leonhard Euler. Digamos que llegas a time ($t$) y mides la población actual ($y$) y la tasa de cambio ($r$). ¿Cuál cree que será la población después de que$\Delta t$ haya transcurrido algún periodo de tiempo?

Si lo dijiste$y + r \Delta t$, ¡enhorabuena! Acabas de inventar el método de Euler.

Esta estimación se basa en el supuesto de que$r$ es constante, pero en general no lo es, por lo que solo esperamos que la estimación sea buena si$r$ cambia lentamente y$\Delta t$ es pequeña.

¿Y si queremos hacer una predicción cuando$\Delta t$ es grande? Una opción es romper$\Delta t$ en pedazos más pequeños, llamados pasos de tiempo. Entonces podemos usar las siguientes ecuaciones para pasar de un paso de tiempo al siguiente:\[\begin{aligned} T_{i+1} &=& T_i + dt \\ Y_{i+1} &=& Y_i + \frac{df}{dt}(t) \, dt\end{aligned}\]

Aquí,$T_i$ es una secuencia de tiempos donde se estima el valor de$y$, y$Y_i$ es la secuencia de estimaciones. Para cada índice$i$,$Y_i$ es una estimación de$y(T_i)$.

Si la tasa no cambia demasiado rápido y el paso de tiempo no es demasiado grande, el método de Euler es lo suficientemente preciso para la mayoría de los propósitos.

Implementando el Método de Euler

Como ejemplo usaremos el método de Euler para resolver la ecuación desde la página,\[\frac{dy}{dt}(t) = a y(t)\] con la condición inicial$y(0) = 5$ mil millones de celdas y el parámetro de crecimiento$a = 0.2$ por hora.

Como primer paso, cree un archivo llamado euler.m con una función de nivel superior y una función auxiliar:

function res = euler()
    T(1) = 0;
    Y(1) = 5;
    r = rate_func(T(1), Y(1))
end

function res = rate_func(t, y)
   a = 0.2;
   dydt = a * y;
   res = dydt;
end

En euler inicializamos las condiciones iniciales y luego llamamos rate_func, así llamado porque calcula la tasa de crecimiento en la población.

Después de probar estas funciones, podemos agregar código a euler para calcular estas ecuaciones de diferencia:\[\begin{aligned} T_{i+1} &=& T_i + \Delta t \\ Y_{i+1} &=& Y_i + r \Delta t\end{aligned}\] donde$r$ está la tasa de crecimiento calculada por rate_func. El listado [lst:euler_method] tiene el código que necesitamos:

function res = euler()
    T(1) = 0;
    Y(1) = 5;
    dt = 0.1;

    for i=1:40
        r = rate_func(T(i), Y(i));
        T(i+1) = T(i) + dt;
        Y(i+1) = Y(i) + r * dt;
    end
    plot(T, Y)
end

Antes del bucle, creamos dos vectores, T e Y, y establecemos el primer elemento de cada uno con las condiciones iniciales; dt, que es el tamaño de los pasos de tiempo, es de 0.1 horas.

Dentro del bucle, calculamos la tasa de crecimiento con base en el tiempo actual, T (i), y la población, Y (i). Podrías notar que la tasa depende solo de la población, pero pasamos el tiempo como variable de entrada de todos modos, por razones que verás pronto.

Después de calcular la tasa de crecimiento, agregamos un elemento tanto T como Y. Entonces, cuando sale el bucle, trazamos Y como una función de T.

Si ejecutas el código, deberías obtener una gráfica de población a lo largo del tiempo, como se muestra en la Figura [fig:euler].

Solución a una ecuación diferencial simple por el método de Euler

Como se puede ver, la población se duplica en un poco menos de 4 horas.

Resolviendo ODEs con ode45

Una limitación del método de Euler es que asume que la derivada es constante entre pasos de tiempo, y eso no es generalmente cierto. Afortunadamente, existen mejores métodos que estiman la derivada entre pasos de tiempo, y son mucho más precisos.

MATLAB proporciona una función llamada ode45 que implementa uno de estos métodos. En esta sección te explicaré cómo usarlo; puedes leer más sobre cómo funciona en “” en la página.

Para poder usar ode45, hay que escribir una función que evalúe$dy/dt$ como una función de$t$ y$y$. Afortunadamente, ya tenemos uno, llamado rate_func:

function res = rate_func(t, y)
   a = 0.2;
   dydt = a * y;
   res = dydt;
end

Podemos llamar a ode45 desde la Ventana de Comando así:

[T, Y] = ode45(@rate_func, [0, 4], 5);
plot(T, Y)

El primer argumento es un manejador de función, como vimos en Chapter. El segundo argumento es el intervalo de tiempo donde queremos evaluar la solución; en este caso el intervalo es de$t=0$ a$t=4$ horas. El tercer argumento es la población inicial, 5 mil millones de células.

La función ode45 es la primera función que hemos visto que devuelve dos variables de salida. Para almacenarlos, tenemos que asignarlos a dos variables, T e Y. La figura [fig:runge] muestra los resultados.

Soluciones a una ecuación diferencial simple usando el método de Euler y ode45

La línea continua es la estimación que calculamos con el método de Euler; la línea discontinua es la solución a partir de ode45.

Durante las primeras 2—3 horas, las dos soluciones son visualmente indistinguibles. Durante la última hora, divergen ligeramente; a las 4 horas, la diferencia es menor al 1 por ciento.

Para muchos propósitos, la diferencia entre el método de Euler y ode45 es la menor de nuestras preocupaciones. En este ejemplo, probablemente no conocemos la población inicial con perfecta precisión o la constante de crecimiento, a. Además, la suposición de que la tasa de crecimiento sólo depende de la población probablemente no sea cierta. Cualquiera de estos errores de modelado podría ser mayor al 1 por ciento.

No obstante, para algunos problemas, el método de Euler puede estar apagado en mucho más del 1 por ciento. En esos casos ode45 es casi siempre más preciso, por dos razones: primero, calcula la función de tasa varias veces por paso de tiempo; segundo, si el paso de tiempo es demasiado grande, ode45 puede detectar el problema y reducir el paso de tiempo. Para más detalles, consulte “” en la página.

Dependencia del Tiempo

Al mirar rate_func en la sección anterior, podrías notar que toma t como variable de entrada pero no la usa. Eso se debe a que la tasa de crecimiento no depende del tiempo, las bacterias no saben qué hora es.

Pero las ratas sí. O, al menos, saben qué temporada es. Supongamos que la tasa de crecimiento de las ratas depende de la población actual y de la disponibilidad de alimento, que varía a lo largo del año. La ecuación diferencial podría ser algo así como\[\frac{dy}{dt}(t) = a y(t) \left(1 - \cos (\omega t) \right)\] dónde$t$ está el tiempo en días y$y(t)$ es la población en el momento$t$. Debido a que la tasa de crecimiento depende del tiempo, esta ecuación diferencial depende del tiempo.

Las variables$a$ y$\omega$ son parámetros, que son valores que cuantifican un aspecto físico del escenario. Los parámetros suelen ser constantes, pero en algunos modelos varían en el tiempo.

En este ejemplo,$a$ caracteriza la tasa reproductiva por día, y$\omega$ es la frecuencia de una función periódica que describe el efecto de la variación del suministro de alimentos en la reproducción.

Usaremos los valores$a = 0.002$ y$\omega = 2 \pi/365$ (un ciclo por año). La tasa de crecimiento es más baja en$t=0$, el 1 de enero, y más alta en$t=365/2$, el 30 de junio.

Ahora podemos escribir una función que evalúe la tasa de crecimiento:

function res = rate_func(t, y)
    a = 0.002;
    omega = 2*pi / 365;
    res = a * y * (1 - cos(omega * t));
end

Para probar esta función, la puse en un archivo llamado rats.m con una función de nivel superior llamada rats:

function res = rats()
    t = 365/2;
    y = 1000;
    res = rate_func(t, y);
end

La función de nivel superior asume, para fines de prueba, que hay 1000 ratas en$t=365/2$ (30 de junio) y calcula la tasa de crecimiento bajo esas condiciones.

Podemos ejecutar la función de nivel superior así:

>> r = rats

r = 4

En estas condiciones, la tasa de crecimiento es de 4 ratas nuevas por día.

Ahora que hemos probado rate_func, podemos usar ode45 para resolver la ecuación diferencial. A continuación, le indicamos cómo llamarlo desde la función de nivel superior en rats.m:

[T, Y] = ode45(@rate_func, [0, 365], 1000)
plot(T, Y)

El primer argumento es un manejador de función, de nuevo. El segundo argumento es el intervalo que nos interesa, una duración de un año, expresado en unidades de días. El tercer argumento es la población inicial,$y(0) = 1000$.

La figura [fig:ratas] muestra los resultados.

Soluciones a una ecuación diferencial simple por el método de Euler y ode45

La población crece lentamente durante el invierno, rápidamente durante el verano, y luego lentamente nuevamente en el otoño.

Para ver la población al final del año, se puede mostrar el último elemento de Y:

Y(end)
2.0751e+03

Eso es un poco más de 2000 ratas, por lo que la población se duplica aproximadamente en un año.

El índice aquí es fin, que es una palabra especial en MATLAB que significa “el índice del último elemento”. Se puede utilizar en una expresión, por lo que Y (fin - 1) es el penúltimo elemento de Y.

¿Qué podría salir mal?

No olvides el @ en el mango de la función. Si lo dejas fuera, como:

[T, Y] = ode45(rate_func, [0, 365], 1000)

MATLAB trata el primer argumento como una llamada a función y llama a rate_func sin proporcionar argumentos. Entonces recibes un mensaje de error:

Not enough input arguments.

Error in rats>rate_func (line 18)
    res = a * y * (1 - cos(omega * t));

Error in rats (line 6)
    [T, Y] = ode45(rate_func, [0, 365], 1000);

Además, la función de tasa que escribas tiene que tomar dos variables de entrada, t e y, en ese orden, y devolver una variable de salida, res.

Si estás trabajando con una función de tarifa como:

\[\frac{dy}{dt}(t) = a y(t)\]

podrías estar tentado a escribir esto:

function res = rate_func(y)        % WRONG
    a = 0.002;
    res = a * y;
end

Pero eso estaría mal. Así que muy mal. ¿Por qué? Porque cuando ode45 llama a rate_func, proporciona dos argumentos. Si solo tomas una variable de entrada, obtendrás un error. Entonces tienes que escribir una función que tome t como variable de entrada, aunque no la uses:

function res = rate_func(t, y)     % RIGHT
    a = 0.002;
    res = a * y;
end

Otro error común es escribir una función que no haga una asignación a la variable de salida. Si escribes algo como:

function res = rate_func(t, y)
    a = 0.002;
    omega = 2*pi / 365;
    r = a * y * (1 - cos(omega * t));    % WRONG
end

y luego lo llaman desde ode45, se obtiene

Output argument "res" (and maybe others) not assigned during call
to "rate_func".

Espero que estas advertencias te ahorren algo de tiempo depurando.

Ejes de Etiquetado

Las tramas de este capítulo tienen etiquetas en los ejes, y una de ellas tiene una leyenda, pero no te mostré cómo hacerlo. Hagámoslo ahora.

Las funciones para etiquetar los ejes son xlabel e ylabel:

xlabel('Time (hours)')
ylabel('Population (billions of cells)')

La función para generar una leyenda es legend:

legend('euler', 'ode45')

Los argumentos son las etiquetas de las líneas, en el orden en que fueron dibujadas. Por lo general, la leyenda está en la esquina superior derecha, pero puedes moverla proporcionando un argumento opcional llamado Ubicación:

legend('euler', 'ode45', 'Location', 'northwest')

Finalmente, guarda las cifras usando saveas:

saveas(gcf, 'runge.eps', 'epsc')

El primer argumento es la figura que queremos guardar; gcf es un comando de MATLAB que significa “get current figure”, que es la figura que acabamos de sacar. El segundo argumento es el nombre del archivo. La extensión especifica el formato que queremos, que es Encapsulated PostScript (.eps). El tercer argumento le dice a MATLAB qué controlador usar. Los detalles no son importantes, pero epsc genera figuras en color.

Revisión del Capítulo

En este capítulo se introdujeron las ecuaciones diferenciales (DE), que son ecuaciones que describen las derivadas de una función desconocida. En una ecuación diferencial ordinaria (ODE), todas las derivadas se toman con respecto a una misma variable, a diferencia de una ecuación diferencial parcial (PDE), que incluye derivadas con respecto a más de una variable.

Un DE de primer orden incluye solo las primeras derivadas, y una DE lineal no incluye productos ni potencias de la función y sus derivados. Una ecuación diferencial depende del tiempo si la función de tasa depende del tiempo.

Cuando resolvemos una ecuación diferencial numéricamente, el paso de tiempo es el intervalo en el tiempo entre elementos sucesivos de la solución. Un parámetro es un valor que aparece en un modelo para cuantificar algún aspecto físico del escenario que se está modelando.

Hasta ahora solo hemos puesto una función en cada archivo M, pero en este capítulo escribimos una función de nivel superior, que es la primera función en un archivo M, y una función auxiliar, que es cualquier función en un archivo M que no sea la función de nivel superior.

En el siguiente capítulo, resolveremos sistemas de ODEs, los cuales se utilizan para describir sistemas físicos con múltiples partes que interactúan. Pero primero, aquí tienes un ejercicio donde puedes aplicar lo que has aprendido hasta ahora.

Ejercicio

Antes de continuar, es posible que desee trabajar en el siguiente ejercicio.

Supongamos que te dan una taza de café de 8 onzas a 90 °C. Por experiencia amarga has aprendido que el café más caliente que puedes tomar cómodamente es de 60 °C.

Si la temperatura del café desciende 0.7 °C durante el primer minuto, ¿cuánto tiempo tendrás que esperar para tomar tu café?

Puedes responder a esta pregunta con Newton's Law of Cooling (ver https://greenteapress.com/matlab/newton):\[\frac{dy}{dt}(t) = -k (y(t) - e)\] dónde$y(t)$ está la temperatura del café en el momento$t$,$e$ es la temperatura del ambiente, y $k$es un parámetro que caracteriza la velocidad de transferencia de calor del café al ambiente.

Supongamos que$e$ es de 20 °C y constante; es decir, el café no calienta la habitación.

Supongamos también que$k$ es constante. En ese caso, podemos estimarlo con base en la información que tenemos. Si la temperatura baja 0.7 °C durante el primer minuto, cuando el café es de 90 °C, podemos escribir\[-0.7 = -k (90 - 20)\]

Resolver esta ecuación rinde$k = 0.01$.

Aquí hay algunas sugerencias para comenzar:

Crea un archivo llamado coffee.m y escribe una función llamada coffee que no tome variables de entrada. Poner una declaración simple como x = 5 en el cuerpo de la función e invocar café desde la Ventana de Comandos.
Agrega una función auxiliar llamada rate_func que toma t e y y calcula$dy/dt$. En este caso, rate_func en realidad no depende de$t$; sin embargo, tu función tiene que tomar$t$ como primera variable de entrada para poder trabajar con ode45.
Pruebe su función agregando una línea como rate_func (0, 90) al café, luego llame al café desde la Ventana de Comandos. Confirmar que la tasa inicial es de -0.7 °C.
Una vez que consiga que rate_func funcione, modifique el café para usar ode45 para calcular la temperatura del café durante 60 minutos. Confirmar que el café se enfría rápidamente al principio, luego se enfría más lentamente, y alcanza la temperatura ambiente después de aproximadamente una hora.
Trazar los resultados y estimar el tiempo en que la temperatura alcanza los 60 °C.

Sistemas de ODEs

En el capítulo anterior se utilizó el método de Euler y ode45 para resolver una única ecuación diferencial de primer orden. En este capítulo, pasaremos a los sistemas de ODE e implementaremos un modelo de un sistema depredador-presa. Pero primero, tenemos que aprender más sobre las matrices.

Matrices

Una matriz es una versión bidimensional de un vector. Al igual que un vector, contiene elementos que son identificados por índices. La diferencia es que los elementos están dispuestos en filas y columnas, por lo que se necesitan dos índices para identificar un elemento.

Creación de una Matriz

Una forma común de crear una matriz es la función ceros, que devuelve una matriz con el tamaño dado lleno de ceros. En este ejemplo se crea una matriz con dos filas y tres columnas.

>> M = zeros(2, 3)

M =  0     0     0
     0     0     0

Si no conoces el tamaño de una matriz, puedes mostrarla usando whos:

>> whos M
  Name      Size            Bytes  Class     Attributes
  M         2x3                48  double

o la función size, que devuelve un vector:

>> V = size(M)

V = 2    3

El primer elemento es el número de filas; el segundo es el número de columnas.

Para leer un elemento de una matriz, especifique la fila y la columna:

>> M(1,2)

ans = 0

>> M(2,3)

ans = 0

Cuando se trabaja con matrices, se necesita un poco de esfuerzo para recordar qué índice viene primero, fila o columna. Me parece útil repetirme “fila, columna”, como un mantra. También podría resultarle útil recordar “abajo, a través” o la abreviatura RC como en “radio control” o RC Cola.

Otra forma de crear una matriz es encerrar los elementos entre paréntesis, con punto y coma entre filas:

>> D = [1,2,3 ; 4,5,6]

D =  1     2     3
     4     5     6

>> size(D)

ans = 2     3

Vectores de fila y columna

Aunque es útil pensar en términos de números, vectores y matrices, desde el punto de vista de MATLAB todo es una matriz. Un número es solo una matriz que pasa a tener una fila y una columna:

>> x = 5;
>> size(x)

ans = 1     1

Y un vector es una matriz con una sola fila:

>> R = 1:5;
>> size(R)

ans = 1     5

Bueno, algunos vectores tienen sólo una fila, de todos modos. En realidad, hay dos tipos de vectores. Los que hemos visto hasta ahora se llaman vectores de fila, porque los elementos están dispuestos en una fila; el otro tipo son vectores de columna, donde los elementos están en una sola columna.

Una forma de crear un vector de columna es crear una matriz con un solo elemento por fila:

>> C = [1;2;3]

C =

     1
     2
     3

>> size(C)

ans = 3     1

La diferencia entre los vectores de fila y columna es importante en álgebra lineal, pero para la mayoría de las operaciones vectoriales básicas, no importa. Por ejemplo, cuando indexas los elementos de un vector, no tienes que saber de qué tipo es:

>> R(2)

ans = 2

>> C(2)

ans = 2

El Operador de Transposición

El operador de transposición, que se parece notablemente a un apóstrofo, calcula la transposición de una matriz, que es una nueva matriz que tiene todos los elementos del original, pero con cada fila transformada en una columna (o se puede pensar en ello al revés).

En este ejemplo D tiene dos filas:

>> D = [1,2,3 ; 4,5,6]

D =  1     2     3
     4     5     6

por lo que su transposición tiene dos columnas:

>> Dt = D'

Dt = 1     4
     2     5
     3     6

¿Qué efecto tiene el operador de transposición en los vectores de fila, vectores de columna y números?

Resolviendo Sistemas de ODEs

Ahora que hemos visto los fundamentos de las matrices, veamos cómo podemos utilizarlas para resolver sistemas de ecuaciones diferenciales.

Lotka-Volterra

El modelo Lotka-Volterra describe las interacciones entre dos especies en un ecosistema, un depredador y su presa. Como ejemplo, consideraremos zorros y conejos.

El modelo se rige por el siguiente sistema de ecuaciones diferenciales:

\[\begin{aligned} \frac{dx}{dt} &=& a x - b x y \\ \frac{dy}{dt} &=& -c y + d x y\end{aligned}\]

donde$x$ y$y$ son las poblaciones de conejos y zorros, y$a$,$b$$c$, y$d$ son parámetros que cuantifican las interacciones entre las dos especies (ver https://greenteapress.com/matlab/lotka).

A primera vista, podrías pensar que podrías resolver estas ecuaciones llamando a ode45 una vez para resolver por$x$ y una vez para resolver por$y$. El problema es que cada ecuación involucra ambas variables, que es lo que hace de este un sistema de ecuaciones y no sólo una lista de ecuaciones no relacionadas. Para resolver un sistema, hay que resolver las ecuaciones simultáneamente.

Afortunadamente, ode45 puede manejar sistemas de ecuaciones. La diferencia es que la condición inicial es un vector que contiene los valores iniciales$x(0)$ y$y(0)$, y la salida es una matriz que contiene una columna para$x$ y una para$y$.

Listing [lst:lotka_volterra] muestra la función rate con los parámetros$a = 0.1$,$b = 0.01$,$c = 0.1$, y$d = 0.002$:

function res = rate_func(t, V)
    % unpack the elements of V
    x = V(1);
    y = V(2);

    % set the parameters
    a = 0.1;
    b = 0.01;
    c = 0.1;
    d = 0.002;

    % compute the derivatives
    dxdt = a*x - b*x*y;
    dydt = -c*y + d*x*y;

    % pack the derivatives into a vector
    res = [dxdt; dydt];
end

La primera variable de entrada, t, es tiempo. A pesar de que la variable de tiempo no se utiliza en esta función de tasa, tiene que estar ahí para que esta función funcione con ode45. La segunda variable de entrada, V, es un vector con dos elementos,$x(t)$ y$y(t)$. El cuerpo de la función incluye cuatro secciones, cada una explicada por un comentario.

La primera sección desempaqueta el vector copiando los elementos en variables. Esto no es necesario, pero dar nombres a estos valores te ayudará a recordar qué es qué. También hace que la tercera sección, donde calculamos las derivadas, se asemejen a las ecuaciones matemáticas que nos dieron, lo que ayuda a prevenir errores.

La segunda sección establece los parámetros que describen las tasas reproductivas de conejos y zorros, y las características de sus interacciones. Si estuviéramos estudiando un sistema real, estos valores vendrían de observaciones de animales reales, pero para este ejemplo elegí valores que arrojan resultados interesantes.

La tercera sección calcula las derivadas de$x$ y$y$, usando las ecuaciones que nos dieron.

La última sección empaqueta las derivadas calculadas de nuevo en un vector. Cuando ode45 llama a esta función, proporciona un vector como entrada y espera obtener un vector como salida.

Los lectores de ojos agudos notarán algo diferente en esta línea:

    res = [drdt; dfdt];

El punto y coma entre los elementos del vector no es un error. Es necesario en este caso porque ode45 requiere que el resultado de esta función sea un vector de columna.

Como siempre, es una buena idea probar tu función tarifaria antes de llamar a ode45. Cree un archivo llamado lotka.m con la siguiente función de nivel superior:

function res = lotka()
    t = 0;
    V_init = [80, 20];
    rate_func(t, V_init)
end

V_init es un vector que representa la condición inicial, 80 conejos y 20 zorros. El resultado de rate_func es

-8.0000
 1.2000

lo que significa que con estas condiciones iniciales, esperamos que la población de conejos disminuya inicialmente a una tasa de 8 por semana y que la población de zorros aumente en 1.2 por semana.

Ahora podemos ejecutar ode45 así:

tspan = [0, 200]
[T, M] = ode45(@rate_func, tspan, V_init)

El primer argumento es el manejador de función para la función de tasa. El segundo argumento es el lapso de tiempo, de 0 a 200 semanas. El tercer argumento es la condición inicial.

Matrices de Salida

La función ode45 devuelve dos valores: T, un vector, y M, una matriz.

>> size(T)
ans = 101     1

>> size(M)
ans = 101     2

T tiene 101 filas y 1 columna, por lo que es un vector de columna con una fila para cada paso de tiempo.

M tiene 101 filas, una por cada paso de tiempo, y 2 columnas, una para cada variable,$x$ y$y$.

Esta estructura, una columna por variable, es una forma común de usar matrices. Y la trama entiende esta estructura, así que si haces lo siguiente:

>> plot(T, M)

MATLAB entiende que debe trazar cada columna de M frente a T.

Puedes copiar las columnas de M en otras variables como esta:

>> R = M(:, 1);
>> F = M(:, 2);

En este contexto, los dos puntos representan el rango del 1 al final, por lo que M (:, 1) significa “todas las filas, columna 1" y M (:, 2) significa “todas las filas, columna 2".

>> size(R)
ans = 101     1

>> size(F)
ans = 101     1

Entonces R y F son vectores de columna.

Ahora podemos trazar estos vectores por separado, lo que facilita darles diferentes cadenas de estilo:

>> plot(T, R, '-')
>> plot(T, F, '--')

La figura [fig:lotka] muestra los resultados. El eje x es tiempo en semanas; el eje y es población. La curva superior muestra la población de conejos; la curva inferior muestra zorros.

Solución para el modelo Lotka-Volterra

Inicialmente, hay demasiados zorros, por lo que la población de conejos disminuye. Entonces no hay suficientes conejos, y la población de zorros disminuye. Eso permite que la población de conejos se recupere, y el patrón se repite.

Este ciclo de “boom y busto” es típico del modelo Lotka-Volterra.

Parcela de fases

En lugar de trazar las dos poblaciones a lo largo del tiempo, a veces es útil trazarlas una contra la otra:

>> plot(R, F)

La figura [fig:phase] muestra el resultado, que se denomina diagrama de fase. Cada punto de esta parcela representa un cierto número de conejos (en el eje x) y un cierto número de zorros (en el eje y). Dado que estas son las únicas dos variables en el sistema, cada punto en este plano describe el estado completo del sistema, es decir, los valores de las variables para las que estamos resolviendo.

Parcela de fases del modelo Lotka-Volterra

Con el tiempo, el estado se mueve alrededor del avión. La figura [fig:phase] muestra la trayectoria trazada por el estado a lo largo del tiempo; esta trayectoria se denomina trayectoria.

Dado que el comportamiento de este sistema es periódico, la trayectoria es un bucle.

Si hay tres variables en el sistema, necesitamos tres dimensiones para mostrar el estado del sistema, por lo que la trayectoria es una curva 3D. Puedes usar plot3 para trazar trayectorias en tres dimensiones, pero para cuatro o más variables, estás solo.

¿Qué podría salir mal?

El vector de salida de la función de tasa tiene que ser un vector de columna, de lo contrario se obtiene un error:

Error using odearguments (line 93)
RATE_FUNC must return a column vector.

que es un mensaje de error bastante bueno. No está claro por qué tiene que ser un vector de columna, pero ese no es nuestro problema.

Otro posible error es invertir el orden de los elementos en las condiciones iniciales o los vectores dentro de lotka. MATLAB no sabe lo que se supone que significan los elementos, por lo que no puede detectar errores como este; solo producirá resultados incorrectos.

El orden de los elementos (conejos y zorros) depende de ti, pero tienes que ser consistente. Es decir, el orden de las condiciones iniciales que proporcionas cuando llamas a ode45 tiene que ser el mismo que el orden dentro de rate_func donde desempaquetas el vector de entrada y el mismo que el orden de las derivadas en el vector de salida.

Revisión del Capítulo

En este capítulo, se utilizó ode45 para resolver un sistema de ecuaciones diferenciales de primer orden. Como ejercicio, tendrás la oportunidad de resolver las famosas ecuaciones de Lorenz, uno de los primeros ejemplos de un sistema caótico.

Aquí están los términos de este capítulo que quizás quieras recordar.

Un vector de fila es una matriz que tiene solo una fila, y un vector de columna es una matriz que tiene solo una columna. La operación de transposición transforma las filas de una matriz en columnas (o al revés, si lo prefieres).

Un sistema de ecuaciones es una colección de ecuaciones escritas en términos del mismo conjunto de variables.

En una función de tasa, a menudo tenemos que desempaquetar la variable de entrada, copiando los elementos de un vector en un conjunto de variables. Entonces tenemos que empaquetar los resultados en un vector como una variable de salida.

El estado de un sistema es un conjunto de variables que cuantifican la condición del sistema a medida que cambia con el tiempo.

Cuando resolvemos un sistema de ecuaciones diferenciales, podemos visualizar los resultados con una gráfica de fases, que muestra el estado de un sistema como punto en el espacio de estados posibles. Una trayectoria es una trayectoria en una gráfica de fases que muestra cómo cambia el estado de un sistema a lo largo del tiempo.

En el siguiente capítulo, pasaremos a los sistemas de segundo orden, que utilizamos para describir sistemas con objetos que se mueven en el espacio, gobernados por las leyes del movimiento de Newton.

Ejercicios

Antes de continuar, es posible que desee trabajar en el siguiente ejercicio.

A partir de los ejemplos que hemos visto hasta ahora, se pensaría que todas las ODE describen a la población en función del tiempo, pero eso no es cierto.

Por ejemplo, el sistema Lorenz es un sistema de ecuaciones diferenciales basado en un modelo de dinámica de fluidos en la atmósfera (ver https://greenteapress.com/matlab/lorenz). Resulta interesante en parte porque sus soluciones son caóticas; es decir, pequeños cambios en las condiciones iniciales producen grandes diferencias en las soluciones.

El sistema es descrito por estas ecuaciones diferenciales: Los valores\[\begin{aligned} x_t &=& \sigma (y - x) \\ y_t &=& x (r - z) - y \\ z_t &=& xy - b z\end{aligned}\] comunes para los parámetros son$\sigma = 10$,$b = 8/3$, y$r=28$.

Utilice ode45 para estimar una solución a este sistema de ecuaciones.

Cree un archivo llamado lorenz.m con una función de nivel superior llamada lorenz y una función auxiliar llamada rate_func.
La función de tasa debe tomar t y V como variables de entrada, donde se entiende que los componentes de V son los valores actuales de x, y y z. Debe computar las derivadas correspondientes y devolverlas en un solo vector de columna.
Pruebe la función llamándola desde la función de nivel superior con valores como$t=0$$x=1$,$y=2$, y$z=3$. Una vez que consigas que tu función funcione, deberías convertirla en una función silenciosa antes de llamar a ode45.
Utilice ode45 para estimar la solución para el lapso de tiempo$[0, 30]$ con la condición inicial$x=1$,$y=2$, y$z=3$.
Trazar los resultados como una serie temporal, es decir, cada una de las variables en función del tiempo.
Utilice plot3 para trazar la trayectoria de$x$,$y$, y$z$.

Sistemas de segundo orden

Hasta ahora hemos visto ecuaciones diferenciales de primer orden y sistemas de ODEs de primer orden. En este capítulo, presentaremos sistemas de segundo orden, que son particularmente útiles para modelar el movimiento newtoniano.

Movimiento Newtoniano

La segunda ley del movimiento de Newton a menudo se escribe así:

\[F = m a\]

donde$F$ está la fuerza neta que actúa sobre un objeto,$m$ es la masa del objeto, y$a$ es la aceleración del objeto.

Esta ecuación sugiere que si sabes$m$ y$a$, puedes calcular la fuerza. Y eso es cierto, pero en la mayoría de las simulaciones físicas es al revés: basado en un modelo físico, ya sabes$F$ y$m$, y calculas$a$.

Entonces, si conocemos la aceleración en función del tiempo, ¿cómo encontramos la posición del objeto,$r$? Bueno, sabemos que la aceleración es la segunda derivada de la posición, así podemos escribir la ecuación diferencial

\[\frac{d^2r}{dt^2} = a\]

donde${d^2r}/{dt^2}$ es la derivada por segunda vez de$r$.

Debido a que esta ecuación incluye una segunda derivada, es una ODE de segundo orden. No podemos resolver la ecuación usando ode45 en esta forma, pero introduciendo una nueva variable, v, para velocidad, podemos reescribirla como un sistema de ODEs de primer orden:

\[\begin{aligned} \frac{dr}{dt} &=& v \\ \frac{dv}{dt} &=& a\end{aligned}\]

La primera ecuación dice que la primera derivada de$r$ es$v$; la segunda ecuación dice que la primera derivada de$v$ es$a$.

Caída Libre

Como ejemplo de moción newtoniana, volvamos a la pregunta de la Sección [centavo]:

Si dejas caer un centavo desde lo alto del Empire State Building, ¿cuánto tiempo se tarda en llegar a la acera y qué tan rápido va cuando llega ahí?

Empezaremos sin resistencia al aire; luego agregaremos resistencia al aire al modelo y veremos qué efecto tiene.

Cerca de la superficie de la tierra, la aceleración por gravedad es de -9.8 m, donde el signo menos indica que la gravedad tira hacia abajo. Si el objeto cae recto hacia abajo, podemos describir su posición con un valor escalar$y$, que representa la altitud.

Listing [lst:penny_rate_func] contiene una función rate que podemos usar con ode45 para resolver este problema:

function res = rate_func(t, X)
    % unpack position and velocity
    y = X(1);
    v = X(2);

    % compute the derivatives
    dydt = v;
    dvdt = -9.8;

    % pack the derivatives into a column vector
    res = [dydt; dvdt];
end

La función rate en Listing [lst:penny_rate_func] toma t y X como variables de entrada, donde se entiende que los elementos de X son la posición y velocidad del penique.

Devuelve un vector de columna que contiene dydt y dvdt, que son velocidad y aceleración, respectivamente. Dado que la velocidad es el segundo elemento de X, simplemente podemos asignar este valor a dydt. Y como la derivada de la velocidad es la aceleración, podemos asignar la aceleración debida a la gravedad a dvdt.

Como siempre, debemos probar la función de tasa antes de llamar a ode45. Aquí está la función de nivel superior que podemos usar para probarla:

function penny()
   t = 0;
   X = [381, 0];
   rate_func(t, X)
end

La condición inicial de X es la posición inicial, que es la altura del Empire State Building, alrededor de 381 m, y la velocidad inicial, que es de 0 m.

El resultado de rate_func es

    0
   -9.8000

que es lo que esperamos.

Ahora podemos ejecutar ode45 con esta función de tasa:

tspan = [0, 10]
[T, M] = ode45(@rate_func, tspan, X)

Como siempre, el primer argumento es el manejador de función, el segundo es el lapso de tiempo (10 segundos), y el tercero es la condición inicial.

El resultado es un vector, T, que contiene los valores de tiempo, y una matriz, M, que contiene dos columnas, una para altitud y otra para velocidad.

Podemos extraer la primera columna y trazarla, así:

Y = M(:, 1)
plot(T, Y)

Altitud versus tiempo para un objeto en caída libre

La figura [fig:penny] muestra el resultado. La altitud baja lentamente al principio y toma velocidad. Entre 8 y 9 segundos, la altitud llega a 0, lo que significa que el centavo golpea la acera. Pero ode45 no sabe dónde está el suelo, así que el centavo sigue pasando por 0 a altitud negativa. Resolveremos ese problema en la siguiente sección.

Eventos de ODE

Normalmente cuando llamas a ode45 especificas una hora de inicio y una hora de finalización. Pero a veces no se sabe con anticipación cuándo debería terminar la simulación. Para resolver este problema podemos definir un evento, algo de interés que ocurre durante una simulación, como el centavo que llega al suelo.

Estos son los pasos:

Primero definimos una función de evento que permite a ode45 averiguar cuándo ocurre un evento. Aquí hay una función de evento para el ejemplo de penny:
```
function [value, isterminal, direction] = event_func(t,X)
    value = X(1);
    isterminal = 1;
    direction = -1;
end
```
La función de evento toma las mismas variables de entrada que la función de tasa y devuelve tres variables de salida: valor determina cuándo puede ocurrir un evento, dirección determina si lo hace y isterminal determina qué sucede. Más específicamente, un evento puede ocurrir cuando el valor pasa por 0. Si la dirección es positiva, el evento solo ocurre si el valor está aumentando. Si la dirección es negativa, el evento solo ocurre si el valor es decreciente. Si la dirección es 0, el evento siempre ocurre. Si isterminal es 1, el evento hace que la simulación termine; si es 0, la simulación continúa.

Esta función de evento utiliza la altitud del centavo como valor para que un evento pueda ocurrir cuando la altitud pasa por 0. Debido a que la dirección es negativa, un evento ocurre solo cuando la altitud está disminuyendo, y debido a que isterminal es 1, la simulación termina si ocurre un evento.
A continuación, usamos odeset para crear un objeto llamado options:
```
options = odeset('Events', @event_func);
```
El nombre de la opción es Eventos y el valor es el manejador de la función de evento.
Finalmente, pasamos opciones como cuarto argumento a ode45:
```
[T, M] = ode45(@rate_func, tspan, X, options);
```
Cuando se ejecuta ode45, invoca event_func después de cada paso de tiempo. Si la función de evento indica que ocurrió un evento de terminal, ode45 detiene la simulación.

Veamos los resultados del centavo ejemplo:

>> T(end)
8.8179

>> M(end, :)
0.0000  -86.4153

El último valor de T es de aproximadamente 8.8, que es el número de segundos que tarda el centavo en llegar a la acera.

La última fila de M indica que la altitud final es 0, que es lo que queríamos, y la velocidad final es de aproximadamente -86 m.

Resistencia al aire

Para que esta simulación sea más realista, podemos agregar resistencia al aire. Para objetos grandes que se mueven rápidamente a través del aire, la fuerza debida a la resistencia del aire, llamada arrastre, es proporcional a la velocidad al cuadrado. Para un objeto que cae, el arrastre se dirige hacia arriba, por lo que si la velocidad es negativa, la fuerza de arrastre es positiva.

Así es como podemos calcular la fuerza de arrastre en función de la velocidad en una dimensión:

\[f_\mathrm{d} = -\mathrm{sign}(v) b v^2\]

donde$v$ es la velocidad y$b$ es una constante de arrastre que depende de la densidad del aire, el área de sección transversal del objeto y la forma del objeto.

La función signo o signum devuelve el valor$1$ para valores positivos de$v$ y$-1$ para valores negativos. Así$f_\mathrm{d}$ es siempre en la dirección opuesta a$v$.

Para convertir de fuerza a aceleración tenemos que saber masa, pero eso es fácil de encontrar: la masa de un centavo es de aproximadamente 2.5 g No es tan fácil encontrar la constante de arrastre, pero a partir de informes de que la velocidad terminal de un centavo es de aproximadamente 18 m, he estimado que es de aproximadamente 75e-6 m.

Listing [lst:acceleration_drag] define una función que toma t y X como variables de entrada y devuelve la aceleración total del centavo debido a la gravedad y arrastre:

function res = acceleration(t, X)
    b = 75e-6;                 % drag constant in kg/m
    v = X(2);                  % velocity in m/s
    f_d = -sign(v) * b * v^2;  % drag force in N

    m = 2.5e-3;               % mass in kg
    a_d = f_d / m;            % drag acceleration in m/s^2

    a_g = -9.8;               % acceleration of gravity in m/s^2
    res = a_g + a_d;          % total acceleration
end

Primero, calculamos la fuerza debida al arrastre. Después calculamos la aceleración por arrastre. Por último, calculamos la aceleración total por arrastre y.

Tenga cuidado cuando esté trabajando con fuerzas y aceleraciones; asegúrese de agregar solo fuerzas a las fuerzas o aceleraciones a las aceleraciones. En mi código, utilizo comentarios para recordarme qué unidades tienen los valores. Eso me ayuda a evitar errores como sumar fuerzas a las aceleraciones.

Para usar esta función, hacemos un pequeño cambio en rate_func:

function res = rate_func(t, X)
    y = X(1);
    v = X(2);

    dydt = v;
    dvdt = acceleration(t, X);   % this line has changed

    res = [dydt; dvdt];
end

En la versión anterior, dvdt es siempre -9.8, la aceleración debida a la gravedad. En esta versión, llamamos aceleración para calcular la aceleración total por gravedad y arrastre.

Altitud versus tiempo por un centavo en caída libre con resistencia al aire

Todo lo demás es igual. La figura [fig:penny2] muestra el resultado.

¡La resistencia al aire hace una gran diferencia! La velocidad aumenta hasta que la aceleración por arrastre equivale a la aceleración debida a la gravedad; después de eso, la velocidad es constante y la posición disminuye linealmente (y mucho más lentamente de lo que lo haría en el vacío).

Con resistencia al aire, el tiempo hasta que el centavo golpea la acera es de 22.4, sustancialmente más largo que antes (8.8).

Y la velocidad final es de 18.1 m, sustancialmente más lenta que antes (86 m).

Revisión del Capítulo

En este capítulo, utilizamos las leyes del movimiento de Newton para escribir una ecuación diferencial que describe el movimiento de un centavo que cae.

Reescribimos esa ecuación como un sistema de ecuaciones diferenciales de primer orden para que pudiéramos usar ode45 para resolverla. Después ejecutamos simulaciones de un centavo que cae con y sin resistencia al aire, también conocido como drag.

Definimos un evento como algo de interés que ocurre durante una simulación, como una colisión entre objetos en movimiento, y escribimos una función de evento, que permite a ode45 averiguar cuándo ocurre un evento.

En el siguiente capítulo, extendemos el movimiento newtoniano a dos dimensiones y modelizamos el vuelo de una pelota de béisbol.

Ejercicios

Antes de continuar, es posible que desee trabajar en el siguiente ejercicio.

En este ejercicio modelaremos el descenso de un paracaidista, tomando en cuenta el cambio de arrastre cuando se abre el paracaídas.

Modifica el código de penique de este capítulo para simular el descenso de un paracaidista 75 desde una altitud inicial de 4000 m. La constante de arrastre para un paracaidista sin paracaídas es de aproximadamente 0.2 m. ¿Cuál sería la velocidad del paracaidista al impactar?
Después de abrir su paracaídas, la velocidad del paracaidista disminuye a unos 5 m. Usa tu simulación para encontrar la constante de arrastre que produce una velocidad terminal de 5 m.
Aumentar la masa del paracaidista, y confirmar que la velocidad terminal aumenta. Este fenómeno es la fuente de la intuición de que los objetos pesados caen más rápido; ¡en el aire, lo hacen!
Ahora supongamos que el paracaidista cae libre hasta llegar a una altitud de 1000 m antes de abrir el paracaídas. ¿Cuánto tardarían en llegar al suelo?
¿Cuál es la altitud más baja donde el paracaidista puede abrir el paracaídas y aún aterrizar a menos de 6 m (suponiendo que el paracaídas se abre y se despliega instantáneamente)?

[tierra]

Aquí hay una pregunta del sitio web Pregunte a un astrónomo (ver https://greenteapress.com/matlab/astro):

Si la Tierra de repente dejara de orbitar el Sol, sé que eventualmente sería arrastrada por la gravedad del Sol y la golpearía. ¿Cuánto tardaría la Tierra en golpear al Sol? Me imagino que iría despacio al principio y luego cogería velocidad.

Usa ode45 para responder a esta pregunta. Aquí hay algunas sugerencias sobre cómo proceder:

Busca la Ley de Gravitación Universal y cualquier constante que necesites. Te sugiero que trabajes íntegramente en unidades SI: metros, kilogramos y newtons.
Cuando la distancia entre la Tierra y el Sol se hace pequeña, este sistema se comporta mal, por lo que se debe utilizar una función de evento para detenerse cuando la superficie de la Tierra llegue a la superficie del Sol.
Expresa tu respuesta en días y traza los resultados como millones de kilómetros versus días.

Dos Dimensiones

En el capítulo anterior, resolvimos un problema unidimensional: un centavo que caía del Empire State Building. Ahora resolveremos un problema bidimensional: encontrar la trayectoria de una pelota de béisbol.

Para ello, usaremos vectores espaciales para representar cantidades en dos y tres dimensiones, incluyendo fuerza, aceleración, velocidad y posición.

Vectores espaciales

La palabra vector significa cosas diferentes para diferentes personas. En MATLAB, un vector es una matriz que tiene una fila o una columna. Hasta ahora, hemos utilizado vectores MATLAB para representar lo siguiente:

Secuencias: Una secuencia matemática, como los números de Fibonacci, es un conjunto de valores identificados por índices enteros; en el Capítulo 5, se utilizó un vector MATLAB para almacenar los elementos de una secuencia.
Vectores de estado: Un vector de estado es un conjunto de valores que describe el estado de un sistema físico. Cuando llamas a ode45, le das condiciones iniciales en un vector de estado. Entonces, cuando ode45 llama a tu función de tarifa, te da un vector de estado.
Series de tiempo: Uno de los resultados de ode45 es un vector que representa una secuencia de valores de tiempo.

En este capítulo, veremos otro uso de los vectores MATLAB: representar vectores espaciales. Un vector espacial representa una cantidad física multidimensional como posición, velocidad, aceleración o fuerza.

Por ejemplo, para representar una posición en el espacio bidimensional, podemos usar un vector con dos elementos:

>> P = [3 4]

Para interpretar este vector, tenemos que conocer el sistema de coordenadas en el que se define. Más comúnmente, utilizamos un sistema cartesiano donde el eje x apunta al este y el eje y apunta al norte. En ese caso P representa un punto 3 unidades al este y 4 unidades al norte del origen.

Cuando un vector espacial se representa de esta manera, podemos usarlo para calcular la magnitud y dirección de una cantidad física. Por ejemplo, la magnitud de P es la distancia desde el origen hasta P, que es la hipotenusa del triángulo con lados P (1) y P (2). Podemos calcularlo usando el teorema de Pitágoras:

>> sqrt(P(1)^2 + P(2)^2)
ans = 5

O podemos hacer lo mismo usando la norma de función, que calcula la norma euclidiana de un vector, que es su magnitud:

>> norm(P)
ans = 5

Hay dos formas de obtener la dirección de un vector. Una convención es calcular el ángulo entre el vector y el eje x:

>> atan2(P(2), P(1))
ans = 0.9273

En este ejemplo, el ángulo es de aproximadamente 0.9. Pero para fines computacionales, a menudo representamos dirección con un vector unitario, que es un vector con longitud 1. Para obtener un vector unitario podemos dividir un vector por su longitud:

function res = hat(V)
    res = V / norm(V)
end

Esta función toma un vector, V, y devuelve un vector unitario con la misma dirección que V. Se llama hat porque en notación matemática, los vectores unitarios se escriben con un símbolo “hat”. Por ejemplo, el vector unitario con la misma dirección$\vec{P}$ que se escribiría$\uvec{P}$.

Adición de vectores

Los vectores son útiles para representar cantidades como fuerza y aceleración porque podemos sumarlas sin tener que pensar explícitamente en la dirección.

Como ejemplo, supongamos que tenemos dos vectores que representan fuerzas:

>> A = [2, 4];
>> B = [2, -2];

A representa una fuerza que tira hacia el noreste; B representa una fuerza que tira hacia el sureste, como se muestra en la Figura [fig:vector2]:

La suma de dos fuerzas representadas por vectores

Para calcular la suma de estas fuerzas, todo lo que tenemos que hacer es sumar los vectores:

>> A + B
ans = 4     2

Más adelante en el capítulo, usaremos la adición vectorial para agregar aceleraciones debido a diferentes fuerzas que actúan sobre una pelota de béisbol.

ODEs en Dos Dimensiones

Hasta ahora hemos utilizado ode45 para resolver un sistema de ecuaciones de primer orden y una sola ecuación de segundo orden. Ahora daremos un paso más, resolviendo un sistema de ecuaciones de segundo orden.

Como ejemplo, simularemos el vuelo de una pelota de béisbol. Si no hay viento y no hay giro en la pelota, la pelota viaja en un plano vertical, por lo que podemos pensar en el sistema como bidimensional, con$x$ representar la distancia horizontal recorrida desde el lugar de partida y$y$ representando la altura o altitud.

Listing [lst:baseball_rate_func] muestra una función de tasa que podemos usar para simular este sistema con ode45:

function res = rate_func(t, W)
    P = W(1:2);
    V = W(3:4);

    dPdt = V;
    dVdt = acceleration(t, P, V);

    res = [dPdt; dVdt];
end

function res = acceleration(t, P, V)
    g = 9.8;             % acceleration due to gravity in m/s^2
    a_gravity = [0; -g];
    res = a_gravity;
end

El segundo argumento de rate_func se entiende como un vector, W, con cuatro elementos. Los dos primeros se asignan a P, que representa la posición; los dos últimos se asignan a V, que representa la velocidad. Tanto P como V tienen dos elementos, que representan los$y$ componentes$x$ y.

El objetivo de la función rate es calcular la derivada de W, por lo que la salida tiene que ser un vector con cuatro elementos, donde los dos primeros representan la derivada de P y los dos últimos representan la derivada de V. La derivada de P es la velocidad. No tenemos que calcularlo; nos lo dieron como parte de W. La derivada de V es la aceleración. Para calcularlo, llamamos aceleración, que toma como variables de entrada tiempo, posición y velocidad. En este ejemplo, no usamos ninguna de las variables de entrada, pero pronto lo haremos.

Por ahora ignoraremos la resistencia aérea, por lo que la única fuerza en el beisbol es la gravedad. Representamos la aceleración por gravedad con un vector que tiene magnitud g y dirección a lo largo del eje y negativo.

Supongamos que una pelota es bateada desde una posición inicial a 1 m por encima de la placa inicial, con una velocidad inicial de 40 m en la horizontal y 30 m en dirección vertical.

Así es como podemos llamar a ode45 con estas condiciones iniciales:

    P = [0; 1];       % initial position in m
    V = [40; 30];     % initial velocity in m/s
    W = [P; V];       % initial condition

    tspan = [0 8]
    [T, M] = ode45(@rate_func, tspan, W);

P y V son vectores de columna porque ponemos punto y coma entre los elementos. Entonces W es un vector de columna con cuatro elementos. Y tspan especifica que queremos ejecutar la simulación para 8.

Las variables de salida de ode45 son un vector, T, que contiene valores de tiempo y una matriz, M, con cuatro columnas: las dos primeras son posición; las dos últimas son velocidad.

Así es como podemos trazar la posición en función del tiempo:

    X = M(:, 1);
    Y = M(:, 2);

    plot(T, X)
    plot(T, Y)

X e Y obtienen la primera y segunda columnas de M, que son las coordenadas$x$ y de$y$ posición.

La figura [fig:béisbol1] muestra cómo se ven estas coordenadas en función del tiempo. La coordenada x aumenta linealmente porque la$x$ velocidad es constante. La coordenada y sube y baja, como esperamos.

Vuelo simulado de un béisbol que descuida la fuerza de arrastre

La simulación termina justo antes de que la pelota aterrice, habiendo recorrido casi 250 m. Eso es sustancialmente más lejos de lo que viajaría un béisbol real, porque hemos ignorado la resistencia aérea, o “fuerza de arrastre”.

Fuerza de Arrastre

Un modelo simple para la fuerza de arrastre en una pelota de béisbol es

\[\vec{F}_\mathrm{d} = - \frac{1}{2} \, \rho v C_\mathrm{d} A \uvec{V}\]

donde$\vec{F}_\mathrm{d}$ es un vector que representa la fuerza sobre el béisbol debido al arrastre,$\rho$ es la densidad del aire,$C_\mathrm{d}$ es el coeficiente de arrastre, y$A$ es el área de sección transversal.

$\vec{V}$es el vector de velocidad del béisbol,$v$ es la magnitud de$\vec{V}$, y$\uvec{V}$ es un vector unitario en la misma dirección que$\vec{V}$. El signo menos al principio significa que el resultado está en la dirección opuesta a$\vec{V}$.

La función en Listing [lst:baseball_drag] calcula la fuerza de arrastre en una pelota de béisbol:

 function res = drag_force(V)
    C_d = 0.3;      % dimensionless
    rho = 1.3;      % kg / m^3
    A = 0.0042;     % m^2
    v = norm(V);    % m/s

    res = -1/2 * C_d * rho * A * v * V;
end

El coeficiente de arrastre para una pelota de béisbol es de aproximadamente 0.3. La densidad del aire al nivel del mar es de aproximadamente 1.3 m. El área transversal de una pelota de béisbol es de 0.0042 m.

Ahora tenemos que actualizar la aceleración para tener en cuenta drag:

function res = acceleration(t, P, V)
    g = 9.8;                       % acceleration due to gravity in m/s^2
    a_gravity = [0; -g];

    m = 0.145;                     % mass in kilograms
    a_drag = drag_force(V) / m;
    res = a_gravity + a_drag;
end

Al igual que en Listing [lst:baseball_rate_func], la aceleración representa la aceleración debida a la gravedad con un vector que tiene magnitud g y dirección a lo largo del eje y negativo. Pero ahora también calcula la fuerza de arrastre y se divide por la masa del beisbol para conseguir aceleración por arrastre. Por último, agrega a_gravity y a_drag para obtener la aceleración total del beisbol.

La figura [fig:vector3] muestra gráficamente las siguientes cantidades: (1) aceleración por arrastre$\vec{D}$, que está en dirección opuesta a (2) velocidad,$\vec{V}$; (3) aceleración por gravedad$\vec{G}$, que es recta hacia abajo; y (4) total aceleración,$\vec{A}$, que es la suma de$\vec{D}$ y$\vec{G}$.

Diagrama de velocidad,$\vec{V}$ aceleración debida a la fuerza de arrastre,$\vec{D}$ aceleración por gravedad,$\vec{G}$ y aceleración total,$\vec{A}$

La figura [fig:béisbol2] muestra los resultados de ode45. La pelota aterriza después de aproximadamente 5, habiendo recorrido menos de 150 m, sustancialmente menos de lo que obtuvimos sin resistencia al aire, alrededor de 250 m.

Vuelo simulado de una pelota de béisbol incluyendo fuerza de arrastre

Este resultado sugiere que ignorar la resistencia del aire no es una buena opción para modelar una pelota de béisbol.

¿Qué podría salir mal?

¿Qué podría salir mal? Bueno, vertcat para uno. Para explicar lo que eso significa, comenzaré con la concatenación, que es la operación de unir dos matrices en una matriz más grande. La concatenación vertical une las matrices apilándolas una encima de la otra; la concatenación horizontal las coloca una al lado de la otra.

Aquí hay un ejemplo de concatenación horizontal con vectores de fila:

>> x = 1:3
x = 1     2     3

>> y = 4:5
y = 4     5

>> z = [x, y]
z = 1     2     3     4     5

Dentro de los corchetes, el operador coma realiza la concatenación horizontal. El operador de concatenación vertical es el punto y coma. Aquí hay un ejemplo con matrices:

>> X = zeros(2, 3)

X =  0     0     0
     0     0     0

>> Y = ones(2, 3)

Y =  1     1     1
     1     1     1

>> Z = [X; Y]

Z =  0     0     0
     0     0     0
     1     1     1
     1     1     1

Estas operaciones solo funcionan si las matrices son del mismo tamaño a lo largo de la dimensión donde están pegadas. Si no, obtienes

>> a = 1:3

a = 1     2     3

>> b = a'

b =  1
     2
     3

>> c = [a, b]
Error using horzcat
Dimensions of matrices being concatenated are not consistent.

>> c = [a; b]
Error using vertcat
Dimensions of matrices being concatenated are not consistent.

En este ejemplo, a es un vector de fila y b es un vector de columna, por lo que no se pueden concatenar en ninguna dirección.

Al leer los mensajes de error, se podría adivinar que horzcat es la función que realiza la concatenación horizontal, e igualmente con la concatenación vertcat y vertical. Estarías en lo correcto.

En Listing [lst:baseball_rate_func] usamos concatenación vertical para empaquetar dPdt y dVdt en la variable de salida:

function res = rate_func(t, W)
    P = W(1:2);
    V = W(3:4);

    dPdt = V;
    dVdt = acceleration(t, P, V);

    res = [dPdt; dVdt];
end

Mientras dPDT y dVDT sean vectores de columna, el punto y coma realiza concatenación vertical, y el resultado es un vector de columna con cuatro elementos. Pero si alguno de ellos es un vector de fila, eso es un problema.

La función ode45 espera que el resultado de rate_func sea un vector de columna, así que si estás trabajando con ode45, probablemente sea una buena idea hacer de todo un vector de columna.

En general, si te encuentras con problemas con horzcat y vertcat, usa size para mostrar las dimensiones de los operandos, y asegúrate de tener claro en qué dirección van tus vectores.

Revisión del Capítulo

En este capítulo, simulamos el vuelo de una pelota de béisbol con y sin resistencia aérea y vimos que la diferencia es sustancial. Podemos concluir que es importante modelar la resistencia del aire si queremos hacer predicciones precisas sobre pelotas de béiscos y proyectiles similares.

Aquí hay algunos términos de este capítulo que quizás quieras recordar.

Un vector espacial es un valor que representa una cantidad física multidimensional como posición, velocidad, aceleración o fuerza. Un vector espacial tiene una dirección y una magnitud. La magnitud también se llama la norma del vector. Un vector unitario es un vector con la norma 1, que a menudo se usa para representar una dirección.

La concatenación es la operación de unir dos vectores o matrices de extremo a extremo para formar un nuevo vector o matriz.

En el siguiente capítulo, continuaremos con el ejemplo del béisbol, usando fzero, que vimos en Chapter, y una nueva herramienta de optimización, llamada fminsearch. También veremos una forma sencilla de animar la solución de una ecuación diferencial.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Cuando los Medias Rojas de Boston ganaron la Serie Mundial en 2007, interpretaron a los Rockies de Colorado en su campo local en Denver, Colorado. Encuentre una estimación de la densidad del aire en la ciudad de Mile High. ¿Qué efecto tiene esto en el arrastre? ¿Qué efecto tiene en la distancia que recorre el beisbol?

El arrastre real sobre una pelota de béisbol es más complicado que lo que capta nuestro sencillo modelo. En particular, el coeficiente de arrastre depende de la velocidad. Puedes obtener algunos de los detalles de The Physics of Baseball (Harper Perennial, 2002) de Robert K. Adair; la figura que necesitas se reproduce en https://greenteapress.com/matlab/drag.

Utilice estos datos para especificar un modelo más realista de arrastre y modifique su programa para implementarlo. ¿Qué tan grande es el efecto en la distancia que recorre el beisbol?

[cañón]

Según Wikipedia, la distancia récord para una bala de cañón humana es de 59.05 m (ver https://greenteapress.com/matlab/cannon).

Modifique el ejemplo de este capítulo para simular el vuelo de una bala de cañón humana. Es posible que tengas que investigar un poco para encontrar el coeficiente de arrastre y el área transversal de un humano volador.

Encuentra la velocidad inicial (tanto la magnitud como la dirección) que necesitarías para romper este récord. Es posible que tengas que experimentar para encontrar el ángulo de lanzamiento óptimo.

¿Cuánta aceleración puede soportar un humano sin lesionarse? A esta aceleración máxima, ¿cuánto tiempo tendría que ser el cañón del cañón para alcanzar la velocidad inicial necesaria para romper el récord?

Optimización

En el capítulo anterior se te pidió encontrar el mejor ángulo de lanzamiento para una bala de cañón humana, es decir, el ángulo que maximiza la distancia recorrida antes de aterrizar. Este tipo de problema, encontrar mínimos y máximos, se llama optimización.

En este capítulo, resolveremos un problema similar, encontrando el mejor ángulo de lanzamiento para una pelota de béisbol. Resolveremos el problema de dos maneras, primero ejecutando simulaciones con un rango de valores y trazando los resultados, luego usando una función MATLAB que automatiza el proceso, fminsearch.

Béisbol Óptimo

En el capítulo anterior escribimos funciones para simular el vuelo de una pelota de béisbol con una velocidad inicial conocida. Ahora usaremos ese código para encontrar el ángulo de lanzamiento que maximiza el alcance, es decir, la distancia que recorre la pelota antes de aterrizar.

Primero, necesitamos una función de evento para detener la simulación cuando la pelota aterriza.

function [value, isterminal, direction] = event_func(t, W)
    value = W(2);
    isterminal = 1;
    direction = -1;
end

Esto es similar a la función de evento que vimos en el Capítulo 3, excepto que usa W (2) como el valor del evento, que es la coordenada y. Esta función de evento detiene la simulación cuando la altitud de la pelota es 0 y cayendo.

Ahora podemos llamar a ode45 así:

    P = [0; 1];       % initial position in m
    V = [40; 30];     % initial velocity in m/s
    W = [P; V];       % initial condition

    tspan = [0 10];
    options = odeset('Events', @event_func);
    [T, M] = ode45(@rate_func, tspan, W, options);

La posición inicial de la pelota está a 1 m por encima del plato inicial. La velocidad inicial de la pelota es de 40 m en la$x$ dirección y 30 m en la$y$ dirección.

La duración máxima de la simulación es de 10, pero esperamos que un evento detenga primero la simulación. Podemos obtener los valores finales de la simulación así:

    T(end)
    M(end, :)

El tiempo final es 5.1. La$x$ posición final es de 131 m; la$y$ posición final es 0, como se esperaba.

Traición

Ahora podemos extraer las$y$ posiciones$x$ - y -:

    X = M(:, 1);
    Y = M(:, 2);

En el Capítulo 3 trazamos X e Y por separado como funciones del tiempo. Como alternativa, podemos tramarlos uno contra el otro, así:

    plot(X, Y)

La figura [fig:béisbol3] muestra el resultado, que es la trayectoria del beisbol desde el lanzamiento, a la izquierda, hasta el aterrizaje, a la derecha.

Vuelo simulado de una pelota de béisbol trazada como trayectoria

Rango Versus Ángulo

Ahora simularemos la trayectoria del beisbol con un rango de ángulos de lanzamiento. Primero, tomaremos el código que tenemos y lo envolveremos en una función que toma el ángulo de lanzamiento como variable de entrada, ejecuta la simulación y devuelve la distancia que recorre la pelota (Listing [lst:baseball_range]).

function res = baseball_range(theta)
    P = [0; 1];
    v = 50;
    [vx, vy] = pol2cart(theta, v);

    V = [vx; vy];     % initial velocity in m/s
    W = [P; V];       % initial condition

    tspan = [0 10];
    options = odeset('Events', @event_func);
    [T, M] = ode45(@rate_func, tspan, W, options);

    res = M(end, 1);
end

El ángulo de lanzamiento, theta, está en radianes. La magnitud de la velocidad, v, es siempre de 50 m. Utilizamos pol2cart para convertir el ángulo y la magnitud (coordenadas polares) a componentes cartesianos, vx y vy.

Después de ejecutar la simulación extraemos la$x$ posición final y la devolvemos como una variable de salida.

Podemos ejecutar esta función para un rango de ángulos como este:

    thetas = linspace(0, pi/2);
    for i = 1:length(thetas)
        ranges(i) = baseball_range(thetas(i));
    end

Y luego trazar rangos en función de thetas:

    plot(thetas, ranges)

La figura [fig:béisbol4] muestra el resultado. Como era de esperar, la pelota no viaja lejos si es golpeada casi horizontal o vertical. El pico al parecer está cerca de 0.7.

Vuelo simulado de una pelota de béisbol trazada como trayectoria

Considerando que nuestro modelo sólo es aproximado, este resultado podría ser lo suficientemente bueno. Pero si queremos encontrar el pico con mayor precisión, podemos usar fminsearch.

fminsearch

La función fminsearch es similar a fzero, que vimos en Chapter. Recordemos que fzero toma un manejador de función y una suposición inicial, y devuelve una raíz de la función. Como ejemplo, para encontrar una raíz de la función

function res = error_func(x)
    res = x^2 - 2;
end

podemos llamar a fzero así:

>> x = fzero(@error_func, 1)
ans = 1.4142

El resultado está cerca de la raíz cuadrada de 2. Llamemos a fminsearch con la misma función:

>> x = fminsearch(@error_func, 1)
x = -8.8818e-16

El resultado es cercano a 0, que es donde se minimiza esta función. Opcionalmente, fminsearch devuelve dos valores:

>> [x, fval] = fminsearch(@error_func, 1)
x = -8.8818e-16

fval = -2

La x es la ubicación del mínimo, y fval es el valor de la función evaluada en x.

Si queremos encontrar el máximo de una función, en lugar del mínimo, aún podemos usar fminsearch escribiendo una función corta que niegue la función que queremos maximizar. En el ejemplo de béisbol, la función que queremos maximizar es baseball_range; podemos envolverla en otra función como esta:

function res = min_func(angle)
    res = -baseball_range(angle);
end

Y entonces llamamos a fminsearch así:

>> [x, fval] = fminsearch(@min_func, pi/4)

x = 0.6921

fval = -131.5851

El ángulo óptimo de lanzamiento para el beisbol es de 0.69; lanzado en ese ángulo, la pelota viaja casi 132 m.

Si tienes curiosidad sobre cómo funciona fminsearch, consulta “” en la página.

Animación

La animación es una herramienta útil para verificar los resultados de un modelo físico. Si algo anda mal, la animación puede hacerlo obvio. Hay dos formas de hacer animación en MATLAB. Una es usar getframe para capturar una serie de imágenes y luego usar película para reproducirlas.

La forma más informal es dibujar una serie de parcelas. Listing [lst:animate] es una función que anima los resultados de una simulación de béisbol:

function animate(T,M)
    X = M(:,1);
    Y = M(:,2);

    minmax = [min([X]), max([X]), min([Y]), max([Y])];

    for i=1:length(T)
        clf; hold on
        axis(minmax)
        plot(X(i), Y(i), 'o')
        drawnow;

        if i < length(T)
            dt = T(i+1) - T(i);
            pause(dt);
        end
    end
end

Las variables de entrada son la salida de ode45: T, que contiene los valores de tiempo, y M, que contiene la posición y velocidad de la pelota de béisbol.

Un vector de cuatro elementos, minmax se usa dentro del bucle para establecer los ejes de la figura. Esto es necesario porque de lo contrario MATLAB escala la figura cada vez a través del bucle, por lo que los ejes siguen cambiando, lo que hace que la animación sea difícil de ver.

Cada vez que pasa por el bucle, animate usa clf para borrar la figura y el eje para restablecer los ejes. Después traza un círculo para representar la posición del.

Tenemos que llamar a drawnow tras plot para que MATLAB realmente muestre cada parcela. De lo contrario, espera hasta que termines de dibujar todas las figuras y luego actualiza la pantalla.

Podemos llamar animate así:

    tspan = [0 10];
    W = [0 1 30 40];
    [T, M] = ode45(@rate_func, tspan, W);
    animate(T, M)

Una limitación de este tipo de animación es que la velocidad de la animación depende de la rapidez con la que su computadora pueda generar las tramas. Dado que los resultados de ode45 generalmente no están igualmente espaciados en el tiempo, su animación podría ralentizarse donde ode45 toma pequeños pasos de tiempo y acelerar donde los pasos de tiempo son más grandes.

Una forma de solucionar este problema es cambiar la forma en que especificamos tspan. Aquí hay un ejemplo:

    tspan = 0:0.1:10;

El resultado es un vector que va de 0 a 10 con un tamaño de paso de 0.1. Pasar tspan a ode45 en esta forma no afecta la precisión de los resultados; ode45 todavía usa pasos de tiempo variables para generar las estimaciones, pero luego las interpola antes de devolver los resultados.

Con pasos iguales de tiempo, la animación debe ser más suave.

Otra opción es usar pausa para reproducir la animación en tiempo real. Después de dibujar cada fotograma y llamar a drawnow, puede calcular el tiempo hasta el siguiente fotograma y usar pausa para esperar:

    dt = T(i+1) - T(i);
    pause(dt);

Una limitación de este método es que ignora el tiempo requerido para dibujar la figura, por lo que tiende a correr lento, sobre todo si la figura es compleja o el paso de tiempo es pequeño.

Revisión del Capítulo

En este capítulo se presentaron dos nuevas herramientas, fminsearch y animate. La función fminsearch de MATLAB busca eficientemente el mínimo de una función y también se puede adaptar para buscar el máximo. La función de animación es una que escribí para leer resultados de ode45 y generar una animación; la versión de este capítulo funciona con los resultados de la simulación de béisbol, pero se puede adaptar para otras simulaciones.

En los ejercicios a continuación, tienes la oportunidad de extender el ejemplo de este capítulo y reunir muchas de las herramientas que hemos utilizado hasta ahora.

En el siguiente capítulo, pasamos a un nuevo ejemplo, la mecánica celeste, que describe el movimiento de los planetas y otros cuerpos en el espacio exterior.

Ejercicios

Antes de continuar, es posible que desee trabajar en los siguientes ejercicios.

Manny Ramírez es un ex miembro de los Medias Rojas de Boston que fue famoso por su actitud relajada. El objetivo de este ejercicio es resolver el siguiente problema inspirado en muchos:

¿Cuál es el mínimo esfuerzo requerido para golpear un jonrón en Fenway Park?

Fenway Park es un estadio de béisbol en Boston, Massachusetts. Una de sus características más famosas es el “Monstruo Verde”, que es una pared en el campo izquierdo que está inusualmente cerca del plato principal, a solo 310 pies de distancia. Para compensar la corta distancia, la pared es inusualmente alta, a 37 pies.

Puedes resolver este problema en dos pasos:

Para una velocidad dada, encuentra el ángulo de lanzamiento que maximiza la altura de la pelota cuando llega a la pared. Observe que esto no es exactamente lo mismo que el ángulo que maximiza la distancia que recorre la pelota.
Encuentra la velocidad mínima que despeja la pared, dado que tiene el ángulo de lanzamiento óptimo. Sugerencia: esto es en realidad un problema de búsqueda de raíces, no un problema de optimización.

[golf]

Una pelota de golf golpeada con retroceso genera sustentación, lo que podría aumentar la distancia que recorre, pero la energía que se destina a generar giro probablemente llegue a costa de una menor velocidad inicial.

Escribe una simulación del vuelo de una pelota de golf y úsala para encontrar el ángulo de lanzamiento y la asignación de giro y velocidad inicial (para un presupuesto de energía fijo) que maximice el rango horizontal de la pelota en el aire.

El levantamiento de una bola giratoria se debe a la fuerza Magnus (ver https://greenteapress.com/matlab/magnus), que es perpendicular al eje de giro y a la trayectoria de vuelo. El coeficiente de sustentación es proporcional a la velocidad de giro; para una bola que gira a 3000 rpm es de aproximadamente 0.1. El coeficiente de arrastre de una pelota de golf es de aproximadamente 0.2 siempre y cuando la pelota se mueva más rápido que 20 m.

Resortes y cosas

Las herramientas computacionales que has aprendido hasta ahora conforman un kit de herramientas versátil para modelar sistemas físicos descritos por ecuaciones diferenciales de primer y segundo orden y sistemas de ecuaciones.

Con ode45 se pueden calcular las variables de estado de estos sistemas a medida que cambian con el tiempo. Al variar los parámetros del modelo, se puede ver qué efecto tienen en los resultados. Entonces puedes usar fminsearch y fzero para encontrar mínimos, máximos y lugares donde las salidas pasan por cero.

Estas herramientas son todo lo que necesitas para resolver muchos problemas, por lo que este capítulo no presenta nuevas herramientas computacionales (¡vaya!). En cambio, veremos algunos sistemas físicos diferentes y algunas fuerzas con las que aún no hemos tratado, incluidas las fuerzas de resorte y la gravitación universal.

Los ejemplos de este capítulo son un poco más indefinidos que los anteriores. Presentaré un problema motivador y algunos antecedentes, y tendrás la oportunidad de implementar los modelos como ejercicios.

Puenting

Supongamos que quieres establecer el récord mundial de la “volcada bungee” más alta, que es un truco en el que un puente bungee moja una galleta en una taza de té en el punto más bajo de un salto. Un ejemplo se muestra en este video: https://greenteapress.com/matlab/dunk.

Dado que el registro es de 70 m, diseñemos un salto para 80 m. Comenzaremos con los siguientes parámetros:

Inicialmente, el saltador se para sobre una plataforma 80 m por encima de una taza de té. Un extremo del cordón elástico está conectado a la plataforma, el otro extremo está unido al puente y el medio cuelga hacia abajo.
La masa del saltador es de 75, y están sujetos a una aceleración gravitacional de 9.8 m.
En caída libre el puente tiene un área de sección transversal de 1 m y una velocidad terminal de 60 m.

Para modelar la fuerza del cordón elástico en el saltador, haré las siguientes suposiciones:

Hasta que el cable esté completamente extendido, no aplica ninguna fuerza al puente. Resulta que esta podría no ser una buena suposición; la revisaremos en la siguiente sección.
Después de que el cordón esté completamente extendido, obedece a la Ley de Hooke; es decir, aplica una fuerza al saltador proporcional a la extensión del cordón más allá de su longitud de reposo.

Podemos escribir la Ley de Hooke como$F_\mathrm{s} = -k x$ donde$F_\mathrm{s}$ está la fuerza del resorte (cuerda elástica) en el saltador en newtons,$x$ es la distancia que el resorte se estira en metros, y$k$ es una constante elástica que representa la fuerza del primavera en newtons por metro. El signo menos indica que la dirección de la fuerza del resorte es opuesta a la dirección en la que se estira el resorte.

La Ley de Hooke no es una ley en el sentido de que siempre es cierta; realmente, es un modelo de cómo se comportan algunas cosas bajo algunas condiciones. Casi todo obedece a la Ley de Hooke cuando$x$ es lo suficientemente pequeño, pero para valores grandes todo se desvía de este comportamiento ideal, de una manera u otra.

En realidad, la constante elástica de una cuerda elástica depende del$x$ rango que nos interese, pero como punto de partida asumiré que$k$ es constante.

Escribir una simulación de este escenario, con base en estos parámetros y suposiciones de modelado. Usa tu simulación para elegir la longitud del cable, L, y su constante elástica, k, para que el saltador caiga hasta la taza de té, ¡pero no más lejos!

Se podría comenzar con la longitud 25 m y la constante de resorte 40 m.

Bungee Revisitado

En la sección anterior, modelamos el movimiento de un puente bungee teniendo en cuenta la gravedad, la resistencia al aire y la fuerza elástica del cordón elástico. Pero ignoramos el peso del cordón.

Es tentador decir que el cordón no tiene ningún efecto porque cae junto con el saltador, pero esa intuición es incorrecta. Al caer el cordón, transfiere energía al saltador.

En https://greenteapress.com/matlab/bungee, encontrarás un artículo de Heck, Uylings y Kedzierska titulado “Entendiendo la física del puenting”; explica este fenómeno y deriva la aceleración del saltador,$a$, en función de la posición,$y$, y velocidad,$v$:\[a = g + \frac{\mu v^2/2}{\mu(L+y) + 2L}\] donde$g$ está la aceleración por gravedad,$L$ es la longitud del cordón, y$\mu$ es la relación de la masa del cordón,$m$, a la masa del puente, $M$.

Si no crees que su modelo sea correcto, este video podría convencerte: https://greenteapress.com/matlab/chain.

Modifique su solución al problema anterior para modelar este efecto. ¿Cómo cambia el comportamiento del sistema a medida que variamos la masa del cordón? Cuando la masa del cordón es igual a la masa del saltador, ¿cuál es el efecto neto en el punto más bajo del salto?

Spider-Man

En este ejemplo desarrollaremos un modelo de Spider-Man balanceándose a partir de un cable elástico de cincha unido a la parte superior del Empire State Building. Inicialmente, Spider-Man se encuentra en la parte superior de un edificio cercano, como se muestra en la Figura [hombre araña].

Diagrama del estado inicial para el ejemplo de Spider-Man

El origen, O, se encuentra en la base del Empire State Building. El vector H representa la posición donde la cincha está unida al edificio, con relación a O. El vector P es la posición de Spider-Man con respecto a O. Y L es el vector desde el punto de unión a Spider-Man.

Siguiendo las flechas de O, a lo largo de H, y a lo largo de L, podemos ver que

H + L = P

Entonces podemos calcular L así:

L = P - H

Como ejercicio, simula este sistema y estima los parámetros que maximizan la distancia de los columpios de Spider-Man.

Implementar un modelo de este escenario para predecir la trayectoria de Spider-Man.
Elige el momento adecuado para que Spider-Man suelte las cinchas para maximizar la distancia que recorre antes de aterrizar.
Elige el mejor ángulo para que Spider-Man salte del edificio, y el mejor momento para soltar la cincha, para maximizar el alcance.

Utilice los siguientes parámetros:

Según el Wiki de Spider-Man (https://greenteapress.com/matlab/spider), Spider-Man pesa 76.
Supongamos que su velocidad terminal es de 60 m.
La longitud de la web es de 100 m.
El ángulo inicial de la banda es de 45 a la izquierda de recto hacia abajo.
La constante de resorte de la banda es de 40 m cuando se estira el cordón y de 0 m cuando se comprime.

Mecánica Celestial

La mecánica celestial describe cómo se mueven los objetos en el espacio exterior. Si hiciste Sección [tierra], simulaste que la Tierra es arrastrada hacia el Sol en una dimensión. Ahora simularemos la Tierra orbitando el Sol en dos dimensiones.

Para que las cosas sean simples, consideraremos solo el efecto del Sol en la Tierra e ignoraremos el efecto de la Tierra sobre el Sol. Entonces colocaremos al Sol en el origen y usaremos un vector espacial,$\vec{P}$, para representar la posición de la Tierra en relación con el Sol.

Dada la masa del Sol,$m_{1}$, y la masa de la Tierra,$m_{2}$, la fuerza gravitacional entre ellos es

\[\vec{F}_\mathrm{g} = -G \frac{m_1 m_2}{r^2} \uvec{P}\]

donde$G$ esta la constante gravitacional universal (ver https://greenteapress.com/matlab/gravity),$r$ es la distancia de la Tierra del Sol, y$\uvec{P}$ es un vector unitario en la dirección de$\vec{P}$.

Escribe una simulación de la Tierra orbitando el Sol. Se puede buscar la velocidad orbital de la Tierra o buscar manualmente la velocidad inicial que hace que la Tierra haga una órbita completa en un año. Opcionalmente, use fminsearch para encontrar la velocidad que acerca la Tierra lo más posible al lugar de inicio después de un año.

Conservación de Energía

Una manera útil de verificar la precisión de un solucionador de ODE es ver si conserva energía. Para el movimiento planetario, resulta que ode45 no.

La energía cinética de un cuerpo en movimiento es

\[KE = m v^2 / 2\]

La energía potencial de un sol con masa$m_1$ y un planeta con masa$m_2$ y una distancia$r$ entre ellos es

\[PE = -G \frac{m_1 m_2}{r}\]

Escribe una función llamada energy_func que tome la salida de tu simulación de la Tierra y calcule la energía total (cinética y potencial) del sistema para cada posición y velocidad estimadas.

Trazar el resultado en función del tiempo y comprobar si aumenta o disminuye a lo largo de la simulación.

Puede reducir la tasa de pérdida de energía disminuyendo la opción de tolerancia de ode45 usando odeset (vea “” en la página):

options = odeset('RelTol', 1e-5);
[T, M] = ode45(@rate_func, tspan, W, options);

El nombre de la opción es RelTol para “tolerancia relativa”. El valor predeterminado es 1e-3, o 0.001. Los valores más pequeños hacen que ode45 sea menos “tolerante”, por lo que usa tamaños de paso más pequeños para reducir los errores.

Ejecute ode45 con un rango de valores para RelTOL y confirme que a medida que la tolerancia se hace más pequeña, la tasa de pérdida de energía disminuye.

Junto con ode45, MATLAB proporciona varios otros solucionadores de ODE (ver https://greenteapress.com/matlab/solver). Ejecute su simulación con uno de los otros solucionadores de ODE que proporciona MATLAB y vea si alguno de ellos conserva energía. Podrías encontrar que ode23 funciona sorprendentemente bien (aunque técnicamente tampoco conserva energía).

Revisión del Capítulo

En este capítulo se presentan ejemplos donde se pueden aplicar las herramientas de este libro para resolver problemas más realistas. Algunos de ellos son más serios que otros, pero espero que te hayas divertido un poco con ellos.

Creo que este kit de herramientas es potente y versátil. Si lo usas para resolver un problema interesante, ¡avísame!

Bajo el capó

En este capítulo “abrimos el capó”, observando más de cerca cómo funcionan algunas de las herramientas que hemos utilizado —ode45, fzero y fminsearch —.

Cómo funciona ode45

Según la documentación de MATLAB, ode45 utiliza “una fórmula explícita de Runge-Kutta, la pareja Príncipe Dormido”. Puedes leer sobre ello en https://greenteapress.com/matlab/runge, pero te voy a dar una idea de ello aquí.

La idea clave detrás de todos los métodos de Runge-Kutta es evaluar la función de tasa varias veces en cada paso de tiempo y usar un promedio ponderado de las pendientes calculadas para estimar el valor en el siguiente paso de tiempo. Diferentes métodos evalúan la función de tasa en diferentes lugares y calculan el promedio con diferentes pesos.

Como ejemplo, resolveremos la siguiente ecuación diferencial:

\[\frac{dy}{dt}(t) = y \sin t\]

Dada una ecuación diferencial, generalmente es sencillo escribir una función de tasa:

function res = rate_func(t, y)
    dydt = y * sin(t);
    res = dydt;
end

Y podemos usarlo así:

    y0 = 1;
    tspan=[0 4];
    options = odeset('Refine', 1);
    [T, Y] = ode45(@rate_func, tspan, y0, options);

Para este ejemplo usaremos odeset para establecer la opción Refinar en 1, que le dice a ode45 que devuelva solo los pasos de tiempo que calcula, en lugar de interpolar entre ellos.

Ahora podemos modificar la función rate para trazar los lugares donde se evalúa:

function res = rate_func(t, y)
    dydt = y * sin(t);
    res = dydt;

    plot(t, y, 'ro')
    dt = 0.01;
    ts = [t t+dt];
    ys = [y y+dydt*dt];
    plot(ts, ys, 'r-')
end

Cuando se ejecuta rate_func, traza un círculo rojo en cada ubicación y una línea roja corta que muestra la pendiente calculada.

La figura [fig:odeplot1] muestra el resultado; ode45 calcula 10 pasos de tiempo (sin contar la condición inicial) y evalúa la función de tasa 61 veces.

Puntos donde ode45 evalúa la función de tasa

La figura [fig:odeplot2] muestra la misma gráfica, ampliada en un solo paso de tiempo. Los cuadrados oscuros en$0.8$ para$1.2$ mostrar los valores que se devolvieron como parte de la solución. Los círculos muestran los lugares donde se evaluó la función de tasa.

Puntos donde ode45 evalúa la función de tasa, ampliada

Podemos ver que ode45 evalúa la función de tasa varias veces por paso de tiempo, en varios lugares entre los puntos finales. También podemos ver que la mayoría de los lugares donde ode45 evalúa la función de tasa no forman parte de la solución que devuelve, y no siempre son buenas estimaciones de la solución. Esto es bueno saber cuando estás escribiendo una función de tasa; no debes asumir que el tiempo y el estado que obtengas como variables de entrada serán parte de la solución.

En cierto sentido, la función de tasa está respondiendo a una pregunta hipotética: “Si el estado en un momento determinado tiene estos valores particulares, ¿cuál sería la pendiente?”

En cada paso de tiempo, ode45 realmente calcula dos estimaciones del siguiente valor. Al compararlos, puede estimar la magnitud del error, que utiliza para ajustar el paso de tiempo. Si el error es demasiado grande, usa un paso de tiempo más pequeño; si el error es lo suficientemente pequeño, usa un paso de tiempo más grande. Debido a que ode45 es adaptativo de esta manera, minimiza el número de veces que llama a la función de tasa para lograr un nivel de precisión dado.

Cómo funciona fzero

Según la documentación de MATLAB, fzero utiliza “una combinación de métodos de interpolación cuadrática inversa, secante e inversa”. (Ver https://greenteapress.com/matlab/fzero)

Para entender lo que eso significa, supongamos que estamos tratando de encontrar una raíz de una función de una variable,$f(x)$, y supongamos que hemos evaluado la función en dos lugares,$x_1$ y$x_2$, y encontramos que los resultados tienen signos opuestos. Específicamente, supongamos$f(x_1) > 0$ y$f(x_2) < 0$, como se muestra en la Figura [fig:secante].

Estado inicial de una búsqueda de raíz

Siempre y cuando$f$ sea continuo, debe haber al menos una raíz en este intervalo. En este caso diríamos eso$x_1$ y$x_2$ paréntesis un cero.

Si esto fuera todo lo que sabías$f$, ¿a dónde irías buscando una raíz? Si dijiste “a medio camino entre$x_1$ y”$x_2$, ¡enhorabuena! ¡Acabas de inventar un método numérico llamado bisección!

Si dijeras: “Yo conectaría los puntos con una línea recta y calcularía el cero de la línea”, ¡enhorabuena! ¡Acabas de inventar el método secante!

Y si decías: “Yo evaluaría$f$ en un tercer punto, encontraría la parábola que pasa por los tres puntos, y calcularía los ceros de la parábola”, felicidades, ¡acabas de inventar la interpolación cuadrática inversa!

Así es la mayor parte de cómo funciona fzero. Los detalles de cómo se combinan estos métodos son interesantes, pero más allá del alcance de este libro. Puedes leer más en https://greenteapress.com/matlab/brent.

Cómo funciona fminsearch

Según la documentación de MATLAB, fminsearch utiliza el algoritmo simplex de Nelder-Mead. Puedes leer sobre ello en https://greenteapress.com/matlab/nelder, pero puede que te resulte abrumador.

Para darle una idea de cómo funciona, voy a presentar un algoritmo más sencillo, la búsqueda de sección dorada. Supongamos que estamos tratando de encontrar el mínimo de una función de una sola variable,$f(x)$.

Como punto de partida, supongamos que hemos evaluado la función en tres lugares,,$x_1$, y$x_2$$x_3$, y encontramos que$x_2$ arroja el valor más bajo. La figura [fig:golden1] muestra este estado inicial.

Estado inicial de una búsqueda de sección dorada

Supondremos que$f(x)$ es continuo y unimodal en este rango, lo que significa que hay exactamente un mínimo entre$x_1$ y$x_3$.

El siguiente paso es elegir un cuarto punto,$x_4$, y evaluar$f(x_4)$. Hay dos posibles resultados, dependiendo de si$f(x_4)$ es mayor que$f(x_2)$ o no. La figura [fig:golden2] muestra los dos estados posibles.

Posibles estados de una búsqueda de sección dorada después de evaluar$f(x_4)$

Si$f(x_4)$ es menor que$f(x_2)$ (se muestra a la izquierda), el mínimo debe estar entre$x_2$ y$x_3$, así descartaríamos$x_1$ y procederíamos con el nuevo triple$(x_2, x_4, x_3)$.

Si$f(x_4)$ es mayor que$f(x_2)$ (se muestra a la derecha), el mínimo local debe estar entre$x_1$ y$x_4$, así descartaríamos$x_3$ y procederíamos con el nuevo triple$(x_1, x_2, x_4)$.

De cualquier manera, el rango se hace más pequeño y nuestra estimación del valor óptimo de$x$ mejora.

Este método funciona para casi cualquier valor de$x_4$, pero algunas opciones son mejores que otras. Puede que tengas la tentación de dividir el intervalo entre$x_2$ y$x_3$, pero esa no resulta ser la mejor opción. Puedes leer sobre una mejor opción en https://greenteapress.com/matlab/golden.

Revisión del Capítulo

La información de este capítulo no es estrictamente necesaria; puedes usar estos métodos sin saber mucho sobre cómo funcionan. Pero hay dos razones que quizás quieras saber.

Una razón es la pura curiosidad. Si usas estos métodos, y sobre todo si llegas a confiar en ellos, puede resultarte insatisfactorio tratarlos como “cajas negras”. A riesgo de mezclar metáforas, espero que hayan disfrutado abriendo el capó.

La otra razón es que estos métodos no son infalibles; a veces las cosas salen mal. Si sabes cómo funcionan, al menos en un sentido general, tal vez te resulte más fácil depurarlos.

Con eso, has llegado al final del libro, ¡así que enhorabuena! Espero que lo hayas disfrutado y hayas aprendido mucho. Creo que las herramientas de este libro son útiles, y las formas de pensar son importantes, no solo en ingeniería y ciencia, sino en prácticamente todos los campos de investigación.

Los modelos son las herramientas que utilizamos para entender el mundo: si construyes buenos modelos, es más probable que hagas las cosas bien. ¡Buena suerte!

Este ejercicio está adaptado de C. F. Gerald y P. O. Wheatley, Análisis Numérico Aplicado, 4ª edición (Boston: Addison-Wesley, 1989). ↩