What's going wrong in the vectorization of my ode function?

조회 수: 3(최근 30일)
Hi All,
I'm trying to setting `vectorization` = 'on' in the ode settings and I clearly don't understand what's going wrong
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', [], '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 = zeros('like', x)
% size(f)
f(1,:) = 0;
f(2:9,:) = mat1*x + mat2*x;
f(10,:) = 2*(x(end-1) - x(end));
% df = [f(1, :); f(2:9, :); f(10, :)];
end
I've tried to vectorize f in fun and also set vectvars=2 to vectorize the jacobian calculation. 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);
Suggestions on how to fix this error will be of great help. Could someone please have a look ?

채택된 답변

Walter Roberson
Walter Roberson 2021년 2월 19일
Put back the
% f = zeros('like', x);
  댓글 수: 1
Deepa Maheshvare
Deepa Maheshvare 2021년 2월 19일
f = zeros(size(x), 'like', x) ;
This helped! thank you . f = zeros('like', x); didn't work for me.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by