mercredi 15 juin 2022

Incorrect results/output of MATLAB vs C++ using Eigen matrices

I am comparing my MATLAB code to c++ and getting incorrect output for something fairly simple.

My MATLAB code:

Ny = 10; 
ygl = -cos(pi*(0:Ny)/Ny)'; %Gauss-Lobatto chebyshev points 
%Chebyshev matrix: 

VGL = cos(acos(ygl(:))*(0:Ny)); 
dVGL = diag(1./sqrt(1-ygl.^2))*sin(acos(ygl)*(0:Ny))*diag(0:Ny);

C++ code:

    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); //This is correct
    
    }
    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(); //this is correct
    m1 = v1.array().matrix().asDiagonal(); //this is correct

    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); //this is incorrect
            cout << dv(j + ny*i) << "\n";
        }
    }

As you can see finding the diagonal matrix of the two column vectors v and v1 are correct but then when multiplied with an additional term give completely off results. I am not sure if it's as simple as order of operation here or something else.

To give more context for this specific example, the correct output should look like:

 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


I am not sure where is the issue here. Thanks

Aucun commentaire:

Enregistrer un commentaire