# 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

## 1 Answer

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.

More info on vectorization can be found here: https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

• 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   [] 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
• Hi.. f = zeros(size(x), 'like', x); helps! – Natasha Feb 19 at 15:45
• Oh yes, looks like I misread Matlab's documentation ! So does this solve your problem ? – Laurent90 Feb 19 at 18:35