Taller de Programación Racket para Matemáticas

[Seminario de matemáticas, CEMATi]

Invitación

Programación Racket para Matemáticas

[viernes 19 y 26/sept. 2014. Unidad Otay: 5-6 pm, Lab. "A distancia"]

A todos los interesados (maestros, estudiantes, visitantes) se les invita a este Seminario de matemáticas. En esta ocasión presentaremos el ambiente de programación DrRacket. El lenguaje Racket (antes PLT/Scheme) viene de la tradición LISP/Scheme. Se cubrirán los elementos esenciales del lenguaje, el módulo Plot, así como el nuevo módulo Math de funciones matemáticas. Se dará cierto énfasis a gráficación en \mathbf{R}^2 y en \mathbf{R}^3.riegraf-0-300x300

Cupo limitado. Para info y registro favor enviar email a: investiga@cemati.org.

Atte. Facilitador, E. Cómer (ITT, Cemati.org)

Figuras en LyX con paquete TikZ

[rev. 2012.11.26]

Rendered by QuickLaTeX.com

Ejemplo (ver [1]) de figura generada con LyX/LaTeX, utilizando el paquete TikZ. (Requisito: Instalar en el preámbulo del documento los paquetes:

  • \usepackage[usenames,dvipsnames,pdftex]{xcolor}
  • \usepackage{tikz,ifthen}

En la sección del documento LyX que desee incluir la figura, puede abrir un bloque de código LaTeX y escribir el código siguiente:

Como ejemplo puede revisar el archivo LyX: tikz1.lyx. Para mayor información se le invita a consultar el manual de Till Tantau (2007) TikZ and PGF (PDF, 405 pp.)

Mucho éxito en sus dibujos en LyX/LaTeX/TikZ

Resortes acoplados y solución en WxMaxima

El siguiente sistema de ecuaciones representa el comportamiento de dos resortes acoplados verticalmente (ver pp. 296-297 [Zill 2009]):

     \begin{align*} \ddot{x_1}+10x_1-4x_2 &= 0 \\ -4x_1+\ddot{x_2}+4x_2 &= 0 \end{align*}

con condiciones iniciales: x_1(0)=0,\; x_1'(0)=1,\; x_2(0)=0,\; x_2'(0)=-1

Como ejemplo del uso de WxMaxima para resolver este tipo de problemas, presentamos  las instrucciones básicas para su solución (ver archivo resuelveode2.wxm)

Evaluando estas instrucciones contenidas en el archivo anterior (con Ctrl-R, en WxMaxima), obtenemos la solución:

     \begin{align*} x_1(t) &= \frac{\sqrt{3}}{5}sen(2 \sqrt{3} t) - \frac{1}{5 \sqrt{2}}sen(\sqrt{2} t) \\ x_2(t) &= -\frac{\sqrt{3}}{10}sen(2 \sqrt{3} t) - \frac{\sqrt{2}}{5}sen(\sqrt{2} t)\end{align*}

Compartimos el archivo visualizaode2.ggb (para GeoGebra) que nos permite visualizar la solución previamente obtenida.

Gráficas de posiciones de resortes acoplados

Le invitamos también a activar el archivo animaode2.ggb que permite de manera básica simular la dinámica de este sistema:

Resortes acoplados y animación básica

Referencia: [Zill 2009] Dennis G. Zill. Ecuaciones Diferenciales con aplicaciones de modelado (9/e) Cengage Learning.


Applet recomendado: Masses & Springs (con audio) PhET Interactive simulations. Univ. of Colorado at Boulder. [acc. 2010.12.11]

Método de Runge-Kutta orden 4 para sistemas y Ecuaciones de Lorenz

El método de Runge-Kutta orden 4 para sistemas es una extensión del método numérico usado para resolver una ecuación de la forma \frac{dy}{dx}=f(x,y) al caso de un sistema de n de tales ecuaciones, es decir permite resolver:

     \begin{align*} \frac{du_1}{dt} &= f_1(t,u_1,\cdots,u_n) \\ \frac{du_2}{dt} &= f_2(t,u_1,\cdots,u_n) \\  \vdots \quad  &= \quad \vdots \\ \frac{du_n}{dt}  &= f_n(t,u_1,\cdots,u_n) \end{align*}

Si \mathbf{u}=(u_1,u_2,\cdots,u_n), entonces el sistema anterior se puede representar vectorialmente como \mathbf{u}' = \mathbf{f}(t,\mathbf{u}), de forma que si se aproxima \mathbf{u}(t_j) con \mathbf{w}_j  (con t_0 \leq t_j \leq t_n),  entonces podemos escribir el método como:

 \mathbf{w}_{j+1} = \mathbf{w}_j+\frac{h}{6}(\mathbf{k}_1+2 \mathbf{k}_2 + 2 \mathbf{k}_3 + \mathbf{k}_4)

donde

  •  \mathbf{k}_1 = \mathbf{f}(t_j, \mathbf{w}_j)
  •  \mathbf{k}_2 = \mathbf{f}(t_j+\frac{h}{2}, \mathbf{w}_j+\frac{h}{2} \mathbf{k}_1)
  •  \mathbf{k}_3 = \mathbf{f}(t_j+\frac{h}{2}, \mathbf{w}_j+\frac{h}{2} \mathbf{k}_2)
  •  \mathbf{k}_4 = \mathbf{f}(t_j+h, \mathbf{w}_j+h \mathbf{k}_3)

para t_j = t_0 + j h y una condición inicial \mathbf{w}(t_0)=\mathbf{w}_0. En el archivo ordinarias.rkt, se implementa en Scheme/DrRacket (versión 5.3)  este método.

Ejemplo: veamos el caso de las ecuaciones clásicas de Lorenz (ver atractor de Lorenz):

     \begin{align*} \frac{dx}{dt} &= \sigma (y-x) \\ \frac{dy}{dt} &= \rho x - y -xz \\  \frac{dz}{dt} &=xy-\beta z \end{align*}

Estas ecuaciones con los valores \sigma =10, \beta = \frac{8}{3} y \rho = 28, generan la trayectoria de solución indicada en la figura, cuando se inicia en el punto (10,7,7).

Le invitamos a explorar el código en lorenz.rkt (y archivos auxiliares: ordinarias.rkt, derivaciones.rkt, vectoriales.rkt, basicas.rkt, para Scheme/DrRacket [versión 5.3]). Las instrucciones básicas utilizadas son:

(define sigma 10)
(define beta 8/3)
(define rho 28)

(define fs (list
            (lambda(t x y z)
              (* 10 (- y x)))
            (lambda(t x y z)
              (- (* 28 x) y (* x z)))
            (lambda(t x y z)
              (- (* x y) (* 8/3 z)))))
(define sols
  (srk4-h fs 0 30 '(10 7 7) 0.005))

Note que (en este caso) se incluye en la lista de puntos solución (j,t_j,x(t_j),y(t_j),z(t_j)) donde t_j=t_0 + j h para h=0.005

Observación: si se desea resolver una ecuación diferencial de orden superior, es suficiente utilizar cambios de variable y convertir dicha ecuación a un sistema equivalente de ecuaciones diferenciales de primer orden. Por ejemplo la ecuación

\frac{d^2x}{dt^2}=-\frac{k}{m}x

asociada a un resorte, se puede convertir a un sistema de dos ED de primer orden, introduciendo la variable u=\frac{dx}{dt}, con lo cual \frac{du}{dt}=\frac{d^2x}{dt^2} y entonces el sistema equivalente es:

     \begin{align*} \frac{du}{dt} &= -\frac{k}{m}x \\ \frac{dx}{dt} &= u \end{align*}

Nota: se corrigió en mordinarias.rkt la función srk4-h (y en consecuencia la gráfica de Lorenz), favor de actualizar a versión 0.4; gracias y una disculpa [2011.12.23]

Quiz de transformaciones lineales (nivel 1)

Quiz AL5 nivel 1

Se le invita a ejercitarse en este importante tema, resolviendo los siguientes problemas:
Start
Return
Shaded items are complete.
123
Return

{Hecho con mTouch Quiz}

Transformaciones lineales y GeoGebra

    \[  \]

Una transformación lineal T: V \rightarrow W del espacio vectorial V con operaciones (+,\cdot) al espacio vectorial W con operaciones (\oplus, \otimes ), es una función que cumple para todo \mathbf{u},\;\mathbf{v} \in V, las siguientes dos propiedades:

  1. T(\mathbf{u}+ \mathbf{v})=T(\mathbf{u}) \oplus T(\mathbf{v})
  2. T(\alpha \cdot \mathbf{v})=\alpha \otimes T(\mathbf{v})

Observación: cuando las operaciones son estándares, utilizamos por convención la notación sencilla de + y \cdot\;\; para la suma de vectores y multiplicación escalar por vector respectivamente, en otros casos es necesario definir claramente las cuatro operaciones implicadas.

Ejemplo: Consideremos a T: \mathbb{R}^2 \rightarrow \mathbb{R}^2 definida por:

    \[ T(x,y)=(x+ay,y) \]

para un valor arbitrario a \in \mathbb{R}. Veamos si cumple las dos condiciones mencionadas:

  1.     \begin{align*} T((x_1,y_1)+(x_2,y_2)) &= T(x_1+x_2,y_1+y_2)  \\ &= (x_1+x_2 + a(y_1+y_2),y_1+y_2) \\ &=(x_1+a y_1 + x_2 + a y_2, y_1 + y_2) \\ &=( x_1+a y_1, y_1) + (x_2 + a y_2, y_2) \\ &= T(x_1,y_1 ) + T(x_2,y_2) \end{align*}

  2.     \begin{align*} T(\alpha (x,y)) &=T(\alpha x, \alpha y) \\ &= (\alpha x + a \alpha y, \alpha y) \\ &= \alpha (x + a y, y) \\ &= \alpha T(x,y) \end{align*}

Ya que cumple ambas propiedades, la transformación T anterior es lineal.

Una forma equivalente de aplicar una transformación T es mediante su matriz de trasformación A_T, que en este caso es: \begin{pmatrix} 1 & a \\ 0 & 1 \end{pmatrix}, es decir:

    \[ \begin{pmatrix} x+ay \\ y \end{pmatrix} = \begin{pmatrix} 1 & a \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]

Compartimos el archivo shear1.ggb (de GeoGebra) que ilustra la transformación lineal “shear” (corte) a lo largo del eje-x.

Ejemplo de “shear” en el eje-x

Para el caso de \mathbb{R}^2, se le invita a explorar en GeoGebra las diferentes transformaciones disponibles en su versión 4.0 : desliza, refleja, rota, corte [homotecia], estira y traslada (en ingl. dilate, reflect, rotate, shear, stretch, translate). Felicitamos a los desarrolladores de GeoGebra, por tan excelente producto.

Solución de un sistema de ED mediante Laplace

Consideremos el PVI definido por:

    \[ \frac{dx}{dt}=1-y, \quad \frac{dy}{dt}=x \]

y las condiciones iniciales x(0)=1, y(0)=0.

Aplicando transformadas de Laplace obtenemos:

  1. \mathcal{L}\{\frac{dx}{dt}=1-y\}\rightarrow sX(s)-1=\frac{1}{s}-Y(s)
  2. \mathcal{L}\{\frac{dy}{dt}=x \} \rightarrow sY(s)=X(s), por lo cual
  3. Y(s)=\frac{1}{s}X(s), que sustituimos en el paso 1 para obtener:
  4. sX(s)-1=\frac{1}{s}-\frac{1}{s}X(s), lo cual factorizamos como X(s)(s+\frac{1}{s})=\frac{1}{s}+1, que se simplifica a X(s)(\frac{s^2+1}{s})=\frac{1+s}{s}, con lo cual: X(s)=\frac{s+1}{s^2+1}=\frac{s}{s^2+1}+\frac{1}{s^2+1}
  5. Podemos ahora aplicar la transformada de Laplace inversa para obtener: x(t)=cost+sent (usando las fórmulas L#7 y L#8)
  6. Por otra parte, sustituyendo X(s) en el paso 2, obtenemos: Y(s)=\frac{1}{s}[\frac{s+1}{s^2+1}]=\frac{1}{s^2+1}+\frac{1}{s(s^2+1)}
  7. Aplicando transformada inversa a la ecuación anterior, obtenemos y(t)=1-cost+sent (usando las fórmulas L#7 y L#30)

Por tanto la solución al PVI original es:

    \[ x(t)=cost+sent, \;y(t)=1-cost+sent \]

Compartimos el archivo sistemaed1.rkt (de Scheme/DrRacket) para visualizar el sistema y su solución:

Campo direccional y curva solución

Como referencia, les invitamos a explorar también el archivo laplacesisej1.wxm (de WxMaxima) que resuelve el mismo sistema.

Método de Runge-Kutta de orden 4

El método de Runge-Kutta de orden 4, para resolver numéricamente una ecuación diferencial de la forma: \displaystyle \frac{dy}{dx}=f(x,y),  (parte de la familia de métodos de Runge-Kutta) nos permite de forma sencilla obtener la sucesión \displaystyle \{ x_i, y_i \}_{0}^{n} dada la condición inicial y(x_0)=y_0, mediante la fórmula:

  • y_{k+1}=y_i+\frac{h}{6}(k_1+2k_2+2k_3+k_4)
  • donde:
  • k_1=f(x_i,y_i)
  • k_2=f(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1 h)
  • k_3=f(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2 h)
  • k_4=f(x_i+h,y_i+k_3 h)

Veamos un ejemplo de su uso resolviendo un circuito eléctrico RC modelado por  R\frac{dq}{dt}+\frac{1}{C}q=E(t) (y resuelto previamente por transformadas de Laplace). Consideraremos el caso cuando: q(0)=0, R=2.5 \Omega, C=0.08\;f, pero ahora definiremos el voltaje como E(t)= (5+sen4t)\;\mathcal{U}(t-2). El siguiente código en Scheme especifica el problema y permite generar la solución numérica.

;; función auxiliar y datos
(define (heaviside t)
  (if (< t 0) 0 1))
(define R 2.5)
(define C 0.08)

;; definicion de E(t)
(define (fem t)
  (* (+ 5 (sin (* 4 t)))
     (heaviside (- t 2))))

;; lado derecho de la ecuacion dq/dt = f(t,q)
(define (rhs t q)
  (/ (- (fem t) (/ q C))
     R))

;; solucion de la ED usando Runge-Kutta orden 4
(define sols
  (rk4-n rhs 0 6. 0 400))

Una vez corrido el programa, podemos solicitar las soluciones para obtener:

>sols
'((0 0 0)
  (1 0.015 0)
  (2 0.03 0)
 ...
  (397 5.955 0.3428757)
  (398 5.97 0.3414624)
  (399 5.985 0.3402598)
  (400 6.0 0.3392721))

En este caso hemos incluido el índice en cada pareja x_i,\;y_i calculada. Compartimos el archivo de aplicación circuitorc1.rkt (y los archivos auxiliares basicas.rkt, vectoriales.rkt,  derivaciones.rkt y ordinarias.rkt). [Nota: estos archivos estan actualizados para DrRacket 5.3]

Carga y corriente en el capacitor, dado E(t)

Se le invita a comparar resultados y métodos con el segmento Ejemplo de PVI con traslación en el tiempo.

Referencia recomendada: Persson, Per-Olof. Runge-Kutta Methods, Department of Mathematics, UC Berkeley [PDF, 8 pp., acc. 2011.12.03].

Nota:  El código utilizado no está optimizado. Se recomienda llamar  rk4-n  con un menor número de subintervalos, por ejemplo (rk4-n rhs 0 6 0 100)

Ejemplo de PVI con traslación en el tiempo

Un circuito eléctrico sencillo RC se puede modelar con la ecuación diferencial: R\frac{dq}{dt}+\frac{1}{C}q=E(t). Resolvamos el PVI siguiente, mediante transformadas de Laplace, el caso cuando: q(0)=0, R=2.5 \Omega, C=0.08\;f, y E(t)=5 \mathcal{U}(t-3). Observe que el voltaje tiene una traslación en el tiempo.

  • Con L#54 para n=2, obtenemos:
  • \mathcal{L}\{ 2.5 \frac{dq}{dt}+\frac{1}{0.08}q\}=2.5[sQ(s)-q(0)]+\frac{100}{8}Q(s), que factorizando se simplifica a Q(s)[\frac{5}{2}(s+5)].
  • Por otra parte \mathcal{L}\{ 5\mathcal{U}(t-3) \}=\frac{5}{s}e^{-3s} (por la fórmula L#51)
  • Igualando ambos resultados obtenemos:
  • Q(s)[\frac{5}{2}(s+5)]=\frac{5}{s}e^{-3s}, de donde obtenemos (con un paso de simplificación) que
  • Q(s)=\displaystyle \frac{2e^{-3s}}{s(s+5)}.
  • Ahora para obtener q(t) aplicamos la transformada inversa (usando la fórmula L#52)
  • q(t)=\mathcal{L}^{-1}\{ \displaystyle \frac{2e^{-3s}}{s(s+5)} \}=2f(t-3)\mathcal{U}(t-3), con f(t-3)=\mathcal{L}^{-1}\{ \frac{1}{s(s+5)}\} |_{t \rightarrow t-3} que (por la fórmula L#28) es igual a \displaystyle \frac{e^0-e^{-5t}}{0+5} |_{t \rightarrow t-3}, es decir f(t-3)=\frac{1}{5}-\frac{1}{5}e^{-5(t-3)}
  • Y por tanto, nuestra solución es: q(t)=2[\frac{1}{5}-\frac{1}{5}e^{-5(t-3)}] \mathcal{U}(t-3)
  • o bien q(t)=\frac{2}{5}\mathcal{U}(t-3)-\frac{2}{5}e^{-5(t-3)} \mathcal{U}(t-3).

Compartimos el archivo traslaciont1.ggb (de GeoGebra) para visualizar el problema, conforme se indica en la siguiente gráfica, generada por el archivo anterior:

Visualización de carga y corriente resultante

Nota: Las referencias L#<num> corresponden a los números de las fórmulas en la tabla de transformadas de Laplace en el libro de Zill 9/e.

Problema Ci-004

Problema Ci-004. [CSP] Sea f(x)=\displaystyle \frac {\alpha}{1+x^4} cos^2(x^2). Determinar (con al menos cuatro cifras decimales) el valor de \alpha para que f(x) cumpla con ser una función de densidad de probabilidad, en el intervalo [0,\infty).

Se le invita a enviar sus propuestas de solución mediante comentarios a este post.

To see the Aside click here.To hide the Aside click here.

[2011.11.29] Sugerencia: utilizar métodos numéricos

Mucho éxito en sus exploraciones matemáticas

P.S. Próximamente incorporaremos los problemas Ci-001 a Ci-002 (con sus respectivas soluciones) antes disponibles en el sitio cemati.com [actualmente en transición]