## Métodos de diferencias finitas para problemas de valores a la frontera
 
Es claro que el método de shooting es **muy lento**, pues el tener que resolver toda la ecuación diferencial para evaluar $F$ y poder resolver $F(z^*) = 0$ 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. 

Así resulta conveniente explorar otro tipo de métodos de solución. Supongamos que en nuetro problema de valores a la frontera estamos trabajando con una  **ecuación diferencial lineal de orden dos**, es decir, una problema de la 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.

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 (i.e. $g(t_i) = g_i$), 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)$ y $1\leq i \leq n-1$

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$, y ya que conocemos quienes son las funciones $a(t)$,$b(t)$, $c(t)$ y $f(t)$ así como los valores $h$, $\alpha$ y $\beta$, entonces **las únicas incógnitas de las ecuaciones son las $y_0,y_1,\ldots,y_n$**

Como todos los demás términos son conocidos, 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  que vimos en las clases de Álgebra Lineal:  eliminación Gaussiana, Jacobi, Gauss-Siedel. A esta forma de resolver ecuacioned diferenciales de segundo orden aproximando las derivadas con diferencias finitas y obteniendo un problema de ecuaciones linealesse les conoce como **métodos de diferencias finitas**

### Ejercicio 1

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

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

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 2

Utiliza tu función `matriz` y alguno de los métodos de las clases de Álgebra lineal para resolver las ecuaciones del ejercicio 3 de la clase 17.2. 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.

## Problemas tipo Sturm-Liouville

Un problema parecido pero fundamentalmente **distinto** a los problemas tradicionales de valores a la frontera son los problemas tipo **Sturm-Liouville**. En dichos problemas, tenemos una ecuación diferencial de forma similar a la ecuación diferencial lineal de orden dos:

$$
\begin{align}
a(t) \ddot{y}(t) + b(t) \dot{y}(t) &+  c(t) y(t) = \lambda y(t) \\
\end{align}
$$

Con $\lambda \in \mathbb{R}$ una constante **desconocida**. Con condiciones iniciales que son de una forma más complicada:

$$
\gamma_1 y(t_0) + \gamma_2 \dot{y}(t_0)  = 0 \\
\delta_1 y(t_f) + \delta_2 \dot{y}(t_f)  = 0 \\
$$

La [teoría de Sturm-Liouville](https://en.wikipedia.org/wiki/Sturm%E2%80%93Liouville_theory) nos garantiza que, bajo ciertas condiciones, existe una infinidad de funciones $y_{(k)}$ y  constantes $\lambda_k$ asociadas a cada una, con $k \in \mathbb{N}$, que satisfacen la ecuación diferencial mostrada, es decir

$$
\begin{align}
a(t) \ddot{y}_{(k)}(t) + b(t) \dot{y}_{(k)}(t) &+  c(t) y_{(k)}(t) = \lambda_k y_{(k)}(t) \\
\end{align}
$$

Por simplicidad, supongamos que las condiciones de frontera son

$$
y(t_0) = 0 \\
y(t_f) = 0 \\
$$


Si volvemos a utilizar diferencias finitas para discretizar las derivadas, de forma totalmente análoga a la sección anterior, podremos entonces expresar la ecuación diferencial 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} = \lambda y_i \\
& \\
y_0 &= 0\\
y_n &= 0 \\
\end{align}
$$

Utilizando los valores ya conocidos de $y_0 = 0$ y de $y_n = 0$, podemos reescribir las ecuaciones para $i=1,n$. Así, las ecuaciones para cada $y_i$ son


$$
\begin{gather}
\left(c_1 - \frac{2a_1}{h^2} \right) y_{1} +  \left(\frac{a_1}{h^2} + \frac{b_1}{2h}\right) y_{2} = \lambda y_1 \\
 \\
\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} = \lambda y_i \quad i = 2,\ldots,n-2 \\
 \\
\left(\frac{a_{n-1}}{h^2} - \frac{b_{n-1}}{2h}\right) y_{n-2} + \left(c_{n-1} - \frac{2a_{n-1}}{h^2} \right) y_{n-1} = \lambda y_{n-1} \\
\end{gather}
$$

Estas ecuaciones se pueden escribir de forma compacta como un conjunto de **$n-1$ ecuaciones lineales** de la forma

$$
\mathbf{M y } = \lambda \mathbf{y}
$$

Con $\mathbf{M}$ otra matriz tridiagonal. Es claro que los vectores $\mathbf{y}$ y las constantes $\lambda$ deben de ser los **valores y vectores propios** de la matriz $\mathbf{M}$. Así, al aproximar las derivadas con diferencias finitas convertimos el problema de Sturm-Liouville en un problema de valores propios para la matriz $\mathbf{M}$. Los vectores propios de la matriz corresponderán a la discretización de las funcions $y_{(k)}$ mientras que sus valores propios asociados serán los $\lambda_k$

### Ejercicio 3

Escribe una función `sturmLiouville(a,b,c,t_0,t_f,n)`. Los argumentos de la función deben de ser

* `a`,`b`,`c` las funciones **escalares** $a(t)$, $b(t)$ y $c(t)$ del problema de Sturm-Liouville

* `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 $\mathbf{M}$, de $n-1 \times n-1$ que correspondan a la discretización del problema de Sturm-Liouville expuesto anteriormente.

**Sugerencia:** puedes ver a la matriz $\mathbf{M}$ como un caso particular de la matriz $\mathbf{A}$. Así, podrías simplemente usar la función `matriz` del ejercicio 1 como base.

### Ejercicio 4

Utiliza la función `sturmLiouville` y las funciones `eigvals` y  `eigvecs` de la librería LinearAlgebra para resolver los siguientes problemas:

* $$\ddot{y} + 3\dot{y} + 2y = -\lambda y \\ y(0) = y(1) = 0 $$

Solución: $y_{(k)}(t) = e^{-3t/2}\sin{(k \pi t)}$, $\lambda_k = (1 + 4k^2 \pi^2)/4$

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

Solución: $y_{(k)}(t) = \sin{\left(\frac{k \pi}{\ln{2}} \ln{(x)} \right)}$, $\lambda_k = \left(\frac{k \pi}{\ln{2}} \right)^2$

Compara, en la misma figura, la función obtenida numéricamente y la solución analítica para valores bajos de $|\lambda_k|$ (i.e., grafica las funciones asociadas a los valores propios **más pequeños en valor absoluto** ). ¿Qué ocurre para valores altos de $|\lambda_k|$? ¿Los valores propios numéricos y analaíticos se parecen?

**Nota:** toma en cuenta que los signos de los valores propios $\lambda$ aquí mostrados son **negativos**, por lo que deberás tener cuidado a la hora de comparar.

## Aplicación en 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 de 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

Resuelve el problema tipo Sturm-Liouville asociado a la particula en el pozo infinito de potencial. Compara las soluciones obtenidas con las soluciones analíticas. ¿Las funciones obtenidas numéricamente se pueden comparar bien con la solución analítica? Justifica tu respuesta

**Nota:** las soluciones analíticas de este problema **están normalizadas**, es decir, cumplen que 

$$
\int_{0}^1 | \psi_k (x) |^2 \: dx = 1
$$

Lo más probable es que la soluciones que obtenas numéricamente **no estén normalizadas**, por lo que tendrás que encontrar una manera equivalente de normalizarlas para que se puedan comparar.

## 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 del sistema.

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

```

b8e5837f6abc7352398e6029f8719828d9904c13
ipynb
Jupyter Notebook
files/fiscomp_2021-1/material/clase18.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2021-1/material/clase18.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2021-1/material/clase18.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
Qwen/Qwen-72B
1. YES
2. YES