En el capítulo de valores iniciales en ecuaciones diferenciales ordinarias dada una ecuación diferencial del tipo \(y'=f(t,y(t))\) donde \(a\leq t\leq b\), hallábamos una solución aproximada al problema anterior dada una condición inicial del tipo \(y(a)=y_0\).

Incluso generalizábamos el problema anterior a ecuaciones diferenciales de órden superior, es decir, dada una ecuación diferencial de orden \(m\), \(y^{(m)}=f(t,y,y',\ldots,y^{(m-1)})\), donde \(a\leq t\leq b\), hallábamos una solución aproximada al problema anterior dada las condiciones iniciales del tipo \(y(a)=y_0,\ y'(a)=y'_0,\ y''(a)=y''_0,\ldots, y^{(m-1)}(a)=y^{(m-1)}_0\).

Fijémonos que las condiciones iniciales se imponen sólo en el extremo izquierdo del intervalo \([a,b]\) de definición de la variable independiente \(t\). No hay ninguna condición sobre el extremo derecho de dicho intervalo \([a,b]\). De ahí que el problema se denomine valores iniciales en ecuaciones diferenciales ordinarias.

Cuando las condiciones iniciales de la ecuación diferencial incluyen los dos extremos del intervalo \([a,b]\) de definición de la variable independiente \(t\), el problema a resolver se denomina problemas de frontera para ecuaciones diferenciales ordinarias.

En este capítulo, hallaremos soluciones aproximadas de problemas de frontera donde la ecuación diferencial considerada tiene orden \(2\), del tipo \(y''=f(t,y,y')\) y las condiciones iniciales son de tipo frontera: \(y(a)=y_0,\ y(b)=y_1\).

Como ya hicimos en el capítulo de valores iniciales en ecuaciones diferenciales ordinarias, dar la solución aproximada significa dar unos valores discretos en \(n\) valores del mallado del intervalo \([a,b]\), es decir, fijado un paso \(h\), daremos la solución aproximada en los valores \(a+ih\), donde \(h=\frac{b-a}{n}\), \(i=0,\ldots,n\), \(\hat{y}(a+ih)\).

Este tipo de problemas aparecen en física e ingeniería.

El siguiente resultado nos asegura que el problema frontera que vamos a resolver tiene solución única:

Teorema: existencia y unicidad del problema frontera.

Sea \(f\) una función continua en el conjunto: \[ \mathbf{D}=\{(t,y,y')\ a\leq t\leq b,\ -\infty <y<\infty<\ -\infty< y'<\infty\}, \] y que sus derivadas parciales \(\frac{\partial f}{\partial y}\), \(\frac{\partial f}{\partial y'}\) son también continuas en el conjunto \(\mathbf{D}\).

Teorema: existencia y unicidad del problema frontera. (continuación)

Si:

1. \(\frac{\partial f}{\partial y}(t,y,y')>0\), para todo \((t,y,y')\in\mathbf{D}\) y,
2. existe una constante \(M\) tal que \[ \left|\frac{\partial f}{\partial y'}\right|\leq M, \mbox{ para todo }(t,y,y')\in\mathbf{D}, \] el problema de hallar una solución de la ecuación diferencial \(y''=f(t,y,y')\), \(a\leq t\leq b\), con condiciones iniciales tipo frontera \(y(a)=y_0,\ y(b)=y_1\) tiene solución única.

Para introducir el método del tiro, vamos a hallar una aproximación a problemas lineales con condiciones frontera donde la función \(f\) tiene la expresión siguiente: \[ f(t,y,y')=p(t)y'+q(t)y+r(t), \] con \(p(t)\), \(q(t)\) y \(r(t)\) funciones que dependen de la variable independiente \(t\).

Es decir, el problema que nos planteamos es el siguiente:

Problema lineal con condiciones frontera.

Dada una ecuación diferencial de segundo orden de la forma \(y''=p(t)y'+q(t)y+r(t)\), \(a\leq t\leq b\) con \(y(a)=y_0,\ y(b)=y_1\), hallar una solución aproximada \(\hat{y}(t)\) de la ecuación diferencial anterior.

En primer lugar, veamos que dicho problema tiene solución única aplicando el Teorema anterior:

Corolario: existencia y unicidad de problema lineal tipo frontera.

Sea \(f(t,y,y')=p(t)y'+q(t)y+r(t)\), con \(a\leq t\leq b\) y \(p(t)\), \(q(t)\) y \(r(t)\) funciones que dependen de la variable independiente \(t\).

Supongamos que se verifica lo siguiente:

1. \(p(t)\), \(q(t)\) y \(r(t)\) son continuas en \([a,b]\),
2. \(q(t)>0\) en \([a,b]\).

Entonces el problema lineal con condición frontera siguiente: hallar una solución \(y(t)\) de la ecuación diferencial \(y''=p(t)y'+q(t)y+r(t)\), \(a\leq t\leq b\) con \(y(a)=y_0,\ y(b)=y_1\) tiene solución única.

Observación.

La razón que se denominen problemas lineales es que la función \(f(t,y,y')\) es lineal en las variables \(y\) e \(y'\): \(f(t,y,y')=p(t)y'+q(t)y+r(t)\).

Veamos el siguiente resultado que dice cómo hallar la solución del problema de condiciones frontera a partir de dos soluciones de problemas de valores iniciales.

Proposición. Expresión de la solución del problema lineal con condiciones frontera.

Sea \(y''=p(t)y'+q(t)y+r(t)\), \(a\leq t\leq b\), \(y(a)=y_0,\ y(b)=y_1\) un problema lineal con condiciones frontera, es decir, se trata de hallar una solución \(y(t)\) de la ecuación diferencial anterior que cumpla las condiciones frontera indicadas.

Supongamos que las funciones \(p(t),q(t)\) y \(r(t)\) verifican las hipótesis del corolario sobre la existencia y unicidad de soluciones de problemas tipo frontera.

Sea \(y^{(1)}(t)\) la solución del problema de condiciones iniciales siguiente: \[ {y^{(1)}}''=p(t){y^{(1)}}'+q(t)y^{(1)}+r(t),\ a\leq t\leq b,\ y(a)=y_0,\ y'(a)=0. \] Sea \(y^{(2)}(t)\) la solución del problema de condiciones iniciales siguiente: \[ {y^{(2)}}''=p(t){y^{(2)}}'+q(t)y^{(2)},\ a\leq t\leq b,\ y(a)=0,\ y'(a)=1. \]

Proposición. Expresión de la solución del problema lineal con condiciones frontera. (continuación)

Entonces la solución \(y(t)\) del problema con condiciones frontera es el siguiente: \[ y(t)=y^{(1)}(t)+\frac{y_1-y^{(1)}(b)}{y^{(2)}(b)}y^{(2)}(t). \]

Entonces, para hallar la solución aproximada \(\hat{y}(t_i)\) en el mallado \(t_i=a+ih\), \(i=0,\ldots,n\), \(h=\frac{b-a}{n}\) que verifique las condiciones frontera,

• primero tenemos que calcular las aproximaciones \({\hat{y}}^{(1)}\) y \({\hat{y}}^{(2)}\) de los problemas de valores iniciales siguientes usando un método de resolución visto en el capítulo de Problemas de valores iniciales en ecuaciones diferenciales ordinarias: \[ \begin{align*} {y^{(1)}}'' = & p(t){y^{(1)}}'+q(t)y^{(1)}+r(t), &\ a\leq t\leq b,\ y^{(1)}(a)=y_0,\ {y^{(1)}}'(a)=0,\\ {y^{(2)}}'' = & p(t){y^{(2)}}'+q(t)y^{(2)}, &\ a\leq t\leq b,\ y^{(2)}(a)=0,\ {y^{(2)}}'(a)=1. \end{align*} \]
• a continuación, aplicar la expresión dada en la proposición. Entonces, la solución aproximada en el mallado anterior será: \[ \hat{y}(t_i)={\hat{y}}^{(1)}(t_i)+\frac{y_1-{\hat{y}}^{(1)}(b)}{{\hat{y}}^{(2)}(b)}{\hat{y}}^{(2)}(t_i),\ i=0,\ldots,n. \]

Vamos a dar el pseudocódigo donde se usa el método de Runge Kutta 4 para resolver los problemas de valores iniciales correspondientes a las dos ecuaciones diferenciales de segundo orden.

Suponemos pues que tenemos definida una rutina usando el pseudocódigo del capítulo de Problemas de valores iniciales en ecuaciones diferenciales que implementa el método de Runge-Kutta 4. Llamaremos a dicha rutina RK4-sistemas(a,b,h,m,y00,g) donde los valores a y b corresponden a los valores inicial y final respecto la variable \(t\), h es el paso del mallado, m es el número de ecuaciones, y00 sería el vector de condiciones iniciales y g la función que tiene \(m\) componentes respecto las variables \(t,y_1,\ldots,y_m\).

Es decir, el método da una aproximación a la resolución del problema de valores iniciales siguiente:

\[ \begin{align*} y_1'(t) & = g_1(t,y_1(t),\ldots,y_m(t)),\\ y_2'(t) & = g_2(t,y_1(t),\ldots,y_m(t)),\\ & \vdots\\ y_m'(t) & = g_m(t,y_1(t),\ldots,y_m(t)), \end{align*} \] para \(t\in [a,b]\) y condiciones iniciales: \[ y_1(a)=y_{01},\ y_2(a)=y_{02},\ldots,y_m(a)=y_{0m}. \] Concretamente, da como resultado una matriz \(\mathbf{\hat{Y}}\), de dimensiones \(m\times (n+1)\) donde la componente \(i,j\) de dicha matriz sería la aproximación de la componente \(i\)-ésima en el valor \(j\)-ésimo del mallado: \[ y_{ij}=\hat{y}_i(t_j)=\hat{y}_i(a+(j-1)h),\ i=1,\ldots,m,\ j=1,\ldots,n+1. \]

Los sistemas de ecuaciones diferenciales correspondientes a los problemas de valores iniciales considerados son los siguientes:

• Solución \(y^{(1)}\): \[ \begin{align*} y_1'(t) & = y_2,\\ y_2'(t) & = p(t)y_2+q(t)y_1+r(t), \end{align*} \] con \(a\leq t\leq b\) y condiciones iniciales \(y_1(a)=y_0\), \(y_2(a)=0\).

• Solución \(y^{(2)}\): \[ \begin{align*} y_1'(t) & = y_2,\\ y_2'(t) & = p(t)y_2+q(t)y_1, \end{align*} \] con \(a\leq t\leq b\) y condiciones iniciales \(y_1(a)=0\), \(y_2(a)=1\).

Supondremos que la salida de la rutina es una matriz \(\mathbf{Y}\) de dimensiones \(m\times (n+1)\) donde la componente \((i,j)\), \(\mathbf{Y}_{ij}\) sería \(\mathbf{y}_i(t_j)=\mathbf{y}_i(a+hj)\), es decir, el valor de la componente \(\hat{y}_i\) en el punto del mallado \(t_j=a+hj\).

• INPUT a,b, y0, y1, n, p, q, r, RK4-sistemas (nos dan los valores inicial y final del intervalo \([a,b]\), los valores \(y(a)=y_0\) e \(y(b)=y_1\), el número de valores del mallado \(n\), las funciones \(p\), \(q\) y \(r\) y la rutina de resolución de los problemas de valores iniciales, en nuestro caso, el método de Runge-Kutta 4, RK4-sistemas)
• Set \(h=\frac{b-a}{n}\). (calculamos el paso del mallado)

• Rutina g1(t,y1,y2) (Definimos la función \(\mathbf{g}_1\) para hallar la aproximación \(\hat{y}^{(1)}\))
  • Return \((y_2,p(t)\cdot y_2+q(t)\cdot y_1+r(t))\)
• Set \(\mathbf{Y}_1 = RK4-sistemas(a,b,h,2,(y_0,0),\mathbf{g}_1)\) (Calculamos la aproximación \(\hat{y}^{(1)}\))
• Rutina g2(t,y1,y2) (Definimos la función \(\mathbf{g}_2\) para hallar la aproximación \(\hat{y}^{(2)}\))
  • Return \((y_2,p(t)\cdot y_2+q(t)\cdot y_1)\)
• Set \(\mathbf{Y}_2 = RK4-sistemas(a,b,h,2,(0,1),\mathbf{g}_2)\) (Calculamos la aproximación \(\hat{y}^{(2)}\))

• Set \(\mathbf{Y}=\mathbf{Y}_1[1,]+\frac{(y_1-\mathbf{Y}_1[1,n+1])}{\mathbf{Y}_2[1,n+1]}\cdot\mathbf{Y}_2[1,]\) (Calculamos la aproximación del problema de valores frontera usando la expresión vista en la proposición)
• Print \(\mathbf{Y}\).

Fijémonos que en la última expresión para calcular \(\mathbf{Y}\), la aproximación del problema de valores frontera, hemos considerado la primera componente de las aproximaciones \(\mathbf{Y}_1\) e \(\mathbf{Y}_2\) que es la componente que aproxima \(y^{(1)}\) e \(y^{(2)}\) respectivamente.

Vamos a generalizar el método del tiro a problemas de valores frontera para ecuaciones diferenciales de segundo orden no lineales: \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y(b)=y_1. \]

En este caso, no podemos escribir la solución aproximada \(\hat{y}(t)\) del problema de valores frontera como combinación lineal de las soluciones aproximadas de dos problemas de valores iniciales tal como hicimos en el caso lineal.

Podemos interpretar el problema a resolver de la forma siguiente:

Problema de valores frontera.

Hallar la condición inicial \(y_0'\) del problema de valor inicial \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y'(a)=y_0', \] tal que \(y(b)=y_1.\)

Dicho de otra forma, se trata de determinar cuál es la pendiente inicial \(y_0'\) con la que debemos partir para que al final lleguemos a la condición frontera \(y(b)=y_1\).

El gráfico siguiente muestra lo que se intenta solucionar: hallar la pendiente \(y_0'\) tal que la solución aproximada alcance el punto rojo en \(t=b\).

Para solucionar el problema del tiro para sistemas no lineales, escribamos la solución del problema de valor inicial \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y'(a)=y_0', \] de la forma \(y(t,y_0,y_0')\) indicando que depende del tiempo \(t\) y de las condiciones iniciales \(y_0\) e \(y_0'\).

Entonces, el problema de valores frontera se reduce a hallar un cero de la función siguiente: \[ g(y_0')=y(b,y_0,y_0')-y_1=0, \] donde los valores de \(b\), \(y_0\) e \(y_1\) son fijos y el valor que vamos variando es la pendiente inicial \(y_0'\).

El el tema de ceros del curso Métodos numéricos con Python. Cálculo Numérico vimos todo un conjunto de métodos para hallar ceros de funciones de una variable real.

Entre ellos, estaba el método de la secante que consistía en hallar una sucesión \((y_{0,k}')_{k\geq 0}\) de la forma siguiente: \[ \begin{align*} y_{0,k}' = & y_{0,k-1}'-\frac{g(y_{0,k-1}')(y_{0,k-1}'-y_{0,k-2}')}{g(y_{0,k-1}')-g(y_{0,k-2}')}\\ = & y_{0,k-1}'-\frac{(y(b,y_0,y_{0,k-1}')-y_1)(y_{0,k-1}'-y_{0,k-2}')}{y(b,y_0,y_{0,k-1}')-y(b,y_0,y_{0,k-2}')}, \end{align*} \] eligiendo dos valores iniciales \(y_{0,0}'\) e \(y_{0,1}'\).

Por último, una vez hallada la aproximación \(y_{0,k}'\) con la tolerancia permitida se resuelve el siguiente problema de valores iniciales: \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y'(a)=y_{0,k}'. \]

• INPUT a,b, y0, y1, n, f, RK4-sistemas (nos dan los valores inicial y final del intervalo \([a,b]\), los valores \(y(a)=y_0\) e \(y(b)=y_1\), el número de valores del mallado \(n\), la función \(f(t,y,y')\), la rutina de resolución de los problemas de valores iniciales, en nuestro caso, el método de Runge-Kutta 4, RK4-sistemas), los valores iniciales \(y_{0,0}'\) e \(y_{0,1}'\), la tolerancia \(TOL\) y el número máximo de iteraciones \(Nmax\).
• Rutina F(t,y1,y2) (Definimos la función \(\mathbf{F}\) para hallar la aproximación \(\hat{y}\))
  • Return \((y_2,f(t,y_1,y_2))\)
• Rutina y(y0') (Rutina que da el valor \(y(b,y_0,y_0')\) dada una pendiente inicial \(y_0'\))
  • Return \(RK4-sistemas(a,b,h,2,(y_0,y_0'),\mathbf{F})\)

• Set \(h=\frac{b-a}{n}\). (calculamos el paso del mallado)
• Set \(k=1\) (inicializamos el contador que da el número de iteraciones)
• While \(k\leq Nmax\)
  • Set \(y_0'=y_{0,1}'-\frac{(y(b,y_0,y_{0,1}')-y_1)(y_{0,1}'-y_{0,0}')}{y(b,y_0,y_{0,1}')-y(b,y_0,y_{0,0}')}\).
  • If \(|y_0'-y_{0,1}'|\leq TOL\)
    • Print: El valor de la pendiente inicial vale
    • Print \(y_0'\).
    • Print: La solución aproximada vale
    • Print: \(RK4-sistemas(a,b,h,2,(y_0,y_0'),\mathbf{F})\)
    • STOP.

• • Set \(y_{0,0}'=y_{0,1}'\).
  • Set \(y_{0,1}'=y_0'\).
  • Set \(k=k+1\).
• Print: El método no converge al cabo de Nmax iteraciones.

Otro método que vimos en el tema de ceros del curso Métodos numéricos con Python. Cálculo Numérico es el método de Newton-Raphson.

Es uno de los métodos más usados sobre todo por su rápida convergencia ya que si el cero es simple, tiene convergencia cuadrática.

Recordemos que tenemos que hallar un cero de la función \(g(y_0')=y(b,y_0,y_0')-y_1=0\). La sucesión que da el método de Newton-Raphson verifica la recurrencia siguiente: \[ y_{0,k}' = y_{0,k-1}'-\frac{g(y_{0,k-1}')}{g'(y_{0,k-1}')}= y_{0,k-1}'-\frac{(y(b,y_0,y_{0,k-1}')-y_1)}{g'(y_{0,k-1}')}, \] a partir de un valor inicial \(y_{0,0}'\).

El valor de \(g'(y_0')\) en un valor \(y_0'\) cualquiera será: \[ g'(y_0')=\frac{\partial{y(b,y_0,y_0')}}{\partial y_0'}. \] La cuestión que se nos plantea es cómo hallamos la parcial anterior \(\frac{\partial{y(b,y_0,y_0')}}{\partial y_0'}\) de cara a poder aplicar el método de Newton-Raphson.

Para resolver la cuestión anterior, definimos la función siguiente \(z(t,y_0,y_0')=\frac{\partial{y(t,y_0,y_0')}}{\partial y_0'}\). De esta manera, nuestro problema se reduce a hallar \(z(b,y_0,y_0')\).

Para hallar \(z(b,y_0,y_0')\) vamos a calcular la ecuación diferencial que verifica y que integraremos numéricamente desde \(t=a\) hasta \(t=b\).

Recordemos que \(y=y(t,y_0,y_0')\) es la solución del problema de valores iniciales: \[ \begin{align*} & y''(t,y_0,y_0')= f(t,y(t,y_0,y_0'),y'(t,y_0,y_0')),\\ & a\leq t\leq b,\ y(a,y_0,y_0')=y_0,\ y(b,y_0,y_0')=y_1, \end{align*} \] donde hemos indicado explícitamente la dependencia de la solución de dicho problema respecto el tiempo \(t\) y respecto las condiciones iniciales \(y_0\) e \(y_0'\).

Si calculamos la parcial con respecto \(y_0'\) en la ecuación diferencial anterior obtenemos: \[ \frac{\partial y''(t,y_0,y_0')}{\partial y_0'}=\frac{\partial f(t,y(t,y_0,y_0'),y'(t,y_0,y_0'))}{\partial y_0'}. \]

Usando el Teorema de Schwartz o Teorema de Clairaut o Teorema de Young (https://es.wikipedia.org/wiki/Teorema_de_Clairaut) que dice que podemos permutar la parcial anterior con la derivada segunda respecto el tiempo \(t\) y usando la regla de la cadena podemos hallar la ecuación diferencial de segundo orden que verifica la función \(z(t,y_0,y_0')\): \[ \begin{align*} \frac{\partial y''(t,y_0,y_0')}{\partial y_0'}= & z''(t,y_0,y_0')\\ = & \frac{\partial f(t,y(t,y_0,y_0'),y'(t,y_0,y_0'))}{\partial y}\cdot\frac{\partial y(t,y_0,y_0')}{\partial y_0'}\\ & \quad +\frac{\partial f(t,y(t,y_0,y_0'),y'(t,y_0,y_0'))}{\partial y'}\cdot\frac{\partial y'(t,y_0,y_0')}{\partial y_0'} \\ = & \frac{\partial f(t,y,y')}{\partial y}\cdot z(t)+\frac{\partial f(t,y,y')}{\partial y'}\cdot z'(t). \end{align*} \]

En resumen, la función \(z(t)\) verifica la ecuación diferencial de segundo orden siguiente: \[ z''(t)=\frac{\partial f(t,y,y')}{\partial y}\cdot z(t)+\frac{\partial f(t,y,y')}{\partial y'}\cdot z'(t). \] Por último, falta hallar las condiciones iniciales para la función \(z(t,y_0,y_0')\).

Usando que \(y(a,y_0,y_0')=y_0\), realizando la parcial respecto \(y_0'\), obtenemos que: \[ z(a)=\frac{\partial y(a,y_0,y_0')}{\partial y_0'}=0, \] y usando que \(y'(a,y_0,y_0')=y_0'\), realizando la parcial respecto \(y_0'\), obtenemos que: \[ z'(a)=\frac{\partial y'(a,y_0,y_0')}{\partial y_0'}=1. \]

La función \(z(t)\) es la solución del siguiente problema de valores iniciales: \[ z''=\frac{\partial f(t,y,y')}{\partial y}\cdot z(t)+\frac{\partial f(t,y,y')}{\partial y'}\cdot z'(t),\ a\leq t\leq b,\ z(a)=0,\ z'(a)=1. \]

Como dicho sistema depende de \(y(t)\) e \(y'(t)\), tenemos que integrar dicha ecuación diferencial junto con la ecuación diferencial que verifica \(y(t)\) resultando el siguiente sistema de ecuaciones diferenciales: \[ \begin{align*} y''= & f(t,y,y'),\\ z''=& \frac{\partial f(t,y,y')}{\partial y}\cdot z(t)+\frac{\partial f(t,y,y')}{\partial y'}\cdot z'(t), \end{align*} \] con \(a\leq t\leq b\) y condiciones iniciales \(y(a)=y_0,\ y'(a)=y_0',\ z(a)=0,\ z'(a)=1\) obteniéndose \(y(b)\) y \(z(b)\).

Vamos a escribir el sistema anterior como un sistema de ecuaciones diferenciales de primer orden de \(4\) ecuaciones con \(4\) funciones a determinar.

Para ello, definimos \(y_1(t)=y(t),\ y_2(t)=y'(t),\ y_3(t)=z(t),\ y_4(t)=z'(t)\). El sistema diferencial de primer orden que verifican las funciones \(y_i(t)\), \(i=1,2,3,4\) es el siguiente: \[ \left. \begin{align*} y_1'= & y_2,\\ y_2'= & f(t,y_1,y_2),\\ y_3'= & y_4,\\ y_4'= & \frac{\partial f(t,y_1,y_2)}{\partial y_1}\cdot y_3 +\frac{\partial f(t,y_1,y_2)}{\partial y_2}\cdot y_4, \end{align*} \right\} \] con condiciones iniciales \(y_1(a)=y_0,\ y_2(a)=y_0',\ y_3(a)=0,\ y_4(a)=1\).

Entonces para aplicar el método de Newton-Raphson, seguimos los pasos siguientes:

• Paso 1: elegimos un valor inicial \(y_{0,0}'\). Integramos el sistema diferencial anterior con condiciones iniciales \(y_1(a)=y_0,\ y_2(a)=y_{0,0}',\ y_3(a)=0,\ y_4(a)=1\). El valor \(y_{0,1}'\) de la sucesión será: \[ y_{0,1}'=y_{0,0}'-\frac{y_1(b)-y_1}{y_3(b)}. \]

• \(\ldots\)

• Paso \(k\): supongamos calculados \(y_{0,0}'\ldots,y_{0,k-1}'\). Integramos el sistema diferencial anterior con condiciones iniciales \(y_1(a)=y_0,\ y_2(a)=y_{0,k-1}',\ y_3(a)=0,\ y_4(a)=1\). El valor \(y_{0,k}'\) de la sucesión será: \[ y_{0,k}'=y_{0,k-1}'-\frac{y_1(b)-y_1}{y_3(b)}. \]

• Y así sucesivamente hasta que \[ |y_{0,k}'-y_{0,k-1}'|\leq \epsilon, \] donde \(\epsilon\) es la tolerancia permitida.

• Por último, para hallar la solución aproximada, se resuelve numéricamente usando el método de Runge-Kutta 4 el siguiente problema de valores iniciales: \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y'(a)=y_{0,k}'. \]

• INPUT a,b, y0, y1, n, f, \(\frac{\partial f}{\partial y}\), \(\frac{\partial f}{\partial y'}\), RK4-sistemas (nos dan los valores inicial y final del intervalo \([a,b]\), los valores \(y(a)=y_0\) e \(y(b)=y_1\), el número de valores del mallado \(n\), las funciones \(f(t,y,y'),\ \frac{\partial f}{\partial y},\ \frac{\partial f}{\partial y'}\), la rutina de resolución de los problemas de valores iniciales, en nuestro caso, el método de Runge-Kutta 4, RK4-sistemas), el valor inicial \(y_{0,0}'\), la tolerancia \(TOL\) y el número máximo de iteraciones \(Nmax\).
• Rutina F(t,y1,y2,y3,y4) (Definimos la función \(\mathbf{F}\) para hallar las aproximación \(\hat{y}\), \(\hat{z}\))
  • Return \(\left(y_2,f(t,y_1,y_2),y_4,\frac{\partial f}{\partial y}(t,y_1,y_2)y_3+\frac{\partial f}{\partial y'}(t,y_1,y_2)y_4\right)\)

• Rutina yz(y0') (Rutina que da los valores \(y(b,y_0,y_0')\) y \(z(b,y_0,y_0')\) dada una pendiente inicial \(y_0'\))
  • Return \(RK4-sistemas(a,b,h,4,(y_0,y_0',0,1),\mathbf{F})\)
• Set \(h=\frac{b-a}{n}\). (calculamos el paso del mallado)
• Set \(k=1\) (inicializamos el contador que da el número de iteraciones)

• While \(k\leq Nmax\)
  • Set \(y_0'=y_{0,0}'-\frac{(yz(b,y_0,y_{0,0}')[1]-y_1)}{yz(b,y_0,y_{0,0}')[3]}\).
  • If \(|y_0'-y_{0,0}'|\leq TOL\)
    • Print: El valor de la pendiente inicial vale
    • Print \(y_0'\).
    • Print: La solución aproximada vale:
    • Print \(RK4-sistemas(a,b,h,4,(y_0,y_0',0,1),\mathbf{F})\).
    • STOP.
  • Set \(y_{0,0}'=y_0'\).
  • Set \(k=k+1\).
• Print: El método no converge al cabo de Nmax iteraciones.

El método del tiro explicado en la sección anterior puede presentar en algunos casos problemas de inestabilidad numérica.

Vamos a introducir un nuevo método basado en las diferencias finitas que tiene mejor comportamiento con respecto a la estabilidad numérica.

La idea fundamental es considerar un mallado regular en el intervalo de integración \([a,b]\), es decir, considerar \(N+1\) puntos equiespaciados de la forma \(t_i=a+ih\), \(h=\frac{b-a}{N}\), \(i=0,1,\ldots,N\).

El objetivo es hallar la solución de la ecuación diferencial en los puntos del mallado, \(y(t_i)=y_i\), \(i=0,\ldots,N\).

Para ello, consideraremos las siguientes aproximaciones de \(y''(t_i)\) y \(y'(t_i)\) vistas en el tema de Derivación e Integración Numérica en el curso Métodos numéricos con Python. Cálculo Numérico: \[ \begin{align*} y'(t_i)=& \frac{y(t_i+h)-y(t_i-h)}{2h}-\frac{h^2}{6}f'''(\xi_{i,1})=\frac{y_{i+1}-y_{i-1}}{2h}-\frac{h^2}{6}f'''(\xi_{i,1}),\\ y''(t_i)= & \frac{y(t_i+h)-2y(t_i)+y(t_i-h)}{h^2}-\frac{h^2}{12}f^{(4)}(\xi_{i,2})\\ = & \frac{y_{i+1}-2y_i +y_{i-1}}{h^2}-\frac{h^2}{12}f^{(4)}(\xi_{i,2}), \end{align*} \] donde \(\xi_{i,1},\xi_{i,2}\in (t_{i-1},t_{i+1})\).

Recordemos que el problemas de condiciones frontera para problemas lineales era hallar una aproximación \(\hat{y}(t)\) de la ecuación diferencial siguiente: \[ y''=p(t)y'+q(t)y+r(t),\ a\leq t\leq b,\ y(a)=y_0,\ y(b)=y_1. \] Entonces, para resolver el problema anterior, dado un entero positivo \(N\), consideramos el mallado siguiente en el intervalo \([a,b]\), \(t_i=a+ih\), \(h=\frac{b-a}{N}\), \(i=0,1,\ldots,n\).

A continuación, aproximamos \(y''\) e \(y'\) por las expresiones vistas anteriormente: \[ \begin{align*} & y''(t_i) = p(t_i)y'(t_i)+q(t_i)y(t_i)+r(t_i),\ \Rightarrow\\ & \frac{\hat{y}_{i+1}-2\hat{y}_i +\hat{y}_{i-1}}{h^2} = p(t_i)\left(\frac{\hat{y}_{i+1}-\hat{y}_{i-1}}{2h}\right)+q(t_i)\hat{y}_i+r(t_i),\ \Rightarrow \\ & \hat{y}_{i+1}-2\hat{y}_i +\hat{y}_{i-1} = \frac{h}{2}p(t_i)(\hat{y}_{i+1}-\hat{y}_{i-1})+h^2q(t_i)\hat{y}_i+h^2r(t_i),\ \Rightarrow\\ & -\left(1-\frac{h}{2}p(t_i)\right)\hat{y}_{i+1}+(2+h^2q(t_i))\hat{y}_i-\left(1+\frac{h}{2}p(t_i)\right)y_{i-1} = -h^2 r(t_i), \end{align*} \] \(i=1,\ldots,N-1\), ya que \(\hat{y}_0 =\hat{y}(a)=y_0\) y \(\hat{y}_N = \hat{y}(t_N)=\hat{y}(b)=y_1\).

Entonces, las aproximaciones \(\hat{y}(t_i)=\hat{y}_i\) del mallado cumplen el siguiente sistema de ecuaciones lineal \(\mathbf{A}\mathbf{\hat{y}}=\mathbf{b}\), donde \(\mathbf{A}\) es una matriz \((N-1)\times (N-1)\) y \(\mathbf{b}\) un vector \((N-1)\times 1\):

\[ \begin{align*} \mathbf{A}= & \begin{bmatrix}2+h^2 q(t_1)&-1+\frac{h}{2}p(t_1)&0&0&\ldots&0 \\-1-\frac{h}{2}p(t_2)&2+h^2 q(t_2)&-1+\frac{h}{2}p(t_2)&0&\ldots&0 \\0&-1-\frac{h}{2}p(t_3)&2+h^2 q(t_3)&-1+\frac{h}{2}p(t_2)&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&-1-\frac{h}{2}p(t_{N-2})&2+h^2 q(t_{N-2})&-1+\frac{h}{2}p(t_{N-2}) \\0&\ldots&\ldots&0&-1-\frac{h}{2}p(t_{N-1})&2+h^2 q(t_{N-1}) \\\end{bmatrix},\\ \mathbf{\hat{y}}= & \begin{bmatrix}\hat{y}_1\\\hat{y}_2\\\vdots\\\hat{y}_{N-1}\end{bmatrix},\quad\mathbf{b}=\begin{bmatrix}-h^2r(t_1)+\left(1+\frac{h}{2}p(t_1)\right)\hat{y}_0 \\ -h^2r(t_2)\\\vdots\\ -h^2 r(t_{N-2})\\ -h^2 r(t_{N-1})+\left(1-\frac{h}{2}p(t_{N-1})\right)\hat{y}_N\end{bmatrix}. \end{align*} \]

El resultado siguiente analiza cuando el sistema anterior tiene solución única:

Proposición: existencia y unicidad de solución para problemas lineales.

Sean \(p,q\) y \(r\) funciones continuas en el intervalo \([a,b]\). Si \(q(t)> 0\) para \(t\in [a,b]\) y \(h<\frac{2}{L}\) siendo \(L=\max_{t\in [a,b]}|p(t)|\), el sistema de ecuaciones anterior tiene solución única, es decir, la matriz \(\mathbf{A}\) es no singular o equivalentemente tiene rango máximo \(N-1\).

Observación.

Tal como se comenta en el tema de Derivación e Integración Numérica en el curso Métodos numéricos con Python. Cálculo Numérico, las expresiones usadas en la derivación numérica son inestables en el sentido que no podemos usar valores \(h\) arbitrariamente pequeños o ir aumentando \(N\) de forma incontrolada.

Por tanto, debemos usar valores de \(h\) moderadamente pequeños pero no demasiado pequeños.

• INPUT a,b, y0, y1, N, p, q, r (nos dan los valores inicial y final del intervalo \([a,b]\), los valores \(y(a)=y_0\) e \(y(b)=y_1\), el valor \(N\), y las funciones \(p\), \(q\) y \(r\) )
• Set \(h=\frac{b-a}{N}\). (calculamos el paso del mallado)
• Set \(\mathbf{A}=\mathbf{0}\) (\((N-1)\times (N-1)\)) (Vamos a definir la matriz \(\mathbf{A}\) y el vector \(\mathbf{b}\))
• Set \(\mathbf{b}=\mathbf{0}\) (\((N-1)\times 1\))

• For i=2,...,N-2
  • Set \(t=a+ih\).
  • Set \(a_{ii}=2+h^2\cdot q(t)\).
  • Set \(a_{i,i+1}=-1+\frac{h}{2}p(t)\).
  • Set \(a_{i,i-1}=-1-\frac{h}{2}p(t)\).
  • Set \(b_i=-h^2\cdot r(t)\).
• Set \(a_{11}=2+h^2\cdot q(a+h)\).
• Set \(a_{N-1,N-1}=2+h^2\cdot q(a+(N-1)h)\).
• Set \(a_{12}=-1+\frac{h}{2}p(a+h)\).
• Set \(a_{N-1,N-2}=-1-\frac{h}{2}p(a+(N-1)h)\).

• Set \(b_1=-h^2r(a+h)+\left(1+\frac{h}{2}p(a+h)\right)y_0\).
• Set \(b_{N-1}=-h^2r(a+(N-1)h)+\left(1-\frac{h}{2}p(a+(N-1)h)\right)y_1\).
• Solve \(\mathbf{A}\mathbf{\hat{y}}=\mathbf{b}\). (Resolver el sistema correspondiente usando un método directo)
• Print \(\mathbf{\hat{y}}\). (Damos la aproximación)

Vamos a aplicar el método de las diferencias finitas a un problema general de valores frontera: \[ y''=f(t,y,y'),\ a\leq t\leq b,\ y(a)=y_0,\ y(b)=y_1. \] Para poder aplicar dicho método supondremos que la función \(f\) verifica las condiciones siguientes:

• la función \(f\) y sus derivadas parciales \(\frac{\partial f}{\partial y}\), \(\frac{\partial f}{\partial y'}\) son continuas en el dominio: \[ \mathbf{D}=\{(t,y,y'),\ |\ a\leq t\leq b,\ -\infty <y<\infty,\ -\infty <y'<\infty\}. \]
• existe un valor \(\delta >0\) tal que \(\frac{\partial f}{\partial y}\geq \delta\).
• existen constantes \(k\) y \(L\) tal que: \[ k=\max_{(t,y,y')\in\mathbf{D}}\left|\frac{\partial f}{\partial y}\right|,\quad L=\max_{(t,y,y')\in\mathbf{D}}\left|\frac{\partial f}{\partial y'}\right|. \]

Entonces, usando el Teorema general de existencia y unicidad de problemas frontera, podemos afirmar que nuestro problema tiene solución única.

A continuación, tal como hicimos con los problemas lineales consideramos

• un mallado regular en el intervalo de integración \([a,b]\): \(t_i=a+ih\), \(h=\frac{b-a}{N}\), \(i=0,1,\ldots,N\),
• las siguientes aproximaciones de \(y''(t_i)\) y \(y'(t_i)\): \[ \begin{align*} y'(t_i)=& \frac{y(t_i+h)-y(t_i-h)}{2h}-\frac{h^2}{6}f'''(\xi_{i,1})\\ =& \frac{y_{i+1}-y_{i-1}}{2h}-\frac{h^2}{6}f'''(\xi_{i,1}),\\ y''(t_i)= & \frac{y(t_i+h)-2y(t_i)+y(t_i-h)}{h^2}-\frac{h^2}{12}f^{(4)}(\xi_{i,2})\\ = & \frac{y_{i+1}-2y_i +y_{i-1}}{h^2}-\frac{h^2}{12}f^{(4)}(\xi_{i,2}), \end{align*} \]

obteniéndose el siguiente sistema no lineal cuyas incógnitas son los valores aproximados \(\hat{y}_i\) en los puntos del mallado: \[ \begin{align*} & \frac{\hat{y}_{i+1}-2\hat{y}_i +\hat{y}_{i-1}}{h^2} = f\left(t_i,\hat{y}_i,\frac{\hat{y}_{i+1}-\hat{y}_{i-1}}{2h}\right),\ i=1,\ldots,N-1,\ \Rightarrow \\ & \hat{y}_{i+1}-2\hat{y}_i +\hat{y}_{i-1}-h^2f\left(t_i,\hat{y}_i,\frac{\hat{y}_{i+1}-\hat{y}_{i-1}}{2h}\right) = 0,\ i=1,\ldots,N-1, \end{align*} \] con \(\hat{y}_0=y_0\) y \(\hat{y}_N=y_1\).

Si escribimos el sistema anterior en componentes, obtenemos:

\[ \begin{align*} \hat{y}_{2}-2\hat{y}_1 -h^2f\left(t_1,\hat{y}_1,\frac{\hat{y}_{2}-y_0}{2h}\right)+y_0 =& 0,\\ \hat{y}_{3}-2\hat{y}_2+\hat{y}_1 -h^2f\left(t_2,\hat{y}_2,\frac{\hat{y}_{3}-\hat{y}_{1}}{2h}\right) =& 0,\\ & \vdots\\ \hat{y}_{N-1}-2\hat{y}_{N-2}+\hat{y}_{N-3} -h^2f\left(t_{N-2},\hat{y}_{N-2},\frac{\hat{y}_{N-1}-\hat{y}_{N-3}}{2h}\right) =& 0,\\ -2\hat{y}_{N-1}+\hat{y}_{N-2} -h^2f\left(t_{N-1},\hat{y}_{N-1},\frac{y_1 -\hat{y}_{N-2}}{2h}\right)+y_1 =& 0. \end{align*} \]

Para resolver el sistema anterior, podemos usar el método de Newton-Raphson visto en el tema de Aproximación a las soluciones de sistemas de ecuaciones no lineales.

En primer lugar, si pensamos que el sistema a resolver es de la forma siguiente: \[ \mathbf{F}(\hat{y}_1,\ldots,\hat{y}_{N-1})=\mathbf{0}, \] con \[ \mathbf{F}(\hat{y}_1,\ldots,\hat{y}_{N-1})=\begin{bmatrix}\hat{y}_{2}-2\hat{y}_1 -h^2f\left(t_1,\hat{y}_1,\frac{\hat{y}_{2}-y_0}{2h}\right)+y_0\\ \hat{y}_{3}-2\hat{y}_2+\hat{y}_1 -h^2f\left(t_2,\hat{y}_2,\frac{\hat{y}_{3}-\hat{y}_{1}}{2h}\right)\\ \vdots\\ \hat{y}_{N-1}-2\hat{y}_{N-2}+\hat{y}_{N-3} -h^2f\left(t_{N-2},\hat{y}_{N-2},\frac{\hat{y}_{N-1}-\hat{y}_{N-3}}{2h}\right)\\ -2\hat{y}_{N-1}+\hat{y}_{N-2} -h^2f\left(t_{N-1},\hat{y}_{N-1},\frac{y_1 -\hat{y}_{N-2}}{2h}\right)+y_1, \end{bmatrix} \]

necesitamos calcular la matriz jacobiana de la función anterior \(\mathbf{F}\): \(\mathbf{J}=(\mathbf{J}_{ij})_{i,j=1,\ldots,N-1}\).

Usando que \[ F_i=\begin{cases}\hat{y}_{2}-2\hat{y}_1 -h^2f\left(t_1,\hat{y}_1,\frac{\hat{y}_{2}-y_0}{2h}\right)+y_0, & \mbox{si }i=1,\\ \hat{y}_{i+1}-2\hat{y}_i+\hat{y}_{i-1}-h^2 f\left(t_i,\hat{y}_i,\frac{\hat{y}_{i+1}-\hat{y}_{i-1}}{2h}\right), & \mbox{si }i=2,\ldots,N-2,\\ -2\hat{y}_{N-1}+\hat{y}_{N-2} -h^2f\left(t_{N-1},\hat{y}_{N-1},\frac{y_1 -\hat{y}_{N-2}}{2h}\right)+y_1, & \mbox{si }i=N-1, \end{cases} \] la matriz jacobiana será la siguiente:

\[ \mathbf{J}_{ij}=\frac{\partial F_i}{\partial \hat{y}_j}=\begin{cases} 1+\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'},&\mbox{si }j=i-1,\\ -2-h^2\frac{\partial f}{\partial\hat{y}_i},&\mbox{si }j=i,\\ 1-\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'},&\mbox{si }j=i+1,\\ 0,& \mbox{en los demás casos,} \end{cases} \] para \(i,j=1,\ldots,N-1\).

El método de Newton-Raphson consiste en definir una sucesión de aproximaciones \((\mathbf{\hat{y}}^{(k)})_{k\geq 0}\) usando la siguiente recurrencia: \[ \mathbf{\hat{y}}^{(k)}=\mathbf{\hat{y}}^{(k-1)}-{\mathbf{J}^{(k)}}^{-1}\mathbf{F}(\mathbf{\hat{y}}^{(k-1)}), \] a partir de una aproximación inicial \(\mathbf{\hat{y}}^{(0)}\).

En la práctica, se procede de la forma siguiente:

Supongamos calculados \(\mathbf{\hat{y}}^{(0)},\ldots,\mathbf{\hat{y}}^{(k-1)}\). Para calcular \(\mathbf{\hat{y}}^{(k)}\),

• en primer lugar, resolvemos el sistema lineal: \[ \mathbf{J}^{(k)}\mathbf{v}=-\mathbf{F}(\mathbf{\hat{y}}^{(k-1)}), \] usando un método directo visto en el tema de Métodos Directos de Resolución de sistemas lineales del curso Métodos Numéricos con Python. Álgebra lineal Numérica
• una vez calculada la solución \(\mathbf{v}\), la nueva aproximación \(\mathbf{\hat{y}}^{(k)}\) se calcula de la forma siguiente: \[ \mathbf{\hat{y}}{(k)}=\mathbf{\hat{y}}{(k-1)}+\mathbf{v}. \]

• repetimos el proceso anterior hasta que \(\|\mathbf{\hat{y}}^{(k)}-\mathbf{\hat{y}}^{(k-1)}\|<\epsilon\), donde \(\epsilon\) es la tolerancia permitida.

• INPUT a,b, y0, y1, n, f, \(\frac{\partial f}{\partial y}\), \(\frac{\partial f}{\partial y'}\), los valores \(y(a)=y_0\) e \(y(b)=y_1\), el número de valores del mallado \(n\), las funciones \(f(t,y,y'),\ \frac{\partial f}{\partial y},\ \frac{\partial f}{\partial y'}\), el valor inicial de la sucesión de aproximaciones \(\mathbf{\hat{y}}^{(0)}\), la tolerancia \(TOL\) y el número máximo de iteraciones \(Nmax\).
• Set \(h=\frac{b-a}{n}\). (Definimos el paso \(h\))
• Set \(k=1\) (Inicializamos el contador de iteraciones)

• While \(k\leq Nmax\)
  • Set \(\mathbf{J}=\mathbf{0}\) (Definimos la matriz jacobiana y el valor \(\mathbf{F}(\mathbf{\hat{y}}^{(0)})\)
  • Set \(\mathbf{F}=\mathbf{0}\)
  • Set \(t=a+h\).
  • Set \(x=\frac{\hat{y}_2^{(0)}-y_0}{2h}\).
  • Set \(\mathbf{F}_1=\hat{y}_2^{(0)}-2\hat{y}_1^{0}-h^2 f(t,\hat{y}_1^{(0)},x)+y_0\).
  • Set \(\mathbf{J}_{11}=-2-h^2\frac{\partial f}{\partial\hat{y}_i}(t,\hat{y}^{(0)}_1,x)\).
  • Set \(\mathbf{J}_{12}=1-\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'}(t,\hat{y}^{(0)}_1,x)\).

• • For i=2,...,N-2
    • Set \(t=a+i\cdot h\)
    • Set \(x=\frac{\hat{y}_{i+1}^{(0)}-\hat{y}_{i-1}^{(0)}}{2h}\).
    • Set \(\mathbf{F}_i=\hat{y}_{i+1}^{(0)}-2\hat{y}_i^{0}+\hat{y}_{i-1}^{(0)}-h^2 f(t,\hat{y}_i^{(0)},x)\).
    • Set \(\mathbf{J}_{i,i-1}=1+\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'}(t,\hat{y}^{(0)}_i,x)\).
    • Set \(\mathbf{J}_{ii}=-2-h^2\frac{\partial f}{\partial\hat{y}_i}(t,\hat{y}^{(0)}_i,x)\).
    • Set \(\mathbf{J}_{i,i+1}=1-\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'}(t,\hat{y}^{(0)}_i,x)\).

• • Set \(t=a+(N-1)h\).
  • Set \(x=\frac{y_1-\hat{y}_{N-2}^{(0)}}{2h}\).
  • Set \(\mathbf{F}_{N-1}=-2\hat{y}_{N-1}^{0}+\hat{y}_{N-2}^{(0)}-h^2 f(t,\hat{y}_{N-1}^{(0)},x)+y_1\).
  • Set \(\mathbf{J}_{N-1,N-1}=-2-h^2\frac{\partial f}{\partial\hat{y}_i}(t,\hat{y}^{(0)}_{N-1},x)\).
  • Set \(\mathbf{J}_{N-1,N-2}=1+\frac{h}{2}\frac{\partial f}{\partial\hat{y}_i'}(t,\hat{y}^{(0)}_{N-1},x)\).
  • Solve \(\mathbf{J}\mathbf{v}=-\mathbf{F}\) (Resolvemos el sistema lineal correspondiente)

• • If \(\|\mathbf{v}\|_2\leq TOL\) (Miramos si la normal euclídea de la diferencia entre dos valores consecutivos de la sucesión \((\mathbf{\hat{y}}^{(k)})_{k\geq 0})\) es menor que la tolerancia)
    • Print: solución aproximada.
    • Print: \(\mathbf{\hat{y}}^{(0)}\).
    • STOP.
  • Set \(\mathbf{\hat{y}}^{(0)}=\mathbf{\hat{y}}^{(0)}+\mathbf{v}\).
  • Set \(k=k+1\).
• Print: el método no converge al cabo de \(Nmax\) iteraciones.
• STOP.