Two Matlab ODE solvers, two different results

Question

I am solving a system of ODEs using Matlab. One particular set of parameters caused the solver to fail, so I worked my way through the different solvers Matlab provides. I was surprised to find that two of the solvers (ode23s and ode23tb) produced completely different, yet reasonable, results. Every other solver failed.

My question: Which one (if any) should I trust and why?

I'm not entirely sure if this is a mathematics question or a programming question, but I suspect it comes to a difference in the numerical methods, hence posting it here.

ode23s:

ode23tb:

For reference the ODEs are: $$\frac{dS}{dt}=-\beta IS+\xi(1-S)+3 \kappa L\\ \frac{dI}{dt}=\beta IS+\alpha IL-\xi I-\gamma I\\ \frac{dR}{dt}=\nu \beta IW+\gamma I - \xi R - 3\kappa R \\ \frac{dL}{dt}=-\alpha IL+3 \kappa W-\xi L -3 \kappa L\\ \frac{dW}{dt}=-\nu \beta IW + 3 \kappa R -\xi W - 3 \kappa W\\$$

With parameters: $$\alpha=26, \beta=260, \kappa=.1, \gamma=17, \nu=5, \xi=0.0125$$

And initial conditions:

$$S_0=.99, I_0=1-S_0, R_0=0, L_0=0, W_0=0$$

Answer 1

Looks a lot like you have (at least one) bifurcation point.

You can perform a stability analysis based on linerization which shows that actually both equilibria have exactly one unstable direction corresponding to an eigenvalue with negative real part.

Thus, if the trajectories from the different solvers are somewhat sufficiently aligned with those unstable directions, this causes them to miss that equilibrium / steady state. Conversely, if they progress in directions more or less aligned with the stable ones, they approach the corresponding equilibrium. In the context of your system, ode23s "finds" thus the stable directions of the $(1, 0, 0, 0, 0)$ equilibrium or is repelled by the unstable direction of the $\sim (0.0589, 0.0021, 0.7937, 0.0655, 0.0799)$ steady state (or mixture of both). Converse holds for the ode23tb solver.

I attach some Matlab files I used throughout the analysis, including a animation of the trajectories in reduced ($I, L, R + L +W$) space, showing that somewhere around $I \approx L \approx 0, R + L + W \approx 1$ the bifurcation happends.

    function f = ODE_System(~, SIRLW)

    % Parameters
    alpha = 26;
    beta  = 260;
    kappa = 0.1;
    k     = 3 * kappa; % Kappa only used after multiplication with 3
    gamma = 17;
    nu    = 5;
    xi    = 0.0125;
    
    S = SIRLW(1); I = SIRLW(2); R = SIRLW(3); L = SIRLW(4); W = SIRLW(5);
    
    f = zeros(5, 1);
    f(1) = -beta * I * S + xi * (1 - S) + k * L;
    f(2) =  beta * I * S + alpha * I * L - xi * I - gamma * I;
    f(3) = nu * beta * I * W + gamma * I - xi * R - k * R;
    f(4) = -alpha * I * L + k * W - xi * L - k * L;
    f(5) = - nu * beta * I * W + k * R - xi * W - k * W;
    end

    function Df_Dx = Eval_Jacob(SIRLW_eq)
    % Evaluates Jacobian at given point
    
    % Parameters
    alpha = 26;
    beta  = 260;
    kappa = 0.1;
    k     = 3 * kappa; % Kappa only used after multiplication with 3
    gamma = 17;
    nu    = 5;
    xi    = 0.0125;
    
    S_eq = SIRLW_eq(1); I_eq = SIRLW_eq(2); R_eq = SIRLW_eq(3); L_eq = SIRLW_eq(4); W_eq = SIRLW_eq(5);
    
    Df_Dx = [-beta * I_eq - xi, -beta * S_eq,                            0,        k,                      0;
             beta * I_eq,       beta * S_eq + alpha * L_eq - xi - gamma, 0,        alpha * I_eq,           0; 
             0,                 nu * beta * W_eq + gamma,                -xi - k,  0,                      nu * beta * I_eq;
             0,                 -alpha * L_eq,                           0,        -alpha * I_eq - xi - k, 0;
             0,                 -nu * beta * W_eq,                       k,        0,                      -xi - k;         ];
    end

    clear;
    close all;
    clc;
    
    % Parameters
    alpha = 26;
    beta  = 260;
    kappa = 0.1;
    k     = 3 * kappa; % Kappa only used after multiplication with 3
    gamma = 17;
    nu    = 5;
    xi    = 0.0125;
    
    %% Equilibrium of upper / first plot
    SILRW_Eq_1 = [1; 0; 0; 0; 0];
    Df_Dx      = Eval_Jacob(SILRW_Eq_1);  
    [EVec_Eq_1, EVal_Eq_1] = eigs(Df_Dx); % Shows one eigenvalue (whose real part) > 0 => Unstable equilibrium
    
    %% Equilibrium of lower / second plot
    tspan   = [0, 100];
    S_0     = 0.99;
    SILRW_0 = [S_0; 1 - S_0; 0; 0; 0];
    [t_tb, SILRW_tb] = ode23tb(@ODE_System, tspan, SILRW_0);
    [t_s,  SILRW_s]  = ode23s( @ODE_System, tspan, SILRW_0);
     
    % Show that this values lead indeed to (more or less) an equilibrium:
    SILRW_Eq_2 = SILRW_tb(end, :)';
    t_dummy = 42;
    disp(ODE_System(t_dummy, SILRW_Eq_2)); % Not really 0 yet
    
    % Other way of obtaining this equilibrium: Solve eqs. numerically with ode solver outcome as start value
    options = optimoptions('fsolve', 'Display', 'iter');
    fun = @(x) ODE_System(t_dummy, x);
    SILRW_Eq_2 = fsolve(fun, SILRW_Eq_2, options);
    % Check values of RHS again, they are now machine precision:
    disp(ODE_System(t_dummy, SILRW_Eq_2));
    
    Df_Dx     = Eval_Jacob(SILRW_Eq_2);  
    [EVec_Eq_2, EVal_Eq_2] = eigs(Df_Dx); % Shows one eigenvalue with real part > 0 => Unstable equilibrium
    
    %% Investigate trajectories
    % Bundle RLW into one variable (leaves one with 3 variables which can be plotted)
    SI_Bundled_tb = [SILRW_tb(:, 1), SILRW_tb(:, 2), SILRW_tb(:, 3) + SILRW_tb(:, 4) + SILRW_tb(:, 5)];
    SI_Bundled_s  = [SILRW_s(:, 1),  SILRW_s(:, 2),  SILRW_s(:, 3)  + SILRW_s(:, 4)  + SILRW_s(:, 5)];
    
    for i = 1:length(t_s)
      plot3(SI_Bundled_s(i, 1), SI_Bundled_s(i, 2), SI_Bundled_s(i, 3), '-rx');
      xlabel('S');
      ylabel('I');
      zlabel('R + L + W');
      title('Trajectories. Red: ode23 S, Blue: ode23 TB');
      hold all;
      pause(0.1);
    end
    
    for i = 1:length(t_tb) / 3 % All the action done by then
      plot3(SI_Bundled_tb(i, 1), SI_Bundled_tb(i, 2), SI_Bundled_tb(i, 3), '-bo', 'Displayname', 'ode23_TB');
      hold all;
      pause(0.1);
    end