Método iterativo lineal para la resolución de ecuaciones ordinarias no lineales

Elena V. Carreras Casanova^1*

¹Universidad Rey Juan Carlos, Madrid, España

^*Autor de correspondencia: e.carreras.2017@alumnos.urjc.es

Resumen
En este artículo se presenta un método iterativo lineal para la resolución de ecuaciones diferenciales no lineales. El método linealiza el problema diferencial para resolverlo iterativamente mediante el cálculo simbólico y/o numérico en el entorno de cálculo de MATLAB ®. Nos permite resolver sistemas diferenciales no lineales que son complejos de solucionar analíticamente. Se aplica exitosamente a la resolución de problemas de valor inicial y sistemas de ecuaciones diferenciales, incluyendo en este último grupo, sistemas que modelizan procesos epidemiológicos. Concretamente, simulamos el comportamiento de enfermedades de transmisión a través de la resolución de un modelo clásico y además, una epidemia real, la COVID-19.

Palabras clave
método iterativo lineal (MIL), métodos numéricos, cálculo simbólico

Introducción

Las ecuaciones diferenciales no lineales juegan un importante papel en la modelización de numerosos fenómenos o problemas de la vida real. Sistemas formados por estas ecuaciones a menudo son difíciles de resolver, pero si estudiamos el problema en detalle, la dificultad surge de la no linealidad de las ecuaciones. Así pues, se describe aquí el método iterativo lineal (MIL) para la resolución de ecuaciones y sistemas diferenciales no lineales.

El método plantea un original enfoque que linealiza el problema a través de una resolución iterativa mediante cálculos simbólicos y/o numéricos. Este método fue propuesto por Temimi y Ansari [1].

Por último, nos centramos en la resolución del modelo SIR que aplicamos a una pandemia real, la Covid-19. Se aplican de nuevo el método MIL y los métodos numéricos de Euler explícito (Euler) y Runge-Kutta de orden 4 (RK4) para predecir su comportamiento y se toman datos reales para contrastar los resultados obtenidos.

1. Método iterativo lineal (MIL)

La idea fundamental del método es la siguiente. Podemos reescribir cualquier ecuación diferencial ordinaria [1] como:

$L (u (t)) + N (u (t)) + g (t) = 0, A (u, u^{'}, u^{''}, u^{'''}, \dots), (1)$

donde t es la variable independiente y u la variable dependiente. Para ello, hay que definir:

◾ Operador lineal L(t): términos lineales de la ecuación diferencial.

◾ Operador N(t): términos no lineales de la ecuación diferencial.

◾ Función g(t): términos no dependientes de la variable u.

◾ Operador auxiliar A: condiciones iniciales o de contorno.

El procedimiento iterativo se basa en la resolución de sistemas sucesivos que denotaremos DS_n [2], donde n ∈ {0, 1, ..., } indica el número de iteración actual, siendo M el número total de iteraciones realizadas. Empezamos suponiendo que u₀ (t) es solución de la ecuación (1) y además, lo consideramos como estimacion inicial de la solucion. El sistema DS₀ se construye sin considerar los términos no lineales de la ecuación diferencial, es decir, sin incluir el operador N(t):

$D S_{0} := {\begin{cases} L (u_{0} (t)) + g (t) = 0 \\ A (u_{0}). \end{cases}, (2)$

Una vez resuelto, el proceso sigue para n = {1, ..., M} incorporando los términos no lineales evaluados con la solución obtenida en la iteración anterior n − 1:

$D S_{n} := {\begin{cases} L (u_{n} (t)) + g (t) + N (u_{n - 1} (t)) = 0 \\ A (u_{n}), \end{cases}, (3)$

siendo u_n(t) la solución calculada en la iteración n. De esta forma, cada vez que se resuelve un sistema, se está resolviendo un problema diferencial lineal. El sistema DS₀ por definición no incluye términos no lineales y los sucesivos sistemas DS_n , pese a añadirlos, los consideran conocidos. Por tanto, no suponen una complejidad añadida. El conjunto de sistemas DS₀ y DS_n puede resolverse mediante el cálculo simbólico y/o el cálculo numérico. En primer lugar, se opta por el cálculo simbólico ya que permite obtener una expresión analítica de la solución. Sin embargo, a medida que aumenta el número de iteraciones, la resolución simbólica es más compleja hasta llegar al punto de saturar las capacidades del entorno de cálculo. Cuando un problema diferencial no puede ser resuelto simbólicamente o la resolución presente limitaciones entonces, recurrimos al cálculo numérico, capaz de hallar soluciones aproximadas.

Se procede mediante el siguiente problema diferencial a exponer el funcionamiento del método MIL y comentar la calidad de sus resultados. Este se va a resolver simbólicamente (MIL-Sim) y numéricamente con los métodos de Euler explícito (MIL-Euler) y Runge-Kutta de orden 4 (MIL-RK4). Además, se establece un criterio de parada que define el número total de iteraciones M a realizar para hallar la solución con la exactitud deseada.

1.1 Problema ilustrativo del método MIL

El problema diferencial propuesto está compuesto por una ecuación diferencial no lineal de segundo orden complementado con dos condiciones iniciales.

${\begin{cases} Ω u^{''} + B u^{'} + k u^{2} = 0, 0 < t < 5, \\ u (0) = 0, \\ u^{'} (0) = 1, \end{cases} (4)$

siendo Ω, B, k constantes conocidas. El primer paso es definir los operadores propios del método MIL:

$L (u) = Ω u^{''} + B u^{'}, N (u) = k u^{2}, g (t) = 0, A (u, u^{'}) = [u (0) = 0, u^{'} (0) = 1] . (5)$

Planteamos los sistemas a resolver en cada iteración siguiendo los esquemas (2) y (3):

$D S_{0} := {\begin{cases} L (u_{0} (t)) + g (t) = Ω u_{0}^{''} + B u_{0}^{'} = 0, 0 < t < 5, \\ A (u_{0}, u_{0}^{'}) = [u_{0} (0) = 0, u_{0}^{'} (0) = 1] . \end{cases} (6)$

$D S_{n} := {\begin{cases} L (u_{n} (t)) + N (u_{n - 1} (t)) + g (t) = Ω u_{n}^{''} + B u_{n}^{'} + k u_{n - 1}^{2} = 0, 0 < t < 5, \\ A (u_{n}, u_{n}^{'}) = [u_{n} (0) = 0, u_{n}^{'} (0) = 1] . \end{cases} (7)$

1.1.1 Resolución simbólica con MIL (MIL-Sim)

En primer lugar, se trata de resolver el problema (4) simbólicamente, sin embargo, MATLAB ® no es capaz de hallar una expresión analítica de la solución debido a la no linealidad de uno de los términos de la ecuación. El método MIL propone linealizar el problema para esquivar así el obstáculo encontrado y poder emplear de nuevo el cálculo simbólico de MATLAB ® en cada iteración n.

Criterio de parada. El criterio de parada consiste en fijar un número máximo de iteraciones debido al elevado tiempo de cómputo empleado por MATLAB ® al resolver los sistemas DS_n. Para el sistema (4) se determina alcanzar un máximo de M = 5 iteraciones, es decir, resolver los problemas DS₀, DS₁ , DS₂, DS₃, DS₄ y DS₅. El cálculo simbólico falla para n = 6 debido a la complejidad de la ecuación generada (su segundo miembro) a pesar de resolver una ecuación lineal.

Resultados por MIL-Sim Las curvas de las soluciones pueden verse en la Figura 1, que muestra las soluciones simbólicas u₅ y v₅ . La gran ventaja de esta implementación consiste en obtener para cada iteración n una expresión simbólica de la solución del sistema. Es decir, tras resolver los sistemas DS_n conocemos la expresión de la solución analítica del sistema (4). Por contra, el principal inconveniente es la limitación de iteraciones del método debido a la complejidad de la solución analítica tras algunas iteraciones. Este factor no nos permite calcular la solución simbólica hasta alcanzar la convergencia.

Una primera conclusión es que MIL-Sim ha permitido obtener una solución cuando el cálculo simbólico no puede de forma directa. Por esta razón, el método MIL se posiciona como un recurso de gran potencial que nos permite obtener soluciones analíticas. Esta resolución podría ser suficiente pero por falta de capacidad del software, nos hemos encontrado con limitaciones en cuanto al número de iteraciones. Se espera poder solventar este problema en el futuro con el aumento de las capacidades informáticas de nuestros equipos.

1.1.2 Resolución numérica con MIL (MIL-Euler, MIL-RK4)

Sea [0, 5] el intervalo donde queremos conocer la solución, el primer paso es discretizar el dominio obteniendo un conjunto discreto de puntos t₁ , t₂, ..., t_I ∈ [0, 5] en los que se calcula una solución aproximada de los sistemas (6) y (7). Para ello se fija un paso de discretización y así se obtiene un total de $I = \frac{5 - 0}{h} + 1$ puntos del soporte ti . Seguidamente, para aplicar los esquemas numéricos de Euler explícito y Runge-Kutta se debe reducir el orden de DS₀ y DS_n para conseguir un sistema de dos ecuaciones diferenciales de primer orden, que son lineales tanto en u₀ y u_n. Realizamos el cambio de variable u′(t) = v(t). Los sistemas anteriores (6) y (7) se reescriben como, 0 < t < 5:

$D S_{0} := {\begin{cases} u_{0}^{'} (t) = v_{0} (t) = f_{u_{0}} (t, u_{0}, v_{0}), \\ u_{0} (0) = 0, \\ v_{0}^{'} (t) = f_{v_{0}} (t, u_{0}, v_{0}) = \frac{- B}{Ω} v_{0} (t), \\ v_{0} (0) = 1 . \end{cases} D S_{n} := {\begin{cases} u_{n}^{'} (t) = v_{n} (t) = f_{u_{n}} (t, u_{n}, v_{n}), \\ u_{n} (0) = 0, \\ v_{n}^{'} (t) = f_{v_{n}} (t, u_{n}, v_{n}) = \frac{- B v_{n} (t) - k u_{n - 1}^{2}}{Ω}, \\ v_{n} (0) = 1 . \end{cases} (8)$

Estos sistemas diferenciales pueden escribirse de forma vectorial. Definimos:

$V_{n} (t) = {(u_{n} (t), v_{n} (t))}^{t}, F_{n} (t, V_{n}) = {(f_{u_{n}} (t, u_{n}, v_{n}), f_{v_{n}} (t, u_{n}, v_{n}))}^{t}, A = (0, 1)^{t}, (9)$

por tanto,

${\begin{cases} V_{0}^{'} (t) = F_{0} (t, V_{0}), \\ V_{0} (0) = A . \end{cases} {\begin{cases} V_{n}^{'} (t) = F_{n} (t, V_{n}), & 0 \leq n \leq M \\ V_{n} (0) = A . \end{cases} (10)$

Además, el proceso de discretización conduce a la siguiente formulación:

${\begin{cases} V_{0, i}^{'} (t) = F_{0} (t_{i}, V_{0, i}), \\ V_{0, 1} (0) = A . \end{cases} {\begin{cases} V_{n, i}^{'} (t) = F_{n} (t_{i}, V_{n, i}) \\ V_{n, 1} (0) = A . & 1 \leq i \leq I - 1 . \end{cases} (11)$

Criterio de parada. Para el cáluculo numérico, el criterio de parada viene dado por una tolerancia ϵ. El proceso termina cuando se cumplan a la vez las dos expresiones siguientes:

$E_{v_{n}}^{L_{2}} = \sqrt{\sum_{i = 1}^{I} {(v_{n, i} - v_{n - 1, i})}^{2}} \leq ϵ, E_{u_{n}}^{L_{2}} = \sqrt{\sum_{i = 1}^{I} {(u_{n, i} - u_{n - 1, i})}^{2}} \leq ϵ . (12)$

donde el error cometido en cada iteración aplicando la norma L₂ entre dos soluciones consecutivas es el siguiente, por una parte entre v_n,i y v_n_−1,_i y por otra entre u_n,i y u_n_−1,_i en cada punto i de los puntos de discretización t_i para i = 1, ..., I y para todo n = 1, ..., M [2].

Resultados por MIL-Euler y MIL-RK4 Utilizando un paso de discretización h = 0,01 se ha resuelto el problema diferencial para diferentes tolerancias ϵ. Veánse las soluciones halladas en la Figura 1, que representa las soluciones obtenidas en la última iteración relizada utilizando una tolerancia de ϵ = 1 · 10⁻¹⁵ . En el cuadro 1 se incluye el número de iteraciones M requeridas para cada método. Como se puede observar, al reducir el valor de tolerancia el

Cuadro 1. Resultados MIL-Euler y MIL-RK4

ϵ	MIL-Euler	MIL-RK4
1·10⁻⁵	M=13	M=3
1 · 10⁻¹⁵	M=21	M=3

número de iteraciones no varía para MIL-RK4. Por el contrario, para MIL-Euler, al endurecer el criterio de parada es necesario iterar más veces. En definitiva, el método MIL-RK4 es más rápido y converge con más velocidad, lo cuál es lógico debido a que éste último método es de orden más elevado.

Ahora bien, ¿qué ganamos respecto a los métodos numéricos directos? De haber resuelto un problema diferencial complejo, podríamos haber topado con el mismo obstáculo que con la resolución simbólica. Es decir, la complejidad de las ecuaciones no lineales del sistema no permite hallar una solución aproximada del problema mediante los métodos de Euler y RK4. En tal caso, disponemos del recurso del método MIL. Permite obtener soluciones aproximadas habiendo linealizado el problema por lo que su resolución no se ve amenazada por las no linealidades del mismo.

1.2 Solución simbólica vs. Solución numérica

MIL-Sim, MIL-Euler y MIL-RK4 hallan soluciones cercanas. Estas pueden compararse en la Figura 1. Para MIL-Sim y MIL-RK4 se observa la mayor coincidencia, no obstante, pese a no presentar limitaciones computacionales, este último es capaz de hallar en un número menor de iteraciones la solución.

De tratar un problema diferencial más complejo, MIL-Sim se vería condicionado por las capacidades del software y se limitaría a un número más bajo de iteraciones sin posibilidad de calcular una solución simbólica hasta alcanzar el grado de exactitud deseado. Por contra, la resolución numérica aplicada al MIL no se vería afectada por estas restricciones.

2. Simulación de una epidemia real mediante el método MIL

Es posible simular y comparar la evolución de una epidemia mediante el modelo SIR [3]. El modelo tiene dos grandes utilidades. En primer lugar, permite simular la evolución de una determinada enfermedad infecciosa y comparar diversos escenarios modificando los valores de los parámetros: tasa de infección β, tasa de recuperación γ, tamaño de la población P y número de infectados en el inicio I₀. En segundo lugar, es posible aplicar el modelo SIR a datos reales procedentes de un determinado país o zona, lo que permite obtener estimaciones sobre los anteriores parámetros de gran ayuda para los epidemiológos.

2.1 Resolución numérica de un brote epidémico

Centrémonos en el segundo objetivo y empleemos los métodos de Euler, RK4 y el método MIL para la resolución del modelo SIR:

${\begin{cases} \frac{d S (t)}{d t} = - β S (t) \frac{I (t)}{P}, \\ \frac{d I (t)}{d t} = β S (t) \frac{I (t)}{P} - γ I (t), t_{0} < t < T \\ \frac{d R (t)}{d t} = γ I (t) . \end{cases} (13)$

Para ello, se han utilizado los datos diarios disponibles sobre el número de infectados por la COVID-19 en España¹ de los días comprendidos entre el 01-03-2020 y el 30-06-2020. Consultamos la tabla de datos reales el día 01-03-2020 (t = 0) para conocer el tamaño de la población española P = 46937060 y el número de infectados I₀ = 28. Las condiciones iniciales, por tanto, son:

$S (0) = S_{0} = P - I_{0} = 46937032, I (0) = I_{0} = 28, R (0) = 0 . (14)$

La evolución de la pandemia modelada con los parámetros (14) y el paso de discretización h = 0,1 en el intervalo [t₀, T] = [0, 120] es tal y como se muestra en la Figuras 2 y 3. Tanto MIL-Euler como MIL-RK4 han hallado en M = 54 iteraciones las soluciones anteriores. Pese a ser un número elevado de iteraciones, el método es rápido y el tiempo de ejecución es muy bajo, de orden 1 · 10⁻¹ segundos. Además, la diferencia de tiempo de ejecución respecto a Euler y RK4 es de orden 1 · 10⁻² segundos. Los resultados son obtenidos ejecutando los códigos implementados en un ordenador de 8 GB RAM con un sistema operativo Windows de 64 bits trabajando con un Intel ® Core(TM) i5-8250U CPU @1.60GHz 1.80 GHz.

Comparemos ahora con los datos reales registrados, en primer lugar con el número de nuevos casos confirmados por día. En la Figura 4a se representa el número de nuevos contagios registrados en el mismo periodo de tiempo. El máximo número de nuevas infecciones confirmadas se produce en t = 27, día 27 de marzo, y es igual a 9181. El método de Euler y MIL-Euler calculan que el mayor número de casos nuevos es igual a 2708598,6021, un valor muy elevado. No obstante, este hecho se produce también en t = 27. Por otra parte, el método de Runge-Kutta y MIL-RK4 estiman un total de 2683923,0293 nuevos casos para t = 28, es decir, un día más tarde.

A continuación, estudiamos en detalle el comportamiento de la enfermedad las dos primeras semanas en las que no había restricciones sanitarias. En la Figura 4b se representan en el intervalo [t₀, T] = [0, 14] los datos reales y las soluciones numéricas obtenidas por el conjunto de métodos. Se observa que inicialmente los resultados de la simulación se acercan a los reales. Sin embargo, a partir del décimo día (t = 10) el número de infectados en la simulación se dispara mientras que la curva de los datos reales sigue creciendo moderadamente.

En definitiva, el modelo predice un crecimiento de la enfermedad más exponencial que el realmente producido. En primer lugar, esta diferencia podría deberse a que los datos empleados únicamente considera los casos confirmados de Covid-19 mediante prueba PCR. De hecho, se estima que el número real de infectados fuera sumamente mayor que el registrado. En segundo lugar, podría deberse a las distintas medidas de confinamiento y contención que se produjeron en España desde el 14 de marzo, mesuras no contempladas por el modelo utilizado.

3. Conclusiones

Se ha presentado un método iterativo lineal para la resolución de ecuaciones diferenciales no lineales. Su gran ventaja es que se aplica sin restricciones de términos no lineales. Esto permite emplear tanto el cálculo simbólico como el cálculo numérico en la resolución. Mediante el cálculo simbólico, una vez linealizado el problema, es posible resolver la ecuación o sistema analíticamente. Sin embargo, a medida que el número de iteraciones aumenta, la resolución simbólica es más compleja hasta llegar al punto de saturar las capacidades del cálculo simbólico. De incrementar el número de iteraciones, el método MIL permitiría obtener una solución correcta en forma analítica. Por otra parte, los resultados numéricos del método MIL muestran que el método es efectivo, preciso y convergente hacia la solución en un número bajo de iteraciones. En caso de no conocer métodos analíticos para la resolución de problemas diferenciales o de ser estos más complejos debido a la no linealidad de sus términos, el método MIL puede hacer frente a estas dificultades.

En conclusión, el método MIL, capaz de obtener resultados fiables y precisos, se caracteriza por su sencillez de aplicación y es útil en la resolución de ecuaciones complejas que no admiten una solución analítica y que además puedan presentar problemas en una resolución numérica directa. Como futura linea de trabajo se desea profundizar en la resolución simbólica del método MIL.

Agradecimientos

Me gustaría agradecer a mi tutor Christian Vanhille la oportunidad brindada con esta publicación y su ayuda durante todo el proyecto.

Referencias

¹ H. Temimi y A. Ansari, “A computational iterative method for solving nonlinear ordinary differential equations”, LMS Journal of Computation and mathematics 18, 730–753 (2015).

² C. Vanhille, “Linear Iterative Procedure to Solve a Rayleigh–Plesset Equation”, Acoustics 3, 212–220 (2021).

³ V. Martın Barroso, “Una breve introducción al modelo SIR aplicado al caso del Covid-19”, Instituto Complutense de Estudios Internacionales (2020).

_______________________________

¹ Los datos están disponibles en https://www.ecdc.europa.eu/en/publications-data

Forum Docentis

Buscador