solving differential equations with jacobian pattern

Question

I'm trying to compare the simulation time for solving a system of differential equations with and without jacobian pattern for a toy model using ode15s in MATLAB.

global mat1 mat2
mat1=[ 
       1    -2     1     0     0     0     0     0     0     0;
       0     1    -2     1     0     0     0     0     0     0;
       0     0     1    -2     1     0     0     0     0     0;
       0     0     0     1    -2     1     0     0     0     0;
       0     0     0     0     1    -2     1     0     0     0;
       0     0     0     0     0     1    -2     1     0     0;
       0     0     0     0     0     0     1    -2     1     0;
       0     0     0     0     0     0     0     1    -2     1;
       ];

mat2 = [
        1    -1     0     0     0     0     0     0     0     0;
        0     1    -1     0     0     0     0     0     0     0;
        0     0     1    -1     0     0     0     0     0     0;
        0     0     0     1    -1     0     0     0     0     0;
        0     0     0     0     1    -1     0     0     0     0;
        0     0     0     0     0     1    -1     0     0     0;
        0     0     0     0     0     0     1    -1     0     0;
        0     0     0     0     0     0     0     1    -1     0;
        ];



x0 = [1 0 0 0 0 0 0 0 0 0]';
tspan = 0:0.01:5;

y = sym('y', [10 1]);
s = mat1*y + mat2*y;
J = jacobian(s, y);
[jr, jc] = size(J);
jpattern = sparse(jr, jc);
jpattern(find(J~=0)) = 1;
jac = sparse(10,10);
jac(2:9, :) = jpattern;
jac(10,:) = [0 0 0 0 0 0 0 0 1 1];

options = odeset('Stats', 'on', 'JPattern', jac);

ttic = tic();
[t, sol]  =  ode15s(@(t,x) fun(t,x), tspan , x0, options);
ttoc = toc(ttic)
fprintf('runtime %f seconds ...\n', ttoc)
plot(t, sol)

function f = fun(t,x)
global mat1 mat2
f(1,1) = 0;     
f(2:9,1) = mat1*x + mat2*x;
f(10,1) = 2*(x(end-1) - x(end));
end

I expected that the simulation time will decrease when the jacobian pattern is specified [t, sol] = ode15s(@(t,x) fun(t,x), tspan , x0, options) . However, I find that the simulation takes longer when jpattern is specified.

Clarifications on why the simulation time increases when jpattern is given will be helpful.

Answer 1

This problem is too small to actually be sparse. Sparse handling has a big overhead because the indexing is not "direct", i.e. you don't necessarily know where the next value will be without branch checking. So you need it to be "sparse enough" that the O(n^3) dense LU-factorization cost shrinking to the purely non-zero terms overcomes the very large constant cost of using a sparse array. This could be overcome by using structured arrays in type-based languages like Julia or Haskell, but if you want to stick to MATLAB, then just make sure to keep sparse operations for very large problems. A good rule of thumb is <0.01% non-zeros (that might be a little conservative, but that depends on the operation).