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

Posted on 2011/12/09 by Matematikai

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

>sols
((0 0 (10 7 7))
 (1 0.005 (9.8783782 7.9989081 7.2774597))
 (2 0.01 (9.8106654 8.9683908 7.5958386))
 (3 0.015 (9.7928907 9.9121323 7.9547719))  
...
 (5997 29.985 (11.0238911 2.493518 37.9423877))
 (5998 29.99 (10.5944435 1.9552066 37.5591314))
 (5999 29.995 (10.1607752 1.4611159 37.1495811))
 (6000 30.0 (9.725236 1.0106386 36.7184511)))

>sols

((0 0 (10 7 7))

(1 0.005 (9.8783782 7.9989081 7.2774597))

(2 0.01 (9.8106654 8.9683908 7.5958386))

(3 0.015 (9.7928907 9.9121323 7.9547719))

...

(5997 29.985 (11.0238911 2.493518 37.9423877))

(5998 29.99 (10.5944435 1.9552066 37.5591314))

(5999 29.995 (10.1607752 1.4611159 37.1495811))

(6000 30.0 (9.725236 1.0106386 36.7184511)))

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]

Solución de un sistema de ED mediante Laplace

Posted on 2011/12/06 by Matematikai

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:

$\mathcal{L}\{\frac{dx}{dt}=1-y\}\rightarrow sX(s)-1=\frac{1}{s}-Y(s)$
$\mathcal{L}\{\frac{dy}{dt}=x \} \rightarrow sY(s)=X(s)$ , por lo cual
$Y(s)=\frac{1}{s}X(s)$ , que sustituimos en el paso 1 para obtener:
$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}$
Podemos ahora aplicar la transformada de Laplace inversa para obtener: $x(t)=cost+sent$ (usando las fórmulas L#7 y L#8)
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)}$
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

Posted on 2011/12/03 by Matematikai

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

(define (fxy x y)

(+ (expt (cos (+ (* x x)

(* y y)))

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.

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

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

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

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

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

Polinomio interpolador de Lagrange

Posted on 2011/11/07 by matematikaizen

El polinomio interpolador de Lagrange, nos permite construir una función polinomial de grado $n$ garantizando que dicha función pasará por los puntos $(x_0,f(x_0)),(x_1,f(x_1),\cdots ,(x_n,f(x_n))$ . Para la definición y más detalles, le invitamos a revisar el artículo Interpolación polinómica de Lagrange [Wikipedia]

Compartimos los archivos lagrange.scm y lagrangeplot.rkt (archivo para DrRacket, ver [1]) donde encontrarán ejemplos de uso, entre ellos:

>(poli-lagrange '(2 5/2 4) '(1/2 2/5 1/4))

El resultado previo indica que el polinomio interpolador asociado a los puntos $(2,1/2), (5/2,2/5),(4,1/4)$ es: $p(x)=\frac{1}{20}x^2 -\frac{17}{40}x+1 \frac{3}{20}$ .

Las instrucciones siguientes ilustran la forma de graficar en Scheme (incluyendo Plot), el polinomio interpolador de Lagrange asociado a otro conjunto de puntos:

>(plot (line (lagrangef '(-3 -1 0 2 3) '(1 2 1 -2 -1)))
 #:title "(lagrangef '(-3 -1 0 2 3) '(1 2 1 -2 -1))")

Mucho éxito en sus exploraciones matemáticas

Nota [1]. Los archivos lagrange.scm y lagrangeplot.rkt son casi idénticos, excepto que langrangeplot.rkt corre bajo DrRacket e incluye la funcionalidad de graficación elemental con plot. El primer archivo lagrange.scm, es compatible con R5RS (en un futuro próximo, esperamos realizarlo bajo el estándar R7RS). Gracias

matikai.com

Matemáticas en pequeños segmentos para disfrutar

Tag Archives: scheme