Se sabe que las ecuaciones diferenciales ordinarias permiten describir una gran variedad de problemas y fenómenos que surgen en las matemáticas, la física, la química, la economía, etc. Por lo cual hay un gran interés por encontrar soluciones a este tipo de problemas. Por lo tanto el objetivo de esta entrada es mostrar como encontrar soluciones aproximadas al problema de valor de frontera de segundo orden (PVF2) mediante diferencias finitas.
El problema que vamos a abordar en esta ocasión es el PVF2, definido como: \begin{equation}\left\{\begin{array}{l} y''=f(x,y,y'),\\ y(a)=\alpha,\\ y(b)=\beta. \end{array}\right. \end{equation} Lo anterior también se puede expresar de la siguiente forma \begin{equation}\label{eq02} y''(x)=p(x)y'(x)+q(x)y(x)+r(x),\,a\leq x\leq b,\, y(a)=\alpha,\, y(b)=\beta, \end{equation} en este caso se dice que el PVF2 es de tipo lineal, en caso contrario, se dice que el PVF2 es del tipo no lineal. Este tipo de problemas son muy comunes en la ciencias básicas o en ingeniería. Sin embargo, su principal dificultad radica en que la solución a dichos problemas no es sencilla. Debido a esto, resulta importante definir métodos numéricos para dar estimaciones numéricas de las soluciones. Existe una gran diversidad de métodos para hacer esto: el método del disparo, el método de elementos finitos, el método de Galerkin discontinuo, el método de diferencias finitas, etc. Aquí solo nos vamos a centrar en el método de diferencias finitas.
El método de diferencias finitas se basa en la idea de aproximar las derivadas que aparecen en el PVF2 mediante una partición del intervalo $[a, b]$. A menudo a estas particiones se le denominan mallas o grillas. Una vez que se hace estas aproximaciones, se sustituyen en la ecuación diferencial del PVF2 y posteriormente se resuelve un sistema de ecuaciones lineales o no lineales, esto depende de si el PVF2 es lineal o no.
El caso lineal
En el caso del problema $\eqref{eq02}$, se procede a discretizar el intervalo $[a,b]$ en una malla o grilla de puntos $a=x_{_{0}} < x_{_{1}} < \cdots < x_{_{n+1}}=b$, en la cual $x_{_{i+1}}-x_{_i}=h$, para todo $i=0,\dots, n$. En general los puntos no necesitan estar distribuidos de esta manera.
En los puntos de la forma $x_{i}=a+ih$ con $i=1,2,\dots,n$. Así la ecuación diferencial a aproximar es \begin{equation}\label{eq:03} y''(x_{_i})= p(x_{_{i}})y'(x_{_i})+q(x_{_i})y(x_{_i})+r(x_{_i}). \end{equation} Ahora empleando el polinomio de Taylor de orden 3 para $y(x)$, alrededor de $x_{_i}$, y se evalúa en $x_{_{i+1}}$ y $x_{_{i-1}}$, entonces, suponiendo que $y\in C^{4}[x_{_{i-1}},x_{_{i+1}}]$, se tiene \begin{equation}\label{eq:04}y(x_{_{i+1}})=y(x_{_i}+h)=y(x_{_i})+hy^{(1)}(x_{_i})+\frac{h^2}{2}y^{(2)}(x_{_i})+\frac{h^3}{6}y_{_i}^{(3)}(x_{_i})+\frac{h^4}{24}y^{(4)}(\xi_{_i}^{+}) \end{equation} para algún $\xi_{_i}^{+}$ en $(x_{_i},x_{_{i+1}})$, y \begin{equation}\label{eq:05}y(x_{_{i-1}})=y(x_{_i}-h)=y(x_{_i})-hy^{(1)}(x_{_i})+\frac{h^2}{2}y^{(2)}(x_{_i})-\frac{h^3}{6}y_{_i}^{(3)}(x_{_i})+\frac{h^4}{24}y^{(4)}(\xi_{_i}^{^-}) \end{equation} para algún $\xi^{^-}_{_i}$ en $(x_{_{i-1}},x_{_i})$. Si se suman adecuadamente $\eqref{eq:04}$ y $\eqref{eq:05}$ se obtiene \[y(x_{_{i+1}})+y(x_{_{i-1}})=2y(x_{_i})+h^2y^{(2)}(x_{_i})+\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right),\] y al despejar $y^{(2)}(x_{_i}),$ \[y^{(2)}(x_{_i})=\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)-\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right).\] Ahora, aplicando el teorema de valor intermedio, se tiene que \begin{equation}\label{eq:06} y^{(2)}(x_{_i})=\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)-\frac{h^4}{12}y^{(4)}(\xi_{_i}), \end{equation} para algún $\xi_{_i}\in (x_{_{i-1}},x_{_{i+1}})$. A este esquema se le llama fórmula de diferencias centradas para $y^{(2)}(x_{_i})$. Ahora si a $\eqref{eq:04}$ se le resta $\eqref{eq:05}$. De manera semejante se obtiene una fórmula de este tipo para $y^{(1)}(x_{_i})$, la cual viene dada por \begin{equation}\label{eq:07} y^{(1)}(x_{_i})=\frac{1}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{h^3}{3}y^{(3)}(\eta_{_i}) \end{equation} para algún $\eta_{_i}$ en $(x_{_{i-1}},x_{_{i+1}})$.
Ahora si se sustituyen $\eqref{eq:06}$ y $\eqref{eq:07}$ en $\eqref{eq:03}$, entonces se obtiene la igualdad \[\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)=\frac{p(x_{_i})}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{p(x{_i})h^3}{3}y^{(3)}(\eta_{_i})+\frac{h^4}{24}\left(y^{(4)}(\xi_{_i}^{+})+y^{(4)}(\xi_{_i}^{^-})\right)\] Por lo tanto hay un error del orden $\mathcal{O}(h^3)$, con $y_{_i}=y(x_{_i})$, $y'_{_i}=y'(x_{_i})$, $y''_{_i}=y'(x_{_i})$, $p_{_i}=p(x_{_i})$, $q_{_i}=q(x_{_i})$ y $r_{_i}=r(x_{_i})$, podemos escribir simplemente \[\frac{1}{h^2}\left(y_{_{i+1}}-2y_{_i}+y_{_{i-1}}\right)=\frac{p_{_i}}{2h}\left(y_{_{i+1}}-y_{_{i-1}}\right)+q_{_i}y_{_i}+r_{_i}+\mathcal{O}(h^{3})\] y por lo tanto de esta manera obtenemos el conjunto de ecuaciones \[\frac{1}{h^2}\left(y_{_{i+1}}-2y_{_i}+y_{_{i-1}}\right)-\frac{p_{_i}}{2h}\left(y_{_{i+1}}+y_{_{i-1}}\right)-q_{_i}y_{_i}=r_{_i},\] para cada $i=1,\dots, n$, junto con $y_{_0}=\alpha$ y $y_{_{n+1}}=\beta$. Reescribiendo la ecuación anterior, se tiene que \[-\left(1+\frac{h}{2}p_{_i}\right)y_{_{i-1}}+\left(2+h^2q_{_i}\right)y_{_i}-\left(1-\frac{h}{2}p_{_i}\right)y_{_{i+1}}=-h^2r_{_i}\] El sistema de ecuaciones resultante se puede expresar en la forma $Ay=b$, donde, \[A=\left[\begin{array}{cccccc} \gamma_{_1} & \delta_{_1} & 0 & \cdots & \cdots & 0 \\ \epsilon_{_2} & \gamma_{_2} & \delta_{_2} & \ddots & & \vdots \\ 0 & \epsilon_{_3} & \gamma_{_3} & \delta_{_3} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & & \ddots & \ddots & \ddots & \delta_{_{n-1}} \\ 0 & \dots & \dots & 0 & \epsilon_{_n} & \gamma_{_n} \\ \end{array}\right],\,\,\begin{array}{ll} \gamma_{_i}=2+h^2q_{_i}, & 1\leq i\leq n,\\ \delta_{_i}=-1+\frac{h}{2}p_{_i}, & 1\leq i\leq n-1,\\ \epsilon_{_i}=-1-\frac{h}{2}p_{_i}, & 2\leq i\leq n,\end{array} \] \[\hat{y}=\left[\begin{array}{c} y_{_1} \\ y_{_2} \\ \vdots \\ y_{_n} \end{array}\right] \mbox{ e } b=\left[\begin{array}{c} -h^2r_{_i}+(1+\frac{h}{2}p_{_i})y_{_0} \\ -h^2r_{_{2}} \\ -h^2r_{_{3}} \\ \vdots \\ -h^2r_{_{n-1}}\\ -h^2r_{_i}+(1+\frac{h}{2}p_{_i})y_{_{n+1}} \end{array}\right].\]
Para revolver este sistema lineal de ecuaciones recomendamos leer factorización LU y para la reconstrucción aproximada de la curva solución recomendamos leer sobre interpolación polinómica de Lagrange.
¿Bajo qué condiciones el sistema lineal $Ay=b$ tiene única solución? La respuesta nos la da el siguiente teorema.
Teorema. Suponga que $p(x)$, $q(x)$ y $r(x)$ son continuas en $[a, b]$. Si $q(x)\geq 0$ en $[a,b]$, entonces el sistema triagonal $Ay=b$ tiene solución única siempre y cuando $h\leq \frac{2}{L}$, donde $L=\max_{a\leq x\leq b}|p(x)|$.
En cuando a los errores el siguiente teorema nos ayuda a deducir las cotas para el error de la aproximación de las ecuaciones diferenciales mediante este método de diferencias finitas de segundo orden.
Teorema Sea $Q^*$ tal que $q(x)>Q^*>0$, para $a\leq x\leq b$, $L=\max_{a\leq x\leq b}|p(x)|$, $M_{_3}=\max_{a\leq x\leq b} |y^{(3)}(x)|$ y $M_{_4}=\max_{a\leq x\leq b}|y^{(4)}(x)|$. Entonces, para $i=0,1,\dots, n+1$, se cumple: \[|\hat{y}_{_i}-y(x_{_i})|\leq h^{3}\left(\frac{M_{_4}+2LM_{_3}}{12Q^*}\right).\]
Este teorema nos dice que la solución de diferencias finitas converge a la solución exacta cuando $h\to 0$ y el error es $\mathcal{O}(h^{3})$. Cuando en la ecuación diferencial lineal, $p(x)=0$, entonces se puede demostrar que el error cometido es del orden de $\mathcal{O}(h^4)$.
El caso no lineal
Para el problema no lineal general con valor de frontera: \[y''=f(x, y,y'),\,a\leq x\leq b,\, y(a)=\alpha,\, y(b)=\beta,\] el método de diferencias finitas es similar al caso lineal. Sin embargo, el sistema de ecuaciones no será lineal y, por lo tanto, se requiere un proceso iterativo para resolverlo. Nuevamente, se procede a discretizar el intervalo $[a,b]$ en una partición de puntos $a=x_{_0}\leq x_{_1}\leq\cdots\leq x_{_{n+1}}=b$, tal que $x_{_{i+1}}-x_{_i}=h$.
En los puntos interiores de la partición $x_{_i}=a+ih$, $i=1,2,\dots,n$, la ecuación diferencial a aproximar es: \begin{equation}\label{eq:08} y''(x_{_i})=f(x_{_i},y(x_{_i}),y'(x_{_i})). \end{equation} Luego sustituyendo $\eqref{eq:04}$ y $\eqref{eq:05}$ en $\eqref{eq:08}$ se obtiene: \[\frac{1}{h^2}\left(y(x_{_{i+1}})-2y(x_{_i})+y(x_{_{i-1}})\right)=f\left(x,y(x_{_i}),\frac{1}{2h}\left(y(x_{_{i+1}})-y(x_{_{i-1}})\right)-\frac{h^3}{3}y^{(3)}(\eta_{_i})\right)+\frac{h^4}{12}y^{(4)}(\xi_{_i})\] para algún $\xi_{_i}$ y $\eta_{_i}$ en el intervalo $(x_{_{i-1}},x_{_{i+1}})$, con $i=1,2,\dots,n$. El método de diferencias finitas se obtiene empleando esta ecuación junto con las condiciones finitas de frontera $y(a)=\alpha$ y $y(b)=\beta$. Por lo tanto si nuevamente se hace $y_{_i}=y(x_{_i})$, $y'_{_i}=y'(x_{_i})$, $y''_{_i}=y'(x_{_i})$, $p_{_i}=p(x_{_i})$, $q_{_i}=q(x_{_i})$ y $r_{_i}=r(x_{_i})$ se obtiene el sistema de ecuaciones \[\frac{y_{_{i+1}}-2y_{_i}+y_{_{i-1}}}{h^2}-f\left(x_{_i},y_{_i},\frac{y_{_{i+1}}-y_{_{i-1}}}{2h}\right)=0,\] para cada $i=1,2,\dots,n$ junto con $y_{_0}=\alpha$ y $y_{_{n+1}}=\beta$. El sistema no lineal de $n\times n$ obtenido con este método $F_{_i}(y_{_1},\dots,y_{_n})=0$ con \[F_{_i}(y_{_1},\dots,y_{_n})=y_{_{i+1}}-2y_{_i}+y_{_{i-1}}-h^2f\left(x_{_i},y_{_i},\frac{y_{_{i+1}}-y_{_{i-1}}}{2h}\right)\] tiene solución única si $h<\frac{2}{L}$, donde $L$ es tal que $\left|\frac{\partial f}{\partial y'}(x,y,y')\right|\leq L$.
Para resolver este sistema de ecuaciones usaremos el método de Newton. Es consiste es suponer una solución $y^{(0)}=(y^{(0)}_{_1},\dots, y^{(0)}_{_{n}})^{\top}$, la iteración con este método es: \[y^{(j+1)}=y^{(j)}-\pmb{J}^{-1}(y^{(j)})F(y^{(j)})\] donde $\pmb{J}(y^{(j)})$ es la matriz jacobiana. En la práctica es más conveniente a efectos computacionales iterar resolviendo las ecuaciones: \[\pmb{J}(y^{(j)})(y^{(j+1)}-y^{(j)})=-F(y^{(j)}).\] En cuanto al jacobiano $\pmb{J}_{rs}(y^{(j)})=\frac{\partial F_{_r}}{\partial y_{_s}}$, se tiene: \[J_{_{kk}}=-2-h^{2}\frac{\partial f}{\partial y}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] \[J_{_{k(k+1)}}=1-\frac{h}{2}\frac{\partial f}{\partial y'}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] \[J_{_{(k+1)k}}=1+\frac{h}{2}\frac{\partial f}{\partial y'}\left(x_{_k},y_{_k},\frac{y_{_{k+1}}-y_{_{k-1}}}{2h}\right),\] con $k=1,\dots, n$. El resto de las entradas de la matris $\pmb{J}$ son nulos. En particular, como todo los elementos $|r-s|>1$ son nunos, entonces $J$ es tridiagonal. En cada iteración del sistema lineal: \[\pmb{J}(y^{(j)})v=-F(y^{(j)}),\] con $v=(v_{_1},\dots,v_{_n})^{\top}$, se requiere obtener $v$, porque \[y^{(j+1)}=y^{(j)}+v.\] Puesto que $\pmb{J}$ es triagonal, se puede utilizar métodos para resolver ecuaciones lineales. Nosotros, para efecto de implementación, en un próximo post usaremos la factorización LU para resolver este tipo de problemas. Para finalizar cabe mencionar que nuevamente la solución que se encuentra con este método se aproxima a la solución exacta con un orden de aproximación del orden $\mathcal{O}(h^{3})$.
Referencias
- Carl D. Meyer. Matrix Analysis and Aplied Linear Algebra. SIAM, 2000.
- Gabriel Nagy . Linear Algebra. Mathematics Department, Michigan State University, 2012.
- Gilbert Strang. MIT Open Course Ware | Free Online Course Materials.
- Richard L. Burden y J. Douglas Faires. Análisis numérico. Thomson Learning, 2002.
- D. N. Arnold. And interior penalty finite element method with discontinuoous elements. SIAM Journal on Numerial Analysis, 1982.
0 comentarios:
Publicar un comentario