计算科学 - 通过在 MATLAB 中设置矢量化“on”来求解微分方程 - 吾爱随笔录

我已经在 MATLAB 中建立了一个 ode 系统，我正在尝试对代码进行矢量化以提高计算速度。

以下是我的玩具模型的代码。

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

(18,1) 中 z 的大小经过几次迭代后变为 (18,18)，因为出现以下错误。

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;

我不确定这里出了什么问题。有人可以看看吗？当“矢量化”、“关闭”时，代码工作正常。