3.2: Técnicas numéricas

Las carreras de ingeniería son la mayoría de los estudiantes en el tipo de curso de física para el que está diseñado este libro, por lo que lo más probable es que caigas en esa categoría. Aunque seguramente reconoces que la física es una parte importante de tu formación, si has tenido alguna exposición a cómo funcionan realmente los ingenieros, probablemente seas escéptico sobre el sabor de la resolución de problemas que se imparte en la mayoría de los cursos de ciencias. Te das cuenta de que no muchos cálculos prácticos de ingeniería caen dentro de la estrecha gama de problemas para los que se puede calcular una solución exacta con una hoja de papel y un lápiz afilado. Los problemas de la vida real suelen ser complicados y, por lo general, deben resolverse mediante el procesamiento de números en una computadora, aunque a menudo podemos obtener información trabajando aproximaciones simples que tienen soluciones algebraicas. No solo es más útil la resolución de problemas numéricos en la vida real, también es educativa; como estudiante principiante de física, realmente solo sentí que entendía el movimiento del proyectil después de haberlo trabajado en ambos sentidos, usando álgebra y luego un programa de computadora. (Esto fue en los días en que 64 kilobytes de memoria se consideraban mucho).

En esta sección, comenzaremos viendo cómo aplicar técnicas numéricas a algunos problemas simples para los que conocemos la respuesta en “forma cerrada”, es decir, una sola expresión algebraica sin ningún cálculo ni sumas infinitas. Después de eso, resolveremos un problema que te habría hecho mundialmente famoso si pudieras haberlo hecho en el siglo XVII ¡usando papel y una pluma! Antes de continuar, deberías leer el Apéndice 1 en la página 908 que te introduce en el lenguaje de programación Python.

Primero resolvamos el trivial problema de encontrar cuánto tiempo le toma a un objeto moverse a velocidad v recorrer una distancia en línea recta dist. Esta respuesta de forma cerrada es, por supuesto, dist/v, pero el punto es introducir las técnicas que podemos utilizar para resolver otros problemas de este tipo. La idea básica es dividir la distancia en n partes iguales, y sumar los tiempos requeridos para recorrer todas las partes. La siguiente función Python hace el trabajo. Tenga en cuenta que no debe escribir los números de línea de la izquierda, y tampoco es necesario escribir los comentarios. He omitido las indicaciones >>> y... con el fin de ahorrar espacio.

import math
def time1(dist,v,n):
  x=0 # Initialize the position.
  dx = dist/n # Divide dist into n equal parts.
  t=0 # Initialize the time.
  for i in range(n):
    x = x+dx # Change x.
    dt=dx/v # time=distance/speed
    t=t+dt # Keep track of elapsed time.
  return t

¿Cuánto tiempo se tarda en mover 1 metro a una velocidad constante de 1 m/s? Si hacemos esto,

>>> print(time1(1.0,1.0,10)) # dist, v, n
0.99999999999999989

Python produce la respuesta esperada dividiendo la distancia en diez segmentos iguales de 0.1 metros, y sumando los diez tiempos de 0.1 segundos requeridos para atravesar cada uno. Dado que el objeto se mueve a velocidad constante, ni siquiera importa si establecemos n en 10, 1 o un millón:

>>> print(time1(1.0,1.0,1)) # dist, v, n
1.0

Ahora hagamos un ejemplo donde la respuesta no sea obvia para las personas que no conocen el cálculo: ¿cuánto tiempo tarda un objeto en caer por una altura h, empezando por el descanso? Sabemos por el ejemplo 8 de la página 83 que la respuesta exacta, encontrada usando cálculo, es\(\sqrt{2h/g}\). Veamos si podemos reproducir esa respuesta numéricamente. La principal diferencia entre este programa y el anterior es que ahora la velocidad no es constante, por lo que necesitamos actualizarlo a medida que avanzamos. La conservación de la energía da\(mgh=(1/2)mv^2+mgy\) por la velocidad\(v\) a la altura\(y\), así\(v=-\sqrt{2g(h-y)}\). (Elegimos la raíz negativa porque el objeto se mueve hacia abajo, y nuestro sistema de coordenadas tiene el\(y\) eje positivo apuntando hacia arriba).

import math
def time2(h,n):
  g=9.8 # gravitational field
  y=h # Initialize the height.
  v=0 # Initialize the velocity.
  dy = -h/n # Divide h into n equal parts.
  t=0 # Initialize the time.
  for i in range(n):
    y = y+dy # Change y. (Note dy<0.)
    v = -math.sqrt(2*g*(h-y)) # from cons. of energy
    dt=dy/v # dy and v are <0, so dt is >0
    t=t+dt # Keep track of elapsed time.
  return t

Para h =1.0 m, el resultado de forma cerrada es\(\sqrt{2\cdot1.0\ \text{m}/9.8\ \text{m}/\text{s}^2}=0.45\ \text{s}\). Con la caída dividida en solo 10 intervalos de altura iguales, la técnica numérica proporciona una aproximación bastante pésima:

>>> print(time2(1.0,10)) # h, n
0.35864270709233342

Pero al aumentar n a diez mil, obtenemos una respuesta que es lo más cercana que necesitamos, dada la precisión limitada de los datos brutos:

>>> print(time2(1.0,10000)) # h, n
0.44846664060793945

Un punto sutil es que cambiamos y en la línea 9, y luego en la línea 10 calculamos v, que depende de y. Como y solo está cambiando en una diecésima milésima de metro con cada paso, podrías pensar que esto no haría mucha diferencia, y estarías casi en lo cierto, excepto por un pequeño problema: si intercambiamos las líneas 9 y 10, entonces la primera vez a través del bucle, tendríamos v =0, que sería producir un error de división por cero cuando calculamos dt! En realidad, lo que tendría más sentido sería calcular la velocidad a la altura y y la velocidad a la altura y+dy (recordando que dy es negativo), promediarlos juntos, y usar ese valor de y para calcular la mejor estimación de la velocidad entre esos dos puntos. Dado que la aceleración es constante en el presente ejemplo, esta modificación da como resultado un programa que da un resultado exacto incluso para n =1:

import math
def time3(h,n):
  g=9.8
  y=h
  v=0
  dy = -h/n
  t=0
  for i in range(n):
    y_old = y
    y = y+dy
    v_old = math.sqrt(2*g*(h-y_old))
    v = math.sqrt(2*g*(h-y))
    v_avg = -(v_old+v)/2.
    dt=dy/v_avg
    t=t+dt
  return t

>>> print(time3(1.0,1)) # h, n
0.45175395145262565

a/Aproximaciones a la curva de brachistocrón usando un polinomio de tercer orden (línea continua) y un polinomio de séptimo orden (discontinuo). Este último sólo mejora el tiempo en cuatro milisegundos.

Ahora estamos listos para atacar un problema que desafiaba a las mejores mentes de Europa allá por los días en que no había computadoras. En 1696, el matemático Johann Bernoulli planteó la siguiente famosa pregunta. A partir del reposo, un objeto se desliza sin fricción sobre una curva que une el punto\((a,b)\) con el punto\((0,0)\). De todas las formas posibles que tal curva podría tener, ¿cuál lleva el objeto a su destino en el menor tiempo posible y cuánto tiempo lleva? La curva óptima se llama brachistochrone, del griego “corto tiempo”. La solución al problema de la brachistocrona evadió al propio Bernoulli, así como a Leibniz, quien había sido uno de los inventores del cálculo. El físico inglés Isaac Newton, sin embargo, se quedó hasta tarde una noche después de un día de trabajo manejando la menta real, y, según la leyenda, produjo una solución algebraica a las cuatro de la mañana. Luego lo publicó de forma anónima, pero se dice que Bernoulli comentó que cuando lo leyó, supo instantáneamente por el estilo que era Newton —podía “distinguir al león por la marca de su garra”.

En lugar de intentar una solución algebraica exacta, como lo hizo Newton, produciremos un resultado numérico para la forma de la curva y el tiempo mínimo, en el caso especial de\(a\) =1.0 m y\(b\) =1.0 m. Intuitivamente, queremos comenzar con una caída bastante pronunciada, ya que cualquier velocidad que podamos acumular al inicio nos ayudará durante el resto de la moción. Por otro lado, es posible ir demasiado lejos con esta idea: si bajamos recto por toda la distancia vertical, y luego hacemos un giro en ángulo recto para cubrir la distancia horizontal, el tiempo resultante de 0.68 s es bastante más largo que el resultado óptimo, la razón es que el camino es innecesariamente largo . Hay infinitamente muchas curvas posibles para las que podríamos calcular el tiempo, pero veamos polinomios de tercer orden,

\[\begin{equation*} y = c_1x+c_2x^2+c_3x^3 , \end{equation*}\]

donde\(c_3=(b-c_1a-c_2a^2)/a^3\) requerimos para que la curva pase por el punto\((a,b)\). El programa Python, a continuación, no es muy diferente de lo que hemos hecho antes. La función sólo pide\(c_1\) y\(c_2\), y calcula\(c_3\) internamente en la línea 4. Dado que el movimiento es bidimensional, tenemos que calcular la distancia entre un punto y el siguiente usando el teorema de Pitágoras, en la línea 16.

import math
def timeb(a,b,c1,c2,n):
  g=9.8
  c3 = (b-c1*a-c2*a**2)/(a**3)
  x=a
  y=b
  dx = -a/n
  t=0
  for i in range(n):
    y_old = y
    x = x+dx
    y = c1*x+c2*x**2+c3*x**3
    dy = y-y_old
    v_old = math.sqrt(2*g*(b-y_old))
    v = math.sqrt(2*g*(b-y))
    v_avg = (v_old+v)/2.
    ds = math.sqrt(dx**2+dy**2) # Pythagorean thm.
    dt=ds/v_avg
    t=t+dt
  return t

Como primera suposición, podríamos intentar una línea diagonal recta,\(y=x\), que corresponde a establecer\(c_1=1\), y todos los demás coeficientes a cero. El resultado es un tiempo bastante largo:

>>> a=1.
>>> b=1.
>>> n=10000
>>> c1=1.
>>> c2=0.
>>> print(timeb(a,b,c1,c2,n))
0.63887656499994161

Lo que realmente necesitamos es una curva que sea muy empinada a la derecha, y más plana a la izquierda, por lo que en realidad tendría más sentido intentar\(y=x^3\):

>>> c1=0.
>>> c2=0.
>>> print(timeb(a,b,c1,c2,n))
0.59458339947087069

Esta es una mejora significativa, ¡y resulta ser solo una centésima de segundo del menor tiempo posible! Es posible, aunque no muy educativo ni entretenido, encontrar mejores aproximaciones a la curva de la brachistocrona jugueteando con los coeficientes del polinomio a mano. El verdadero punto de esta discusión fue dar un ejemplo de un problema no trivial que puede ser atacado exitosamente con técnicas numéricas. Encontré la primera aproximación que se muestra en la figura a,

\[\begin{equation*} y = (0.62)x+(-0.93)x^2+(1.31)x^3 \end{equation*}\]

utilizando el programa que figura en el apéndice 2 de la página 911 para buscar automáticamente la curva óptima. La aproximación de séptimo orden que se muestra en la figura provino de una extensión directa del mismo programa.

Colaboradores y Atribuciones

Template:ContribCrowell