Author Archives:

Ejemplo de integración numérica doble

Posted on by Matematikai

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:

(define (fxy x y)
  (+ (expt (cos (+ (* x x)
                   (* y y)))
           2)
     2))

(define (vol n)
  (int2g 0 1
         (λ(x) 0)
         (λ(x) (+ (- x) 1))
         fxy
         n n))

(map vol '(2 3 4 5 6 7 8))

(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

Posted on by Matematikai

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

Posted on by Matematikai

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) = 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

Posted on by Matematikai

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

Posted on by Matematikai

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:

>(map  (lambda(n)
         (intg (lambda(x)
                  (* 5 x (exp (* -5 x))))
                0 1/2 n))
      '(2 3 4 5 6 7 8))

(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: Buchanan, J. Robert (2006 ) Gaussian Quadrature (pdf, 21 pp.)

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

Posted on by Matematikai

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

Posted on by Matematikai

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:

>(romberg (lambda(x)
            (+ (/ 1 (- 1 x))
               (/ 5 (log x))
               (* 1/8 (sin (* 10 x)))))
          2 5 1e-8)

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

Potencias y raíces de números complejos

Posted on by Matematikai

Dado un número complejo en su forma polar z=re^{i\theta}, podemos calcular la potencia y raices n-ésimas mediante:

z^n=r^ne^{i\theta n}

z^{1/n}=r^{1/n}exp\{i(\frac{\theta}{n}+\frac{2\pi k}{n})\} para k=0,1,\dots ,n-1

Bajo GeoGebra podemos definir los ángulos asociados a las raíces mediante el comando iterationList[] con los argumentos x+\frac{2 \pi}{n}, \frac{\theta}{n} y n-1. La siguiente figura ilustra las potencias de un número complejo (ver archivo potencias.ggb)

Figura generada por potencias.ggb

Visualizamos también las raíces mediante el archivo raices.ggb , que genera la siguiente figura:

Figura generada con raices.ggb

Referencias selectas: