# EDOs de segundo orden: problemas de valores a la frontera

Todos los métodos vistos en las tres clases anteriores pretendían resolver **problemas de valor inicial**  en los que conocíamos el valor de la función $\mathbf{x}(t_0)$ y buscábamos conocer $\mathbf{x}(t)$ para $t \geq t_0$.

Sin embargo, muchas veces nos enfrentamos a problemas en los que no se conocen en todas las variables en el mismo tiempo, si no en en dos tiempos: los extremos del intervalo $[t_0,t_f]$. Nos interesa conocer $\mathbf{x}(t)$ en ese intervalo. En ese caso, decimos que estamos tratando con un problema de **valores a la frontera**.

## Preludio Matemático: Problema con condiciones en varios puntos

Existe todo un formalismo para manejar esta clase de problemas y para describirlos. Por simplicidad, nosotros nos concentraremos en problemas de sistemas de dos ecuaciones diferenciales de primer orden, obtenidos a partir de una **ecuación diferencial lineal de segundo orden**, que tiene la siguiente forma

$$
\begin{align}
a(t) \ddot{y}(t) + b(t) \dot{y}(t) &+  c(t) y(t) = f(t) \\
& \\
y(t_0) &= \alpha \\
y(t_f) &= \beta \\
\end{align}
$$

Con $a(t) \neq 0$, $b(t)$, $c
(t)$ y $f(t)$ funciones conocidas, $\alpha$ y $\beta$ constantes.

## Métodos de Shooting

### Planteamiento

Usando el truco para convertir la ecuación de segundo orden en un sistema de ecuaciones diferenciales ordinarias, sea $x_1(t) = y(t)$, $x_2(t) = \dot{y}(t)$. Así, el problema se convierte en resolver el sistema:


$$
\begin{align}
\dot{x_1} &= \dot{y} =  x_2 \\
\dot{x_2} &= \ddot{y} =  \frac{ f(t) -  b(t) x_2(t) -  c(t) x_1(t)  }{a(t)} \\
\end{align}
$$

Y las condiciones se vuelven

$$
\begin{align}
x_1(t_0) &= \alpha \\
x_1(t_f) &= \beta \\
\end{align}
$$

Para simplificar la notación, denotemos 

$$
g(x_1,x_2,t) = \frac{ f(t) -  b(t) x_2(t) -  c(t) x_1(t)  }{a_2(t)} \\
$$

El sistema obtenido es muy parecido a los problemas de valores iniciales que ya sabemos resolver. La única diferencia es que no conocemos $x_2(t_0)$, si no $x_1(t_f)$. Sin embargo, **supogamos que si lo conozco**: sea $x_2(t_0)=z$ con $z$ una constante arbitraria. Entonces obtengo el sistema:

$$
\begin{align}
\dot{x_1} &= x_2 \\
\dot{x_2} &= g(x_1,x_2,t)  \\
& \\
x_1(t_0) &= \alpha \\
x_2(t_0) &= z \\
\end{align}
$$

El sistema ya cumple una de las condiciones de frontera requeridas, que $x_1(t_0) = \alpha$. Sin embargo, no sabemos si la otra condición original ($x_1(t_f)=\beta$) se cumple. Si ahora utilizamos uno de los métodos (no adaptativos, con paso de tiempo constante) para problemas de valores iniciales a dicho sistema, desde el tiempo $t_0$ hasta el tiempo $t_f$, usando un paso $h$ que **divide excatamente** al intervalo $[t_0,t_f]$, es decir, que cumple

$$
\frac{t_f - t_0}{h} = n \; \in \mathbb{N}
$$

los últimos valores obtenidos por el método corresponderán a  los valores $x_1(t_f)$ y $x_2(t_f)$. 

Es claro que dichos valores dependerán de la condición incial $x_2(t_0)=z$. Así, denotemos al valor de $x_1(t_f)$ obtenido utilizando la condición incial $x_2(t_0)=z$ como $x_1(t_f ; z)$. Sabemos que dicho valor **segurmante no va a coincidir** con el valor real $\beta$. Así, podemos definir una función de diferencia $F(z)$, definida como

$$
F(z) = x_1(t_f ; z) - \beta
$$

Es claro que si existe un valor $z^*$ tal que $F(z^*) =0$ , entonces para ese $z^*$ la solución al problema de valor incial cumplirá la condición de frontera que queremos, por lo que para ese valor $z$ el problema de valor inicial se vuelve completamente equivalente al problema de valores a la frontera, por lo que **sus soluciones $x_{1}(t)$, $x_2(t)$ corresponden a las mismas funciones**

### ¿Cómo encontramos el valor $z^*$?

Nuevamente, notemos que el problema de encontrar $z^*$ tal que $F(z^*) =0$ corresponde a **resolver una ecuación trascendental** $F(z^*) =0$. Ya estudiamos previamente los métodos para resolver dicha ecuación, por lo que podemos utilizar cualquiera de ellos, pero utilizaremos en este caso el de Newton-Rhapson.

Así, los **métodos de Shooting** (*shooting* por que, al igual que con un cañon, vamos haciendo pruebas desde un punto inicial para ver si acertamos en el punto final) consisten en resolver la ecuación trascendental $F(z^*) =0$ y luego utilizar el valor $z^*$ para resolver el problema de valores iniciales que es equivalente al problema de valores a la frontera. 

Noten que, por la definición de $F(z)$, para calcular el valor de la función en un punto **es necesario resolver numéricamente todo el sistema de ecuaciones diferenciales** con condición inicial $x_2(t_0)=z$. Esto va a volver  muy lento el procedimiento de resolver la ecuación trascendental.

### Ejercicio 1 y 2

Implementa una función `shootingRK4(alpha,beta,g,t_0,t_f,n,epsilon)` que tome como argumentos:

* `alpha`, `beta`: las condiciones a la frontera $\alpha$, $\beta$ del problema

* `g` la función **vectorial** $g(x_1,x_2,t)$

* `t_0`,`t_f` los tiempos inicial y final, `n` $>2$ un natural correspondiente el número de subintervalos que se van a tomar entre `t_0` y `t_f`. Eso **define** el paso de tiempo $h$ del sistema a resolver como

$$
\frac{t_f - t_0}{n} = h
$$

* `epsilon` $>0$, un  número real arbitrario que será nuestra tolerancia para el método de Newton-Rhapson.

La función va a implementar el método de Shooting resolviendo el sistema con la función `RK4` de la clase 16. La función debería proceder de la siguiente forma:

1. Primero definimos la función $F(z)$ que utilizando la función `RK4`, obtiene $x_1(t_f;z)$ y calcula $x_1(t_f;z) - \beta$

2. Ya definida $F(z)$, usar la función `miNewtonTol` de la clase 08 para aplicarle el método de Newton-Rhapson a la función $F$, utilizando una diferencia finita centrada y una tolerancia `epsilon`. Esto nos dará el valor $z^*$

3. Ya obtenido el valor de $z^*$, se debe de volver a utilizar `RK4` con la condición inicial $x_2(t_0) = z^*$ para resolver el sistema de ecuaciones.

Tu función debe de regresar el arreglo de tiempos `[t_0,...,t_n]` y las soluciones `[x_0,..,x_n]`.

### Ejercicio 3

Usa tu función `shootingRK4` para resolver los siguientes problemas de valores a la frontera

1. $\ddot{y} + 4y =0$, $y(0) = -2$, $y(\pi/4) = 10$

Respuesta: $y(t) = -2 \cos{(2t)} + 10 \sin{(2t)}$

2.  $t^2 \ddot{y} - 2 t \dot{y} + 2y =0$, $y(1) = -1$, $y(2)=1$

Respuesta: $y(t) = -\frac{5}{2}x + \frac{3}{2} x^2$

3. $\ddot{y} + 9y = \cos{x}$, $y(0) = 1$, $y(\pi/2)=-5/3$

Respuesta: $y(t) = \frac{1}{8}\cos{t} + \frac{7}{8}\cos{(3t)} + \frac{5}{3} \sin{(3t)}$

En una misma figura, compara tus soluciones con la solución analítica.

**Sugerencia** Ya que el método de Newton-Raphson necesita un valor inicial para $z$, puedes ayudarte de la solución analítica para proponer un valor de $z$ inicial cercano al real para que tu método de Newton no necesite demasiadas iteraciones.

## Métodos de diferencias finitas

Es claro que el método de shooting es **muy lento**, pues el tener que resolver toda la ecuación diferencial para evaluar $F$ puede ser muy tardado, además de que la solución obtenida siempre es aproximada, por lo que puede haber errores ahí que nos alejen de la condición verdadera y que nos den una mala solución.

Una alternativa para resolver problemas de valores a la frontera es utilizar diferencias finitas para aproximar las derivadas $\ddot{y}$ y $\dot{y}$. Utilizando una diferencia centrada para ambas derivadas, llegamos a la siguiente aproximación:

$$
\begin{align}
a(t) \frac{y(t+h) + y(t-h) - 2y(t)}{h^2} &+ b(t) \frac{y(t+h) - y(t-h)}{2h} +  c(t) y(t) = f(t) \\
& \\
y(t_0) &= \alpha \\
y(t_f) &= \beta \\
\end{align}
$$

Así, para un a partición regular de orden $n+1$ (i.e con $n$ subintervalos) del intervalo $[t_0,t_f]$, obtenemos los tiempos $t_0,\ldots,t_n$ que corresponden a un paso  $h = (t_f - t_0)/n$. Utilizando la notación de índices para el tiempo, el problema se convierte en:

$$
\begin{align}
a_i \frac{y_{i+1} + y_{i-1} - 2y_i}{h^2} &+ b_i \frac{y_{i+1} - y_{i-1}}{2h} +  c_i y_i = f_i \\
& \\
y_0 &= \alpha \\
y_n &= \beta \\
\end{align}
$$

Con $a_i = a(t_i)$, $b_i=b(t_i)$, $c_i = c(t_i)$, $f_i = f(t_i)$, $y_{i} = y(t_i)$

Podemos reordenar las ecuaciones para que quede de la siguiente forma:

$$
\begin{align}
\left(\frac{a_i}{h^2} - \frac{b_i}{2h}\right) y_{i-1} + \left(c_i - \frac{2a_i}{h^2} \right)& y_{i} +  \left(\frac{a_i}{h^2} + \frac{b_i}{2h}\right) y_{i+1} = f_i \\
& \\
y_0 &= \alpha \\
y_n &= \beta \\
\end{align}
$$

Si consideramos que la primera ecuación solo se cumple para $1 \leq i \leq n-1$, entonces el conjunto de ecuaciones representa un **sistema de $n+1$ ecuaciones lineales** de la forma $\mathbf{Ay} = \mathbf{b}$, para el que la solución $\mathbf{y} = (y_0,y_1,..,y_n)$ es la discretización de la función solución $y(t)$ a nuestra ecuación diferencial.

Dicho sistema puede resolverse utilizando los métodos (eliminación Gaussiana, Jacobi, Gauss-Siedel) que vimos en las clases de Álgebra Lineal.

### Ejercicio 4

Implementa una función `matriz(a,b,c,f,alpha,beta,t_0,t_f,n)`. Los argumentos de la función deben de ser

* `a`,`b`,`c`,`f` las funciones **escalares** $a(t)$,$b(t)$,$c(t)$ y $f(t)$

* `alpha`, `beta`: las condiciones a la frontera $\alpha$, $\beta$ del problema

* `t_0`,`t_f` los tiempos inicial y final, `n` $>2$ un natural correspondiente el número de subintervalos que se van a tomar. Eso define el paso de tiempo $h$ del sistema a resolver como

$$
\frac{t_f - t_0}{n} = h
$$

La función debe de regresar la matriz $A$, de $n+1 \times n+1$ y un vector $b$ de longitud $n+1$ que correspondan a la discretización del problema de valores a la frontera expuesto.

**Sugerencia:** la matriz $A$ debe de ser **tridiagonal**, es decir, los únicos elementos distintos de cero son elementos de la forma $A_{i,i}$, $A_{i+1,i}$ y $A_{i-1,i}$

### Ejercicio 5

Utiliza tu función `matriz` y alguno de los métodos de las clases de Álgebra lineal para resolver las ecuaciones del ejercicio 3. Compara tus soluciones con la solución analítica y con las soluciones obtenidas por el método de Shooting. ¿Es mejor este método o el método de Shooting?

**Sugerencia** si tienes dudas sobre si tu método de álgebra lineal converge para esa matriz, usa el operador `\` para apoyarte.

## Mecánica cuántica: la ecuación de Schrödinger independiente del tiempo

Un problema de valores a la frontera que surge comúnmente en física es la solución de la ecuación de Schrödinger dependiente del tiempo

La evolución en el tiempo de un sistema cuántico está definida por la ecuación de Schrödinger dependiente del tiempo:

$$
i \hbar \frac{\partial \psi}{\partial t} = H \psi
$$

Con $i$ la unidad imaginaria, $\hbar$ la constante reducida de Planck, $\psi$ la amplitud de probabildad del sistema y $H$ el hamiltoniano. La solución general a la ecuación es

$$
\psi = \sum_{k=1}^{\infty}A_k e^{-E_k t} \psi_k 
$$

Con $\psi_k$ y $E_k$ las funciones propias y valores propios, respectivamente, del operador hamiltoniano, es decir, las que cumplen con la propiedad

$$
H \psi_k = E_k \psi_k
$$

A esa ecuación se le llama la **ecuación de Schrödinger independiente del tiempo**. Trabajando un sistema d una partícula en una dimensión, cuya energía potencial es un pozo infinito de la forma

$$
V(x) = \begin{cases}
0 & \text{si } 0 \leq x \leq 1 \\
\infty & \text{en otro caso} 
\end{cases}
$$

Escribiendo el hamiltoniano de dicho sistema cuántico y haciendo una adimensionalización, la ecuación de Schrödinger independiente del tiempo se reduce al siguiente problema de valores iniciales

$$
\frac{1}{2} \frac{d^2 \psi_k}{dx^2}(x) + \epsilon_k \psi_k (x) = 0
 \\
\psi_{k} (0) = 0 \\
\psi_{k} (1) = 0 \\
$$

La solución analítica está dada por 
$$
\psi_k(x) = \sqrt{2} \sin{(k\pi x)} \\
\epsilon_k = \frac{k^2 \pi^2}{2} \\
$$

### Ejercicio 6

Utiliza los valores de la solución analítica para $\epsilon_k$ y resuelve, utilizando el método de tu preferencia (Shooting o diferencias finitas), el problema de valor incial para obtener $\psi_k$. Compara tu solución con la solución analítica. 

## La visión general

Resolver ecuaciones diferenciales ordinarias en la computadora es un problema fundamental de la física computacional por su relación íntima con mecánica clásica. No todos los métodos pueden resolver cualquier sistema, por lo que es importante buscar métodos con alta precisión para no tener errores como el del péndulo simple con el método de Euler

Existe un area dedicada a modelar los átomos de manera clásica. Así, resolver un sistema **termodinámico** (moleculas de aire en un cuarto, un fluido ideal) o **biológicos** (átomos en una proteína, moléculas en una membrana,etc) se convierte en resolver las ecuaciones de movimiento de todos los átomos del sistema. En 3 dimensiones, si tenemos $n$ átomos, habrá un sistema de $6n$ ecuaciones difernciales de primer orden que debemos de resolver para saber la evolución.

A esas simulaciones se les llama **dinámica molecular**. La dinámica molecular es una herramienta **revolucionaria** de la ciencia, pues permite resolver muchisísimos problemas. Hoy en día, se utlizan simulaciones de dinámica molecular para probar la eficiencia de fármacos y segurmante se están utilizando dichas simulaciones para probar y crear vacunas para el COVID19


```julia

```
