Programación Racket para Matemáticas

[Taller de herramientas matemáticas, CEMATi [2]]

Invitación

Programación Racket para Matemáticas

[viernes 24/oct [1]. 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]

Problema Ci-003

Problema Ci-003 [CSP]. Calcular el área sobre el eje-x acotada superiormente por la polilínea formada por los puntos A_0, A_1, A_2,\dots, A_n,\dots; donde A_0=(0,0) y A_n=(x_n,y_n) para n \ge 1,  con x_n=\alpha +\alpha^2+\alpha^3+\cdots +\alpha^n, y_n=\beta +\beta^2+\beta^3+\cdots +\beta^n, para los valores \alpha =\frac{1}{3} y \beta =2.

Solución: [Presentada el 28 de julio, 2011, al problema originalmente publicado el 27 de julio, 2011] por Damian Moreno [Estudiante, Los Alamos, Chile] (Felicitaciones ver Aside para solución)

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

[2011.08.03] Solución por Damian Moreno (2011.07.28). Ya que tanto x_n como y_n corresponden a series geométricas, definimos x_n=\frac{1}{3}\frac{1-(1/3)^n}{1-\frac{1}{3}}=\frac{1}{2}(1-(\frac{1}{3})^n), así mismo, y_n=2\frac{1-2^n}{1-2}=(-2)(1-2^n)=2(2^n-1). Con lo cual podemos definir \Delta x_n=x_n-x_{n-1}=\alpha^n=(\frac{1}{3})^n, y \Delta y_n=y_n-y_{n-1}=\beta^n=2^n. Si consideramos A_n como el área en el intervalo (x_{n-1},x_n), entonces el área pedida será A=\Sigma_{n=1}^\infty A_n. Ya que A_n=\frac{1}{2}\Delta x_n \cdot \Delta y_n+\Delta x_n \cdot y_{n-1}, entonces podemos escribir A_n=\frac{1}{2}(\frac{1}{3})^n 2^n +(\frac{1}{3})^n\cdot 2 (2^{n-1}-1), que una vez simplificado queda como A_n=(\frac{2}{3})^{n-1}-2(\frac{1}{3})^n. Ahora, para encontrar A aplicamos sencillamente sumatoria, desde 1 hasta \infty [a las series geométricas implicadas] para obtener A=\frac{1}{1-(2/3)}-2\frac{1/3}{1-(1/3)}=3-2\frac{1}{2}=3-1=2\;\star

Mucho éxito en sus reflexiones matemáticas

Cambio de base en la representación de un vector

Supongamos que tenemos el vector \mathbf{u}_w=\binom{2}{2} representado en la base B_1=\{ \mathbf{w}_1,\mathbf{w}_2\} donde \mathbf{w}_1 = \binom{1}{0} y \mathbf{w}_2 = \binom{0}{1}. Para representar dicho vector \mathbf{u}_w en una nueva base, por ejemplo B_2=\{ \mathbf{v}_1,\mathbf{v}_2\} donde \mathbf{v}_1 = \binom{-1}{-1} y \mathbf{v}_2 = \binom{-1}{1}, lo que hacemos es crear una matriz de transición T cuyas columnas corresponden a la representación de los vectores de la base B_1 en términos de la nueva base B_2. Es decir, determinamos \alpha_1, \beta_1, \alpha_2, \beta_2 tal que:

\mathbf{w}_1 = \alpha_1 \mathbf{v}_1 + \beta_1 \mathbf{v}_2

\mathbf{w}_2 = \alpha_2 \mathbf{v}_1 + \beta_2 \mathbf{v}_2

En este caso (resolviendo los sistemas lineales resultantes) obtenemos: \alpha_1 = -\frac{1}{2}, \beta_1 = -\frac{1}{2}, \alpha_2 = -\frac{1}{2} y \beta_2 = \frac{1}{2}, es decir: T=\binom{\alpha_1 \quad \alpha_2}{\beta_1 \quad \beta_2}=\binom{-0.5 \quad -0.5}{-0.5 \quad 0.5}.

Para obtener la nueva representación, sencillamente hacemos el producto T \mathbf{u}_w para obtener \mathbf{u}_v=\binom{-2}{0}.

Compartimos el archivo transicion.ggb (de GeoGebra) para experimentar en 2D este tipo de conversión.

Figura generada por transicion.ggb

Ejemplo de integración numérica doble

Los métodos de integración numérica de funciones de una variable, se pueden extender a integrales múltiples. En particular el método de cuadratura gaussiana es conveniente por el menor número de operaciones requeridas, en comparación (p. ej.) con la regla compuesta de Simpson. A continuación presentamos un caso de integración doble [1]

Por ejemplo usando una generalización del método de cuadratura gaussiana a dos dimensiones, podemos calcular V=\int_0^1 \int_0^{-x+1} cos^2(x^2+y^2)dydx mediante el siguiente código:

(1.4270959233395155
  1.4307179095080604
  1.4304856648923563
  1.4304734682663138
  1.430474413405907
  1.4304744206786033
  1.4304744197597712)

Es decir, obtenemos V=1.4304744. La función map permite evaluar el volúmen para diferentes valores de nodos (en este caso con igual número en ambos ejes). Compartimos el archivo ejemplo1int2d.rkt [requiere: basicas.rkt e integraciones.rkt] que incluye el código anterior y la generación de la siguiente gráfica asociada

Figura generada por ejemplo1int2d.rkt

Observación: La figura generada en el ambiente de DrRacket puede ser reorientada con ayuda del mouse. Estas gráficas se pueden generar a partir de la versión 5.2 de DrRacket.

Como referencia compartimos también el archivo ejemplo1int2d.sce (de Scilab) que ilustra cómo realizar este mismo cálculo. La  “restricción” del comando int2d en Scilab es que es necesario definir la región de integración en términos de un conjunto de triángulos en el plano. Observación: esto no representa realmente una restricción, sino que señala (de manera indirecta) la importancia de las triangulaciones en geometría computacional, en particular en sus aplicaciones a ingeniería (ver p. ej. notas de mesh generation, del curso Applied Parallel Computating, MIT).

Nota [1] Cuando el número de dimensiones aumenta a más de tres, hay que recurrir a otros métodos, entre ellos el llamado método de Monte Carlo.

Bases y conjunto de soluciones para un SEL homogéneo

Una base \mathbf{B} para un espacio vectorial \mathbf{V} cumple las siguientes condiciones:

  1. El conjunto \mathbf{B} está formado por vectores linealmente independientes.
  2. Los vectores en \mathbf{B} generan el espacio vectorial \mathbf{V}.

Una de las aplicaciones fundamentales de este concepto, es el poder expresar el conjunto de soluciones de un SEL (sistema de ecuaciones lineales) homogéneo de n variables, mediante una base para el conjunto de sus soluciones como vectores en \mathbb{R}^n.

Por ejemplo, el sistema:

\begin{matrix} x+z+3w & = 0 \\ 3x-y+2z+w &=0 \end{matrix}

Se puede transformar a su forma escalonada reducida por renglones:

1 0 1 3
0 1 1 8

Lo cual es equivalente a las siguientes dos ecuaciones:

\begin{matrix} x+z+3w & = 0 \\ y + z + 8w &=0 \end{matrix}

Tomando los parámetros s=z y t=w, podemos despejar x y y para obtener el conjunto de soluciones:

S=\{ (-s-3t,-s-8t,s,t)\in \mathbb{R}^4\; | \;s,t \in \mathbb{R} \}

Utilizando ahora el concepto de combinación lineal, con \mathbf{v}_1=(-1,-1,1,0)^{T} y \mathbf{v}_2=(-3,-8,0,1)^{T}  podemos escribir:

S=\{ s\mathbf{v}_1+t\mathbf{v}_2\; | \;s,t \in \mathbb{R} \}

Decimos entonces que \mathbf{B}=\{ (-1,-1,1,0)^{T}, (-3,-8,0,1)^{T} \} es una base para el conjunto de soluciones del sistema lineal en cuestión.

Compartimos el archivo espaciosol.ggb con el que puede visualizar (hasta cierto punto) el espacio solución indicado:

Figura generada con espaciosol.ggb

Actividad: ¿Podría verificar que los elementos de la base \mathbf{B} anterior, son linealmente independientes?

Propiedades de traslación para transformadas de Laplace

Los siguientes teoremas son muy útiles para resolver ciertos PVI mediante transformada de Laplace:

  1. [traslación en el eje s] Si \mathcal{L}\{ f(t) \}=F(s)  y a \in \mathbb{R}, entonces \mathcal{L}\{e^{at} f(t) \}=F(s-a). Para este caso usamos la notación: \mathcal{L}\{e^{at} f(t) \}=\mathcal{L}\{ f(t) \} |_{s \rightarrow s-a}
  2. [forma inversa de la traslación en el eje s] \mathcal{L}^{-1}\{ F(s-a) \}=\mathcal{L}^{-1}\{ F(s) \} |_{s \rightarrow s-a}=e^{at} f(t)
  3. [traslación en el eje t] Si F(s)=\mathcal{L} \{ f(t) \} y a>0, entonces \mathcal{L} \{ f(t-a) \;\mathcal{U}(t-a) \}=e^{-as}F(s). También podemos escribir de manera alternativa: \mathcal{L} \{g(t) \;\mathcal{U}(t-a) \}=e^{-as}\mathcal{L} \{g(t+a) \}
  4. [forma inversa de la traslación en el eje t] Si f(t)=\mathcal{L}^{-1} \{ F(s) \} y a>0, entonces \mathcal{L}^{-1} \{ e^{-as}F(s) \}=f(t-a) \;\mathcal{U}(t-a)

Veamos un ejemplo para cada uno de los casos anteriores:

  1. \mathcal{L}\{ e^{-3t} t^2 \} = (L#50) \mathcal{L}\{ t^2 \} |_{s \rightarrow s+3}= (L#3) \frac{2}{s^3} |_{s \rightarrow s+3} = \frac{2}{(s+3)^3} \star
  2. \mathcal{L}^{-1}\{ \frac{s-2}{(s-2)^2+1} \}=\mathcal{L}^{-1} \{ \frac{s}{s^2+1}\} |_{s \rightarrow s-2} = (L#50) e^{2t} \mathcal{L}^{-1} \{ \frac{s}{s^2+1}\} = (L#8) e^{2t} cost \star
  3. \mathcal{L} \{ e^{-2t} \mathcal{U}(t-5) \} = (L#53) e^{-5s} \mathcal{L} \{ e^{-2(t+5)} \} = e^{-5s} \mathcal{L} \{ e^{-2t} e^{-10}\} = e^{-5s-10} \mathcal{L} \{ e^{-2t} \} \\ \quad =(L#11)e^{-5s-10}\frac{1}{s+2} = e^{-10} \displaystyle \frac{e^{-5s}}{s+2}\star
  4. \mathcal{L}^{-1}\{ e^{-5s}\frac{2}{s^2+4} \}= (L#52) \mathcal{L}^{-1}\{ \frac{2}{s^2+4} \} |_{t \rightarrow t-5} \mathcal{U}(t-5) = (L#7) sen(2t) |_{t \rightarrow t-5}\mathcal{U}(t-5) \\ \quad = sen(2(t-5)) \mathcal{U}(t-5) = sen(2t-10) \mathcal{U}(t-5)\star

Compartimos el archivo traslacioneslap.wxm (de WxMaxima) como ejemplo del potencial de este excelente software de cómputo numérico-simbólico.

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.

Combinación lineal y espacio generado

Dados los vectores \mathbf{v}_1,\mathbf{v}_2 de un espacio vectorial V, llamamos combinación lineal de \mathbf{v}_1,\mathbf{v}_2 a cualquier vector de la forma: \mathbf{v}=c_1\mathbf{v}_1 +c_2\mathbf{v}_2, donde c_1,c_2\in \mathbf{F} (i.e. son escalares, por ejemplo, en \mathbb{R}). Lamamos espacio generado por \mathbf{v}_1,\mathbf{v}_2 al conjunto:

gen(\{\mathbf{v}_1,\mathbf{v}_2 \})=\{c_1\mathbf{v}_1 +c_2\mathbf{v}_2 \; | c_1,c_2 \in \mathbf{F})

Lo anterior se generaliza de forma natural a n vectores.

Ejemplo: Si \mathbf{v}_1=x-1 y \mathbf{v}_2=x^2, ambos considerados como vectores en \mathbf{P}_n, y \mathbf{F} = \mathbb{R} , entonces tenemos que \mathbf{v}=x^2+x-1 es una combinación lineal de dichos vectores y que

\mathbf{G}=gen(\{x-1,x^2 \})=\{c_1 x - c_1 + c_2 x^2 \; | c_1,c_2 \in \mathbb{R})

Note que dicho espacio generado es menor al espacio \mathbf{P}_2. Por ejemplo x-x^2 \not\in \mathbf{G}.

Figura generada con contraejemploclin.ggb

Puede explorar el archivo contraejemploclin.ggb (de GeoGebra) y experimentar con otros ejemplos.

Método de cuadratura gaussiana

Uno de los métodos más eficientes para integración numérica es la cuadratura de Gauss (también llamada cuadratura de Gauss-Legendre).

Para aproximar \int_a^b f(x)dx primero consideramos la ecuación de cuadratura \int_{-1}^1 f(x)dx \approx \displaystyle\sum_{i=1}^n w_i f(x_i) donde las x_i son las raices del polinomio de Legendre de grado n [ver nota al pie] y los w_i son los pesos en la cuadratura de Gauss que se calculan como:

w_i = \int_{-1}^1 \displaystyle \prod_{j=1,j\neq i}^n \frac{x-x_j}{x_i-x_j}dx

Finalmente para calcular la integral original \int_a^b f(x)dx realizamos un cambio de variable x=\frac{b-a}{2}t + \frac{a+b}{2} para obtener:

\int_a^b f(x)dx = \frac{b-a}{2}\int_{-1}^1 f(\frac{b-a}{2}t + \frac{a+b}{2})dx

Compartimos el archivo en Scheme cuadratura.scm, el cual implementa la cuadratura de Gauss hasta 8 nodos, e incluye casos para integración doble y triple (que veremos en un próximo segmento).

Como ejemplo de la evaluación en integración simple, podemos aproximar \int_0^{1/2} 5 x e^{-5x}dx, para diferentes números de nodos, con el código siguiente:

(0.14649853795928922
 0.14262865093540064
 0.14254137312775894
 0.14254050586883563
 0.14254050097880808
 0.14254050096928556
 0.14254050097726964)

Figura generada con ejemploint2.ggb

El valor calculado usando integración por partes es: 0.142540500963271, como puede verificar con el archivo ejemploint2.ggb (de GeoGebra). Este valor coincide en 11 cifras decimales, cuando se utilizan 7 nodos en el método de cuadratura gaussiana.

Referencia recomendada: González, Purificación (c. 2012) Integración numérica (pdf, 21 pp., sección 7.5)

Nota: Los polinomios de Legendre se pueden generar fácilmente utilizando la fórmula recursiva de Bonnet [ver Legendre polynomials en Wikipedia]

P_0(x)=1, \quad P_1(x)=x

(n+1) P_{n+1}(x)=(2n+1) x P_n (x) - n P_{n-1} (x)

Subespacios de un espacio vectorial

Dado un subconjunto H de un espacio vectorial V con operaciones asociadas, podemos determinar si H es en sí mismo un espacio vectorial (es decir un subespacio vectorial) de V si es que se cumplen los dos axiomas de cerradura.

  • [cerradura bajo la suma de vectores] \forall x,y \in H, se cumple x \oplus y \in H
  • [cerradura bajo la multiplicación por escalar] Si x \in H, y \alpha es un escalar, entonces \alpha \otimes x \in H.

Ejemplo 1: Sea V un plano en \mathbb{R}^3 que pasa por el origen, con las operaciones normales de vectores en \mathbb{R}^3, entonces cualquier recta H que pasa por el origen y que es subconjunto de dicho plano, es un subespacio vectorial de V. ¿Podría justificar lo anterior aplicando la suma normal de vectores en \mathbb{R}^3 y la multiplicación por escalares del campo \mathbb{R}? (compartimos el archivo subespacios.rkt que genera la siguiente figura, que bajo DrRacket puede mover con el mouse)

Figura generada con Plot (de DrRacket) versión 5.2

Ejemplo 2: Si V=\mathbb{R}^2 y H={(x,y): x^2+y^3<1}. ¿Se cumple que H es un subespacio vectorial de V?. Resp.  Podemos argumentar que no es un subespacio vectorial porque no cumple los axiomas de cerradura, en particular: sea x=\frac{1}{2} e y=\frac{1}{2}, se cumple que x^2+y^3=\frac{1}{2}^2+\frac{1}{2}^3<1, pero 2(x,y)=(2x,2y)=(1,1) no cumple la segunda cerradura, ya que 1^2+1^3\not < 1. Puede explorar el archivo nosubespacio.ggb que genera la siguiente figura.

Figura generada con GeoGebra 4

Método de Romberg

Uno de los métodos más interesantes de integración numérica es el método de Romberg. En este caso no se requiere indicar como parámetro el número de los subintervalos o la longitud h  de cada subintervalo. Este método opera de manera recursiva realizando tantas subdivisiones como sea necesario del subintervalo original [a,b] hasta obtener la exactitud \epsilon buscada. La estrategia consiste en aplicar una combinación de la regla del trapecio y la extrapolación de Richardson, a particiones sucesivas del intervalo [a,b]

Compartimos el archivo romberg.scm que implementa de manera básica (sin optimizaciones) este método. Por ejemplo, para calcular \int_2^5 (\frac{1}{1-x}+\frac{5}{lnx}+\frac{1}{8}sin(10x))dx usando \epsilon = 10^{-8}, usamos el siguiente código:

(r (1 : 1) 13.727206765306258)
(r (2 : 2) 11.651045922540662)
(r (3 : 3) 11.53126127609862)
(r (4 : 4) 11.584473332128743)
(r (5 : 5) 11.55074252580743)
(r (6 : 6) 11.553939909441697)
(r (7 : 7) 11.553866831518363)
11.553867238872222

Figura generada por ejemplo-intf.ggb

Si gusta puede comparar este resultado, con el obtenido en el archivo ejemplo-int1f.ggb (de GeoGebra) que coincide en nueve cifras decimales.

Solución de un PVI mediante Laplace

Problema de Valor Inicial: y''-6y'+9y=t, y(0)=0, y'(0)=1 (ver #25, p. 219 Zill 9/e)

Solución mediante transformada de Laplace [ver artículo en Wikipedia]:

  1. Por las propiedades de linealidad: \mathcal{L}\{y''\}-6\mathcal{L}\{y'\}+9Y(s)=\frac{1}{s^2}
  2. Por la transformada de Laplace de derivadas: s^2Y(s)-sy(0)-y'(0)-6\left[ sY(s)-y(0)\right] +9Y(s)=\frac{1}{s^2}
  3. Aplicamos las condiciones iniciales: Y(s)\left[ s^2-6S+9\right] -1=\frac{1}{s^2}
  4. De donde Y(s)=\frac{1+1/s^2}{s^2-6S+9}=\frac{s^2+1}{s^2(s-3)^2}
  5. Descomponemos en fracciones parciales: Y(s)=\frac{As+B}{s^2}+\frac{C}{s-3}+\frac{D}{\left( s-3 \right) ^2}
  6. Por lo cual: s^2+1=\left( As+B \right) \left( s-3 \right) ^2+Cs^2 \left( s-3 \right) +Ds^2
  7. Para s=0 obtenemos: 1=9B de donde B=\frac{1}{9}
  8. Para s=3 obtenemos: 10=9D por lo cual D=\frac{10}{9}
  9. Con s=1 y s=-1 obtenemos dos ecuaciones \frac{1}{9}=A-\frac{1}{2}C y -\frac{1}{18}=-A-\frac{1}{4}C las cuales resolvemos fácilmente por eliminación, para obtener: C=-\frac{2}{27} y A=\frac{2}{27}
  10. Por lo anterior Y(s)=\left( \frac{2}{27}s+\frac{1}{9}\right)\frac{1}{s^2}-\frac{2}{27}\frac{1}{s-3}+\frac{10}{9}\frac{1}{\left( s-3 \right)^2}
  11. Y entonces \mathcal{L}^{-1} \{ Y(s) \} =\frac{2}{27}\mathcal{L}^{-1} \{ \frac{1}{s} \} +\frac{1}{9}\mathcal{L}^{-1} \{ \frac{1}{s^2}\} -\frac{2}{27}\mathcal{L}^{-1} \{ \frac{1}{s} \| _{s \to s-3} \} \\ \quad+\frac{10}{9}\mathcal{L}^{-1} \{ \frac{1}{s^2} \| _{s \to s-3} \}
  12. Por lo cual: y(t)=\frac{2}{27}\cdot 1+\frac{1}{9}\cdot t-\frac{2}{27}e^{3t}\cdot 1+\frac{10}{9}e^{3t}\cdot t
  13. Que reescribimos como: y(t)=\frac{2}{27}+\frac{t}{9}-\frac{2}{27}e^{3t}+\frac{10}{9}te^{3t}

Se recomienda reconstruir a mano por su cuenta los pasos anteriores, después utilizar WxMaxima para ejercitarse con dicha herramienta (y automatizar parte del proceso) para finalmente visualizar en GeoGebra que la solución cumple con el PVI dado.