# solving differential equations with jacobian pattern

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.

• Your problem is very small, and your jacobian is constant as your problem is linear, so it's only evaluated once. Try a nonlinear problem instead, with more variables, and will see a benefit from this technique if your Jacobian structure is "sparse enough". – Laurent90 Feb 17 at 18:15

## 1 Answer

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