# Solving differential equation by specifying jacobian pattern

This is a follow up to my previous question posted here

I'm trying to construct the sparsity pattern of the jacobian matrix to speed up the computation of a large system of odes. The following is the code in which I am trying to set up the jpattern in odeset for a toy model 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;

f0 = fun(0, x0);
joptions = struct('diffvar', 2, 'vectvars', 1, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions);
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'JPattern', sparsity_pattern);

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

fprintf('size of x: %d %d\n', size(x))

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


When I run the code, I notice that the size of the vector x in fun changes for the second iteration because of which an error occurs.

size of x: 10 1
size of x: 10 10


Error:

Unable to perform assignment because the size of the left side is 8-by-1 and the size of the right side is 8-by-10.

Error in Untitled>fun (line 47)
f(2:9,1) = mat1*x + mat2*x;

Error in odenumjac (line 143)
Fdel = feval(F,Fargs_expanded{:});

Error in Untitled (line 31)
J = odenumjac(@fun,{0 x0}, f0, joptions);


Suggestions on how to fix this error will be really helpful.

EDIT:

I've tried to vectorize f in fun and also set vectvars=2 to vectorize the jacobian calculation.

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;

f0 = fun(0, x0);
joptions = struct('diffvar', 2, 'vectvars', 2, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions);
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'JPattern', sparsity_pattern, 'Vectorized', 'on');

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,:) = 0;
f(2:9,:) = mat1*x + mat2*x;
f(10,:) = 2*(x(end-1) - x(end));
end


However, there is a problem again when vectvars = 2 in joptions and/or 'Vectorized', 'on' in options defined for ode15s.

Unable to perform assignment because the size of the left side is 8-by-1 and the size of the right side is 8-by-10.

Error in cse_11_5_20>fun (line 44)
f(2:9,:) = mat1*x + mat2*x;

Error in odenumjac (line 143)
Fdel = feval(F,Fargs_expanded{:});

Error in cse_11_5_20 (line 31)
J = odenumjac(@fun,{0 x0}, f0, joptions);

• This question seems to be more about programming than computational science. – nicoguaro Feb 18 at 15:28
• @nicoguaro Do you suggest I migrate this to stackoverflow? If this is not the right place to post this question, could one of the moderators consider migrating my post to stackoverflow? – Natasha Feb 18 at 15:29
• I am not sure if it is suitable there. – nicoguaro Feb 18 at 15:36
• Regarding your code, mat1*x + mat2*x is a vector with length 10 and you are assigning it to one with length 8, f(2:9, 1). – nicoguaro Feb 18 at 15:37
• @nicoguaro Thank you. mat1 is of size (8, 10) and x is of size (10,1) so the resulting vector is of size 8 x 1. This works fine in the 1st iteration . In the second iteration, x is of size (10, 10). This results in a matrix of size 8 x 10 (here is the issue and this happens while using odenumjac). I'm not sure why the size of x changes. – Natasha Feb 18 at 15:40

odenumjac calls your function in a vectorized manner it seems, and your function is not vectorized. You can easily change that by changing the second index of f in your function to : instead of 1, for instance: f(10,:) = 2*(x(end-1,:) - x(end,:));
I thought the setting joptions.vectvars=1 would not allow the vectorised call (see one of your other questions). I realize that this is actually not the case, you should instead set it to [] (empty). If you want to vectorize the Jacobian calculation, set it to 2. You can type open odenumjacin the Matlab console to access the function file and read its documentation.
• Thanks a lot again. Setting vectvars = [ ] solved my problem and for my real system the computation time reducted from ~2 hours to ~7 minutes . But I am still facing problem in using the vectorized form of both jacobian and the function f in fun . Could you please have a look edit 2? I looked the function file of odenumjac, unfortunately I couldn't understand the documentation provided for vectvars. – Natasha Feb 19 at 12:33
• @Natasha this may come from the first line f(1,:) = 0 which may result in the second dimension of f being initialised as 1. I would try initialising f properly at the start of fun, for example with f = zeros('like',x). Maybe you also need to play with the value of vectvars. For that you need to look at the odenumjac function to see how it works. – Laurent90 Feb 19 at 13:22
• f = zeros('like', x) gives size(f) = [1,1]. I tried all three options [1] [2] [] for vectvars. The simulation fails when I set 'Vectorized', 'on' when vectvars=[]. I couldn't understand the details given for vectvars in odenumjac function :/ – Natasha Feb 19 at 13:43