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