Resolución numérica de ecuaciones diferenciales con RK.

[Ghiret]

Hola a todos, escribo esto por que ando intentando resolver una ecuación de segundo orden de forma numérica pero no consigo entender el método de Runge-Kutta. La ecuación es:

[Error LaTeX: Compilación LaTeX fallida]

donde todo es constante menos $x$ y $\dot{x}$. Ésta se puede reducir a un sistema de ecuaciones diferenciales de primer orden con el cambio $\dot{x}=v$ y $\dot{v}=g(x,v)$, siendo $g(x,v)$ el chorizo de ahí arriba.

Si uso un Runge-Kutta de cuarto orden, se supone que tengo que hacer:

$<br /> k_{1,1} = hg(x_{n},v_{n})<br /> \\<br /> k_{1,2}=hg\left(x_{n}+\frac{h}{2},v_{n}+\frac{k_{1}}{2}\right)<br /> \\<br /> k_{1,3}=hg\left(x_{n}+\frac{h}{2},v_{n}+\frac{k_{2}}{2}\right)<br /> \\<br /> k_{1,4}=hg\left(x_{n}+h,v_{n}+k3\right)<br /> \\<br /> v_{n+1} = v_{n}+\frac{1}{6}\left( k_{1}+2k_{2}+2k_{3}+k_{4}\right)<br />$

y luego calcular $\dot{x}$ usando un esquema parecido y así reiteradamente. Mi problema es que no sé que función $f(t,x)$ tengo que usar para resolver la ecuación $\dot{x}=v$; quizás es que me he liado, pero no sé por donde meterle mano a la segunda ecuación del sistema y en google no he encontrado nada(o nada suficientemente simple como para que yo lo entendiera).

[pod]

Si tienes:

$\ddot x = g(x, \dot x) \ .$

y defines $v = \dot x$, entonces te quedaan dos ecuaciones

$\dot v = g(x, v) \ , \quad \dot x = v = f(x ,v) \ .$

Pues esa es la función $f(x,v)$ que tienes que utilizar.

Yo debo tener un programa en fortran que resuelve ese tipo de ecuaciones de segundo orden, de cuando hacía la carrera. Si quieres podria enviártelos, aunque no es que sean un ejemplo de elegancia. Estan programados para detectar secciones de Poincaré, o sea que no dan toda la trayectoria. Pero eso sería una modificación trivial Además en mi caso tenia dos funciones, lo que en tu caso vendría a ser $x_1$ y $x_2$.

[Ghiret]

Si tienes:
$<br /> \ddot x = g(x, \dot x) \ .<br />$
y defines v = \dot x, entonces te quedaan dos ecuaciones

$<br /> \dot v = g(x, v) \ , \quad \dot x = v = f(x ,v) \ .<br />$
Pues esa es la función $f(x,v)$ que tienes que utilizar.

Ahí está mi problema, que no sé cuál es la forma particular(para la ecuación dada) que tengo que tomar para $f(x,v)$ ni tampoco como calcular las $k$'s para hacerle Runge-Kutta a esa expresión.

Yo debo tener un programa en fortran que resuelve ese tipo de ecuaciones de segundo orden, de cuando hacía la carrera. Si quieres podria enviártelos, aunque no es que sean un ejemplo de elegancia. Estan programados para detectar secciones de Poincaré, o sea que no dan toda la trayectoria. Pero eso sería una modificación trivial Wink Además en mi caso tenia dos funciones, lo que en tu caso vendría a ser $x_1$ y $x_2$

Te lo agradecería mucho, a ver si así me aclaro como funciona el método. Tienes mi dirección de correo por ahí, ¿verdad?

[pod]

Si tienes:
$<br /> \ddot x = g(x, \dot x) \ .<br />$
y defines v = \dot x, entonces te quedaan dos ecuaciones

$<br /> \dot v = g(x, v) \ , \quad \dot x = v = f(x ,v) \ .<br />$
Pues esa es la función $f(x,v)$ que tienes que utilizar.

Ahí está mi problema, que no sé cuál es la forma particular(para la ecuación dada) que tengo que tomar para $f(x,v)$ ni tampoco como calcular las $k$'s para hacerle Runge-Kutta a esa expresión.

Pero si te lo he puesto clarisimo:

$f(x,v) = v \ .$

Yo debo tener un programa en fortran que resuelve ese tipo de ecuaciones de segundo orden, de cuando hacía la carrera. Si quieres podria enviártelos, aunque no es que sean un ejemplo de elegancia. Estan programados para detectar secciones de Poincaré, o sea que no dan toda la trayectoria. Pero eso sería una modificación trivial Wink Además en mi caso tenia dos funciones, lo que en tu caso vendría a ser $x_1$ y $x_2$

Te lo agradecería mucho, a ver si así me aclaro como funciona el método. Tienes mi dirección de correo por ahí, ¿verdad?

Pues no lo tengo.

La parte importante es:

c              * Mètode RK4 d'evolució dinàmica
c                * 1r càlcul a t
                 q01  = q1(i)
                 q02  = q2(i)
                 dq01 = p1(i)
                 dq02 = p2(i)
                 dp01 = dp1(i, q01, q02)
                 dp02 = dp2(i, q01, q02)

c                * 2n càlcul a t+dt/2
                 q11  = q01 + dq01*dt/dfloat(2)
                 q12  = q02 + dq02*dt/dfloat(2)
                 dq11 = dq01+ dp01*dt/dfloat(2)
                 dq12 = dq02+ dp02*dt/dfloat(2)
                 dp11 = dp1(i, q11, q12)
                 dp12 = dp2(i, q11, q12)

c                * 3r càlcul a t+dt/2
                 q21  = q01 + dq11*dt/dfloat(2)
                 q22  = q02 + dq12*dt/dfloat(2)
                 dq21 = dq01+ dp11*dt/dfloat(2)
                 dq22 = dq02+ dp12*dt/dfloat(2)
                 dp21 = dp1(i, q21, q22)
                 dp22 = dp2(i, q21, q22)

c                * 4t càlcul a t+dt
                 q31  = q01 + dq21*dt
                 q32  = q02 + dq22*dt
                 dq31 = dq01+ dp21*dt
                 dq32 = dq02+ dp22*dt
                 dp31 = dp1(i, q31, q32)
                 dp32 = dp2(i, q31, q32)

c                * Sumació final 
                 q1n(i) = q1(i) + dt*( ( dq01 + dq31 )/dfloat(6) 
     +                + ( dq11 + dq21 )/dfloat(3) )

                 q2n(i) = q2(i) + dt*( ( dq02 + dq32 )/dfloat(6) 
     +                + ( dq12 + dq22 )/dfloat(3) )

                 p1n(i) = p1(i) + dt*( ( dp01 + dp31 )/dfloat(6)
     +                + ( dp11 + dp21 )/dfloat(3) )

                 p2n(i) = p2(i) + dt*( ( dp02 + dp32 )/dfloat(6)
     +                + ( dp12 + dp22 )/dfloat(3) )

Ignora el índice i que aparece varias veces, se utilizaba para simular dos sistemas distnitos. La notación es algo distinta, es la que se utiliza en mecanica de Hamilton. Las ecuaciones son:

$\displaystyle \frac{d \mathrm{q}_j}{d t} = p_j \ , \qquad \frac{d \mathrm{p}_j}{d t} = g(p_j, q_j) \ .$

ya que, por suerte, ninguna ecuación depende del tiempo (variable independiente).
En tu caso, como sólo tienes una función incógnita, tendrás la mitad de lineas, deberas eliminar una sí, una no. En el programa, la función "g" está implementada a parte, como puedes ver en llamadas a sendas funciones "dp1" y "dp2".

[Ghiret]

Pero si te lo he puesto clarisimo:
$<br /> f(x,v) = v \ .<br />$

A eso me refiero, a que $v$ no lo conozco de forma analítica como para calcular las $k$'s ni nada, sólo conozco los puntos $v_{n}$ que me da Runge-Kutta aplicado a $v'=g(x,v)$. Mi duda es como usar esos valores de $v$ para calcular las $k$'s y usar Runge-Kutta para calcular losvalores de $x(t)=x_{n}$.

Gracias por el programa, a ver si puedo aclararme mirándolo.

[pod]

Pero si te lo he puesto clarisimo:
$<br /> f(x,v) = v \ .<br />$

A eso me refiero, a que $v$ no lo conozco de forma analítica como para calcular las $k$'s ni nada, sólo conozco los puntos $v_{n}$ que me da Runge-Kutta aplicado a $v'=g(x,v)$. Mi duda es como usar esos valores de $v$ para calcular las $k$'s y usar Runge-Kutta para calcular losvalores de $x(t)=x_{n}$.

Gracias por el programa, a ver si puedo aclararme mirándolo.

$v$ es la forma analítica de $v$. En primer lugar, no es $g(v,x)$, sino $g(t,x)$. Para evalularlo, que utilizas? Los puntos que te da runge-kutta, no? Pues lo mismo. Las k's para $v$ son similares a lo que tienes tu escrito, pero cambiando $x \to t$ (el primer arguento siempre es el tiempo) y $v \to x$ (te importa variar los valores de la variable que integras en ese momento). Las k's para x:

$k_{1,1} = h f(t_n, v_n) = h v_n \ ,$
$k_ {1,2} = h f\left(t_n + \frac h2, \, v_n + \frac{k_1}2\right) = h\cdot\left(v_n + \frac{k_1}2\right) \ ,$
y etcétera...

Pareces tener un poco de lío... [url=http://en.wikipedia.org/wiki/Runge_kutta]repasa la teoría[/tex].

[pod]

Si tienes:
$<br /> \ddot x = g(x, \dot x) \ .<br />$
y defines v = \dot x, entonces te quedaan dos ecuaciones

$<br /> \dot v = g(x, v) \ , \quad \dot x = v = f(x ,v) \ .<br />$
Pues esa es la función $f(x,v)$ que tienes que utilizar.

Ahí está mi problema, que no sé cuál es la forma particular(para la ecuación dada) que tengo que tomar para $f(x,v)$ ni tampoco como calcular las $k$'s para hacerle Runge-Kutta a esa expresión.

Pero si te lo he puesto clarisimo:

$f(t,v) = v \ .$

Yo debo tener un programa en fortran que resuelve ese tipo de ecuaciones de segundo orden, de cuando hacía la carrera. Si quieres podria enviártelos, aunque no es que sean un ejemplo de elegancia. Estan programados para detectar secciones de Poincaré, o sea que no dan toda la trayectoria. Pero eso sería una modificación trivial Wink Además en mi caso tenia dos funciones, lo que en tu caso vendría a ser $x_1$ y $x_2$

Te lo agradecería mucho, a ver si así me aclaro como funciona el método. Tienes mi dirección de correo por ahí, ¿verdad?

Pues no lo tengo.

La parte importante es:

c              * Mètode RK4 d'evolució dinàmica
c                * 1r càlcul a t
                 q01  = q1(i)
                 q02  = q2(i)
                 dq01 = p1(i)
                 dq02 = p2(i)
                 dp01 = dp1(i, q01, q02)
                 dp02 = dp2(i, q01, q02)

c                * 2n càlcul a t+dt/2
                 q11  = q01 + dq01*dt/dfloat(2)
                 q12  = q02 + dq02*dt/dfloat(2)
                 dq11 = dq01+ dp01*dt/dfloat(2)
                 dq12 = dq02+ dp02*dt/dfloat(2)
                 dp11 = dp1(i, q11, q12)
                 dp12 = dp2(i, q11, q12)

c                * 3r càlcul a t+dt/2
                 q21  = q01 + dq11*dt/dfloat(2)
                 q22  = q02 + dq12*dt/dfloat(2)
                 dq21 = dq01+ dp11*dt/dfloat(2)
                 dq22 = dq02+ dp12*dt/dfloat(2)
                 dp21 = dp1(i, q21, q22)
                 dp22 = dp2(i, q21, q22)

c                * 4t càlcul a t+dt
                 q31  = q01 + dq21*dt
                 q32  = q02 + dq22*dt
                 dq31 = dq01+ dp21*dt
                 dq32 = dq02+ dp22*dt
                 dp31 = dp1(i, q31, q32)
                 dp32 = dp2(i, q31, q32)

c                * Sumació final
                 q1n(i) = q1(i) + dt*( ( dq01 + dq31 )/dfloat(6)
     +                + ( dq11 + dq21 )/dfloat(3) )

                 q2n(i) = q2(i) + dt*( ( dq02 + dq32 )/dfloat(6)
     +                + ( dq12 + dq22 )/dfloat(3) )

                 p1n(i) = p1(i) + dt*( ( dp01 + dp31 )/dfloat(6)
     +                + ( dp11 + dp21 )/dfloat(3) )

                 p2n(i) = p2(i) + dt*( ( dp02 + dp32 )/dfloat(6)
     +                + ( dp12 + dp22 )/dfloat(3) )

Ignora el índice i que aparece varias veces, se utilizaba para simular dos sistemas distnitos. La notación es algo distinta, es la que se utiliza en mecanica de Hamilton. Las ecuaciones son:

$\displaystyle \frac{d \mathrm{q}_j}{d t} = p_j \ , \qquad \frac{d \mathrm{p}_j}{d t} = g(p_j, q_j) \ ,$

ya que, por suerte, ninguna ecuación depende del tiempo (variable independiente).
En tu caso, como sólo tienes una función incógnita, tendrás la mitad de lineas, deberas eliminar una sí, una no. En el programa, la función "g" está implementada a parte, como puedes ver en llamadas a sendas funciones "dp1" y "dp2".

diez caracteres sin quote

[pod]

Si tienes:

$\ddot x = g(x, \dot x) \ .$

y defines $v = \dot x$, entonces te quedaan dos ecuaciones

$\dot v = g(t, x, v) \ , \quad \dot x = v = f(t, x ,v) \ .$

Pues esa es la función $f(x,v)$ que tienes que utilizar.

Yo debo tener un programa en fortran que resuelve ese tipo de ecuaciones de segundo orden, de cuando hacía la carrera. Si quieres podria enviártelos, aunque no es que sean un ejemplo de elegancia. Estan programados para detectar secciones de Poincaré, o sea que no dan toda la trayectoria. Pero eso sería una modificación trivial Además en mi caso tenia dos funciones, lo que en tu caso vendría a ser $x_1$ y $x_2$.

diez caracteres sin quote