(811) - Catapultas, fundíbulos y trabuquetes. Ingeniería y PDE medievales

Supongamos que estás al servicio de un señor feudal entre los siglo XII y XIV, y en medio de una campaña llega la hora de asediar una ciudad amurallada. Para atacar la ciudadela lo mejor sería artillería, pero los cañones, si hay, son muy rudimentarios, y muy poco fiables. Aunque tu señor feudal tenga una clara superioridad numérica y logística, el sitio puede durar de semanas a meses, tiempo en el que el enemigo puede venir a levantar el asedio. La forma más efectiva de tomar la ciudad es de alguna manera hacer un boquete en la muralla enemiga y tomar la ciudad a la fuerza, pero, ¿cómo conseguirlo?

Ahí es donde entra en juego las catapultas, pero sobre todo los fundíbulos y los trabuquetes (hay unos vídeos muy interesantes en YouTube, sobre todo si se busca trechubet).

Los fundíbulos desplazan y rotan una piedra respecto de su eje de sujeción hasta soltarla por los aires para arrasar lo que encuentren por su paso. En esto presentan tres fases:

1. La primera desde que se suelta el soporte activándose el fundíbulo.
2. La segunda empieza cuando justo la cuerda unida a la piedra se tensa y comienza a deslizarse esta a lo largo del riel.
3. La última comienza cuando la piedra, habiéndose desplazado por el riel, tiene suficiente momento, rota un ángulo de 120º-140º (dependiendo del "modelo"), instante en el que con suficiente impulso se manda una piedra rotando sobre so centro de masa, de varios cientos de quilogramos con una velocidad de varios quilómetros por hora.

Los fundíbulos presentan en cada una de las tres etapas un sistema no-lineal y no-homogéneo de ecuaciones diferenciales en derivadas parciales. Resolver analíticamente es prácticamente imposible, mientras que numéricamente puede ser un horror plantearlo y resolverlo computacionalmente. Horror, pero viable.
Aun así, es sorprente que en la Edad Media, sin calculadoras, derivadas, y sin algoritmos rápidos para la multiplicación, como la prostaféresis o los logaritmos (pincha en los enlaces para ver otros artículos al respecto), y solo con un ábaco, no solo tuviesen un conocimiento cualitativo bastante bueno, sino también uno cuantitativo.
Como curiosidad divulgativa, dejo un vídeo de cómo era la vide de un constructor y operador de fundíbulos: https://www.youtube.com/watch?v=HG8wt9alyag

Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(809) - ¿Qué es realmente la mediana?

Vamos a ver el concepto de mediana a distintos niveles de profundidad.
Media, moda y mediana. Esto nos enseñaban en el tema de estadística de matemáticas antes de entrar al instituto. Media es sumarlos todos y dividir entre el número de sumandos. Moda es el [elemento] que más se repite. Mediana es el [elemento] de en medio. Sin embargo, ¿qué significa esto? Tanto media y moda como mediana están entre el mínimo y el máximo de los elementos. Asimismo el de en medio puede ser poco intuitivo, por ejemplo $\{1,3,5,5,5\}$ tiene mediana $5$ , no $3$ como uno podría pensar al ser el de en medio en $\{1,3,5\}$ . Esto es un error recurrente: no se pueden eliminar los elementos repetidos.

La mediana de una lista $\{x_i\}_{i=1}^n$ se define como el ecuador de dicha lista, es decir, como el elemento que parte la lista en dos sublistas ordenadas de igual longitud: una con todos sus elementos menores o iguales que la mediana, mientras que la otra tiene todos sus elementos mayores o igual que la mediana. Para una distribución simétrica con una única moda (unimodal), la mediana coincide con la moda. Para calcular la mediana de una lista $\{x_i\}_{i=1}^n$ realmente se calcula la mediana de una lista equivalente: se reordenan los elementos de menor a mayor, que se denota son el subíndice entre paréntesis, $\{x_{(j)}\}_{j=1}^n$ $$ \overset{n}{\underset{i=1}{\operatorname{med}}}(x_i) = \overset{n}{\underset{j=1}{\operatorname{med}}}(x_{(j)}) \overset{\mathrm{def}}{=} \begin{cases} x_{(\ell+1)} & n = 2\ell+1 \\[2ex] \displaystyle \frac{x_{(\ell)}+x_{(\ell+1)}}{2} & n = 2\ell \\ \end{cases} $$ Si hay un número impar de elementos, $2\ell +1$ , se toma el elemento $(\ell+1)-$ésimo tal que hay $\ell$ elementos menores o igual que ese, y otros $\ell$ elementos mayores o igual que ese.
Si hay un número par de elementos, $2\ell $ , se tiene que tomar el elemento tal que haya $\ell$ elementos menores o igual que ese, y los otros $\ell$ elementos mayores o igual que ese, por lo que se hace la semisuma de los elementos $\ell-$ésimo y $(\ell+1)-$ésimo.

Desde un punto de vista analítico se puede ver como el argumento que minimiza la función $\rho$ , es decir, $\operatorname{med}(x_i) = \operatorname{arg\, min}(\rho)$ , siendo dicha función: $$ \rho(x) = \sum_{i=1}^n | x-x_i | = \sum_{j=1}^n \big| x-x_{(j)} \big| $$ También se puede ver que la mediana $m$ de una variable aleatoria $X$ es el valor de $c$ (argumento) que minimiza la función $\operatorname{E}(|X-c|)$ .

La mediana es un L-estimador (un estimador que es una combinación lineal de estadísticos de orden $k$ siendo estos el $k-$ésimo valor más pequeño de la muestra), por lo que se tiene la siguiente relación de afinidad: $$ \operatorname{med}(a\, x_i + b) = a\, \operatorname{med}(x_i) + b$$ Sin embargo , la mediana de la suma dos listas (o dos vectores) no es la suma de las medianas. Un ejemplo sería $(1,0,0)$ , con mediana $0$ , $(0,2,0)$ , con mediana también $0$ , mientras que la suma $(1,2,0)$ , tiene mediana $1$ , no satisfacen la desigualdad triangular ya que $ 1 \not\leqslant 0 + 0 $. Mientras que $(4,-8,6)$ , con mediana $4$ , $(-9,4,5)$ , con mediana también $4$ , mientras que la suma $(-5,-4,11)$ , tiene mediana $-4$ : se tiene $-4 \leqslant 4 + 4$ , es decir: $$ \operatorname{med}(x_i+y_i) \lesseqgtr \operatorname{med}(x_i) + \operatorname{med}(y_i)$$ La mediana es un caso muy particular de los cuantiles. Para un $p\in [0,1]$ se llama $p-$cuantiles a los valores que parten la muestra en sublistas de igual longitud que representan $p$ frente a la longitud total de la lista. Así pues se tiene (en frecuencia relativa) que el primer $p-$cuantil es mayor o igual que $p$ elementos respecto al total, y menor o igual que $(1-p)$ elementos respecto al total, el segundo $p-$cuantil es mayor o igual que $2p$ elementos respecto al total (y por ende mayor que el primer $p-$cuantil ), y menor o igual que $(1-2p)$ elementos respecto al total...
Cuando $p=0.25$ se conocen como cuartiles (de cuarto) y hay $3$ , mientras que el cuantil$-0.5$ es la mediana.

Cuando los elementos aparecen repetidos, o no todos los $x_i$ tienen la misma importancia, sino que cada uno tiene un peso de $\omega_i$ (y denotamos por $\omega_{(j)}$ es peso asociado a $x_{(j)}$ ) se tiene que la mediana es el elemento $x_k$ que satisface las siguientes desigualdades (aunque solo importan la primera y la última): $$ \sum_{j=1}^{k-1} \omega_{(j)} \leqslant \frac{1}{2} \qquad \sum_{j=1}^k \omega_{(j)} \geqslant \frac{1}{2} \qquad\qquad \sum_{j=k}^{n} \omega_{(j)} \geqslant \frac{1}{2} \qquad \sum_{j=k+1}^{n} \omega_{(j)} \leqslant \frac{1}{2} $$ Para una función de densidad de probabilidad $\operatorname{pdf}-X(x)$ de una variable aleatoria continua $X$ es el valor $m$ que satisface. $$ \int\limits_{(-\infty,m)}\!\! \operatorname{pdf}_X(x)\;\mathrm{d}x \geqslant \frac{1}{2} \qquad \int\limits_{(m,\infty)}\!\! \operatorname{pdf}_X(x)\;\mathrm{d}x \geqslant \frac{1}{2} $$ Volviendo al tema de media, moda y mediana, si denotamos como $\mu$ la media, $M$ la moda y $m$ la mediana con $\sigma$ la desviación típica tenemos la desigualdad: $ |\mu-m| \leqslant \sigma$ , y si la distribución es unimodal se tiene una desigualdad más estricta $\displaystyle |\mu-m| \leqslant \sqrt{\frac{3}{5}\;}\sigma$ , y además $|m-M| \leqslant \sqrt{3\;}\sigma$ . Además para distribuiones unimodales moderamente asimetricas se tiene la aproximación $\displaystyle \mu-M \approx 3(\mu-m) \iff m \approx \frac{2}{3}\mu + \frac{1}{3}M $ .

¿Existe una desigualdad tipo Jensen para la mediana? Sí, pero es ligeramente diferente. Tomemos $ f : A \to B $ con $A,B \subset \mathbb{R}$ y decimos que es función tipo $\mathcal{C}$ , $f\in\mathcal{C}(A,B)$ , si $\forall t\in B$ se tiene que $f^{[-1]}\big( \,(-\infty, t]\, \big) = \big\{ x \in \mathbb{R} \mid f(x) \leqslant t \big\}\subset A$ es un conjunto cerrado. Entonces se cumple $f\big(\operatorname{med}(X)\big) \leqslant \operatorname{med}\!\big(f(X)\big)$ . Un ejemplo de función tipo $\mathcal{C}$ son las convexas.

La mediana es muy útil, pero a veces muy costosa de calcular computacionalmente ya que hay que ordenar el conjunto para calcularla. Sin embargo, el matemático Tukey propuso en 1978 un estimador de la mediana que solo necesita hacer cuatro comprobaciones de tres elementos cada una: ninther ("novenero") . Sin embargo este estimador depende de cómo sea el conjunto $A$ , ya que dará otra estimación si permutan algunos elementos. Para un conjunto $A$ se subdivide en tres de igual longitud, en cada uno se hace la mediana del primer, el del medio y el último elemento, y el ninther es la mediana de estos tres, es decir: $$ a_1 = \operatorname{med}\!\big(A[1],A[n/6],A[n/3]\big) \quad a_2 = \operatorname{med}\!\big(A[n/3+1],A[n/2],A[2n/3]\big) \quad a_3 = \operatorname{med}\!\big(A[2n/3+1],A[5n/6],A[n]\big) \quad \operatorname{ninther}(A)=\operatorname{med}(a_1,a_2,a_3)$$ La mediana a veces también se utiliza en regresión: en mínimos cuadrados ordinarios (OLS por sus siglas en inglés) se minimiza la suma de los cuadrados de los residuales, mientras que en mínima mediana de cuadrados (LMS por sus siglas en inglés) se miniza la mediana de los cuadrados de los residuales. $$ \mathcal{Z}_{OLS} = \sum_{i=1}^n {e_i}^2 = \sum_{j=1}^n {e_{(j)}}^2 \qquad \mathcal{Z}_{LMS} = \overset{n}{\underset{i=1}{\operatorname{med}}}({e_i}^2) = \overset{n}{\underset{j=1}{\operatorname{med}}}({e_{(j)}}^2)$$ Sin embargo, como ya hemos comentado, a diferencia de OLS, está función no es una norma para $n\geqslant 3$.

Puntos $(x,y,z)\in\mathbb{R}^3$ tales que satisfacen $\operatorname{med}(x^2,y^2,z^2)=\operatorname{med}(|x|,|y|,|z|)=1$

 
Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(797) - Sobre el número π (pi)

Hoy \textit{no} es el día del número Pi. Pero me la refanfinfla. uwu $\pi \approx 3.14$ se puede definir de muchos modos. El más geométrico, consiste en decir que $\pi$ es el ratio de proporcionalidad entre el diámetro de una circunferencia y su perímetro. Pero yo prefiero definirlo de otra forma. Sea $\exp$ la única solución (holomorfa) de la siguiente ecuación diferencial: \[ \begin{cases} f \colon \mathbb{C} \to \mathbb{C}\\ f'(z)=z,\\ f(0)=1, \end{cases} \] Entonces, \[ \pi=-2\text{í} \min \{\ |z|\ \colon z\in\mathbb{C}\setminus\{0\},\ \exp(z) = 0 \} \] En la misma vena, la función $\operatorname{arc\,tg}$ permite definir $\pi$ como suma de una serie (la serie de Madhava-Gregory-Leibniz): \[ \pi=4 \sum_{n=1}^\infty \frac{(-1)^n}{2n+1} = 4-\frac{4}{3}+\frac{4}{5}-\ldots \] Esta serie converge bastante lento, lo cuál queda en evidencia por el lento decrecimiento del valor absoluto de sus sumandos. Converge más rápido esta otra: \[ \pi= 2\sum_{k=0}^\infty \frac{2^kk!^2}{(2k+1)!} = 2+\frac{2}{3}+\frac{4}{15}+\ldots \] La derivación de esta fórmula aparece en este enlace. Dejamos al lector la tarea de probar que esta serie converge estrictamente más rápido que la anterior (se debe precisar que es ``estrictamente más rápido'' y utilizar la fórmula de Stirling para probarlo). Comentar también que es una lectura muy interesante para cualquiera interesado en la aproximación computacional de $\pi$. owo Terminamos esta muy breve entrada con una conocida pero igualmente genial anécdota: en 1897, el Congreso del Estado de Indiana, EEUU votó unánimamente a favor de declarar por ley que: ``$\pi$ tiene dos valores: $3,2$ y $4$ '' Por suerte, el Senado del Estado de Indiana tumbó semejante estupidez. Para la historia completa, míralo aquí.

(787) - Sobre el número e

El lector conocerá la definición de número irracional: \[ i \in \mathbb{R}\text{ es irracional } \Leftrightarrow \nexists p,q \in \mathbb{Z}\text{ tal que }i=\frac{p}{q}. \] Es trivial probar que $\sqrt{2}$ es irracional (de ahí su estatus de ejemplo estándar). Pero cuando nos salimos de los números algebraicos, la cosa ya no es tan sencilla. \phantom{meter un espacio así} Aquí nos fijaremos en el número $e$. Recordemos que \[ e=\sum_{k=0}^\infty \frac{1}{k!} = 1+1+\frac{1}{2!}+\frac{1}{3!}+\ldots \approx 2,718. \] Supongamos que $e$ fuera racional, $e=p/q$ donde $p$ y $q$ no comparten factores primos. Llamemos $\displaystyle S_n=\sum_{k=0}^n \frac{1}{k!}$ a las sumas parciales. Se pueden dar dos casos:

• que $e=S_n$ para cierto $n$.
• que $e \neq S_n$ para ningún $n$.

El primer caso es obvio que no va a ser. Para estudiar el segundo caso, \mbox{consideremos} el resto $e-S_n$: \begin{align*} e-\sum_{k=0}^n \frac{1}{k!} &= \frac{1}{(n+1)!}\left( 1 + \frac{1}{n+2} + \frac{1}{(n+2)(n+3)}+\ldots \right)\\ &\leqslant \frac{1}{(n+1)!}\left( 1 + \frac{1}{n+2} + \frac{1}{(n+2)^2}+\ldots \right)\\ &\leqslant \frac{1}{(n+1)!}\left( 1 + \frac{1}{n+1} + \frac{1}{(n+1)^2}+\ldots \right) = (\#) \end{align*} Ahora bien, $\displaystyle \left(1+\frac{1}{n+1}+\frac{1}{(n+1)^2}+\ldots\right)$ es una serie geométrica. Por tanto \[ (\#)=\frac{1}{(n+1)!}\frac{1}{1-\frac{1}{n+1}}=\frac{1}{(n+1)!}\frac{n+1}{n}=\frac{1}{n}\frac{1}{n!}. \] Si $e$ fuera racional, entonces $r=e-S_n$ también sería racional. Pero la cuenta anterior nos dice que \[ e-S_n \leqslant \frac{1}{n}\frac{1}{n!} \xrightarrow[\substack{n\to\infty}]{} 0, \] lo que implica que $r=0$, absurdo. está mal hecho, pero son las 5 de la mañana
Probar que $e$ es trascendente (es decir, no algebraico) es poco trivial (enlace). Un poco más fácil es probar que $\pi$ es irracional (enlace). Por último, probar que $\pi$ es trascendente es aún más interesante (léase, complicado).

(773) - Caos deterministíco. No podemos predecir el futuro, ni saber el pasado

En Occidente tenemos la noción de representarnos mirando hacia el futuro, y de espaldas al pasado, pues tenemos la percepción de que el pasado lo hemos dejado atrás y estamos yendo hacia el futuro. Sin embargo, en algunas culturas como la aymara en los Andes pasado y delante tienen la misma palabra (nayra), y futuro y espaldas comparten otra (quipa), ya que según su cosmovisión el pasado, al haberlo vivido y experimentado lo conocemos bien, mientras que el futuro al ser el porvenir, no. Los antiguos germanos, los romanos tenían parte de esta filosofía en sus lenguas (en especial en algunos prefijos temporales que originalmente eran locales), así como en vasco.

En $1814$ Laplace publica «Ensayo filosófico sobre las probabilidades» (Essai philosophique sur les probabilités) donde propone la idea de lo que acabaría conocienco como el demonio de Laplace, un intelecto superior. El término demonio quizá no es el más acertado, originalmente es del griego antiguo pero con otro significado: δαίμων - daímōn «divinidad, genio, espíritu protector». Laplace argumentaba que dado que las ecuaciones que rigen la mecánica clásica con puramente deterministas, dadas unas condiciones iniciales, es posible saber el estado del sistema en cualquier momento a futuro. Es más, como a su vez las ecuaciones se pueden utilizar a pasado, es decir, no hay preferencia temporal, dadas las mismas condiciones iniciales es posible saber de qué estado proviene el sistema.
Esto es bastante sorprendente, porque si bien una instantánea de las posiciones de un sistema nos dice momentáneamente cómo es dicho sistema, podemos averiguar de qué estado venía y cómo va evolucionar.

Sin embargo, a pesar de que Laplace no tuvo en cuenta que su idea no es siempre aplicable por efectos cuánticos - que él desconocía - no tuvo en cuenta el caos. Llamamos caos en matemáticas al grado de sensibilidad de las condiciones iniciales, es decir, cómo varía el sistema si comienza con unas condiciones iniciales similares pero no idénticas (tras una perturbación diferencial). Muchos sistemas dinámicos son caóticos como el problema de los $N$ cuerpos - cuyo mejor ejemplo es el Sistema Solar - , el movimiento del péndulo doble entre otros. En ambos casos, dado una situación inicial es posible saber en todo momento cómo evoluciona el sistema, pero si se varían ligeramente las condiciones iniciales, para un mismo instante ambos sistemas con categóricamente diferentes.
La incertidumbre en el conocimiento "profundo del sistema" (la precisión de los datos) en un instante inicial (presente), nos imposibilita asegurar a futuro con plena certeza el estado del sistema en cualquier momento, en especial cuanto más separado esté del presente. Por ello a la par que se va resolviendo las ecuaciones diferenciales que presenta el sistema hay que ir comprobando empíricamente cómo está el sistema.

Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(769) - Expresiones para soluciones de ecuaciones diferenciales. ¿Existen o se crean?

Poco antes de verano estaba con un amigo del doble grado y redactor de este blog hablando de ecuaciones diferenciales y salieron dos ejemplos en particular: la del péndulo realista, $y^{\prime\prime}+\beta\,y^\prime+{\omega_0}^2\sin(y)=0$ (con amortiguamiento/fricción y sin aproximaciones de ángulo pequeño), y el sistema de Lotka-Volterra: $\displaystyle \begin{cases} x^\prime & = & x(\alpha-\beta y) \\ y^\prime & = & -y(\gamma-\delta x) \end{cases}$ . Él me preguntó que si acaso tenían solución, y yo le dije que sí, que se podían construir. Él realmente no me inquirió en si eran resolubles (ya fueren analítica o numéricamente), sino si existía alguna función conocida, que pudiese ser composición y/o combinación de otras funciones elementales o no (por ejemplo) tal que fuese solución suya, es decir, si se podía expresar su solución en términos de funciones conocidas. Él se quejó de mi respuesta, a lo que yo le dije que realmente eso poco importaba a la hora de resolver ecuaciones diferenciales y al estudiarlas. A veces la propia solución puede ser tan engorrosa analíticamente que es reconfortante pensar que siempre nos quedará el cálculo numérico...

Por ejemplo, tomemos la función error gaussiano, una función no-elemental, que se define como: $\displaystyle \operatorname{erf}(x)\overset{\text{def}}{=}\frac{2}{\sqrt{\pi\;}}\int_0^x e^{-t^2}\;\text{d}t$ . Resolvamos la ecuación diferencial ordinaria de I orden $y^\prime = e^{-x^2}$ . Ahora bien, ¿la solución de esta ecuación diferencial es debida a dicha definición anterior, tal que $\displaystyle y\triangleq\frac{\sqrt{\pi\;}}{2}\operatorname{erf}(x)+y_0$ , o bien la definición anterior se dio para que satisficiere la ecuación diferencial $y^\prime = e^{-x^2}$ ?, es decir, ¿qué vino antes: la definición (y después la ecuación diferencial que resulta ser que la tiene como solución salvo constantes) o la ecuación diferencial (y la definición se dio posteriormente como su solución a este problema)?

Sin la definición de $\operatorname{erf}(x)$, el problema es irresoluble analíticamente, pero con dicha definición es trivial, ya que la propia definición es la solución. Es como si hubiésemos definido que esta función sea su solución. Es equivalente a construir la solución, y darla como definición, por lo que el problema ya está resuelto.

Esto puede parecer una tontería, o un poco enrevesado, pero si se quiere resolver $y^\prime = 3x^2$ , que tiene como solución $y=x^3+y_0$ , no se requiere de ninguna definición previa que no esté contemplada ya.

Uno siempre puede decir: Sea $y\overset{\text{def}}{=}\theta(t)=\theta(t|\beta,\omega_0;y_0,{y_0}^\prime)$ la solución a la ecuación diferencial $y^{\prime\prime}+\beta\,y^\prime+{\omega_0}^2\sin(y)=0$ en la variable $t$ donde $\beta,\omega_0\geqslant 0$ son parámetros que dependen de la ecuación diferencial, con $y_0,{y_0}^\prime$ como condiciones iniciales (hay soluciones analíticas si se linealiza o si al menos uno de los parámetros se anula, aunque con $\beta=0$ es enrevesada). Después mediante cálculo y análisis numéricos, con un conocimiento previo de ecuaciones diferenciales, uno puede estudiar la solución numéricamente, y a veces a lo sumo, hasta cierto punto analíticamente.

Otro ejemplo puede ser el oscilador de Van der Pol: $y^{\prime\prime}-\mu(1-y^2)y^\prime+{\omega_0}^2y=0$ (nótese que el caso $\mu=0$ es el oscilador armónico simple, por lo que sí tiene una solución analítica en este caso), otra vez ahí definir $y\overset{\text{def}}{=}\operatorname{vdp}(t)=\operatorname{vdp}(t|\mu,\omega_0;y_0,{y_0}^\prime)$ como una solución obtenida numéricamente y a partir de ella estudiar la solución "analítica".

Lo que estoy argumentando es equivalente a si un matemático en el $\text{siglo XVI-XVII}$, desconociendo tanto el número $e$ como las funciones exponencial $e^x$ y logaritmo $\ln(x)$ , decidiese analizar la ecuación diferencial ordinaria de primer orden $y^\prime=k\, y$ definiendo su solución como $y\overset{\text{def}}{=}\operatorname{sol}(t) = \operatorname{sol}(t|k;y_0)$ y estudiarla a partir de su solución construida. Hoy sabemos que tiene una forma cerrada, $y=y_0\, e^{k(t-t_0)}$ .

Uno puede incluso argumentar que estos ejemplos están sesgados y que no tienen validez, pero no es así, hay muchas ecuaciones diferenciales que no tienen necesariamente una solución analítica y para entenderlas es necesario su estudio con cálculo numérico y definir sus soluciones. En especial, en física, como dijo Steven Strogatz ($1959-$): "Desde Newton, la humanidad se ha dado cuenta que las leyes de la física siempre se expresan en la lengua de ecuaciones diferenciales." Veamos algunos ejemplos:

El sistema de $2$ ecuaciones diferenciales ordinarias de II orden que rigen el péndulo doble $$ \begin{cases} \ddot{\theta}_1 & = & \displaystyle \frac{-g (2m_1+m_2)\sin(\theta_1)-m_2g\sin(\theta_1-2\theta_2)-2\sin(\theta_1-\theta_2)m_2\big({\dot{\theta_2}}^2 l_2+{\dot{\theta_1}}^2 l_1\cos(\theta_1-\theta_2)\big)} {l_1\big(2m_1+m_2-m_2\cos(2\theta_1-2\theta_2)\big)} \\ \ddot{\theta}_2 & = & \displaystyle \frac{2 \sin(\theta_1-\theta_2) \big({\dot{\theta}_1}^2 l_1 (m_1 + m_2)+ g(m_1 + m_2) \cos(\theta_1) + {\dot{\theta}_2}^2 l_2 m_2 \cos(\theta_1 - \theta_2)\big)} {l_2 \big(2 m_1 + m_2 - m_2 \cos(2 \theta_1 - 2 \theta_2)\big)} \end{cases} $$ El sistema de $2$ ecuaciones diferenciales ordinarias de II orden que rigen el péndulo elástico $$ \begin{cases} \ddot{r} & = & \displaystyle (r+l_0)\dot{\theta}^2 -\frac{k}{m}r + g\cos(\theta) \\ \ddot{\theta} & = & \displaystyle -\frac{g\sin(\theta)+2\dot{r}\dot{\theta}}{r+l_0} \end{cases} $$ La ecuación diferencial ordinaria de II orden que rige el péndulo esférico $$ \ddot{\theta}(t)- \left(\frac{p_\varphi}{ml^2}\right)^2\frac{\cos\big(\theta(t)\big)}{\sin^3\big(\theta(t)\big)}+\frac{g}{l}\sin\big(\theta(t)\big) = 0 $$ El sistema de $6$ ecuaciones diferenciales ordinarias de II orden que rigen el problema de los dos cuerpos $$ \begin{cases} \ddot{\vec{r}}_1 & = &\displaystyle Gm_2 \frac{\vec{r}_2-\vec{r}_1}{{\|\vec{r}_2-\vec{r}_1\|_2}^3}\\ \ddot{\vec{r}}_2 & = &\displaystyle Gm_1 \frac{\vec{r}_1-\vec{r}_2}{{\|\vec{r}_1-\vec{r}_2\|_2}^3} \end{cases} \implies m_1\ddot{\vec{r}}_1 + m_2\ddot{\vec{r}}_2 = \vec{0} $$ El sistema de $3$ ecuaciones diferenciales ordinarias de I orden conocido como el atractor de Lorentz, un ejemplo famoso de la Teoría del Caos: $$ \begin{cases} \dot{x} & = & \sigma (y-x) \\ \dot{y} & = & x (\rho-z)-y \\ \dot{z} & = & xy-\beta z \end{cases} $$
Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(761) - Mínimos cuadrados en machine learning. Least trimmed squares (LTS)

En el día de hoy veamos cómo se utiliza la regresión lineal en machine learning. Recapitulemos un poco cómo funcionan los mínimos cuadrados ordinarios, ordinary least squares u OLS por sus siglas en inglés, el caso más simple: Dadas dos sucesiones de valores, una predictora $\{x_i\}_{i=1}^n$ (abscisas), y otra respuesta $\{y_i\}_{i=1}^n$ (ordenadas), se buscan los argumentos (parámetros) que minimicen la siguiente restricción $Z_{OLS}$ (del alemán Zwang - obligación): $$ Z_{OLS} = \sum_{i=1}^n \big(y_i-(\beta_0+\beta_1x_i)\big)^2 \qquad \hat{\beta}_0,\hat{\beta}_1 = \operatorname{arg\,min} Z_{OLS}$$ Nótese que esta restricción siempre es no-negativa, $Z_{OLS}\geqslant 0$ , y la igualdad (que sería lo deseable) solo se da si el ajuste es perfecto y no hay ninguna desviación.

Dado que queremos hallar la recta que mejor se ajuste a $\{y_i\}_{i=1}^n$ en la expresión aparece $\beta_0+\beta_1x_i$ , sino, si supiésemos o si intuyésemos que los $\{y_i\}_{i=1}^n$ siguen una función $f(x)$ genérica (polinómica, sinusoidal, exponencial,...), entonces se buscaría qué parámetros de $f(x)$ minimizan $\displaystyle Z_{OLS} = \sum_{i=1}^n \big(y_i-f(x_i)\big)^2$ . Se suele denotar por $\hat{y}_i=f(x_i)$ a los valores predichos. El método de mínimos cuadrados busca los argumentos, i.e. parámetros, $\hat{\beta}_0,\hat{\beta}_1,\cdots$ que minimizan la suma total de los cuadrados de todas las $n$ diferencias entre los valores reales $y_i$ , y los valores predichos $\hat{y}_i$ , $Z_{OLS}$ , (esta diferencia es el residual $r_i\overset{\text{def}}{=}y_i-\hat{y}_i$, por eso realmente se suele poner como hallar los argumentos que minimizan $\displaystyle Z_{OLS}=\sum_{i=1}^n {r_i}^2$ ). En nuestro caso para la regresión de una recta si denotamos como $\hat{\beta}_0,\hat{\beta}_1$ a la ordenada en el origen y pendiente predichas según el modelo respetivamente, se tiene que $\hat{y}_i = \hat{\beta}_0+\hat{\beta}_1x_i$ .

Sin embargo, en machine learning no se utiliza esto, al menos no así. Se hace lo que se conoce como mínimos cuadrados recortados, least trimmed squares o LTS por sus siglas en inglés: Primero se preselecciona un número $k$ de puntos que se tendrán en cuenta (y los $(n-k)$ restantes para un estudio posterior), donde $k$ suele estar entre el $75\%-80\%$ del total $n$ de puntos. Empero dichos $k$ puntos que se toman no son al azar, sino que se escogen de una manera muy particular: de entre todos los $n$ puntos $\big\{(x_i,y_i)\big\}_{i=1}^n$ , se busca entre todas las $\displaystyle \binom{n}{k}$ posibles combinaciones (de $k$ elementos tomados de un supraconjunto de $n$ elementos) aquella cuya restricción $Z$ de $k$ sumandos sea la mínima de todas las posibles, es decir, en cualquier otro subconjunto de tamaño $k$ sobre el supraconjunto $n$ se ajustan peor entre sí a mínimos cuadrados (su restricción $Z$ es mayor).

Sea $\mathcal{Z}$ el conjunto de las $\displaystyle |\mathcal{Z}|=\binom{n}{k}$ posibles mínimas sumas parciales de los cuadrados de $k$ residuales, $\displaystyle \mathcal{Z} \overset{\text{def}}{=} \bigg\{ \min\!\sum_{i=1}^k {r_{\sigma(i)}}^2 \;|\; k\leqslant n\, ,\,\sigma\in S_n\bigg\}$ con $k$ fijo , donde $\sigma(\cdot)$ indica una permutación de los $n$ puntos. Buscamos qué argumentos minimizan $Z_{LTS}\overset{\text{def}}{=}\min\{\mathcal{Z}\}\geqslant 0$ . Por esta definición el coeficiente de correlación $R^2$ de $Z_{LTS}$ es el máximo de los del conjunto $\mathcal{Z}$ , es decir, $R^2(Z_{LTS})\triangleq\max\big\{R^2(\mathcal{Z})\big\}$ .

Es más, si denotamos por $\big\{{r_{(i)}}^2\big\}_{i=1}^n$ la sucesión de los cuadrados de los residuales ordenados de menor a mayor, es decir, ${r_{(i)}}^2\leqslant {r_{(i+1)}}^2$ , buscamos qué parámetros minimizan la suma parcial $\displaystyle Z_{LTS}=\sum_{i=1}^k {r_{(i)}}^2$ , que es menor o igual que la suma total $\displaystyle \sum_{i=1}^n {r_{(i)}}^2=\sum_{i=1}^n {r_i}^2 = Z_{OLS}$ . Esta última restricción es como ya hemos dicho, la que se busca minimizar en OLS. La suma parcial será igual a la suma total si y solo si $r_{(i)}=0 \quad i = 1,\cdots, n$, es decir, si todos, los $n$ puntos, están perfectamente alineados.

Para comparar ambos resultados, ya hemos visto que $0\leqslant Z_{LTS}\leqslant Z_{OLS}$ , pero estaría bien comparar sendos errores cuadráticos medios, es decir, $\displaystyle \frac{1}{k} Z_{LTS}$ frente a $\displaystyle \frac{1}{n} Z_{OLS}$ , o también $\displaystyle \frac{1}{k-p} Z_{LTS}$ frente a $\displaystyle \frac{1}{n-p} Z_{OLS}$ (que se pueden entender como estimadores de las varianzas de sendos residuales) donde $p$ es el número de parámetros que se estiman (en una recta por ejemplo son dos: la ordenada en el origen y la pendiente), siendo $(k-p)$ y $(n-p)$ sendos grados de libertad.

A este conjunto de $k$ elementos que minimizan la suma parcial de los $k$ cuadrados de los residuales entre todas las posibles se lo conoce como conjunto de entrenamiento (training data set en inglés), mientras que el conjunto restante de $(n-k)$ elementos se lo conoce como conjunto de validación (test data set en inglés). El porqué de separar estos dos es debido a que el conjunto con los $(n-k)$ elementos restantes puede estar corrupto, contaminado, o ser defectuoso (son los outliers). Por eso este método es robusto: considera el conjunto de $k$ elementos que mejor se ajustan entre sí, mientras que el conjunto de $(n-k)$ elementos se utiliza para comprobar cómo de bueno es nuestro modelo, si las hipótesis son las correctas, y si hay que replantearse algo. Sin embargo, el conjunto de $(n-k)$ elementos no influyen en hallar los argumentos que minimizan la restricción $Z_{LTS}$ .

Cabe resaltar que $\displaystyle \frac{n^k}{k!}\geqslant \binom{n}{k} \geqslant \frac{(n-k+1)^k}{k!}$ siendo equivalentes cuando $n\to\infty$, es decir, que aumentando $n$ (al aumentar también $k$ ) implica que las posibles elecciones tienen un crecimiento exponencial-potencial, por lo que hallar qué conjunto de $k$ elementos es el adecuado puede ser muy costoso computacionalmente, en especial con $n$ medianos o con algoritmos no-óptimos. Por otro lado intentar estimar el conjunto de los $k$ con una simulación tipo Monte Carlo no parece ser muy fructífera, ya que la probabilidad de acertar aleatoriamente es $\displaystyle \frac{1}{\binom{n}{k}}$ .

Como último comentario, también existe un método que generaliza OLS conocido como mínimos cuadrados ponderados, weighted least squares o WLS por sus siglas en inglés, donde en vez de minimizar $\displaystyle Z_{OLS}=\sum_{i=1}^n {r_i}^2$ , se busca minimizar $\displaystyle Z_{WLS}=\sum_{i=1}^n \omega_i\,{r_i}^2$ donde $\{\omega_i\}_{i=1}^n$ es una sucesión de pesos no-negativos con cada $\omega_i$ asociado a su respectivo $r_i$ (realmente a $x_i$ , pero es equivalente) y normalmente normalizados, $\displaystyle \sum_{i=1}^n \omega_i=1$ . Nuestro modelo de LTS es en realidad un caso muy particular de WLS donde $\displaystyle\omega_{\sigma(i)}=\begin{cases} \displaystyle \frac{1}{k} & i\in\{1,\cdots,k\} \\ 0 & i\in\{k+1,\cdots,n\} \end{cases}$ . Es más, se podría generalizar esto aún más a un modelo de mínimos cuadrados ponderados y recortados, weighted least trimmed squares o WLTS por sus siglas en inglés donde se buscan: $$\operatorname{arg\,\min} \sum_{i=1}^k \omega_{(i)}\,{r_{(i)}}^2 $$
Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(757) - Cuando Gauss usó la Estadística para la Física y creó una nueva rama de la Mecánica

Es innegable que Sir Isaac Newton ($1643-1727$) revolucionó la Física con la publicación de su Principia Mathematica en $1687$ donde no solo introdujo la Ley de gravitación universal, sino que además introdujo tres axiomas que hoy día se conocen como las (tres) leyes de Newton, los postulados de la Mecánica newtoniana. Cabe resaltar la Ley Fundamental de la Dinámica (II Ley), que se suele presentar como $\vec{F}=m\vec{a}$ , o al decir que la aceleración es la segunda derivada respecto al tiempo de la posición, $\vec{F}=m\ddot{\vec{r}}$ , o sea, $\displaystyle \ddot{\vec{r}}=\frac{\vec{F}}{m}$ .
Toda la mecánica newtoniana trata sobre fuerzas y cómo afectan a las partículas. Sin embargo, desde finales del $\text{siglo XVIII}$ hay varias reformulaciones de la mecánica (clásica) que trabajan en su mayoría las energías en vez de las fuerzas, por ejemplo de mecánica lagrangiana (de $1788$ , donde se introducen coordenadas generalizadas), la hamiltoniana (de $1833$ , donde introduce además momentos conjugados), la routhiana (de $1855/1877$ , que es una mezcla de la dos últimas), y la que vino con la ecuación de movimiento de Gibbs-Appel (de $1879/1900$).
Sin embargo, Gauss ($1777-1855$) propuso aplicar su procedimiento de mínimos cuadrados (descubierto independientemente por él en $1809$ y por Legendre en $1805$) a una versión de la II Ley de Newton, en su artículo de $1829$ Über ein neues allgemeines Grundgesetz der Mechanik (Sobre una nueva constitución universal de la Mecánica). Lo llamó Prinzip des kleinsten Zwanges, es decir, Principio de la mínima restricción (obligación).
Supongamos que tenemos un sistema de $n$ partículas que sabemos sendas posiciones y velocidades en un instante determinado (el inicial e.g.), $\big\{\vec{r}_i(t_0),\dot{\vec{r}}_i(t_0)\big\}_{i=1}^n$ , sendas masas $m_i$ sobre las que se aplican sendas fuerzas (conocidas) $\vec{F}_i$ . En un sistema sin restricciones en virtud de la II Ley de Newton, se tendrá que $\displaystyle\ddot{\vec{r}}_i=\frac{\vec{F}_i}{m_i}$ . Gauss definía la restricción $Z$ como: $$ Z\Big(\big\{\ddot{\vec{r}}_i\big\}_{i=1}^n\Big) \overset{\text{def}}{=}\frac{1}{2}\sum_{i=1}^n m_i \left(\ddot{\vec{r}}_i-\frac{\vec{F}_i}{m_i}\right)^2 = \sum_{i=1}^n \frac{1}{2m_i} (m_i\ddot{\vec{r}}_i-\vec{F}_i)^2 \geqslant 0$$ Gauss lo describió como el principio dinámico más fundamental, ya que establece que la verdadera aceleración instantánea de cada una de las $n$ partículas del sistema es puntualmente el argumento que minimiza dicha restricción $Z$ , es decir, $\big\{\ddot{\vec{r}}_i\big\}_{i=1}^n\ = \operatorname{arg\,min} Z$ . Nótese que $Z\Big(\big\{\ddot{\vec{r}}_i\big\}_{i=1}^n\Big) = Z\big(\ddot{\vec{r}}_1,\cdots,\ddot{\vec{r}}_n\big)$ es la función restricción de un conjunto de aceleraciones $\big\{\ddot{\vec{r}}_i\big\}_{i=1}^n$ , la cual hay que minimizar . Si el sistema no está sujeto a una restricción, entonces $Z=0$ y evoluciona siguiendo las ecuaciones del movimiento de Newton.

Realmente el Principio de Gauss de la mínima restricción no utiliza mínimos cuadrados ordinarios (OLS por sus siglas en inglés), ya que supondría que todas las $n$ aceleraciones tienen la misma importancia en $Z$ (todos los $n$ sumadnos realmente), sino que hace uso de mínimos cuadrados ponderados (WLS por sus siglas en inglés), cuya ponderación individual es la masa $m_i$ . El caso particular de que todas las masas fuesen iguales sí que sería OLS.

Nótese que la solución minimizante es un total de $n$ vectores tridimensionales, por lo que buscamos un total de $3n$ argumentos que minimizan la restricción $Z$ , por ende puede ser un proceso computacionalmente costoso para sistemas de muchas partículas. Si se quiere resolver numéricamente, sin hacer uso de la solución general de mínimos cuadrados, se puede usar el algortimo de bassinhopping y tomando como iterante inicial la secuencia $\displaystyle \left\{\frac{\vec{F}_i}{m_i}\right\}_{i=1}^n$ , o una pequeña perturbación cercana a ella .

Sin embargo, a veces es preferible hacer un cambio de variables, en especial en un sistema restringido, por lo que consideraremos las nuevas variables $\vec{\rho}_i = \sqrt{m_i\,}\,\vec{r}_i$ , y $\displaystyle \vec{\phi}_i =\frac{1}{\sqrt{m_i\;}}\,\vec{F}_i$ , por ende la II Ley de Newton queda ahora como $\vec{\phi} \triangleq \ddot{\vec{\rho}}$ . Las coordenadas $\{\vec{\rho}_i\}_{i=1}^n$ , posiciones de Jacobi, están relacionadas con la métrica de Jacobi, por lo que este sistema de coordenadas se suele llamar sistema de referencia de Jacobi (Evans and Morriss, $1990$). La restricción $Z\Big(\big\{\ddot{\vec{\rho}}_i\big\}_{i=1}^n\Big)$ es (ahora sí OLS): $$ Z\Big(\big\{\ddot{\vec{\rho}}_i\big\}_{i=1}^n\Big) \triangleq\frac{1}{2}\sum_{i=1}^n \big(\ddot{\vec{\rho}}_i-\vec{\phi}_i\big)^2\geqslant 0$$ El Principio de Gauss de la mínima restricción establece que las trayectorias que realmente se siguen son las que se desvían tan poco como sea posible, en el sentido de mínimos cuadrados (ordinarios o ponderados), respecto de las trayectorias newtonianas sin-restricción. Hertz se refirió a la función $Z$ como el cuadrado de la curvatura.

Realmente el Principio de Gauss de la mínima restricción se puede generalizar: en vez de usar mínimos cuadrados, usar otro tipo de regresión con el mismo fin, hallar las trayectorias con restricción que menos se desvíen en el sentido de dicha regresión respecto de las trayectorias newtonianas sin-restricción.

Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(751) - Racionales en trigonometría de racionales. Teorema de Niven-Hadwiger

La pregunta de hoy es bien simple: ¿Qué ángulos que son un número racional de vueltas tienen como seno, coseno o tangente también un número racional?Este resultado se concoce como Teorema de Niven ($1915-1999$) de $1956$, pero el matemático Hadwiger ($1908-1981$) ya hizo una demostración en $1948$.
Una primera idea descartable es argumentando como las funciones trigonométricas se pueden expresar como series, pero la serie de racionales no es necesariamente racional ( el ejemplo más claro $\displaystyle \frac{1}{n!}\in\mathbb{Q}$ pero $\displaystyle e=\sum_{n=0}^\infty \frac{1}{n!} \not\in\mathbb{Q}$ ).

Otra idea sería considerar los polinomios de Chebyshov de I especie [Чебышёв - Čebyšëv], polinomios de coeficientes enteros, $T_N(x)\in\big(\mathbb{Z}[x]\big)_N $ , que satisfacen la relación: $T_N\big(\cos(\theta)\big)=\cos(N\theta)$ . Sin embargo, este método solo nos dice que $\cos(2\pi\theta)\in\mathbb{Q} \implies \cos(2\pi N\theta)\in\mathbb{Q}$ , es decir, si el coseno es racional, el coseno de múltiplos de ángulo también lo es. No podemos decir lo mismo de los de II especie $U_N(x)\in\big(\mathbb{Z}[x]\big)_N $ , que satisfacen la relación: $\sin(\theta)U_{N-1}\big(\cos(\theta)\big)=\sin(N\theta)$ .

La idea es buscar los conjuntos maximales $\varnothing\subset Q_{0,r},Q_{1,r}\subset\mathbb{Q}$ donde $Q_{1,r}\overset{\text{def}}{=}\operatorname{r}(2\pi Q_{0,r})$ para alguna razón trigonométrica $\operatorname{r}$ , es decir, hallar los ángulos racionales, $\varphi\in\mathbb{Q}$ , tales que alguna razón trigonométrica es racional, $\operatorname{r}(2\pi\varphi)\in\mathbb{Q}$ .
Si $\varphi\in Q_{0,r}\subsetneq\mathbb{Q} \implies \operatorname{r}(2\pi\varphi)\in\mathbb{Q}$ , es decir que si el ángulo $\varphi$ es "racional de Niven", su razón trigonométrica también es racional.
Si $\varphi\in(\mathbb{Q}\setminus Q_{0,r})\subsetneq\mathbb{Q} \implies \operatorname{r}(2\pi\varphi)\in(\mathbb{R}\setminus\mathbb{Q})$, es decir que si el ángulo $\varphi$ es racional pero no "[racional] de Niven", su razón trigonométrica es estrictamente irracional.
Si $\varphi\in (\mathbb{R}\setminus\mathbb{Q})\subsetneq\mathbb{R} \implies \operatorname{r}(2\pi\varphi)\in \big((\mathbb{R}\setminus\mathbb{Q})\bigcup \hspace{ -7pt }\raise-.5ex{\scriptsize | } \hspace{3pt} (\mathbb{Q}\setminus Q_{1,r})\big) \triangleq (\mathbb{R}\setminus Q_{1,r}) $ , es decir que si el ángulo $\varphi$ es irracional, su razón trigonométrica también es o bien irracional o bien racional.

Si para algún ángulo el seno o coseno es $\displaystyle \frac{p}{q}$ donde $0\leqslant |p| \leqslant |q| $ y $p,q\in\mathbb{Z}$ , entonces el otro es $\displaystyle \pm\frac{\sqrt{q^2-p^2\;}}{q}$ . La pregunta es entonces $\sqrt{q^2-p^2\;}\overset{\text{?}}{\in}\mathbb{N}$ . Esta pregunta es equivalente a preguntar si existe una terna pitagórica con hipotenusa $|q|$ y un cateto $|p|$ . Al final resulta ser que los únicos que seno y coseno son racionales son las soluciones triviales.

Veamos los pocos valores que satisfacen la relación en $\displaystyle 0\leqslant \theta\leqslant \frac{\pi}{2}$ (en otros cuadrantes solo hay que tener en cuenta las relaciones de los demás cuadrantes con el primero):

Para el seno se tiene $\displaystyle 0,\frac{\pi}{6},\frac{\pi}{2}$ que valen respectivamente $\displaystyle  0,\frac{1}{2},1$ .

Para el coseno se tiene $\displaystyle 0,\frac{\pi}{3},\frac{\pi}{2}$ que valen respectivamente $\displaystyle  0,\frac{1}{2},1$ .

Para la tangente se tiene $\displaystyle 0,\frac{\pi}{4}$ que valen respectivamente $\displaystyle  0,1$ .

Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(743) - Resolver ecuaciones diferenciales numéricamente en Excel

A finales de septiembre estaba ayudando a un novato de física con un ejercicio de tiro parabólico, pero con amortiguamiento por la fricción con el aire. Como él no sabía programar nada, le hice una pequeña simulación en Excel (sí Excel, ¿qué pasa?) para poder cambiar los parámetros al instante. Quiero exponer aquí cómo lo hice. Aviso que este método no es muy bueno numéricamente, en particular porque Excel no está pensado para resolver ecuaciones diferenciales. Esto es rápido, para tener una idea preliminar, no para tener exactitud. Para ver cómo funciona, hay que repasar algunos conceptos:
Recordemos la definición de derivada [analítica]: $$ f^{(1)}(x) \overset{\text{def}}{=} \lim_{h\to0} \frac{f(x+h)-f(x)}{h}$$ Si en vez de tomar el límite, se toma un $h$ suficientemente pequeño, podemos calcular lo que se llama en computación derivada numérica (una aproximación numérica de la derivada analítica). Esto nos permite aproximar la función $f$ en un punto $x+h$ si sabemos $f(x)$ y $f^{(1)}(x)$ : $$ f^{(1)}(x) \approx \frac{f(x+h)-f(x)}{h} \implies f(x+h) \approx f(x) + f^{(1)}(x)\cdot h$$ Es decir, podemos aproximar (casi) cualquier función $f(t)$ por una recta tangente en $t=a$ , $T_1\big(f,a\big)(t)$ , con el error que decrece hasta hacerse nulo en $t=a$ . $$ f(t) = T_1\big(f,a\big)(t)+{\scriptstyle \mathcal{ O } }\big( (t-a)^1\big) \qquad T_1\big(f,a\big)(t) = f(a)+ f^{(1)}(a)\cdot (t-a)$$ En la resolución numérica se suele discretizar la función (aunque no sea discreta) en intervalos de longitud pequeña, $[x_{n-1},x_n]$ . Vamos a tomar la longitud de cada uno de estos intervalo igual (nodos equiespaciados) de longitud $h$ . Cuanto más cercano a $0$ sea $h$ , mejor la aproximación, pero más cálculos son necesarios. Llamemos $y_n$ a la función $f(x)$ evaluada en el $n-$ésimo nodo $x_n$ , $f(x_n)$ , y denotemos por ${y_n}^{(k)}$ a la $k-$ésima derivada de la función $f(x)$ evaluada en el $n-$ésimo nodo $x_n$ , $f^{(k)}(x_n)$ . Es decir: $$ \begin{matrix} h & = & x_n-x_{n-1} & = & \Delta x_n \\ y_n & = & f(x_n) \\ {y_n}^{(k)} & = & f^{(k)}(x_n) & & k\in\mathbb{N} \end{matrix} $$ Es decir, cuando solo se tiene la primera derivada: $$ y_{n+1} = y_n + {y_n}^{(1)} h $$ Si tenemos una función que conocemos ambas derivadas (como por ejemplo la posición $r$ con su velocidad inicial $v_0$ y su aceleración $a$ que $\displaystyle r = r_0 + v_0t+\frac{a}{2}t^2 $ ) se tiene que: $$y_{n+1} = y_n +{y_n}^{(1)} h + \frac{{y_n}^{(2)}}{2}h^2 = y_n + h\left( {y_n}^{(1)} + \frac{{y_n}^{(2)}}{2}h\right) $$ Al final lo que estamos haciendo es aproximar las sucesivas derivadas $\displaystyle \frac{\text{d}^{(n)}y}{\text{d}x^{(n)}}$ como diferencia finitas $\displaystyle \frac{\Delta^{[n]} y}{\Delta x^n}$ , hacer una expansión de Taylor localmente en cada punto, $\displaystyle \sum_{n=0}^\infty \frac{1}{n!}\frac{\text{d}^{(n)}y}{\text{d}x^{(n)}}\Delta x^n$ , y parar tras $(N+1)$ sumandos cuando consideramos que nuestro error $\varepsilon_N$ es suficientemente pequeño (disminuye cuanto menor sea $\Delta x = h$ ): $$ \frac{\text{d}y}{\text{d}x} \bumpeq \frac{\Delta y}{\Delta x} \implies \frac{\text{d}^{(n)}y}{\text{d}x^{(n)}} \bumpeq \frac{\Delta^{[n]} y}{\Delta x^n} \implies \frac{1}{n!}\Delta^{[n]} y \bumpeq \frac{1}{n!}\frac{\text{d}^{(n)}y}{\text{d}x^{(n)}}\Delta x^n \\ \sum_{n=0}^\infty \frac{1}{n!}\Delta^{[n]}y = \sum_{k=0}^N \frac{1}{k!}\Delta^{[k]}y +\varepsilon_N$$ Por ejemplo para resolver el péndulo simple, de ecuación $\ddot{\theta} + {\omega_0}^2 \sin(\theta)=0 \implies \ddot{\theta} =-{\omega_0}^2 \sin(\theta) $ . En la primera columna ponemos nuestra variable indepeniente, $t$ , con $A1=0$ y $A2=A1+h$ donde $h$ lo hemos definido antes y es un valor fijo. Para que sea más útil, como Excel permite hacer gráficas de $32.000$ puntos, arrastramos $A2$ hasta $A32000$ . Necesitamos una condiciones iniciales de posición y velocidad que esas las ponemos en $B1$ y $C1$ respectivamente. Repitiendo la notación de antes se tiene que (con $Bn=y_n$ , $Cn={y_n}^{(1)}$ , $Dn={y_n}^{(2)}$ ) : $$ \begin{matrix} {y_n}^{(2)} = -{\omega_0}^2 \sin(y_n) \\ {y_{n+1}}^{(1)} = {y_n}^{(1)} + {y_n}^{(2)}h \\ \displaystyle y_{n+1} = y_n + h\Big( {y_n}^{(1)} + {y_n}^{(2)} \frac{h}{2}\Big) \\ \end{matrix}$$ Este método, aunque efectivo, es poco práctico. Se podría haber usado la integración de Verlet, que es mucho más preciso y no hace falta calcular la primera derivada en todo punto: $$ \begin{matrix} \displaystyle y_1 = y_0+{y_n}^{(1)}h+\frac{{y_n}^{(2)}}{2}h^2 \\ y_{n+1} = 2y_n-y_{n-1}+{y_n}^{(2)}h^2 \qquad n\geqslant 1\\ \end{matrix}$$

Péndulo simple, $\ddot{\theta} + {\omega_0}^2 \sin(\theta)=0$ , y péndulo amortiguado $\ddot{\theta} + \gamma \dot{\theta} + {\omega_0}^2 \sin(\theta)=0$ con mismas condiciones iniciales $\theta_0=2,6$ y $\dot{\theta}_0=0,5$ y con $\gamma=0,25$

Este método es útil a la hora de dar propiedades cualitativas, y no tanto cuantitativas, de la solución y para entender qué pasa con la misma (en el ejemplo del péndulo amortiguado que cuanto mayor sea $\gamma$ , más deprisa tiende a ser idénticamente nula).

  Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.