issue with vectorization of ode

1 view (last 30 days)
Deepa Maheshvare
Deepa Maheshvare on 21 Feb 2021
Commented: Deepa Maheshvare on 21 Feb 2021
Hi All,
This is a follow up to my previous question here.
I'm trying run the code below.
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', 'on'); %, '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' is set 'off'.
  2 Comments
Walter Roberson
Walter Roberson on 21 Feb 2021
f(:,:) = fun_f(z(1:10), z(11:end));
g(:,:) = fun_g(z(11:end));
Should be
f = fun_f(z(1:10,:), z(11:end,:));
g = fun_g(z(11:end,:));

Sign in to comment.

Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!