Численная ошибка в вычислениях направления вектора - PullRequest
Новая
       71

Численная ошибка в вычислениях направления вектора

0 голосов
/ 12 апреля 2020

Я пишу код, который будет вводить вектор ri и вычислять углы, которые он образует с осями x, y и z. Таким образом я узнаю направление вектора и смогу рассчитать составляющие силы (которые я вычисляю отдельно) в этом направлении. alfa - это угол между вектором и осью z, beta - это угол между вектором и осью x, gamma - это угол между вектором и осью y.

Вот код, который делает это: 

long double length_of_vector, alfa, betta, gamma, x_force, y_force, z_force;

vector<double> ri = {0, -3.85443e-07, 5.98761e-06} ;

length_of_vector = sqrt(pow(ri[0],2)+pow(ri[1],2)+pow(ri[2],2)) ;

    alfa = acos(ri[2]/length_of_vector);  //angle with z axis

    betta = acos ( ri[0] / (length_of_vector * sin(alfa)) ) ; //angle with x axis

    gamma = acos ( ri[1] / (length_of_vector * sin(alfa)) ) ; //angle with y axis

    cout<< "alfa is: " << alfa * 180/M_PI << endl ;
    cout<< "betta is: " << betta * 180/M_PI<< endl ;
    cout<< "gamma is: " << gamma* 180/M_PI << endl ;
    cout<< "length is: " << length_of_vector << endl ;
    //------------------------------------------------END OF ANGLE CALCULATIONS

    if ((ri[1] == 0)&&(ri[2]==0)) {
    x_forces += length_of_vector * sin(alfa) * cos(betta) ;
    }
    else if ((ri[0] == 0)&&(ri[2]==0)) {
        y_forces += length_of_vector * sin(alfa) * cos(gamma) ;
    }
    else if ((ri[0] == 0)&&(ri[1]==0)) {
        z_forces += length_of_vector * cos(alfa) ;
    }
    else if (ri[0] == 0) {
        y_forces += length_of_vector * sin(alfa) * cos(gamma) ;
        z_forces += length_of_vector * cos(alfa) ;
    }
    else if (ri[1] == 0) {
        x_forces += length_of_vector * sin(alfa) * cos(betta) ;
        z_forces += length_of_vector * cos(alfa) ;
    }
    else if (ri[2] == 0) {
        x_forces += length_of_vector * sin(alfa) * cos(betta) ;
        y_forces += length_of_vector * sin(alfa) * cos(gamma) ;
    }
    else {
        x_forces += length_of_vector * sin(alfa) * cos(betta) ;
        y_forces += length_of_vector * sin(alfa) * cos(gamma) ;
        z_forces += length_of_vector * cos(alfa) ;
    }

Как вы можете видеть здесь, вместо length_of_vector в операторах if я должен вставить величину силы, которая будет суммироваться в x_force, y_force и z_force. Тем не менее, для целей проверки я поставил length_of_vector там. Таким образом, ответ, который я получаю, должен представлять компоненты вектора, который я первоначально ввел.

Я ввожу ri (что является vector<double>) как {0, -3.85443e-07, 5.98761e-06}. Вот что я получаю: 

x_force = 0
y_force = nan
z_force = 5.98761e-06

Как видно, я получаю, что y_force вычисляется как nan. Я проследил код, я понял, что это nan, потому что значение gamma становится nan, пытаясь вычислить этот угол.

Если я изменю значение компонента y ri на -3.85443e-05, я получу правильный ответ, но все, что меньше этого значения, даст мне gamma = nan. Я не понимаю, почему это происходит вообще. Значит ли это, что вектор {0, -3.85443e-07, 5.98761e-06} не существует? Пожалуйста, помогите!

...