MATLAB使用ode45后如何批量计算3x3三维矩阵的误差

方案1：R2020b及以上版本最优解 #

MATLAB R2020b推出的pagemtimes函数原生支持三维矩阵的逐页矩阵乘法，是当前效率最高的实现方式，配合pagetranspose实现逐页转置，代码如下：

% 逐页计算 C_i * C_i^T ，再批量减去单位阵
err = pagemtimes(C, pagetranspose(C)) - repmat(eye(3), 1, 1, size(C,3));

单步即可完成所有切片的误差计算，运算效率远高于循环和其他向量化方案。

方案2：低版本兼容实现 #

如果使用的MATLAB版本低于R2020b，没有内置的页乘法函数，可以选择以下两种实现：

循环实现 #

MATLAB新版本对循环的JIT优化已经非常成熟，3945次循环几乎无延迟，代码可读性最高、不易出错：

slice_num = size(C, 3);
err = zeros(3, 3, slice_num);
for k = 1:slice_num
    err(:,:,k) = C(:,:,k) * C(:,:,k).' - eye(3);
end

arrayfun向量化实现 #

适合偏好单行代码的场景，可读性稍弱：

err_cell = arrayfun(@(k) C(:,:,k)*C(:,:,k).' - eye(3), 1:size(C,3), 'UniformOutput', false);
err = cat(3, err_cell{:});

注意事项 #

你之前直接写C(:,:,:)*C(:,:,:)无法得到预期结果，是因为MATLAB的*运算符默认执行标准二维矩阵乘法，三维矩阵直接调用*会将前两个维度之外的维度当做批次做广播乘法，不是逐页独立计算的逻辑，因此不符合需求。

内容的提问来源于stack exchange，提问作者biscuit2120