Category Archive

You are currently browsing the category archive for the ‘ODE’ category.

Diferencias finitas para la segunda derivada con puntos separados por distancias que siguen una serie geometrica respecto del punto central (I)

noviembre 21, 2013 in Matemáticas, Métodos Numéricos, ODE, PDE | Tags: diferencias finitas, distancia serie geometrica, no equidistancia | 1 comentario

El caso mas sencillo es cuando tenemos cinco puntos:

$(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3,y_3), (x_4,y_4)$ ,

de manera que:

$x_1-x_0 = 2(x_2-x_1) = 2(x_3-x_2) = x_4-x_3$ .

Tal y como escribimos en el anterior post, el polinomio general para cinco puntos quedaria:

$L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x) + y_3 l_3(x) + y_4 l_4(x)$ ,

donde:

$l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4}$ ,

$l_1(x) = \frac{x-x_0}{x_1-x_0} \frac{x-x_2}{x_1-x_2} \frac{x-x_3}{x_1-x_3} \frac{x-x_4}{x_1-x_4}$ ,

$l_2(x) = \frac{x-x_0}{x_2-x_0} \frac{x-x_1}{x_2-x_1} \frac{x-x_3}{x_2-x_3} \frac{x-x_4}{x_2-x_4}$ ,

$l_3(x) = \frac{x-x_0}{x_3-x_0} \frac{x-x_1}{x_3-x_1} \frac{x-x_2}{x_3-x_2} \frac{x-x_4}{x_3-x_4}$ ,

$l_4(x) = \frac{x-x_0}{x_4-x_0} \frac{x-x_1}{x_4-x_1} \frac{x-x_2}{x_4-x_2} \frac{x-x_3}{x_4-x_3}$ .

Dado que tratamos de aproximar la segunda derivada, tenemos, en el caso de equidistancia quedaría:

$\frac{d^2}{dx^2}u_i = \frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12 h^2}$ .

Vamos a ver que queda ahora en nuestro caso. Reescribimos los $l_i(x)$ en función de $x_0$ de manera que, por ejemplo, tenemos:

$l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4} =$

$=\frac{x-x_0-2h}{x_0-x_0-2h} \frac{x-x_0-3h}{x_0-x_0-3h} \frac{x-x_0-4h}{x_0-x_0-4h} \frac{x-x_0-6h}{x_0-x_0-6h} =$

$=\frac{(x-x_0-2h)(x-x_0-3h)(x-x_0-4h)(x-x_0-6h)}{144h^4}$ ,

que para $x=x_2=x_0+3h (2h+h)$ , la segunda derivada queda:

$\frac{d^2}{dx^2} l_0(x) |_{x=x_2} = -\frac{1}{72h^2}$ .

Haciendo lo mismo para todas las derivadas segundas de todos los $l_i(x)$ , obtenemos la siguiente formula en diferencias finitas:

$\frac{d^2}{dx^2}u_i = \frac{-u_{i-3} + 81 u_{u-1} - 160 u_i + 81 u_{i+1} - u_{i+3}}{72h^2}$ ,

donde los índices hacen referencia a una malla inicial equiespaciada.

Para el mismo denominador, antes teniamos:

$\frac{d^2}{dx^2}u_i = \frac{-6u_{i-2}+96u_{i-1}-180u_i+96u_{i+1}-6u_{i+2}}{72 h^2}$ .

Podemos hacer lo mismo para cada uno de los puntos $x=x_0$ , $x=x_1=x_0+2h$ , $x=x_3=x_0+4h$ y $x=x_4=x_0+6h$ :

$\frac{d^2}{dx^2}u_{i-3} = \frac{40 u_{i-3} -243 u_{i-1} + 352 u_i -162 u_{i+1} + 13 u_{i+3}}{36 h^2}$

$\frac{d^2}{dx^2}u_{i-2} = \frac{7 u_{i-3} - 32 u_i + 27 u_{i+1} - 2 u_{i+3}}{36 h^2}$

$\frac{d^2}{dx^2}u_{i+1} =\frac{-2 u_{i-3} + 27 u_{i-1} - 32 u_i -+7 u_{i+3}}{36 h^2}$

$\frac{d^2}{dx^2}u_{i+3} =\frac{13 u_{i-3} -162 u_{i-1} + 352 u_i -243 u_{i+1} + 40 u_{i+3}}{36 h^2}$

En el caso de $x=x_0+4h (3h+h)$ , tenemos:

$\frac{d^2}{dx^2}u_i = \frac{-u_{i-4} + 256_{u-1} - 510 u_i + 256 u_{i+1} - u_{i+4}}{140h^2}$ .

Interpolación de Lagrange y formula en diferencias finitas para la aproximación de las derivadas

noviembre 20, 2013 in Matemáticas, Métodos Numéricos, ODE, PDE | Tags: diferencias finitas, interpolación Lagrange, polinomio Lagrange | 1 comentario

Dados $n+1$ puntos:

$(x_0,y_0), \ldots, (x_n,y_n)$ ,

definimos el polinomio interpolador de Lagrange:

$L(x) := \Sigma_{i=0}^n y_i l_i(x)$

donde:

$l_i(x) := \Pi_{0 \leq j \leq n,j \neq i} \frac{x-x_j}{x_i-x_j}$ .

Con esto, tenemos:

$\frac{d^n}{dx^n}L(x) = \Sigma_{i=0}^n y_i \frac{d^n}{dx^n}l_i(x)$

Supongamos que tenemos $n=2$ :

$(x_0,y_0), (x_1,y_1)$ .

En este caso, tenemos:

$L(x) = y_0 l_0(x) + y_1 l_1(x)$

con:

$l_0(x) = \frac{x-x_1}{x_0-x_1}$

$l_1(x) = \frac{x-x_0}{x_1-x_0}$

Como tenemos dos puntos, únicamente podemos aproximar la primera derivada:

$\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x)$

con:

$\frac{d}{dx} l_0(x) = \frac{1}{x_0-x_1}$

$\frac{d}{dx} l_1(x) = \frac{1}{x_1-x_0}$ ,

de manera que:

$\frac{d}{dx}L(x) = \frac{y_1 - y_0}{x_1 - x_0}$ .

Con esto, tenemos las aproximaciones de primer orden:

$u_x(x_0) = u_x(x_1) = \frac{y_1 - y_0}{x_1-x_0}$ ,

que, en índices, queda:

$\frac{d}{dx}u_i \approx \frac{u_{i+1}-u_i}{h_x}$ ,

$\frac{d}{dx}u_{i+1} \approx \frac{u_{i+1}-u_i}{h_x}$ ,

donde $h_x = x_1 - x_0$ (y que es lógico, ya que la derivada será una constante).

Supongamos ahora que tenemos tres puntos:

$(x_0,y_0),(x_1,y_1),(x_2,y_2)$ .

En este caso tenemos:

$L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x)$ ,

$\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x) + y_2 \frac{d}{dx} l_2(x)$ ,

$\frac{d^2}{dx^2} L(x) = y_0 \frac{d^2}{dx^2} l_0(x) + y_1 \frac{d^2}{dx^2} l_1(x) + y_2 \frac{d^2}{dx^2} l_2(x)$ .

Ahora tenemos:

$l_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{x^2 + (-x_1-x_2)x + x_1x_2}{(x_0-x_1)(x_0-x_2)}$ ,

$l_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{x^2 + (-x_0-x_2)x + x_0x_2}{(x_1-x_0)(x_1-x_2)}$ ,

$l_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{x^2 + (-x_0-x_1)x + x_0x_1}{(x_2-x_0)(x_2-x_1)}$ .

$\frac{d}{dx} l_0(x) = \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)}$ ,

$\frac{d}{dx} l_1(x) = \frac{2x+(-x_0-x_2)}{(x_1-x_0)(x_1-x_2)}$ ,

$\frac{d}{dx} l_2(x) = \frac{2x+(-x_0-x_1)}{(x_2-x_0)(x_2-x_1)}$ .

$\frac{d^2}{dx^2} l_0(x) = \frac{2}{(x_0-x_1)(x_0-x_2)}$ ,

$\frac{d^2}{dx^2} l_1(x) = \frac{2}{(x_1-x_0)(x_1-x_2)}$ ,

$\frac{d^2}{dx^2} l_2(x) = \frac{2}{(x_2-x_0)(x_2-x_1)}$ .

De esta manera, por ejemplo, tenemos que:

$\frac{d^2}{dx^2}L(x) = y_0 \frac{2}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2}{(x_2-x_0)(x_2-x_1)}$ ,

que, en el caso de tener los puntos equiespaciados ( $h_x:=x_{i+1}-x_{i}$ con $0 \leq i \leq 1$ ), queda la aproximación de segundo orden:

$\frac{d^2}{dx^2}u(x) = \frac{y_0 - 2 y_1 + y_2}{h_x^2}$ ,

que además, como era de esperar, es independiente de $x$ :

$\frac{d^2}{dx^2}u_{i-1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}$ ,

$\frac{d^2}{dx^2}u_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}$ ,

$\frac{d^2}{dx^2}u_{i+1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}$ .

¿Qué pasa con la primera derivada? En este caso, el resultado si que depende de $x$ (por lo que tendremos un resultado diferente en función de si la $x$ vale $x_0$ , $x_1$ o $x_2$ ) y lo que obtenemos es:

$\frac{d}{dx} L(x) = y_0 \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}$ ,

por lo que:

$\frac{d}{dx} L(x_0) = y_0 \frac{2x_0 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_0 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_0 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}$

$\frac{d}{dx} L(x_1) = y_0 \frac{2x_1 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_1 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_1 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}$

$\frac{d}{dx} L(x_2) = y_0 \frac{2x_2 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_2 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_2 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}$ ,

que, con equiespaciado, quedan:

$L'(x_0) = \frac{3 y_0 -4 y_1 + y_2}{2 h_x}$ ,

$L'(x_1) = \frac{y_0 -2 y_1 + y_2}{2 h_x}$ y

$L'(x_2) = \frac{y_0 -4 y_1 + 3 y_2}{2 h_x}$ .

Reescritura del Laplaciano nd en diferencias finitas mediante el 1d con fuentes de las (n-1)d restantes

noviembre 18, 2013 in Física, Matemáticas, Métodos Numéricos, ODE, PDE | Tags: diferencias finitas, Laplaciano, Laplaciano 1d, Laplaciano 2d, Laplaciano 3d, Laplaciano nd | Deja un comentario

Laplaciano en cartesianas:

$\Delta u = \Sigma_i \frac{\partial^2}{\partial x_i^2}u$

$\frac{u_{i-1}-2u_i+u_{i+1}}{h^2} = f_i$

$\frac{1}{h^2}u_{i-1} + \frac{1}{h^2}u_{i+1} +\frac{-2}{h^2}u_i= f_i$

$\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2} + \frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2} = f_{i,j}$

$i$ fijo:

$\frac{1}{h_y^2}u_{i,j-1} + \frac{1}{h_y^2}u_{i,j+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_x^2}u_{i-1,j} + \frac{-1}{h_x^2}u_{i+1,j})$

$j$ fijo:

$\frac{1}{h_x^2}u_{i-1,j} + \frac{1}{h_x^2}u_{i+1,j} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_y^2}u_{i,j-1} + \frac{-1}{h_y^2}u_{i,j+1})$

$\frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{h_x^2} + \frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{h_y^2} + \frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{h_z^2} = f_{i,j,k}$

$i,j$ fijos:

$\frac{1}{h_z^2}u_{i,j,k-1} + \frac{1}{h_z^2}u_{i,j,k+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}$

$(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k})$

$i,k$ fijos:

$\frac{1}{h_y^2}u_{i,j-1,k} + \frac{1}{h_y^2}u_{i,j+1,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}$

$(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})$

$j,k$ fijos:

$\frac{1}{h_x^2}u_{i-1,j,k} + \frac{1}{h_x^2}u_{i+1,j,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}$

$(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})$

Aplicación exponencial de GL(n,K)

agosto 20, 2013 in Álgebra, Geometría, Matemáticas, ODE | Tags: aplicación exponencial, trayectorias campo | Deja un comentario

Dado el campo vectorial invariante por la izquierda

$u(a) = (L_a)_* u$ del grupo de Lie $G$ ,

llamamos trayectorias del campo a la soluciones de la ODE:

$\frac{d}{dt}a = u(a)$ ,

y que, fijada una condición inicial $a(0) = a$ , la trayectoria queda unívocamente determinada en forma de serie exponencial:

$a(t) = \exp (tu) a$ con $\exp(tu) = I + tu + \frac{t^2}{2}u^2 + \ldots$ ,

donde $u^r(a) = u(a) \circ u^{r-1}(a)$ .

Las trayectorias:

son analíticas y esta definidas para todo valor de $t$ por el carácter anaítico y completo de los campos invariantes.
que pasan por la unidad del grupo de Lie, $a(0) = e$ , son grupos uniparamétricos de G, ya que cumplen la propiedad $a(s+t)=a(s)a(t)$ . En este caso:

$\exp: \mathfrak{g} \longrightarrow G \, / \, u \mapsto \exp u = a(1)$ .

En el caso de $GL(n,\mathbb{K})$ , hemos visto que los campos invariantes a la izquierda tienen la forma:

$U(A) = AU$ con $A \in GL(n,\mathbb{K})$ y $U \in M(n,\mathbb{K})$ .

El subgrupo uniparamétrico $A(t)$ es la solución a:

$\frac{d}{dt}A = AU$ con $A(0)=I$ ,

y que es:

$A(t) = \exp(tU) = \Sigma_{s=0}^{\infty}\frac{t^s}{s!}U^s$ .

De esta manera, tenemos:

$\exp: \mathfrak{g}(n,\mathbb{K}) \longrightarrow GL(n,\mathbb{K}) \, / \, A = e^U$ ,

es decir, la exponencial matricial habitual.

Nueva minisimulación con 1000 péndulos

octubre 24, 2012 in Física, Gráficos por Computador, Informática, Matemáticas, ODE | Tags: Runge Kutta 23, Runge Kutta 4, Sistema ODE, SPH | Deja un comentario

En la siguiente simulación tenemos $1000$ partículas. Cada partícula $a$ representa un péndulo de masa $m_a$ y posición $(x_a(t),y_a(t),z_a(t))$ , con $z_a(t)$ constante, sometido a las fuerzas:

$F_1 = -m_ag$ , la fuerza de la gravedad, con $g$ constante.
$F_2$ , la fuerza elática de la suspensión, de sentido opuesta a la posición del péndulo y de magnitud $k_1d$ , con $k_1 > 0$ y $d$ la variación de la longitud del péndulo respecto de su longitud en reposo $l$ .
$F_3$ , la fuerza de fricción, con sentido opuesta a la velocidad del péndulo y magnitud $k_2v$ , con $v = \sqrt{(x_a')^2+(y_a')^2}$ y $k_2 > 0$ .

La ley de Newton $F = ma$ queda, para cada uno de los péndulos, es decir, para cada partícula, como el sistema de dos ecuaciones diferenciales de segundo orden:

$m_a \begin{bmatrix} x_a''(t) \\ y_a''(t) \end{bmatrix} = F_1 + F_2 + F_3$

que es:

$m_ax_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})x_a - k_2\frac{x'_a}{\sqrt{x_a'^2+y_a'^2}}$

$m_ay_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})y_a - k_2\frac{y'_a}{\sqrt{x_a'^2+y_a'^2}}$

Podemos expresar este sistema como un sistemad de ecuaciones de primer orden introduciendo incongnitas adicionales: $u_a = x_a'$ y $v_a = y_a'$ , con lo que nos queda:

$\left \{ \begin{array}{rcl} x_a' &=& u_a \\ y_a' &=& v_a \\ m_au_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)x_a - k_2\frac{u_a}{\sqrt{u_a^2+v_a^2}} \\ m_av_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)y_a - k_2\frac{v_a}{\sqrt{u_a^2+v_a^2}} - m_ag \\ \end{array} \right .$

En la siguiente simulación tomamos $m_a$ aleatoria para cada partícula $a$ , por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos $l=1$ , $k_1 = 10000$ y $k_2=1$ en todos los péndulos. Cada uno de éstos se distribuye a lo largo de una espilar cilíndrica y las velocidades iniciales son, para todas ellas, $u_{a_0}=0$ , $v_{a_0} = -0.7$ y $v_{z_0} = 0$ . El intervalo de tiempo es desde $t_i = 0$ hasta $t_f=10$ que se divide en $400$ subintervalos.

El resultado es el siguiente:

Métodos numéricos para ODEs

agosto 18, 2012 in Matemáticas, Métodos Numéricos, ODE | 5 comentarios

En [Rosswog 2009] comenta algunos métodos para ODEs, ya que una de las características relevantes de SPH es que transforma las PDEs de la hidrodinámica en un sistema de ODEs y hacer una simulación en SPH no es mas que resolver dicho sistema.

Las ODEs en astrofísica suelen ser de segundo orden pero podemos transformarlas en un sistema de ODEs de primer orden introduciendo como nuevas variables a las velocidades. De esta manera obtenemos $n$ ecuaciones acopladas de la forma:

$\frac{d}{dt}y_i = f_i(t,y_1,\ldots,y_n)$ con $i = 1,\ldots,n$

Tenemos, en este caso, dos tipos de problemas: los problemas de contorno y los problemas de valor inicial. Este último caso, que es el que nos interesa, se trata de determinar el valor de $y_i^n$ en $t^n$ a partir de la condición inicial $y_i^0$ en $t^0$ .

El método de Euler es el método de integración pedagógico. La solución avanza de $t^n$ a $t^{n+1} \equiv t^n + \Delta t$ según:

$\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^n,\vec{y}^n)$

donde $\vec{y}^{n+1} = \vec{y}(t^{n+1})$ y $\Delta t$ es el paso de tiempo escogido. Es un método explícito de primer orden y es asimétrico en tiempo, pues la derivada es correcta al principio del paso de tiempo pero no al final. De ahí el nombre de Euler hacia adelante. En el caso de posiciones, velocidades y aceleraciones, tenemos:

$\vec{x}^{n+1} = \vec{x}^n + \Delta t \vec{v}^n$

$\vec{v}^{n+1} = \vec{v}^n + \Delta t \vec{a}^n$

La posibilidad contraria es:

$\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^{n+1},\vec{y}^{n+1})$

llamado el método de Euler hacia atrás, que es un método implícito. Todos los siguientes métodos que comentaremos son explícitos.

El siguiente método es el Euler modificado que utiliza el valor de la derivada al principio y al final del paso de tiempo:

$\vec{y}^{n+1} = \vec{y}^n + \Delta t \frac{\vec{f}^n + \vec{f}^{n+1,p}}{2}$

que es un método de segundo orden. En el caso de posiciones, velocidades y aceleraciones tenemos:

$\vec{x}^{n+1} = \vec{x}^n + \Delta t \frac{\vec{v}^n + \vec{v}^{n+1,p}}{2}$

$\vec{v}^{n+1} = \vec{v}^n + \Delta t \frac{\vec{a}^n + \vec{a}^{n+1,p}}{2}$

En el método de Runge-Kutta de segundo orden tomamos la derivada en el punto medio de los intervalos de tiempo tomados:

Tomamos como paso $\frac{\Delta t}{2}$ , calculamos $\vec{y}^{mid} = \vec{y}^n + \frac{1}{2} \Delta t \vec{f}(t^n,\vec{y}^n)$ y evaluamos la derivada en ese punto $\vec{f}^{mid} = \vec{f}(t^{\frac{n+1}{2}},\vec{y}^{mid})$
Tomamos el paso entero con esta derivada: $\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}^{mid}$