# Due to numerical inaccuracy, the solution of a boundary value problems becomes negative

I treat a toy example to get my point across. In reality I have to deal with a much more complex model.

Let us consider a one dimensional boundary value problem using the bvp5c solver in Matlab. Two liquids enter at different points, move in opposite directions and react with each other. Liquid 1 enters at $x = 0$ and moves in the positive direction. Liquid 2 enters at $x = 1$ and moves in the negative direction. We are interested in the equilibrium distribution of the two concentrations over the domain $[0,1]$.

Let $c_1(x)$ be the concentration of liquid 1 at the point $x$ and similarly for liquid two. The differential equations that need to be solved are the following

\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}

with $K$ a constant that we can choose. Note that the $c_1(x)$ goes down for increasing $x$ and $c_2(x)$ goes down for decreasing $x$, which is exactly as we expect. The two boundary values are $c_1(0) = 10$ and $c_2(1) = 6$.

In Matlab, one models this as a boundary value problem with two boundaries and two components. By solving this problem for $K = 1$ we obtain a nice solution shown in the following figure.

However, by increasing $K$ to, say, $K = 10^5$ we run into a numerical inaccuracy, see this figure. I also receive a warning message from Matlab, see the bottom of this post.

The concentration $c_2(x)$ drops below zero, while that should not be possible according to the second differential equation.

My intuition is that this happens because Matlab takes discrete sizes of $dx$ and computes the derivative in a point (which is large) and updates the solution, making the concentration negative. One could combat this phenomenon by increasing the maximum allowed mesh size, but this is not a solution that can be incorporated in my actual model due to the computation time and memory issues.

How to make sure that both concentrations are nonnegative?

Matlab code

The script that runs the bvp5c solver

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);


The function with the ODEs. The dependence on $x$ is implied in the input of the function.

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


The function with boundary values. When there are two boundary values, Matlab refers to them as $a$ and $b$, which are $x = 0$ and $x = 1$, respectively.

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.


Error message for $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

• This special problem could be solved by setting $c_2(x)=\exp(u_2(x))$. Consider also that with large K the ODE is stiff, necessitating a solver for stiff problems. – LutzL May 8 '14 at 13:12
• Also, the parameter x = linspace(0,1,100) does not prescribe the used step size in the ODE solver, that is probably dynamically adapted. x only gives the sample points where the y values of the solution are reported. – LutzL May 8 '14 at 13:15
• @LutzL I am aware that is the case. I was referring to actually increasing the 'NMax' option. You mean by introducing a different variable $u_2(x)$ that is allowed to take any value in $\mathbb{R}$? What do you mean by a 'stiff' ODE? Is bvp5c suited for 'stiff' problems? – Ritz May 8 '14 at 13:35
• Have you tried adjusting the absolute and relative integration tolerances using bvpset rather than using the defaults? – horchler May 8 '14 at 14:57
• Things like warning messages and errors should always be mentioned in the question , not left for us to discover if we run your code ourselves. If you can convert this to an initial value problem and use ode45, you could use odeset to specify the 'NonNegative' option. – horchler May 8 '14 at 15:10

The solver uses a collocation method, that is, at each mesh point it treats the discretized differential equation as an algebraic equation (for function value and a number of derivatives) and then tries to find a numerical solution of the BVP by solving this large non-linear system. After that, the error between the mesh points is estimated and the mesh refined towards a more homogeneous error estimate.

The point now is that the discretized differential equation will have more than one solution. Typically, the other solutions are not smooth, have large jumps between them, so that they are easily detectable by the solver. Or more often, the initial approximation is much closer to the smooth solution than to any of the extraneous solutions.

However, for stiff systems that might no longer be the case. "stiff" means that you have fast and slow components. The fast components require a very small step size to be modeled correctly in the discretization, while the many steps that that requires might lead to an undue accumulation of floating point errors in the other components. So it requires extra effort to find a workable balance.

In the collocation model, the large derivatives of the fast components lead to the situation where even the sampled exact solution is visually not smooth. That is, even some of the normal steps of the exact solution look like big jumps. So the true numerical solution is no longer easily distinguishable from the extraneous jumpy solutions of the algebraic system.

In the given system, the exact solution for large $K$ jumps at $x=1$ from essentially $c_2(x)\approx 0$ to $c_2(1)$ in the last part of the interval of length about $1/\sqrt{K}$. The square root is arbitrary, but captures the essence of the situation.

Close to $x=1$ the slow component $c_1$ can be treated as a constant, so that $c_2$ behaves qualitatively like an exponential function $$c_2(x)=c_2(1)\,e^{K\,c_1(1)\,(x-1)}=c_2(1)\,e^{-K\,c_1(1)\,(1-x)}$$ This situation is captured by $e^{-Kt}$ close to $t=0$, i.e., using the substitution $t=c_1(1)(1-x)$. For $K$ large the function $e^{-Kt}$ goes from the appreciable value $1$ at $t=0$ to the very small value $e^{-\sqrt{K}}$ only a very small distance away at $t=1/\sqrt{K}$. Even more, if one wants the steps in the function values of $e^{-K\cdot n\,dt}$ to be small, one would need $K\,dt$ to be small, for instance by choosing $dt=K^{-3/2}$.

The mesh refinement algorithm would have to detect this situation. However, without telling the algorithm about this situation, it would have to create a sizable number, up to $O(K)$, of mesh points with a distance of $K^{-1-\epsilon}$ close to $x=1$ to zero in to the singularity. For $K=10^5$ this would requite mesh points in the order of some thousands or ten thousands, which will likely exceed any moderate size of Nmax.

• Thank you for the elaborate explanation of the collocation method, stiff problems, and what is happening at the boundary $x = 1$. Am I right by stating that, in essence, I should define the concentrations $c_i(x) := f(u_i(x))$ with $f : \mathbb{R} \to [0,+\infty)$ and thus let $u_i(x) \in \mathbb{R}$? – Ritz May 8 '14 at 15:22
• Yes, if possible. Another problem, apart from the almost vertical movement close to $x=1$, is numerical underflow in rest of the interval. Close to $x=0$ the values of $c_2(x)$ are of the size $e^{-10^{4}}$ or smaller, which can only be represented by $0$. In transforming the equations you have to achieve that the symbolic cancellation of these small terms is possible. If, for example, the exponentials in $e^{u(x)}u'(x)=f(e^{u(x)})$ resp. $u'(x)=e^{-u(x)}f(e^{u(x)})$ do not cancel, nothing has been gained. – LutzL May 8 '14 at 16:09