Ejemplo de integración numérica doble

Posted on 2011/11/24 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 2011/11/22 by Matematikai

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

El conjunto $\mathbf{B}$ está formado por vectores linealmente independientes.
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 2011/11/22 by Matematikai

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

[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}$
[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)$
[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) \}$
[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:

$\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$
$\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$
$\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$
$\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 2011/11/18 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 2011/11/17 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 2011/11/15 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 2011/11/15 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.

Solución de un PVI mediante Laplace

Posted on 2011/11/11 by matematikaizen

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]:

Por las propiedades de linealidad: $\mathcal{L}\{y''\}-6\mathcal{L}\{y'\}+9Y(s)=\frac{1}{s^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}$
Aplicamos las condiciones iniciales: $Y(s)\left[ s^2-6S+9\right] -1=\frac{1}{s^2}$
De donde $Y(s)=\frac{1+1/s^2}{s^2-6S+9}=\frac{s^2+1}{s^2(s-3)^2}$
Descomponemos en fracciones parciales: $Y(s)=\frac{As+B}{s^2}+\frac{C}{s-3}+\frac{D}{\left( s-3 \right) ^2}$
Por lo cual: $s^2+1=\left( As+B \right) \left( s-3 \right) ^2+Cs^2 \left( s-3 \right) +Ds^2$
Para $s=0$ obtenemos: $1=9B$ de donde $B=\frac{1}{9}$
Para $s=3$ obtenemos: $10=9D$ por lo cual $D=\frac{10}{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}$
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}$
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} \}$
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$
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.

Regla compuesta de Simpson

Posted on 2011/11/10 by matematikaizen

Para aproximar la integral definida $\int_a^b f(x)dx$ , uno de los métodos clásicos, aunque no el más eficiente, es la regla compuesta de Simpson (ver el artículo Regla de Simpson, en Wikipedia). Si dividimos un intervalo $[a,b]$ en un número par de subintervalos $n$ , podemos definir $h=\frac{b-a}{n}$ y $x_i=a+ih$ , para luego aplicar la fórmula:

$\int_a^b f(x)dx \approx \frac{h}{3}[f(a)+f(b)+2 \displaystyle\sum_{j=1}^{\frac{n}{2}-1} f(x_{2j})+4 \displaystyle\sum_{j=1}^{n/2} f(x_{2j-1})]$

la cual nos aproxima la integral con un error máximo igual a

$(b-a)\frac{h^4}{180}max_{a\leq \gamma \leq b} | f^{(4)}(\gamma)| \quad (*)$

Compartimos los archivos simpson.scm y simpsonplot.rkt (para DrRacket) que incluyen una implementación básica de esta regla. En la implementación, note que si se entrega $h$ como parámetro, entonces podemos determinar el valor de $n$ como sigue: Sea $m=\lceil \frac{b-a}{h} \rceil$ . Si $m$ es par, entonces $n:=m$ en otro caso $n:=m+1$ .

Consideremos la integral $\int_0^5 t cos^2(t) \sqrt{1+t^2}dt$

Como ejemplo del código Scheme, podemos aproximar dicha integral por la regla compuesta de Simpson para diferentes valores de $h$ , como sigue:

>(map (lambda(h)
       (simpson (lambda(t)
                  (* t (expt (cos t) 2)
                     (sqrt (+ 1 (* t t)))))
                0 5 h))
     '(0.1 0.01 0.001 0.0001 0.00001))

(17.423338274907277
17.423337562707726
17.423337562644544
17.42333756264443
17.423337562644658)

La siguiente gráfica nos permite visualizar que el reducir el valor de $h$ en general reduce el error, sin embargo, debemos tener presente, que el error máximo depende también del valor máximo de $f^{(4)}$ en el intervalo bajo integración.

Nota: Para algunas aplicaciones, podemos combinar de forma efectiva los métodos de derivación e integración numérica. Para tal fin, compartimos como referencia los archivos derivaciones.scm y derivacionesplot.rkt. Mucho éxito en sus matemáticas experimentales.

Axiomas de un espacio vectorial

Posted on 2011/11/08 by matematikaizen

El enfoque axiomático es fundamental para definir objetos matemáticos ya que nos ayudan a determinar si se cumplen o no ciertas propiedades. Decimos que un conjunto de objetos matemáticos $V$ , junto con dos operaciones $\oplus$ y $\otimes$ (que llamamos la suma vectorial y la multiplicación de un escalar¹ por un vector, respectivamente) forman un espacio vectorial, si se cumplen los diez axiomas siguientes:

[cerradura bajo la suma de vectores] $\forall x,y \in V$ , se cumple $x \oplus y \in V$
[asociatividad de la suma de vectores] $\forall x,y,z \in V$ , se cumple $(x \oplus y) \oplus z = x \oplus (y \oplus z)$
[elemento identidad de la suma de vectores] $\exists 0 \in V$ tal que $0 \oplus x = x \oplus 0 = x$
[elemento inverso aditivo] Si $x \in V$ , $\exists \; -x \in V$ tal que $x + (-x) = 0$
[conmutatividad de la suma de vectores] $\forall x,y \in V$ , se cumple $x \oplus y = y \oplus x$
[cerradura bajo la multiplicación por escalar] Si $x \in V$ , y $\alpha$ es un escalar, entonces $\alpha \otimes x \in V$ .
[primera ley distributiva] Si $x,y \in V$ y $\alpha$ un escalar, entonces $\alpha \otimes (x \oplus y)=(\alpha \otimes x) \oplus (\alpha \otimes y)$
[segunda ley distributiva] Si $x \in V$ y $\alpha, \beta$ son escalares, entonces $(\alpha + \beta) \otimes x = (\alpha \otimes x) \oplus (\beta \otimes x)$
[asociatividad de multiplicación por escalares] Si $x \in V$ y $\alpha, \beta$ son escalares, entonces $\alpha (\beta \otimes x) = (\alpha \beta) \otimes x$
[identidad de multiplicación por escalar] Para cada vector $x \in V$ , se cumple $1 \otimes x = x$

La belleza de estos axiomas es permiten identificar como espacios vectoriales las estructuras $(V, \oplus, \otimes)$ más inesperadas. Se le invita a revisar el excelente artículo Espacio vectorial [ en Wikipedia].

Compartimos el archivo de GeoGebra vectoressyp.ggb [que por cierto está estrenando la versión 4.0.8], en el que pueden interactuar para visualizar la suma de vectores y producto de escalar con vector, para el caso de $\mathbb{R}^2$ .

Nota: ¹Al conjunto del que se obtienen los escalares le llamamos campo (o cuerpo) [ver artículo Cuerpo (matemáticas) en Wikipedia]. Entre los campos principales se tienen los reales $\mathbb{R}$ y los complejos $\mathbb{C}$ (además de los campos finitos). Como nota adicional, los campos y sus extensiones son claves para entender la famosa Teoría de Galois.

matikai.com

Matemáticas en pequeños segmentos para disfrutar