这是我之前在此处发布的问题的后续行动
我正在尝试构建雅可比矩阵的稀疏模式以加速大型颂歌系统的计算。以下是我试图
odeset在 MATLAB 中为玩具模型设置 jpattern 的代码。
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
当我运行代码时，我注意到向量的大小在第二次迭代x中fun发生了变化，因此发生了错误。
size of x: 10 1
size of x: 10 10
错误：
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);
有关如何修复此错误的建议将非常有帮助。
编辑：
我试图向量化f并fun设置vectvars=2 来向量化雅可比计算。
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
但是，当vectvarsjoptions 中的 = 2 和/或为 ode15s 定义的选项中的 'Vectorized', 'on' 时，再次出现问题。
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);