Examen 5

Enunciado

FinalJulio2007-1.jpg
FinalJulio2007-2.jpg

Solución

Problema 1

Es muy fácil confundirse en este ejercicio ya que no se sabe bien para que exactamente nos da el código del algóritmo. Lo que el ejercicio pide es resolver la IVP con el mismo método descrito en el álgoritmo que en a su vez resuelve otra IVP que realmente no importa hasta el ejercicio Teórico 1.

El algoritmo implementado allí es Runge-Kutta de orden 2 sobre una ODE de orden 2. Uno puede darse cuenta de que hay una ODE de orden 2 al ver que hay 2 variables con sucesiones, $u$ y $v$. Que es RK2 se puede observar al apreciar que cada una de estas sucesiones tiene la influencia de la suma de 2 variables con el mismo peso multiplicadas por $\frac{1}{2}$.

RK2 puede verse como el siguiente ciclo for

for (x = 0, t = 0; t <= 0.5; t += h)
    k1 = h * f(t, x);
    k2 = h * f(t+h, x + k1);
    x += (k1 + k2)/2;

( Advertencia: Al parecer el algoritmo tiene un error, m = delta u, no delta v. Esto esta aclarado en el ejercicio de la teórica. )
Suponiendo que x, k1 y k2 son vectores columna x = (u,v), k1 = (l1, m1) y k2 = (l2, m2), y como la f(t,x) sale de una ODE de primer grado con 2 incógnitas tiene la siguiente expresión en vector columna f(t, x) = (v, -0.9 * sin(u) - 0.24 * v) y evidentemente, f(t+h, x+k1) = (v + m1, -0.9 sin(u+l1) - 0.24 * (v+m1))

De esta manera podemos dejar el código

h = 0.05;
u = 0.0;
v = 0.0;
for (x = 0.0; x < 1.0; x += h)
    m1 = h * v;
    m2 = h * (v + m1);
    l1 = h * (-0.9 * sin(u) - 0.24 * v);
    l2 = h * (-0.9 * sin(u + l1) - 0.24 * (v + m1));
    u = u + (l1 + l2) * 0.5;
    v = v + (m1 + m2) * 0.5;

que es equivalente al del enunciado.

Dado que se identifica RK2 en la implementación, es necesario ahora resolver el IVP propuesto utilizando otro algoritmo de orden 2 (global), por lo que utilizamos el Método de Taylor. Se dispone de $x'$, y para poder aplicarlo debemos obtener la segunda derivada completa de $x$:
Siendo $x' = F(t, x) = \frac{t-x}{2}$, la derivada 2da de x, es entonces $x'' = F'(t,x) = F_t(t,x) t' + F_x(t,x) x'(t) t' = F_t(t,x) + F_x(t,x) F(t,x)$ (siempre respecto de t estamos derivando, ya que x es una función de t). Por lo tanto

(1)
\begin{align} x'' = \frac{1}{4} (2 - t + x) \end{align}

Entonces la sucesión con Taylor queda:

(2)
\begin{align} x_{i+1} = x_i + h \frac{ih - x_i}{2} + \frac{h^2 (2 - ih + x_i)}{8} \end{align}

Aplicando RK2 la sucesión para esta ODE resulta:

(3)
\begin{align} x_{i+1} = x_i +0.5 \Bigg(h \frac{ih - x_i}{2} +h \frac{h (i + 1) - x_i - h \frac{ih - x_i}{2}}{2}\Bigg) \end{align}

Con esas dos sucesiones se obtienen los resultados expuestos a continuación (con $h = \frac{1}{8}$):

$i$ $ih$ $x_i_{Taylor}$ $x_i_{RK2}$
0 0 1.0000 1.0000
1 0.125 0.9470 0.9509
2 0.250 0.9050 0.9124
3 0.375 0.8733 0.8838
4 0.500 0.8512 0.8644

chequear estos valores porque las expresiones estaban mal

Problema 2

Primero, necesitamos 3 ecuaciones para obtener nuestras 3 incógnitas, $y(0.25)$, $y(0.50)$ y $y(0.75)$. Desarrollamos Taylor hasta un orden que nos permita luego obtener una fórmula centrada de orden de error $h^2$:

(4)
\begin{array} {l} y(x+h) = y(x) + hy'(x) + \frac{h^2}{2}y''(x) + \frac{h^3}{6}y'''(x) + O(h^4) \\ y(x-h) = y(x) - hy'(x) +\frac{h^2}{2}y''(x) - \frac{h^3}{6}y'''(x) + O(h^4) \\ \end{array}

Sumamos las ecuaciones para obtener:

(5)
\begin{array} {l} y(x+h) + y(x-h) = 2y(x) + h^2 y''(x) + O(h^4) \\ y''(x) = \frac{y(x-h) - 2y(x) + y(x+h)}{h^2} + O(h^2) \\ \end{array}

Reemplazando en la ecuación dada:

(6)
\begin{array} {l} y"(x) + 4y(x) = x + 1 \\ \frac{y(x-h) - 2y(x) + y(x+h)}{h^2} + 4y(x) = x + 1 \\ \end{array}

Nuestro $h = \frac{1}{4}$, porque pasamos de 0.25 a 0.5 a 0.75. $h^2 = 0.0625$ por lo que la se multiplica por un factor de 16.

(7)
\begin{array} {l} 16y(x-h) - 32y(x) + 16y(x+h) + 4y(x) = x + 1 \\ 16y(x-h) - 28y(x) + 16y(x+h) = x + 1 \\ 16y(x-0.25) - 28y(x) + 16y(x+0.25) = x + 1 \\ \end{array}

Ahora, hago 3 ecuaciones, con $x=0.25$, $x=0.50$ y $x=0.75$:

(8)
\begin{align} \left\{ \begin{array}{l} 16y(0) - 28y(0.25) + 16y(0.50) = 1.25 \\ 16y(0.25) - 28y(0.50) + 16y(0.75) = 1.5 \\ 16y(0.50) - 28y(0.75) + 16y(1) = 1.75 \\ \end{array}\right \end{align}

Sabemos que $y(0) = y(1) = 0$, entonces hacemos una matriz de ecuaciones con nuestras 3 incógnitas

(9)
\begin{align} \left[ \begin{array} {ccc|c} -28 & 16 & 0 & 1.25 \\ 16 & -28 & 16 & 1.5 \\ 0 & 16 & -28 & 1.75 \\ \end{array} \right] \end{align}

Ahora usando "Gauss Elimination Method" (GEM) resolvemos la matriz.

(10)
\begin{align} \left[ \begin{array} {ccc|c} -28 & 16 & 0 & 1.25 \\ 0 & -18.85 & 16 & 2.21 \\ 0 & 16 & -28 & 1.75 \\ \end{array} \right] \end{align}

Pivoteo, porque me quedó un cero…

(11)
\begin{align} \left[ \begin{array} {ccc|c} -28 & 16 & 0 & 1.25 \\ 0 & -18.85 & 16 & 2.21 \\ 0 & 0 & -14.43 & 3.628 \\ \end{array} \right] \end{align}

Ya tenemos una matriz triangular superior, y con 4 cifras decimales significativas queda:

(12)
\begin{align} \left\{ \begin{array}{l} y(0.75) = -0.2516 \\ y(0.50) = -0.3309 \\ y(0.25) = -0.2337 \\ \end{array}\right \end{align}

Teoría 1

Parte a

Corrección al enunciado:

m1 = delta u

Debemos determinar qué ODE de orden 2 resuelve esta implementación. Para ello hay que reconocer la forma de RK2 para cada una de las dos variables que se encuentran aquí, por el proceso de renombramiento.

u = u + 0.5 (l1 + l2)
v = v + 0.5 (m1 + m2)

De allí podemos distinguir que las variables m1, m2, l1 y l2 corresponden a los $k_1$ y $k_2$ del método RK2. En particular, l1 y l2 corresponden a la variable u, y m1 y m2 corresponden a la variable v. Para poder determinar la ecuación diferencial que se resuelve sólo observaremos los $k_1$, cuyo valor siempre es:

(13)
\begin{equation} k_1 = h f(t, x) \end{equation}

La variable m1 es la $k$ asociada a v. Con el siguiente fragmento de código:

m1 = delta u

podemos determinar que delta representa $h$ y u representa la derivada de v, es decir:(14)
\begin{equation} u = v' \end{equation}

La variable l1 es la $k$ asociada a u. Con el siguiente fragmento de código:

l1 = delta (-0.9 sin u - 0.24v)

podemos determinar que delta representa $h$ y $-0.9 \sin(u) - 0.24v$ es la derivada de u, es decir:(15)
\begin{align} u' = -0.9 \sin(u) - 0.24v \end{align}

Por las ecuaciones (14) y (15) podemos decir que la ecuación diferencial es:

(16)
\begin{align} v'' + 0.9 \sin(v') + 0.24v = 0 \end{align}

El paso de integración es delta, 0.05. El método es Runge-Kutta de orden 2.

Parte b

Parte c

Teoría 2

Parte a

Considerando que

(17)
\begin{align} X(\omega)=H(\omega)F(\omega) \end{align}

Tenemos por el enunciado que

(18)
\begin{align} H(\omega) = \frac{1}{6-\omega^2+4\omega i} = \frac{X(\omega)}{F(\omega)} \end{align}

Despejamos $F(\omega)$ y tenemos

(19)
\begin{align} F(\omega)=X(\omega)(6-\omega^2+4\omega i) \end{align}
(20)
\begin{align} F(\omega)=6X(\omega)-\omega^2X(\omega)+4\omega i X(\omega) \end{align}

En el término de $\omega^2$ dividimos y multiplicamos por $i^2$:

(21)
\begin{align} F(\omega)=6X(\omega)-\frac{i^2\omega^2}{i^2}+4\omega i X(\omega) \end{align}

El $i^2$ del denominador lo eliminamos y cambiamos el signo de ese término, y como $i^2\omega^2$ indica la transformada de la derivada segunda (análogamente para la derivada primera con $i\omega$) tenemos:

(22)
\begin{align} F(\omega) = 6X(\omega) + \omega^2 i^2 X(\omega) +4 \omega i X(\omega) \end{align}

Antitransformamos esto, y la ecuación diferencial queda:

(23)
\begin{equation} f(t) = 6x(t)+x''(t)+4x'(t) \end{equation}

Parte b

El código primero transforma los valores de la función de entrada para poder luego hacer la multiplicación de $H(w)$ y $F(w)$ para obtener el espectro de la salida, que luego antitransforma para obtenerlo en el tiempo. Esa multiplicación la efectúa dado que la transformada de la convolución en el tiempo es la multiplicación en la frecuencia. Estas son dos posibles demostraciones de esa propiedad:

(24)
\begin{align} \text{Producto de Convoluci\'on:}\quad \int_{-\infty}^{+\infty} g(z) H(t-z) dz = (g*H)(t)\\ \end{align}

Demostración 1

(25)
\begin{eqnarray} \mathcal{F}[(g*H)(t)] &= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \\ \mathcal{F}\Big[\int_{-\infty}^{+\infty} g(z) H(t-z) dz\Big] &= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \\ \int_{-\infty}^{+\infty}\Big[\int_{-\infty}^{+\infty} g(z) H(t-z) dz\Big] e^{-j\omega t}dt &= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \\ \int_{-\infty}^{+\infty}\Big[\int_{-\infty}^{+\infty} g(z) H(t-z) e^{-j\omega t} dt\Big] dz&= \mathcal{F}[g(t)]\mathcal{F}[H(t)] & \text{por Fubbini}\\ \int_{-\infty}^{+\infty}g(z)\Big[\int_{-\infty}^{+\infty} H(t-z) e^{-j\omega t} dt\Big] dz&= \mathcal{F}[g(t)]\mathcal{F}[H(t)]\\ \int_{-\infty}^{+\infty}g(z) H(\omega) e^{-j\omega z} dz&= \mathcal{F}[g(t)]\mathcal{F}[H(t)] & \text{por Transformada de un desplazamiento} \\ H(\omega) \int_{-\infty}^{+\infty}g(z) e^{-j\omega z} dz&= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \\ H(\omega) g(\omega) &= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \\ \mathcal{F}[g(t)]\mathcal{F}[H(t)] &= \mathcal{F}[g(t)]\mathcal{F}[H(t)] \end{eqnarray}

Demostración 2

Aplicando la definición de transformada de Fourier y de convolución:

(26)
\begin{align} \mathcal{F}\{h(t) * f(t)\} = \frac{1}{2\pi} \int_{-\infty}^{\infty} \left[ \int_{-\infty}^{\infty} h(x) f(t - x) dx \right] e^{-jwt} dt \end{align}

Haciendo que la integral interna sólo contenga todo lo relacionado con $t$, extrayendo a la externa todo lo relacionado con $x$ se obtiene:

(27)
\begin{align} \frac{1}{2\pi} \int_{-\infty}^{\infty} h(x) \left[ \int_{-\infty}^{\infty} f(t - x) e^{-jwt} dt \right] dx \end{align}

Ahora, con un cambio de variable $u = t-x$ y reacomodando el factor $\frac{1}{2\pi}$ se llega a:

(28)
\begin{align} \frac{1}{2\pi} \int_{-\infty}^{\infty} h(x) \left[ \int_{-\infty}^{\infty} f(u) e^{-jw(u+x)} dt \right] dx = \int_{-\infty}^{\infty} h(x) \left[ \frac{1}{2\pi} \int_{-\infty}^{\infty} f(u) e^{-jwu} dt \right] e^{-jwx} dx \end{align}

Allí se puede apreciar la forma de la transformada de fourier de $f$, por lo que:

(29)
\begin{align} F(w) \int_{-\infty}^{\infty} h(x) e^{-jwx} dx \end{align}

Basta ahora con multiplicar por $\frac{2\pi}{2\pi}$ para poder obtener la forma de $H(w)$:

(30)
\begin{align} 2\pi F(w) \frac{1}{2\pi} \int_{-\infty}^{\infty} h(x) e^{-jwx} dx \end{align}

Entonces queda:

(31)
\begin{align} \mathcal{F}\{h(t) * f(t)\} = 2\pi F(w) H(w) \end{align}

Nota: El algoritmo presentado en la figura del enunciado no incluye en su fórmula $\verb#X(k) = F(k) * H(k)#$ un factor de $2\pi$. Esto se debe a que quien lo escribió no considera que la transformada de Fourier incluya un factor $\frac{1}{2\pi}$.

Parte c

Del código se obtiene que:

(32)
\begin{align} f(t) = 0,3 e^{-25 t^2} \cos(2 \pi t) + 4,5 \sin (4 \pi t) \end{align}

Utilizando la expresión del seno y coseno como exponenciales se obtiene que

(33)
\begin{align} f(t) = \frac{0,3}{2} e^{-25 t^2} \left(e^{j 2 \pi t} + e^{-j 2 \pi t}\right)+ \frac{4,5}{2j} \left( e^{j 4 \pi t} - e^{-j 4 \pi t}\right) \end{align}

Para la primera parte (la expansión del coseno) podemos utilizar la ayuda que nos da llamando a $a = 25$

(34)
\begin{align} \mathcal{F}\left\{e^{-at^2}\right\} = \sqrt{\frac{\pi}{a}} e^{-\frac{w^2}{4a}} \end{align}

Para la segunda parte (la del seno) podemos pensar que

(35)
\begin{align} \mathcal{F}\left\{1\right\} = \delta(w) \end{align}

En ambos tenemos que aplicar la propiedad de desplazamiento en la frecuencia, $F(w - w_0) = \mathcal{F} \left[f(t) e^{jw_0 t}\right]$:

(36)
\begin{align} F(w) = \frac{0,3}{2} \left( \mathcal{F}\left\{e^{-25 t^2}\right\}(w - 2 \pi) + \mathcal{F}\left\{e^{-25 t^2}\right\}(w + 2 \pi )\right)+ \frac{4,5}{2j} \left( \mathcal{F}\left\{1\right\}(w - 4 \pi) - \mathcal{F}\left\{1\right\}(w + 4 \pi)\right) \end{align}

Remplazando

(37)
\begin{align} F(w) = \frac{0,3}{2}\sqrt{\frac{\pi}{25}} \left(e^{-\frac{(w-2\pi)^2}{100}} + e^{-\frac{(w+2\pi)^2}{100}}\right)+ \frac{4,5}{2j} \left( \delta(w-4\pi) - \delta(w+4\pi) \right) \end{align}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License