Solving numerically a strongly stiff nonlinear ODE system with ill-conditioned Jacobian

Using Matlab, I am trying to solve numerically the following nonlinear system of ODEs:

\begin{aligned} \dot B &= -\alpha B -\nu BV \\ \dot X &= A-\mu_1 X -c E(B)VX \\ \dot Y &= -\mu_2 Y +c E(B)VX \\ \dot V &= kY - (\gamma B + E_0)V - \delta V \\ \end{aligned}

with

$$E(B)=(\gamma B +E_0)e^{-\beta \gamma B}$$

The system tries to capture the dynamics of how the body of an infant reacts to the presence of the dengue virus.

• $$V$$ is the number of virus.

• $$B$$ the number of antibody.

• $$X$$ the number of healthy cells.

• $$Y$$ the number of infected cells.

• $$E$$ is a function enhancing the effect of antibodies, depending on the number of antibodies.

Implementing this model and trying different solver, I noticed that solvers ode15s and ode23s are performing way better than ode45. Hence, I deduced that my problem was stiff. But with certain set of parameters I had the following warning from Matlab:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.334107e-017.

I worried for a while but I found the malicious ill-conditionned matrix. It's actually that solvers ode23s and ode15s use the Jacobian of the system. I have then computed myself this Jacobian and found the following:

$$DJ(B,X,Y,V)= \begin{pmatrix} -\alpha - \nu V & 0 & 0 & - \nu B \\ -c V X E'(B) & -c V E(B)-\mu_1 & 0 & -c X E(B) \\ c V X E'(B) & c V E(B) & -\mu_2 & c E(B) X\\ -\gamma V & 0 & k & -(\gamma B+E_0+\delta) \\ \end{pmatrix}$$

with

$$E'(B)=e^{-\beta \gamma B}(\gamma-\beta \gamma^2 B-E_0 \beta \gamma)$$

at this point you have to know a bit more about the dynamics of the system. There is two stable points , one where the virus dies , the antibody disappears , the cells stays healthy

$$P_0=(0, \frac{A}{\mu_1},0,0)$$

And another one, depending of the pararameters in a more complex way :

$$P_1 = \left(0,\mu_2 \frac{E_0+\delta}{k c E_0},\bar{Y}, \frac{k}{E_0+\delta} \bar{Y }\right) \text{ with } \bar{Y}= \frac{A}{\mu_2} - \frac{\mu_1(E_0+\delta)}{k c E_0}$$

which means, basically, that some cells become infected and stay infected, the virus population explodes before stabilising, and the antibodies disappear.

This second point is the problem because lots of initial conditions seems to be attracted by it and it means that $$V$$ becomes huge (around $$10^{12}$$) when the total number of cells doesn't exceed $$10^{7}$$ and the population of antibodies stays below $$10^{4}$$.

When $$V$$ becomes huge, I check the Jacobian again, knowing that $$\gamma \approx 1$$. We have this one value in the left down corner becoming really huge compared to the other values, also the parameter $$c$$ is around $$10^{-10}$$, so really it's mostly about this left corner term.

From there I'm not sure what to do. Is there a way to make a time-varying, well-conditioned Jacobian matrix? Or do you know about any other numerical method that could fit my needs?

It doesn't have to be implemented in Matlab. For now, the results I get are mainly garbage.

• Why not compute the Jacobian symbolically? Does the ODE solver allow you to provide the Jacobian? – Rodrigo de Azevedo Apr 11 at 20:01
• hello, yes the solver allows me to do that and that's what I do but I think it replace the symbols with numeric values at computation time before inversion of the Jacobian or whatever it does with the Jacobian. – nathan raynal Apr 22 at 16:40
• Uhmm... are you sure there even is a solution? Nonlinear ODEs can have finite escape time. – Rodrigo de Azevedo Apr 22 at 16:54
• if all initial parameters are above zero , yes it always have a solution . Anyway, I translated my problem in C++ to allow myself to get a better control over the numerical methods used . If I implement myself the implicit method, I should be able to use substitution to reduce the non linear system to a big non linear equation in dimension 1 and could then use an algorithm not involving a Jacobian (a simple newton method could do the trick) . Do you think that's a valid strategy ? I can get this 1d equation with the trapezoidal rule already but didn't try to implement it . – nathan raynal Apr 24 at 9:48
• I haven't touched numerical ODEs in over a decade, so I can't offer a lot. However, I am surprised MATLAB does not get the job done. Their solvers should be sophisticated enough to handle stiff problems. Did you check the documentation thoroughly? – Rodrigo de Azevedo Apr 24 at 21:50