Métodos numéricos generales no estándar para la solución directa de ecuaciones diferenciales no despejadas en formas canónicas.

==Métodos numéricos generales no estándar para la solución directa de ecuaciones diferenciales no despejadas en formas canónicas.==

'''Carlos Oscar Rodríguez Leal 

Universidad de Guadalajara. CUCEI.

e-mail: misifu17@gmail.com'''

==Resumen==

En este trabajo se desarrollan algoritmos numéricos que se pueden aplicar de manera directa a ecuaciones diferenciales de la forma general <math display="inline">f(t,x,x') = 0</math>, sin necesidad de despejar <math display="inline">x'</math>. Mis métodos son algoritmos híbridos entre los métodos estándar de solución de ecuaciones diferenciales y métodos de solución de ecuaciones algebraicas, con los que se despeja numéricamente la variable <math display="inline">x'</math>.

La aplicación de estos métodos va desde ecuaciones diferenciales ordinarias de orden uno, hasta el caso más general de sistemas de <math display="inline">m</math> ecuaciones de orden <math display="inline">n</math>. Dichos algoritmos se aplican a la solución de diversas ecuaciones de la física-matemática.

Por último, se hace el análisis numérico correspondiente de existencia, unicidad, estabilidad, consistencia y convergencia, principalmente para el caso más sencillo de una sola ecuación diferencial ordinaria de primer orden.

'''Palabras clave''' Ecuaciones diferenciales, métodos numéricos no estándar, función derivada implícita, método de Euler-bisección, consistencia, estabilidad, convergencia, ecuaciones de estrella de bosones.

Non-standard general numerical methods for the direct solution of differential equations not cleared in canonical forms

==Abstract==

In this work I develop numerical algorithms that can be applied directly to differential equations of the general form <math display="inline">f (t, x, x ') = 0</math>, without the need to cleared <math display="inline">x'</math>. My methods are hybrid algorithms between standard methods of solving differential equations and methods of solving algebraic equations, with which the variable <math display="inline">x'</math> is numerically cleared.

The application of these methods ranges from the ordinary differential equations of order one, to the more general case of systems of <math display="inline">m</math> equations of order <math display="inline">n</math>. These algorithms are applied to the solution of different physical-mathematical equations.

Finally, the corresponding numerical analysis of existence, uniqueness, stability, consistency and convergence is made, mainly for the simplest case of a single ordinary differential equation of the first order.

'''Keywords''' Differential equations, non-standard numerical methods, implicit derivative function, Euler-bisection method, consistency, stability, convergence, boson star equations.

==1 Introducción==

En este trabajo se desarrollan algoritmos numéricos que se pueden aplicar de manera directa a ecuaciones diferenciales de la forma general <math display="inline">f(t,x,x') = 0</math>, sin necesidad de despejar algebraicamente <math display="inline">x'</math> de la forma <math display="inline">x' = f(t,x)</math> para poder emplear los métodos ordinarios.

Los algoritmos desarrollados son modificaciones de métodos estándar para resolver ecuaciones diferenciales con condiciones iniciales, en particular el método de Euler. La esencia de la modificación consiste en desarrollar algoritmos híbridos entre los métodos estándar de solución de ecuaciones diferenciales y métodos de solución de ecuaciones algebraicas, con los cuales se despejan numéricamente la variable <math display="inline">y'</math>. Los métodos se desarrollan tanto para ecuaciones diferenciales ordinarias de orden uno, como para el caso más general de sistemas de <math display="inline">m</math> ecuaciones con <math display="inline">m</math> incógnitas de orden <math display="inline">n</math>.

La importancia de estos métodos radica en que en la práctica muchas veces es muy difícil hacer los despejes algebraicos de <math display="inline">y'</math> en las ecuaciones diferenciales, si no es que imposible.

Dichos procedimientos directos los aplico a ecuaciones de la física-matemática, entre ellas la ecuación de caída libre con fricción y las ecuaciones de Einstein para las estrellas de bosones.

Por último, realizo el análisis numérico correspondiente de existencia, unicidad, estabilidad, consistencia y convergencia de las soluciones numéricas proporcionadas por mis métodos, principalmente para el caso más sencillo de una sola ecuación diferencial ordinaria de primer orden resuleta por mi método de Euler modificado, dejando para trabajos posteriores el análisis correspondiente de los demás casos.

==2 Definiciones y teoremas previos==

En esta sección se enunciarán algunos teoremas ya establecidos que serán necesarios para la elaboración de este trabajo.

<span id='theorem-teorema5.19Burden'></span>Teorema 1:     Supóngase que el problema de valor inicial

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x' = g(t,x),  t_0 \leq t \leq t_f, x(t_0) = x_0,    </math>
|}
|}

se aproxima por el método de diferencias de un paso de la forma

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{x}_0 = t_0,          \tilde{x}_{n+1} = \tilde{x}_n + h \phi (t_n, \tilde{x}_n; h).    </math>
|}
|}

Considérese también que existe un número <math display="inline">h_0 > 0</math> y <math display="inline">\phi (t,x;h)</math> es continua y satisface una condición de Lipschitz en la variable <math display="inline">x</math> con constante de Lipschitz <math display="inline">L</math> en el conjunto

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>F = \{ (t,x,h | t_0 \leq t \leq t_f, -\infty < x < \infty , 0 \leq h \leq h_0 \} .    </math>
|}
|}

Entonces el método es estable.

<span id='theorem-y_biprima_Burden'></span>Teorema 2:     Si en la ecuación <math display="inline">x' = g(t,x)</math>, <math display="inline">\partial g / \partial t</math> y <math display="inline">\partial g / \partial y</math> existen, entonces

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x''(t) = \frac{dx'}{dt}(t) = \frac{dg}{dt}(t,x(t)) = \frac{\partial g}{\partial t}(t, x(t)) + \frac{\partial g}{\partial x}(t,x(t)) \cdot g(t, x(t)).    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

==3 Método de Euler modificado o de Euler-Leal==

En esta primera sección, comenzando por los casos más simples, se explican tanto la elaboración como el análisis numérico para el método de Euler modificado, combinando el método de Euler con el método algebraico de bisección.

===3.1 Método de Euler-Bisección===

En este primer método, como punto de partida, se emplan los dos métodos más simples, intuitivos y fáciles de analizar, los cuales son el método de Euler, para ecuaciones diferenciales, y el método de la bisección.

====3.1.1 Desarrollo del algoritmo====

Así pues, comenzaremos considerando la ecuación general sin despejar

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(t,x,x') = 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

para la cual calcularemos los distintos valores aproximados de <math display="inline">x'</math> en un intervalo <math display="inline">t_0 \leq t \leq t_f</math>, con un tamaño de paso <math display="inline">h</math> y con la condición inicial <math display="inline">x(t_0) = x_0</math>; la cual, al sustituirla en la función <math display="inline">f</math>, generará la expresión

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(t_0,x_0,x') = 0.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

Entonces, dado que no tenemos despejada de manera explícita la variable <math display="inline">x'</math>, para encontrar su valor, de manera aproximada, usaremos el método de la bisección, donde tomamos dos valores de arranque inamovibles <math display="inline">\alpha _i</math> y <math display="inline">\alpha _d</math>, <math display="inline">\alpha _i < \alpha _d</math>, tales que

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(t_0,x_0,\alpha _i) \cdot f(t_0,x_0,\alpha _d) < 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

Además la diferencia

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _d - \alpha _i,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

deberá ser convenientemente grande en lugar de ser convenientemente pequeña (lo cual se explicará más adelante).

Luego de ello comenzamos las iteraciones del núcleo del algoritmo.

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | '''1.''' <math> {\tilde{x}'}_{0} = Biseccion(t_{0},\tilde{x}_{0},\alpha _i,\alpha _d,\tau ) </math>
 '''2.'''  para  <math> (n = 0 \hbox{ hasta } N-1) </math>
'''a)'''   <math> t_{n+1} = t_n + h </math>'''b)'''   <math> \tilde{x}_{n+1} = \tilde{x}_n + {\tilde{x}'}_n h </math> '''c)'''   <math> \tau = \tau (h) </math>'''d)'''   <math>{\tilde{x}'}_{n+1} = Biseccion(t_{n+1},\tilde{x}_{n+1},\alpha _i,\alpha _d,\tau ) </math> '''3.'''  fin_para
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde en la ecuación del inciso c)

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = \tau (h),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

<math display="inline">\tau </math> es la tolerancia o error para el método de la bisección, de tal manera que la subrutina de bisección se detendrá cuando <math display="inline">{\tilde{x}'}_{{(n+1)}_m} - {\tilde{x}'}_{{(n+1)}_{m-1}} \leq \tau </math>, para alguna <math display="inline">m</math>-ésima iteración del algoritmo anidado “Biseccion” dentro del núcleo del algoritmo de Euler-biseccion en su iteración <math display="inline">n</math>-ésima. El porqué se ha elegido a la tolerancia <math display="inline">\tau </math> como función de <math display="inline">h</math> será aclararado en la sección del análisis de este método.

Por ahora diremos que la elección de la diferencia ([[#eq-5|5]]) debe ser suficientemente amplia, pues si dicha diferencia es pequeña, aunque el teorema del valor intermedio nos garantiza que existe una raiz <math display="inline">x'</math> para <math display="inline">f</math> dados <math display="inline">t_0</math> y <math display="inline">x_0</math> fijos y cumpliéndose ([[#eq-4|4]]) (bajo el supuesto de que la función <math display="inline">f</math> es continua, lo que se analizará más adelante, en la subsección de convergencia), como tal función <math display="inline">f</math> irá variando sus valores <math display="inline">t_{n}</math> y <math display="inline">x_{n}</math> al entrar al núcleo recursivo ([[#eq-6|6]]), para esos nuevos valores fijados (con <math display="inline">n</math> suficientemente grande) dicha función podría ya no tener una raiz para <math display="inline">x'</math> dentro del intervalo pequenño <math display="inline">[\alpha _i, \alpha _d]</math>. Pero con un intervalo adecuadamente amplio y donde <math display="inline">{x'}_{0}</math> esté más o menos centrada, disminuirá considerablemente el riesgo de que la raiz se localice fuera de dicho intervalo y el método aborte dándonos un mensaje de error (pues está programado para tolerar ese fallo).

====3.1.2 Condiciones generales para un problema bien planteado====

Ahora procederemos a hacer el correspondiente análisis de unicidad y buen planteamiento, del método numérico anterior, estableciendo los criterios bajo los cuales se cumplirán ambos requisitos.

Observación 1: En todos los análisis de este artículo siempre tomaremos como punto de partida algún espacio <math display="inline">\mathbb{R}^n</math>, para algún <math display="inline">n \in \mathbb{N}</math> y consideraremos una métrica euclideana, es decir, partiremos de un espacio euclidiano <math display="inline">n</math>-dimensional.

Así pues, emplearemos las condiciones de la sección 3.1.1, con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>t_0 \leq t \leq t_f, </math>
|-
| style="text-align: center;" | <math> - \infty < x < \infty , </math>
|-
| style="text-align: center;" | <math> - \infty < x' < \infty , </math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Ahora consideraremos las siguientes hipótesis de esta subsección que han de cumplirse, y sus correspondientes tesis.

Definición 1: Definiremos primeramente a <math display="inline">F = [t_0,t_f]\times \mathbb{R} \times \mathbb{R}</math>, lo cual significa que <math display="inline">F</math> es el dominio de definición de <math display="inline">f</math>.

<span id='theorem-condiciones'></span>Teorema 3:

<ol>

<li>Si extendemos el conjunto <math display="inline">F</math>, eligiendo <math display="inline">G \supset F</math> definido como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

G = (\beta _i,\beta _d) \times \mathbb{R}^2, </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

de tal forma que <math display="inline">\beta _i < t_0</math> y <math display="inline">\beta _d > t_f</math>, tendremos que <math display="inline">G</math> será un subconjunto abierto de <math display="inline">R^3</math> que contendrá a <math display="inline">F</math>.

<li>Por lo tanto, si ampliamos el dominio de definición de <math display="inline">f</math> a <math display="inline">G</math>, es decir <math display="inline">f:(\beta _i,\beta _d) \rightarrow \mathbb{R}</math>, y se cumple que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\frac{\partial f}{\partial t},  \frac{\partial f}{\partial x}, \textrm{ y } \frac{\partial f}{\partial x'} </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

existen y son continuas en el conjunto abierto convexo <math display="inline">G</math>, entonces, por el teorema de la continuidad de una función solamente si sus derivadas parciales son continuas (Rudin, 1980, p. 236), tendremos que <math display="inline">f \in C'(G)</math> y por lo tanto <math display="inline">f \in C'(F)</math> (definiendo las derivadas por derecha y por izquierda en los respectivos extremos <math display="inline">\alpha _i</math> y <math display="inline">\alpha _d</math>).

<li>Si además, dados <math display="inline">t = t_r \in (\beta _i,\beta _d)</math> y <math display="inline">x = x_r \in \mathbb{R}</math> fijos, se cumple

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\frac{\partial f}{\partial x'}|_{t_r,x_r} \neq 0, </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

para todo <math display="inline">x' \in \mathbb{R}</math>, entonces podemos definir a la función

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

u: _{t_r, x_r}\mathbb{R} \rightarrow \mathbb{R}, x' \mapsto u_{t_r, x_r}(x') := f(t_r,x_r,x'), 

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

siendo con ello la derivada

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\frac{\partial f}{\partial x'}|_{t_r,x_r} = \frac{d u_{t_r,x_r}}{d x'} \neq 0 

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

continua en <math display="inline">\mathbb{R}</math>, i.e., por el teorema del valor intermedio, claramente <math display="inline">{u'}_{t_r,x_r}(x') > 0</math> ó <math display="inline">{u'}_{t_r,x_r}(x') < 0</math> (pero no ambos) para todo <math display="inline">x' \in \mathbb{R}</math>. De lo anterior, usando el teorema de la función estrictamente monótona (Lang, 1990, p. 60) se deduce que <math display="inline">u_{t_r,x_r}</math> es estrictamente monótona en <math display="inline">\mathbb{R}</math> para todo par <math display="inline">(t,x) \in H = (\beta _i, \beta _d) \times \mathbb{R}</math>. En particular <math display="inline">u_{t_r,x_r}</math> es estrictamente monótona en <math display="inline">\mathbb{R}</math> para todo <math display="inline">(t,x) \in I = [t_0, t_f]\times \mathbb{R}</math>.

<li>Y si también se cumple que para <math display="inline">t = t_r \in (\beta _i,\beta _d)</math> y <math display="inline">x = x_r \in \mathbb{R}</math> fijos existen <math display="inline">{x'}_{i_r}</math> y <math display="inline">{x'}_{d_r}</math>, con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\infty {< x'}_{i_r}{ < x'}_{d_r} < \infty , </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

tales que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

{f(t_r,x_r,x'}_{i_r}{) f(t_r,x_r,x'}_{d_r}) < 0, 

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

es decir

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

u_{t_r,x_r}{(x'}_{i_r}) u_{t_r,x_r}{(x'}_{d_r}) < 0, 

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

entonces, usando una vez más el teorema de valor intermedio, se puede demostrar que existe un único <math display="inline">{x'}_r</math>, con <math display="inline">{x'}_{i_r}{ < x'}x'_< x'_{d_r}</math>, que satisface <math display="inline">u_{t_r,x_r}{(x'}_r) = 0</math>, o lo que es lo mismo <math display="inline">{f(t_r,x_r,x'}_r) = 0</math>. En particular lo anterior es cierto para <math display="inline">t \in [t_0,t_f]</math>. 

</ol>

<span id='theorem-T.F.I'></span>Teorema 4: Cosiderando los puntos del teorem anterior, y observando que la ecuación <math display="inline">f(t,x,x') = 0</math> definida en <math display="inline">G</math> determina implícitamente a <math display="inline">x' = g(t,x)</math> (que es la función que no podemos despejar), entonces, de acuerdo al punto 4. de dicho teorema y utilizando el teorema de la función implícita, se demuestra que para cada par <math display="inline">(t_r,x_r) \in H</math> existe un <math display="inline">{x'}_r \in \mathbb{R}</math> que satisface <math display="inline">{f(t_r,x_r,x'}_r) = 0</math>, de tal forma que para dicha triada <math display="inline">{(t_r,x_r,x'}_r)</math> se cumple lo siguiente (es decir, en cada punto <math display="inline">(t_r,x_r)</math> del dominio <math display="inline">H</math> para <math display="inline">g</math>):

<ol>

<li><math display="inline">{x'}_r=g(r_r,y_r)</math> es función monovaluada.    </li>
<li><math display="inline">{x'}_r=g(r_r,y_r)</math> es continua y con primera derivada continua.    </li>
<li><math display="inline">\frac{\partial g}{\partial x}|_{t_r,x_r}  = - \frac{\frac{\partial f}{\partial x}|_{t_r,x_r}}{\frac{\partial f}{\partial x'}|_{t_r,x_r}}</math>. </li>

</ol>

Corolario 1: En particular, el teorema anterior es cierto para el dominio restringido <math display="inline">I</math> para <math display="inline">g</math>.

Teorema 5: Utilizando el teorema de la derivada parcial para la condición de Lipschitz <span id='citeF-1'></span>[[#cite-1|[1]]] y el punto 3. del teorema ([[#theorem-T.F.I|4]]) resulta evidente que si existe una constante <math display="inline">L > 0</math> que satisface

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>| \frac{\partial g}{\partial x} | = | \frac{\frac{\partial f}{\partial x}}{\frac{\partial f}{\partial x'}} | \leq L, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

para toda <math display="inline">(t,x) \in I</math>, entonces <math display="inline">g</math> satisface una condición de Lipschitz sobre <math display="inline">I</math> en la variable <math display="inline">x</math> y con constante de Lipschitz <math display="inline">L</math>.

<span id='theorem-problema_bien_planteado'></span>Teorema 6: Ahora, haciendo uso del teorema del problema de valor inicial bien planteado (Burden y Faires, 1996, p. 242), si la función <math display="inline">g</math> definida implícitamente es continua en <math display="inline">I</math> y  satisface una condición de Lipschitz en la variable <math display="inline">x</math> sobre <math display="inline">I</math>, entonces el problema de valor inicial

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x' = g(t,x),  t_0 \leq t \leq t_f,  x(t_0) = x_0,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

está bien planteado.

====3.1.3 Análisis del error del método de Euler-bisección.====

En esta subsección se hará un análisis particular para el error de <math display="inline">\tilde{x}</math> respecto a <math display="inline">x</math>, su orden y el error óptimo, para el método de Euler-bisección, suponeindo que se cumplen las hipótesis de la subsección anterior.

<span id='theorem-Euler_perturbado'></span>Teorema 7: Dada la tolerancia <math display="inline">\tau (h)</math> definida en ([[#eq-7|7]]), así como los extremos <math display="inline">t_0</math> y <math display="inline">t_f</math>, entonces es claro que la subrutina de bisección se detendrá en cuanto se cumpla que <math display="inline">|\tilde{x}' - x'| \leq \frac{t_f - t_0}{2^{N(h)} } \leq \tau </math>, para algún <math display="inline">N \in \mathbb{N}_0</math>, donde <math display="inline">N(h)</math> no ha de depender de <math display="inline">t</math>, <math display="inline">x</math> o <math display="inline">x'</math> y por lo tanto siempre valdrá lo mismo para un <math display="inline">h</math> dado. Por lo tanto podemos expresar a dicho error como

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0 \leq |\tilde{x}' - x'| < \frac{t_f - t_0}{2^{N(h)} } =: \delta (h) = \tau (h) - \epsilon (h) \leq \tau (h), \hbox{ con } 0 \leq \epsilon < \tau /2,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

donde <math display="inline">\delta (h)</math> tampoco dependerá de <math display="inline">t</math>, <math display="inline">x</math> o <math display="inline">x'</math> y por lo tanto será constante dado un valor fijo para <math display="inline">h</math>. En consecuencia tendremos que <math display="inline">-\delta < \tilde{x}' - x' < \delta </math>, o lo que es lo mismo, <math display="inline">g(t,x) - \delta < \tilde{x}' < g(t,x) + \delta </math>. Por ello, al emplear el método de Euler-biseccion tendremos la expresión

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | \begin{align} \tilde{x}_0 & = x_0, \\
\tilde{x}_{n+1} & = \tilde{x}_n + h {\tilde{x}'}_n        \end{align}     
|}
|}

lo cual equivale a

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}• \tilde{x}_0 = x_0 + (\delta _{-1} = 0), \\        • \tilde{x}_{n+1} = \tilde{x}_n + h g(t_n,\tilde{x}_n) + (h \delta _n := \Delta _n), \hbox{ para } n = 0,1 \ldots N-1, \\        \end{array}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

con <math display="inline">\phi = g + \delta _n</math> y <math display="inline">0 \leq |\delta _n| < \delta \leq \tau </math>, y por ende

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0 \leq |\Delta _n| = h |\delta _n| <  h \delta \leq h \tau := \Delta{.}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

La anterior es una ecuación en diferencias como la que corresponde al método de Euler perturbado <span id='citeF-1'></span>[[#cite-1|[1]]], con una perturbación <math display="inline">h \delta _n = \Delta _n</math>, para <math display="inline">n = 0, \ldots, N-1</math>.

Ahora deduciremos el siguiente teorema.

<span id='theorem-x_biprima_Euler-biseccion'></span>Teorema 8: Empleando la fórmula del teorema ([[#theorem-y_biprima_Burden|2]]) y el teorema ([[#theorem-T.F.I|4]]) podemos expresar a <math display="inline">x''</math> como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x'' = \frac{-\frac{\partial f}{\partial t} -\frac{\partial f}{\partial x} x'}{\frac{\partial f}{\partial x'}}.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

Además, si existe una constante <math display="inline">M</math> con la propiedad <math display="inline">|x''| \leq M</math>, para alguna constate <math display="inline">M > 0</math> y <math display="inline">t \in [t_0,t_f]</math>, entonces

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>|x''| = | \frac{\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x} x'}{\frac{\partial f}{\partial x'}} | \leq M,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

para todo <math display="inline">t \in [t_0,t_f]</math>.

Observación 2: En el teorema anterior, en la expresión ([[#eq-23|23]]) aparece <math display="inline">x'</math>. Sin embargo, si dicha <math display="inline">x'</math> no se puede eliminar de la expresión por simplificación algebraica, entonces, una vez que ya se ha demostrado la unicidad en la solución de la ecuación, si el sistema físico nos da información extra respecto a un acotamiento para la velocidad (<math display="inline">x'</math>), entonces la ecuación modelo debe poseer el mismo acotamiento para <math display="inline">x'</math>, ya que la ecuación modelo y el sistema físico se corresponden en sus valores siempre que la ecuación con condiciones iniciales sea única.

Ahora enunciaremos el teorema de la cota del error para <math display="inline">\tilde{x}</math>

<span id='theorem-cota_error'></span>Teorema 9: Usando el teorema del acotamiento del error global (Burden et al., 1996, p. 248), la fórmula ([[#eq-23|23]]), el teorema ([[#theorem-Euler_perturbado|7]]) y el teorema del error global para el método de Euler perturbado (Burden et al., 1996, p. 251), podemos obtener la siguiente cota para el error, en su iteración n-ésima, de la solución numérica respecto a la exacta

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>|x_n - \tilde{x}_n| \leq \frac{1}{L} \left(\frac{h M}{2} + \frac{\Delta }{h} \right)(e^{L(t_n - t_0)} - 1) = \frac{1}{L} \left(\frac{h M}{2} + \tau (h) \right)(e^{L(t_n - t_0)} - 1),          </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

donde <math display="inline">L</math> es la constante de Lipschitz para <math display="inline">x' = g(t,x)</math>    Además, el valor óptimo para <math display="inline">h</math> (el mínimo error posible) se presenta para un <math display="inline">h</math> aproximado a <math display="inline">h = \sqrt{2 \Delta / M}</math>, o, despejando a <math display="inline">Delta</math>, cuando <math display="inline">\Delta = \frac{M}{2} h^2</math>, es decir, cuando

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = \frac{M}{2} h \approx \delta ,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

tomando entonces la expresión ([[#eq-19|19]]) la forma

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>|x_n - \tilde{x}_n| \leq \frac{M h}{L}  (e^{L(t_n - t_0)} - 1)     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

<span id='theorem-error_de_truncamiento_local'></span>Teorema 10: Usando la definición del error de truncamiento local (Burden et al., 1996, p. 254), el teorema ([[#theorem-Euler_perturbado|7]]) y suponiendo que <math display="inline">x'' \leq M</math> para <math display="inline">t \in [t_0, t_f]</math>, podemos observar que el error de truncamiento local para el método de Euler-bisección es <math display="inline">e_{n+1}(h) = \frac{x_{n+1}-x_{n}}{h} - {\tilde{x}'}_ == \frac{x_{n+1}-x_{n}}{h} - g(t_n,x_n) - \delta _n = \frac{h}{2} x''(\xi ) - \delta _n \leq \frac{h}{2} M - \delta _n</math>, por lo que al aplicar el valor absoluto nos quedará <math display="inline">|e_{n+1}(h)| = |\frac{h}{2} x''(\xi ) - \delta _n| \leq |\frac{h}{2} x''(\xi )| + |\delta _n| < \frac{h}{2} M + \delta \leq \frac{h}{2} M + \tau </math>. En conclusión, obtenemos el siguiente acotamiento para el error de truncamiento local

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0 \leq |e_{n+1}(h)| < \frac{h}{2} M + \tau ,  n = 0,1, \ldots , N-1     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

donde vemos que dicho acotamiento no depende de la iteración <math display="inline">(n+1)</math>-ésima en la que estemos, ni tampoco de <math display="inline">t</math>, <math display="inline">x</math> ni <math display="inline">x'</math>.

<span id='theorem-tau_lineal'></span>Teorema 11: Del teorema anterior vemos que el método de Euler-biseccion es de orden lineal o menor, dependiendo de cómo esté definida la función <math display="inline">\tau (h)</math>. Así, si <math display="inline">\tau </math> es constante, el método es de orden cero, y ello significa que el error <math display="inline">e</math> no cambia, sin importa que hagamos más pequeña a <math display="inline">h</math>. Por ello, lo conveniente es definir a

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = C \cdot h,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

de esa forma el método será por lo menos lineal. Además, es fácil ver que esta definición para <math display="inline">\tau </math> es compatible con la dada en el teorema ([[#theorem-cota_error|9]]) para el error mínimo, donde <math display="inline">C = M/2</math>, i.e.

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = \frac{M}{2} h  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

es el valor óptimo para la tolerancia, el cual nos dará, usando la fórmula ([[#theorem-error_de_truncamiento_local|10]]), un error mínimo de

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0 \leq |e_{n+1}(h)| < M h,  n = 1, 2, \ldots , N-1  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Observación 3: Las fórmulas ([[#eq-29|29]]) y ([[#eq-30|30]]) se corresponden bellamente con las fórmulas ([[#eq-25|25]]) y ([[#eq-26|26]]), respectivamente, si desarrollamos la exponencial en su serie de Maclaurin.

Observación 4: Ahora ya se puede comenzar a dar respuesta a la pregunta de la subsección pirmera del por qué es conveniente definir a <math display="inline">\tau = \tau (h)</math>.

====3.1.4 Análisis de estabilidad, consistencia y convergencia del método de Euler-bisección====

Ahora se hará el análisis de estabilidad, consistencia y convergencia para el método de Euler-bisección suponiendo que se cumplen todas las condiciones de la subsección anterior.

Teorema 12: Empleando la definición de consistencia de un método  de ecuación en diferencias (Burden et al., 1996, p. 312), la fórmula ([[#theorem-error_de_truncamiento_local|10]]), la definición de punto adherido o punto límite, el teorema de compresión (Lang, 1990, p. 41) (teorema del emparedado) y definiendo a <math display="inline">\tau </math> como <math display="inline">\tau = C h</math>, de acuerdo a ([[#theorem-tau_lineal|11]]), concluimos que <math display="inline">\lim _{h \rightarrow 0} |e_{n+1}(h)| = 0</math>, por lo que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lim _{h \rightarrow 0}  \max _{1 \leq n \leq N} |e_n(h)| = 0,  n = 1,2, \dots , N,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

y en consecuencia el método es consistente.

Teorema 13: Ahora, empleando la definición de convergencia para un método de ecuación en diferencias <span id='citeF-1'></span>[[#cite-1|[1]]], y la fórmula ([[#eq-26|26]]), obtendremos, de manera análoga al teorema anterior, <math display="inline">lim_{h \rightarrow 0} |x_n - \tilde{x}_n| = 0</math>, es decir

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lim _{h \rightarrow 0}  \max _{1 \leq n \leq N} |x_n - \tilde{x}_n| = 0,  n = 1,2, \dots , N,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

i.e. el método es convergente.

Teorema 14: Considerando ahora el teorema ([[#theorem-teorema5.19Burden|1]]), la subrutina de bisección y el problema perturbado mostrado en ([[#eq-20|20]]), entonces es obvio que <math display="inline">\phi (t_n,\tilde{x}_n; h) = biseccion(t_n,\tilde{x}_n; h, \alpha _i, \alpha _d) = biseccion(t_n,\tilde{x}_n; h) = g + \delta _n = {\tilde{x}'}_n</math>

<span id='theorem-Estabilidad'></span>Teorema 15:     Retomando el teorema anterior y empleando la condición de Lipschitz para <math display="inline">g</math>, así como la fórmula ([[#eq-21|21]]), y la expresión ([[#theorem-tau_lineal|11]]), tendremos para <math display="inline">\phi (t,x; h)</math> que <math display="inline">|\phi (t,x,h) - \phi (t,y,h)| = |\tilde{x}' - \tilde{y}'| = |g(t,x) + \tilde{\delta }(t,x,h) - g(t,y) - \tilde{\delta }(t,y,h)| \leq |g(t,x) - g(t,y)| +  |\tilde{\delta }(t,x,h) - \tilde{\delta }(t,y,h)| \leq L|x-y| + |\tilde{\delta }(t,x,h)| + |\tilde{\delta }(t,y,h)|< L|x-y| + 2 \tau = L|x-y| + 2 C h</math>, donde hemos hecho la transición <math display="inline">\delta _n \rightarrow \tilde{\delta }</math> para el caso continuo, resultando la siguiente desigualdad

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>|\phi (t,x,h) - \phi (t,y,h)| < L|x-y| + 2 C h,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

donde, del lado derecho de la última ecuación no podemos deshacernos del término <math display="inline">2 C h</math>. Sin embargo, debido a que ya se demostró que el problema está bien planteado si cumple el criterio dado en ([[#theorem-problema_bien_planteado|6]]), además de la consistencia y convergencia, entonces es de esperarse que para un <math display="inline">h</math> razonablemente pequeño, dicho término se puede tomar como equivalente a pequeños errores de redondeo de computadora, y por ende que el método se comporte como el de Euler ligeramente perturbado, siendo entonces un método estable.

<span id='theorem-PBP'></span>Observación 5:     Desafortunadamente, para este método no pude completar la prueba final que demostrara que <math display="inline">\phi </math> es una función que satisface la condición de Lipschitz para <math display="inline">x</math>, y con ello que el método es estable. No obstante dí los argumentos de que es un problema bien planteado para tratar de justificar esta última prueba.

==4 Análisis rápido de los demás métodos híbridos para el caso de una ecuación individual de orden uno==

Ahora abordaremos de manera rápida y práctica el problema de consistencia, convergencia y estabilidad de las demás combinaciones de métodos híbridos.

Consideremos primeramente que el análisis de la sección 3.1.2, referente al problema bien planteado se aplica de manera general para ecuaciones implícitas de la forma <math display="inline">f(t,x,x')</math>, sin importar qué combinación de métodos usemos.

En segundo lugar, en este artículo se omitirá el análisis del error de los demás métodos, dejándolo para futuros trabajos.

Y por último, respecto al análisis de estabilidad, consistencia y convergencia de los otros métodos, solo argumentaré algo análogo a lo dicho en el teorema ([[#theorem-Estabilidad|15]]) y en la observación ([[#theorem-PBP|5]]), en cuanto a que mis métodos híbridos se pueden considerar como pequeñas perturbaciones al problema de valor inicial (para <math display="inline">h</math> suficientemente pequeña) empleando los métodos convencionales de solución de ecuaciones diferenciales, como se observa en la fórmula ([[#eq-20|20]]), (donde la perturbación es <math display="inline">h \delta _n = \Delta </math>) pues solo se calcula de manera aproximada el valor de <math display="inline">x'</math> por métodos numéricos algebraicos, lo que nos va agregando un error extra en cada iteración. Mas como para todos esos métodos convencionales ya se han hecho los respectivos análisis, al tener un problema bien planteado significa que las soluciones dadas por mis métodos tendrán un pequeño error respecto a las soluciones mostradas por tales métodos estándar y por lo tanto también un pequeño error aceptable respecto a la solución exacta si <math display="inline">h</math> es muy pequeño.

==5 Análisis gráficos del método de Euler-bisección==

Ahora se procederá a hacer el análisis gráfico para el método de Euler-bisección, tomando como modelo el problema de caída libre con friccón laminar (bajas velocidades), por lo que la fuerza de fricción es proporcional a la velocidad y en sentido opuesto <span id='citeF-5'></span>[[#cite-5|[5]]].

Cabe señalar que aunque este modelo se puede resolver por métodos estándar e incluso analíticamente, será utilizado justamente para evaluar todo lo afirmado en las secciones anteriores.

La ecuación diferencial con condiciones iniciales ha de ser entonces

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}• f(t,v,v') = m v' + k v + m g= 0\\     • v(0) = v_0,    \end{array}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

donde <math display="inline">v</math> indica la velocidad del cuerpo, <math display="inline">v_0</math> es su velocidad inicial, <math display="inline">m</math> su masa, <math display="inline">g = 9.8 m/ s^2</math> es la aceleración de la gravedad en unidades S.I., <math display="inline">k</math> es la constante de fricción y hemos reducido el orden de la ecuación con el cambio de variable <math display="inline">v = x'</math>.

Para este problema en particular elegimos <math display="inline">v_0 = 0</math>, <math display="inline">m = 1</math>, y <math display="inline">k = .1</math>.

Además, el tiempo irá de <math display="inline">0 \leq t \leq 100</math>. También tomaremos un tamaño de paso de <math display="inline">h = .1</math>.

Ahora, en aras de la brevedad solo diremos que todas las condiciones para un problema bien planteado se cumplen (el lector interesado puede verificarlo), siendo la constante de Lipschitz <math display="inline">L = k/m</math> (para <math display="inline">v</math> en <math display="inline">g(t,v)</math>).

Y en lo referente al análisis del error, de acuerdo al sistema físico que modelamos la aceleración máxima se da al inicio, cuando la velocidad inicial es cero y por lo tanto la fuerza de fricción también es cero. Y dicha aceleración máxima ha de ser igual a <math display="inline">|g| = 9.8</math>.

Así que de acuerdo a la fórmula ([[#eq-23|23]]), podemo ver <math display="inline">|v''| = (k/m) |v'| \leq .98 = M</math>, siendo por lo tanto el valor del error óptimo <math display="inline">\tau = (M/2) h = .049</math> (según ([[#eq-25|25]])).

Con ello ya se pueden hacer los análisis respectivos de cotas del error y ver que el método tiene una convergencia de orden lineal, de acuerdo a ([[#eq-26|26]]).

A continuación se presentan los gráficos correspondientes, donde se puede observar que se cumple todo lo señalado.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Rodríguez_421804131-velocidad_y_aceleracion.png|480px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' 
|}

<span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]]

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Rodríguez_421804131-d_aceleracion_y_error.png|480px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' 
|}

<span id="fnc-2"></span>[[#fn-2|<sup>2</sup>]]

<span id="fn-1"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) El gráfico izquierdo es para la velocidad del cuerpo que cae, el derecho para su aceleración. El eje horizontal corresponde al tiempo. Se muestra en rojo la solución analítica y en azul la numérica; puede verse que prácticamente se traslapan.</span>

<span id="fn-2"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-2|<sup>2</sup>]]) En el gráfico izquierdo aparece la derivada de la aceleración, es decir <math>v''</math>, en el derecho el error global <math>\tilde{v} - v</math></span>

==6 Sistemas de ecuaciones diferenciales de orden <math>n</math>==

===6.1 Explicación breve del algoritmo===

En esta sección se explicará rápidamente cómo se desarrollan los métodos más generales de solución de sistemas de <math display="inline">m</math> ecuaciones diferenciales con <math display="inline">m</math> funciones incógnitas de <math display="inline">n</math>-ésimo orden. El respectivo análisis numérico se dejará para futuros trabajos.

Primeramente habrá de señalarse que todo sistema de ecuaciones diferenciales de orden <math display="inline">n</math> puede transformarse en un sistema equivalente de orden uno. Una vez transformado, se aplican métodos semejantes a los del caso de una sola ecuación de orden uno, como se verá a continuación.

Al considerar las condiciones iniciales, el sistema se transforma en un sistema algebraico para las variables incógnitas derivadas, el cual se puede resolver mediante métodos de resolución de sistemas de ecuaciones algebraicas, tales como punto fijo multivariable o Newton-Raphson multivariable. Una vez obtenidos los valores aproximados para las funciones derivadas, se procede a aplicar cualquier método estándar de solución numérica de sistemas de ecuaciones diferenciales, pero con un pequeño error en el cálculo de las funciones <math display="inline">g</math> para las derivadas (problema perturbado). Una vez aplicados los métodos estándar se obtiene el valor aproximado de las variables dependientes en <math display="inline">t_1</math>. Y el proceso se repite hasta llegar al <math display="inline">t_f</math>.

===6.2 Gráficos del método generalizado aplicado a las ecuaciones de campo relativistas de las estrellas de bosones===

La métrica para una estrella de bosones sin dependencia temporal (objeto compacto, esféricamente simétrico, estático) es <span id='citeF-6'></span>[[#cite-6|[6]]]

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>ds^2 = - \alpha (r)^2 dt^2+ a(r)^2 dr^2 + r^2 d\theta ^2 + r^2 \sin ^2(\theta ) d\varphi ^2,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

de la cual se deducen, usando el formalismo de la geometría diferencial, las siguientes ecuaciones de Einstein <span id='citeF-6'></span>[[#cite-6|[6]]]

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{a_{,r}}{a} + \frac{\alpha _{,r}}{\alpha } = \frac{k_{0}}{2} r \frac{a^2 \omega ^2 \phi _{0}^2}{\alpha ^2} + \frac{k_{0}}{2} r \phi _{0,r}^2,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{a_{,r}}{a} - \frac{\alpha _{,r}}{\alpha } = \frac{k_{0}}{2} r a^2 \phi _{0}^2 m^2 + \frac{k_{0}}{4} r \lambda a^2 \phi _{0}^4 + \frac{1-a^2}{r}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Además, elegimos un potencial de la forma <math display="inline">V = m^2 |\phi |^2 + \frac{\lambda }{2}|\phi |^4</math> <span id='citeF-7'></span><span id='citeF-8'></span>[[#cite-7|[7,8]]], siendo <math display="inline">\phi = \phi _0(r) e^{-i \omega t}</math>. Por otro lado, la ecuación de la conservación de la energía <math display="inline">\nabla _\nu T^\nu _\mu </math> se reduce a la ecuación de Klein-Gordon <math display="inline"> \left(\Box - \frac{dV}{d|\phi |^2} \right)\phi = 0</math>, de lo que se obtiene la tercera ecuación de Einstein

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi _{0,rr} + \phi _{0,r} \left(\frac{2}{r} + \frac{\alpha _{,r}}{\alpha } - \frac{a_{,r}}{a} \right)+ \omega ^2 \phi _{0} \frac{a^2}{\alpha ^2}  - a^2 (m^2 + \lambda \phi _{0}^2)  \phi _{0} = 0 \textrm{,}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

donde <math display="inline">a, \alpha </math> y <math display="inline">\phi _{0}</math> dependen solamente del radio <math display="inline">r</math>.

La solución numérica de este sistema usando mis métodos con las condiciones iniciales <math display="inline">a(0) = 1</math>, <math display="inline">\phi _0(0) = .1</math> y <math display="inline">\alpha (0) = 0.908710</math>, origina los siguientes gráficos

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Rodríguez_421804131-aBS.png|360px|Gráfico para a(r) vs r.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Gráfico para <math>a(r)</math> vs <math>r</math>.
|}

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Rodríguez_421804131-alphaBS.png|360px|Gráfico para α(r) vs r.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Gráfico para <math>\alpha (r)</math> vs <math>r</math>.
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Rodríguez_421804131-Phi0BS.png|360px|Gráfico para ϕ₀(r) vs r.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Gráfico para <math>\phi _0(r)</math> vs <math>r</math>.
|}

<span id="fnc-3"></span>[[#fn-3|<sup>1</sup>]]

<span id="fn-3"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-3|<sup>1</sup>]]) Estos gráficos se corresponden muy bien con los obtenidos usando métodos numéricos tradicionales, como se puede apreciar al compararlos con <span id='citeF-9'></span>[[#cite-9|[9]]]. </span>

==7 Conclusiones y futuros trabajos==

Los métodos numéricos alternativos para resolver ecuaciones diferenciales con condiciones iniciales, desarrollados en este trabajo, son útiles para los casos en que no se puede despejar explícitamente a la funciones derivadas. Además, en general pasaron todas las pruebas del análisis numérico (solo con ciertos detalles para la prueba de estabilidad) y se determinaron cotas del error, las cuales, los ejemplos aquí expuestos las satisfacen muy bien. En particular me enfoqué en el método de Euler-bisección, pero en futuros trabajos se desarrollarán otras combinaciones de métodos, así como los correspondientes métodos numéricos no estándar para ecuaciones diferenciales en derivadas parciales.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  R. L. Burden y J. D. Faires, ``Análisis Numérico", Segunda edición, Grupo Editorial Iberoamérica, Ciudad de México, México, (1996). <div id="cite-2"></div>
'''[2]'''  R. Gómez, D. Escobar, E. A. Gómez, J. L. Guerrero, J. M.M A. Gutiérrez, M. A. Hernández, A. P. Puerto y L. M. Sandoval, ``Elementos de métodos numéricos para ingeniería", Primera edición, Mc Graw Hill, Ciudad de México, México, (2002). <div id="cite-3"></div>
'''[3]'''  W. Rudin, ``Principios de Análisis Matemático", Tercera edición, McGraw Hill, Edo. de México, México, (1980). <div id="cite-4"></div>
'''[4]'''  S. Lang, ``Introducción al Análisis Matemático", Primera edición, Addison-Wesley Iberoaméricana, Ciudad de México, México, (1990). <div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  M. R. Spiegel, “Theoretical Mechanics”, Segunda edición, Schaum's Outline Series, McGraw-Hill, Estados Unidos de América, (1967).

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  C. Rodríguez, ``Solución a las ecuaciones de Einstein para las estrellas de bosones usando el método de álgebras y simetrías”, Tesis de Maestría, Universidad de Guadalajara, Ameca, México, (2017). <div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  F. S. Guzmán, ``The three dinamical fates of Boson Stars”, Revista Mexicana de Física, pp. 321-326, (2009). <div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  F. E. Schunck y D. F. Torres, ``Boson stars with generic self-interactions”, arXiv:gr-qc/9911038v1, (1999). <div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  A. Bernal, ``Estudio Dinámico de Campos Escalares Autogravitantes”, Tesis Doctoral, CINVESTAV, México, (2007).