Debido a la imprecisión numérica, la solución de un problema de valor límite se vuelve negativa

Question

Trato un ejemplo de juguete para entender mi punto de vista. En la realidad tengo que lidiar con un modelo mucho más complejo.

Consideremos un problema de valor límite unidimensional utilizando el solucionador bvp5c en Matlab. Dos líquidos entran por puntos diferentes, se mueven en direcciones opuestas y reaccionan entre sí. El líquido 1 entra por $x = 0$ y se mueve en dirección positiva. El líquido 2 entra en $x = 1$ y se mueve en sentido negativo. Nos interesa la distribución de equilibrio de las dos concentraciones sobre el dominio $[0,1]$ .

Sea $c_1(x)$ sea la concentración del líquido 1 en el punto $x$ y lo mismo para el líquido dos. Las ecuaciones diferenciales que hay que resolver son las siguientes

\begin{align} \frac{dc_1(x)}{dx} &= -(c_1(x) + c_2(x)), \\ \frac{dc_2(x)}{dx} &= K c_1(x) c_2(x), \end{align}

con $K$ una constante que podemos elegir. Obsérvese que la $c_1(x)$ baja por aumentando $x$ et $c_2(x)$ baja por disminuyendo $x$ que es exactamente lo que esperábamos. Los dos valores límite son $c_1(0) = 10$ et $c_2(1) = 6$ .

En Matlab, esto se modela como un problema de valor límite con dos límites y dos componentes. Resolviendo este problema para $K = 1$ obtenemos una bonita solución que se muestra en la siguiente figura.

Sin embargo, al aumentar $K$ digamos, $K = 10^5$ nos encontramos con una imprecisión numérica, véase esta figura. También recibo un mensaje de advertencia de Matlab, ver la parte inferior de este post.

La concentración $c_2(x)$ cae por debajo de cero, mientras que eso no debería ser posible según la segunda ecuación diferencial.

Mi intuición es que esto ocurre porque Matlab toma tamaños discretos de $dx$ y calcula la derivada en un punto (que es grande) y actualiza la solución, haciendo que la concentración sea negativa. Se podría combatir este fenómeno aumentando el tamaño máximo de malla permitido, pero no es una solución que se pueda incorporar a mi modelo real debido a los problemas de tiempo de cálculo y de memoria.

Cómo asegurarse de que ambas concentraciones no negativo ?

---

Código Matlab

El script que ejecuta el solucionador bvp5c

x = linspace(0,1,100); % the mesh
yinit = [ 5 0 ]; % an initial guess for the solution for both liquids
K = 1; % the factor K in the second differential equation

solinit = bvpinit(x,yinit); % setting the mesh and the initial guess 

% Calling the bvp solver. First argument are the ODEs, the second argument 
% are the boundary conditions at x = 0 (for liquid 1) and at x = 1 (for
% liquid 2), the last argument specifies the mesh and initial guess
sol = bvp5c(@(x,c)Concentrations(x,c,K),@(ca,cb)BoundaryCond(ca,cb),solinit);

La función con las ODE. La dependencia de $x$ está implícito en la entrada de la función.

function dcdx = Concentrations(x,c,K)

dcdx(1) = -(c(1) + c(2)); % differential equation for liquid 1
dcdx(2) = K*c(1)*c(2); % differential equation for liquid 2

La función con valores límite. Cuando hay dos valores límite, Matlab se refiere a ellos como $a$ et $b$ que son $x = 0$ et $x = 1$ respectivamente.

function res = BoundaryCond(ca,cb)

res = [ ca(1)-10; cb(2)-6 ]; % the concentration of liquid 1 at boundary a
is 10 and the concentration of liquid 2 at boundary b is 6.

---

Mensaje de error para $K = 10^5$

Warning: Unable to meet the tolerance without using more than 5000 mesh points.
 The last mesh of 892 points and the solution are available in the output argument.
 The maximum error is 536.678, while requested accuracy is 0.001. 
> In bvp5c at 339

Answer 1

El solucionador utiliza un método de colocación, es decir, en cada intervalo de la malla discretiza la ecuación diferencial aproximando la solución mediante un polinomio de grado fijo y añadiendo el número necesario de ecuaciones para igualar el número de parámetros libres exigiendo la igualdad exacta de la EDO en algunos puntos interiores del intervalo. Para la EDO polinómica dada, esto da lugar a ecuaciones algebraicas para los coeficientes del polinomio. El solucionador procede entonces a encontrar una solución numérica del BVP resolviendo este gran sistema no lineal. A continuación, se estima el error entre los puntos de la malla y se refina la malla hacia una estimación de error más homogénea.

La cuestión ahora es que la ecuación diferencial discretizada tendrá más de una solución. Normalmente, las otras soluciones no son suaves, tienen grandes saltos entre ellas, de modo que son fácilmente detectables por el solucionador. O más a menudo, la aproximación inicial es mucho más cercana a la solución suave que a cualquiera de las soluciones extrañas.

Sin embargo, en los sistemas rígidos puede que ya no sea así. "Rígido" significa que hay componentes rápidos y lentos. Los componentes rápidos requieren un tamaño de paso muy pequeño para modelarse correctamente en la discretización, mientras que los muchos pasos que eso requiere podrían dar lugar a una acumulación indebida de errores de coma flotante en los demás componentes. Por tanto, es necesario un esfuerzo adicional para encontrar un equilibrio viable.

En el modelo de colocación, las grandes derivadas de los componentes rápidos conducen a una situación en la que incluso la solución exacta muestreada no es visualmente suave. Es decir, incluso algunos de los pasos normales de la solución exacta parecen grandes saltos. Por lo tanto, la solución numérica real ya no se distingue fácilmente de las soluciones con saltos extraños del sistema algebraico.

---

En el sistema dado, la solución exacta para grandes $K$ salta a $x=1$ de esencialmente $c_2(x)\approx 0$ a $c_2(1)$ en la última parte del intervalo de longitud alrededor de $\Delta t=1/\sqrt{K}$ . La raíz cuadrada es arbitraria para $1\gg\Delta t\gg1/K$ pero capta la esencia de la situación.

Cerca de $x=1$ (no puede haber capas límite en otros puntos del intervalo) la componente lenta $c_1$ puede tratarse como una constante, de modo que $c_2$ se comporta cualitativamente como una función exponencial $$ c_2(x)=c_2(1)\,e^{K\,c_1(1)\,(x-1)}=c_2(1)\,e^{-K\,c_1(1)\,(1-x)} $$ Esta situación se recoge en $e^{-Kt}$ cerca de $t=0$ es decir, utilizando la sustitución $t=c_1(1)(1-x)$ . Para $K$ grande la función $e^{-Kt}$ pasa del valor apreciable $1$ en $t=0$ al valor muy pequeño $e^{-\sqrt{K}}$ a muy poca distancia en $t=1/\sqrt{K}$ . Más aún, si se quieren los pasos en los valores de la función de $e^{-K\cdot n\,dt}$ para ser pequeño, se necesitaría $K\,dt$ ser pequeño, $1\gg K\,dt\gg 1/K$ por ejemplo, eligiendo $dt=K^{-3/2}$ .

---

El algoritmo de refinamiento de la malla tendría que detectar esta situación. Sin embargo, sin informar al algoritmo de esta situación (por ejemplo, estructurando la malla inicial de esta manera), tendría que crear un número considerable, hasta $O(K)$ de puntos de malla con una distancia de $K^{-1-\epsilon}$ cerca de $x=1$ para llegar a la singularidad. En $K=10^5$ esto requeriría puntos de malla del orden de algunos miles o diez miles, lo que probablemente superará cualquier tamaño moderado de Nmax.