Cálculo diferencial 6: una solución numérica para el paracaidista

Última parte del recorrido hacia las ecuaciones diferenciales. Para un índice de los artículos sobre este tema, puede consultarse este enlace.


Aquí estamos llegando al final de esta extraña y para nada rigurosa introducción al cálculo diferencial. Y como lo prometido es deuda, veremos de hacer un cálculo numérico que nos de la velocidad en función del tiempo para la caída libre con una atmósfera que limita la velocidad. Veremos aquí solamente la primer parte de la caída, antes de la apertura del paracaídas, la que inicia con velocidad cero y llega a la velocidad límite de la cual hablamos en la primer entrada: la segunda parte de la caída la dejo como ejercicio para el lector.

[Ya que estamos, Alan Eustace ha superado el record de Baumgartner… 😉 ]

Para esta tarea usaremos resultados vistos en la entrada anterior. ¿Ya repasaron? Bien, allí vamos.

Volvamos a escribir nuestra ecuación diferencial, olvidándonos por un momento de que la aceleración es la derivada de la velocidad (lo recordaremos pronto):

\displaystyle M g - k v^2 = M a

y despejemos de allí la aceleración

\displaystyle a = g -\frac{k}{M} v^2

La aproximación que haremos comienza con separar la caída en pequeños intervalos de tiempo \Delta t, todos iguales entre sí. Nada particular, pensará el lector atento, después de todo el cálculo diferencial (y el integral) no es otra cosa que tomar pequeños intervalos y luego hacer un límite para esos intervalos tendiendo a longitud cero. Y efectivamente no hay nada de particular en esta separación: el paso importante está en aproximar la aceleración en cada uno de esos pequeños intervalos por el valor que esta tendría al principio del mismo. Es decir, en lugar de una función que varía en forma continua diremos que la aceleración va tomando valores constantes en cada uno de los intervalos.

Esto claramente es erróneo, ya sabemos que al aumentar la velocidad de caída aumenta la fuerza de fricción con el aire y por lo tanto disminuye la aceleración. Esto hace que los valores de aceleración y velocidad que calculemos con nuestra aproximación sean siempre mayores que los reales, pero si tomamos intervalos de tiempo suficientemente pequeños (ya veremos qué significa «suficientemente») nos irá bastante bien.

Entonces el procedimiento es el siguiente. En el primer intervalo, donde la velocidad es aún pequeña, podemos despreciar el término de fricción para obtener así una la caída con una aceleración igual a la de la gravedad. Luego, utilizamos las fórmulas de cinemática que encontramos en la entrada anterior para calcular velocidad y posición al final de ese primer intervalo. Y ahora damos el «salto»: usamos el valor de la velocidad al final del primer intervalo para calcular la aceleración que utilizaremos en el segundo, aceleración con la cual calcularemos una nueva velocidad que será utilizada para calcular la aceleración del tercer intervalo y…

Es decir, en el intervalo 1 tenemos

a_1 \approx g;

v_1 \approx g \Delta t;

x_1 \approx \frac{1}{2} g \Delta t ^2.

Para el segundo intervalo, tomaremos por lo tanto los valores finales del primer intervalo como valores iniciales de este segundo, calculando una nueva aceleración:

\displaystyle a_2 \approx g -\frac{k}{M} v_1 ^2;

v_2 \approx v_1 + a_2 \Delta t;

x_2 \approx x_1 + \frac{1}{2} a_2 \Delta t ^2.

Y en general, para el paso n-ésimo

\displaystyle a_n \approx g -\frac{k}{M} v_{n-1} ^2;

v_n \approx v_{n-1} + a_n \Delta t;

x_n \approx x_{n-1} + \frac{1}{2} a_n \Delta t ^2.

Es decir, tenemos una fórmula recursiva que nos permite calcular una aproximación para un momento determinado en función del valor del momento anterior.

A este punto es importante remarcar qué estamos haciendo realmente en este algoritmo: en cada uno de los pequeños intervalos temporales en los que estamos dividiendo el problema, la aproximación realizada no es otra cosa que reemplazar la ecuación diferencial original por otra más simple, donde la aceleración es una constante. Resuelta esta ecuación simple en uno de los intervalos, usamos los resultados como valores iniciales del intervalo siguiente.

Antes de mostrar el guión para Octave, veamos un poco lo que significa realizar la aproximación de este modo y cuáles son sus «riesgos». Por ejemplo, si tomamos \Delta t = 2 s , obtendremos lo siguiente:

caer-2s

Como referencia, la línea roja muestra la solución exacta de la ecuación diferencial, la cual está dada por la función

\displaystyle v(t) = \sqrt{\frac{M g}{k}}\cdot \tanh{\left(\sqrt{\frac{k g}{M}}\cdot t \right)}

Como podemos ver, tomando dos segundos la velocidad dada por nuestro simple modelo numérico crece mucho más rápido de lo que debería, llegando a «pasarse de largo» de la velocidad límite: en lugar de aproximarse lentamente a su valor final nuestra solución numérica comienza a oscilar violentamente.

Claramente, dos segundos no es un intervalo temporal «suficientemente pequeño»…

Si en lugar de dos segundos tomamos \Delta t = 1 s , la situación mejora notablemente… aunque no del todo:

caer-1s

Seguimos teniendo un crecimiento exagerado al principio, pero ya casi no oscila y llega al valor correcto rápidamente. En cambio, con \Delta t = 0,1 s

caer-01s

la solución numérica es ya casi indistinguible de la «real», por lo que finalmente podemos decir que una décima de segundo es un intervalo «suficientemente pequeño».

Como era de esperarse, tomar un valor para \Delta t más pequeño nos da mejores resultados. El problema, claro está, es que más pequeño es el intervalo mayor será el número de pasos necesarios por lo que se requiere cada vez más poder de cálculo. Ciertamente esto no es un inconveniente para el problema aquí propuesto: en mi procesador i3 tomando \Delta t = 0,0001 s octave no tarda más de quince segundos en realizar los más de 400000 pasos necesarios, pero en problemas más complejos como el cálculo de la órbita de una sonda interplanetaria, o las órbitas halo en torno al punto de Lagrange Tierra-Sol que usan los satélites que controlan la actividad solar, o realizar un pronóstico del tiempo, o calcular las energías disponibles para los electrones en los semiconductores que hacen funcionar nuestras computadoras, o encontrar las estructuras de las moléculas utilizadas en la industria farmacéutica… se requiere el uso de computadoras que difícilmente podrían caber en un escritorio.

Además, si bien en este caso de la caída libre es posible encontrar una solución exacta de la ecuación diferencial en otras situaciones la única posibilidad de obtener un resultado son las aproximaciones numéricas: los problemas que tienen una solución exacta son la excepción, no la regla. Pero lo importante es que para todos estos problemas existe una herramienta matemática que resulta fundamental en el desarrollo de un algoritmo que permita encontrar, aunque más no sea aproximadamente, la solución buscada: el cálculo diferencial.

Pues bien, hasta aquí hemos llegado. Antes de dejarte con el guión de Octave, estimado lector, solo me queda por decir que espero que estas seis entradas hayan sido de tu agrado.


Código a colocar en un archivo de extensión .m para ser llamado desde la línea de comando de Octave. El valor de \Delta t se maneja con el parámetro dt. Todo lo que está a la derecha de un % es un comentario que no será ejecutado.

M = 90; % la masa
g = 9.8; % aceleración de la gravedad
k = 0.29; % el coeficiente del término de rozamiento
h1 = 3000; % altura del salto
h2 = 1000; % altura al abrir el paracaídas
fin1 = h1-h2; % recorrido hecho hasta abrir el paracaídas

dt = 0.1; % intervalo de tiempo en segundos

pos = 0; % valor inicial de posición, velocidad aceleración y tiempo
vel = 0;
acc = g;
tem = 0;

j = 1; % un índice

x(1) = pos; % todo irá en vectores
v(1) = vel;
a(1) = acc;
t(1) = tem;

while pos < fin1 % mientras no llegamos al paracaídas, seguimos con el cálculo
j = j+1;
vel = vel + acc*dt;
v(j) = vel;
acc = g - (k/M)*(v(j-1))^2;
pos = pos + vel*dt + (acc*dt^2)/2;
x(j) = pos;
tem = tem + dt;
t(j) = tem;
end

vcorr = sqrt(M*g/k)*tanh(sqrt(k*g/M)*(t)); % la velocidad correcta
plot(t,v,t,vcorr,'r');xlabel('tiempo (s)');ylabel('velocidad (m/s)')
figure;plot(x,v,x,vcorr,'r');xlabel('posición (m)');ylabel('velocidad (m/s)')
Anuncios

, ,

A %d blogueros les gusta esto: