MATLAB与Eigen实现的C++矩阵运算输出结果不一致问题

问题描述 #

我正在对比MATLAB代码与C++实现的运行结果，一段逻辑较为简单的矩阵计算代码出现了输出错误，代码功能为生成Gauss-Lobatto切比雪夫点并构造切比雪夫矩阵。

基准MATLAB代码 #

Ny = 10; 
ygl = -cos(pi*(0:Ny)/Ny)'; % 生成Gauss-Lobatto切比雪夫点
% 构造切比雪夫矩阵
VGL = cos(acos(ygl(:))*(0:Ny)); 
dVGL = diag(1./sqrt(1-ygl.^2))*sin(acos(ygl)*(0:Ny))*diag(0:Ny);

初始C++ Eigen实现 #

static const int ny = 10; 
std::vector<double> ygl(ny+1); 

for (int i = 0; i< ny+1; i++){
    ygl[i] = -1. * cos(((i) * EIGEN_PI )/ny); // 该部分计算正确
}
Eigen::Matrix< double, ny+1, 1> v ;
v.setZero();

Eigen::Matrix< double, ny+1, 1> v1 ;
v1.setZero();

Eigen::Matrix< double, ny+1, ny+1> m; 
m.setZero();

Eigen::Matrix< double, ny+1, ny+1> m1;
m1.setZero();

Eigen::Matrix< double, ny+1, ny+1> dv; 
dv.setZero();

for (int j = 0; j < ny+1; j++){
        v[j] = 1./(sqrt(1-pow(ygl[j],2)));
        v1[j] = 1. * (j);
}
m = v.array().sqrt().matrix().asDiagonal(); // 该部分逻辑正确
m1 = v1.array().matrix().asDiagonal(); // 该部分逻辑正确

for (int i = 0; i < ny+1; i++){
    for (int j = 0; j < ny+1; j++){
        dv(j + ny*i) = m(j) * sin(acos(ygl[j]) * (i)) * m1(j); // 该部分计算错误
        cout << dv(j + ny*i) << "\n";
    }
}

初始版本排查进展 #

• 已确认两个列向量转换为对角矩阵m、m1的逻辑正确
• 对角矩阵与正弦项相乘得到目标矩阵的环节输出完全偏离预期，暂未定位到是运算顺序错误还是其他逻辑问题

正确输出参考 #

NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN
         0    1.0000   -3.8042    7.8541  -12.3107   16.1803  -18.4661   18.3262  -15.2169    9.0000    0.0000
         0    1.0000   -3.2361    4.8541   -4.0000   -0.0000    6.0000  -11.3262   12.9443   -9.0000   -0.0000
         0    1.0000   -2.3511    1.1459    2.9062   -6.1803    4.3593    2.6738   -9.4046    9.0000    0.0000
         0    1.0000   -1.2361   -1.8541    4.0000    0.0000   -6.0000    4.3262    4.9443   -9.0000   -0.0000
         0    1.0000   -0.0000   -3.0000    0.0000    5.0000   -0.0000   -7.0000    0.0000    9.0000   -0.0000
         0    1.0000    1.2361   -1.8541   -4.0000   -0.0000    6.0000    4.3262   -4.9443   -9.0000   -0.0000
         0    1.0000    2.3511    1.1459   -2.9062   -6.1803   -4.3593    2.6738    9.4046    9.0000   -0.0000
         0    1.0000    3.2361    4.8541    4.0000   -0.0000   -6.0000  -11.3262  -12.9443   -9.0000    0.0000
         0    1.0000    3.8042    7.8541   12.3107   16.1803   18.4661   18.3262   15.2169    9.0000   -0.0000
       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN

第一次调整后的代码及问题 #

采纳建议后做了如下调整：

• 将v1改为行向量
• 使用std::iota完成0到ny的序列赋值
• 单独构造正弦项中间矩阵dv后，再通过矩阵乘法计算最终结果dv1

调整后代码如下：

Eigen::Matrix< double, 1, ny+1> v1;// 重写原循环内v1[j] = 1. * (j)的赋值逻辑
v1.setZero();
std::iota(v1.begin(), v1.end(), 0);

for (int j = 0; j < ny+1; j++){
        v[j] = 1./(sqrt(1-pow(ygl[j],2)));
}
m = v.array().matrix().asDiagonal(); 
m1 = v1.array().matrix().asDiagonal(); 

for (int i = 0; i < ny+1; i++){
    for (int j = 0; j < ny+1; j++){
        ygl[j] = -1. * cos(((j) * EIGEN_PI )/ny);// 该部分计算正确
        dv(j + ny*i) =   sin(acos(ygl[j]) * (i));// 该部分计算正确
    }
}
dv1 = (m) * (dv) * (m1);
std::cout << dv1 << "\n"; // 输出结果不符合预期

调整后dv1的输出结果依然与参考结果不符。

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