Solving differential equation by setting vectorization `on` in MATLAB

Question

This is a follow up to my previous question posted here.

I've set up an ode system in MATLAB and I'm trying to vectorize the code to increase the speed of computation.

The follow is the code for my toy model.

global mat1 mat2 k
k = 2;
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]';
y0 = [1 1 1 1 1 1 1 1]';
z0 = [x0; y0];
tspan = 0:0.01:0.1;

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

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

function dz = fun(t,z)
fprintf('size of z %d %d...\n', size(z))
dz = zeros(size(z), 'like', z);
f(:,:) = fun_f(z(1:10), z(11:end));
g(:,:) = fun_g(z(11:end));
dz(:,:) = [f;g];
size(dz)
end

function f = fun_f(x, y)
global mat1 mat2 k
f = zeros(size(x), 'like', x);
fprintf('size of f %d %d...\n', size(f))
% size(x)
% size(y)
f(1,:) = 0;
f(2:9,:) = mat1*x + mat2*x -k*y;
f(10,:) = 2*(x(end-1) - x(end));
end

function g = fun_g(y)
global k
g = zeros(size(y), 'like', y);
fprintf('size of g %d %d...\n', size(g))
g(:,:) = k*y;
end

The size of z in (18,1) for few iterations and later it changes to (18,18) becasue of which the following error occurs.

size of z 18 1...
size of f 10 1...
size of g 8 1...

ans =

    18     1

size of z 18 18...
size of f 1 10...
Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise
multiplication, use '.*'.

Error in cse_2_20_21>fun_f (line 60)
f(2:9,:) = mat1*x + mat2*x -k*y;

I'm not sure what's going wrong here. Could someone please have a look? The code works fine when when 'Vectorized', 'off'.

Answer 1

You have not completely vectorized your code. z(1:10) is the first 10 numbers in the data of z, not the first 10 rows. For that you need to write z(1:10,:) like you have already correctly done on the left side.

Why the flip in the matrix orientation occurs between the two calls, that is why z(1:10) inherits the column-ness in the first case and the row-ness in the second, has certainly a good reason in some application case that is now lost to history.