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}$$
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ʒ.
Hola,
ResponderEliminarInteresante artículo sobre este método "para tener una idea preliminar" (sic).
Siguiendo con la reflexión del autor (y sus conclusiones metodológicas implícitas), me lleva a pensar que los métodos de resolución matemática de los ingenieros (telecos y similares) no son al final tan malos como parecen y por supuesto no tanto como para "desprecialos" (vamos, "aquello" de Euler de e ^ i*pi + 1 = 0 aprox= 3 ^ i*3 + 1 = 0).
En cualquier caso, enhorabuena al autor,
Angel.