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