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ʒ.
