Создание набора ортогональных осей

Question

Мне нужно создать серию векторов для моделирования, которые дают направление распространения и поляризацию света (3 измерения). Из-за того, что направление и поляризация должны быть ортогональны, так что мне нужно что-то вроде этого:

D= [dir1;  P=[pol11; 
    dir1;   pol12; 
    dir2;   pol21; 
    dir2;   pol21; 
    dir3;   pol31; 
    .... ]  ... ] 

Как вы можете видеть каждое направление имеет две поляризации. Здесь важно, что dir1 * pol11 '= 0, dir1 * pol12' = 0 и т. Д. Направления должны охватывать весь телесный угол, в то время как направление поляризации не является строго важным, хотя было бы неплохо, если бы они ортогональны друг другу. Я пробовал два разных подхода, один создает базовый ортогональный базис и вращает его, другой создает матрицу направлений и использует функцию null() для создания поляризации. В обоих случаях я получаю, что если я делаю D * P ', я получаю серию из 0, но некоторые из значений не равны нулю - очень маленькие (например, 1е-17), но тем не менее не равны нулю. Код 1:

bDir=[1,0,0;1,0,0] 
bPol=[0,1,0;0,0,1] 
dir=bDir 
pol=bPol 
for phi=0:pi/5:2*pi 
    for theta=0:pi/5:pi 
     rotatePhi=[cos(phi) -sin(phi) 0;... 
         sin(phi) cos(phi) 0;... 
         0 0 1]; 
     rotateTheta=[cos(theta) 0 sin(theta);... 
           0 1 0;... 
           -sin(theta) 0 cos(theta)]; 
     rDir=bDir*rotateTheta*rotatePhi; 
     rPol=bPol*rotateTheta*rotatePhi; 

     dir=vertcat(dir,rDir); 
     pol=vertcat(pol,rPol); 

    end 
end 

Код 2:

bDir=[1,0,0;1,0,0] 
dir=bDir 
pol=[0,1,0;0,0,1] 
for phi=0:pi/5:2*pi 
    for theta=0:pi/5:pi 
     rotatePhi=[cos(phi) -sin(phi) 0;... 
         sin(phi) cos(phi) 0;... 
         0 0 1]; 
     rotateTheta=[cos(theta) 0 sin(theta);... 
           0 1 0;... 
           -sin(theta) 0 cos(theta)]; 
     rDir=bDir*rotateTheta*rotatePhi; 
     rPol=null(rDir)'; 

     dir=vertcat(dir,rDir); 
     pol=vertcat(pol,rPol); 

    end 
end 

Я считаю, что проблема в том, что MATLAB вводит некоторые точные ошибки, но это может быть не так. Может ли кто-нибудь сказать мне, есть ли ошибка в моем коде, или если есть лучший способ получить 2 матрицы, которые я ищу?

Answer 1

Я думаю, что ваши ожидания относительно floating point arithmetic немного ушли. Арифметика с плавающей точкой не является точной. В MATLAB обычно арифметические операции с плавающей запятой выполняются только с точностью до 16 знаков после запятой (1e-16). Любые значения, меньшие этого, обычно можно рассматривать как 0.

Функция eps в MATLAB может использоваться для определения наименьшего числа с плавающей запятой, которое MATLAB может представлять. В моей системе это производит

>> eps 
ans = 

    2.2204e-16 

Любые числа, меньшие, чем этот MATLAB, не могут точно представлять.

---

А также точность арифметики с плавающей точкой это также, как правило, не так, чтобы сравнить их с целочисленным значением 0. Более точное сравнение было бы

if dir*pol < 1e-16 
    fprintf(1, 'Orthogonal\n'); 
else 
    fprintf(1, 'Not Orthogonal\n'); 
end 

В большинстве случаев при использовании арифметики с плавающей точкой, и, как правило, в MATLAB дали выше вывод eps любое значение меньше, чем 2.2204e-16 эффективно 0.

---

Если вы хотите переменную точность арифметических операций можно использовать vpa. Тем не менее, это, как правило, не нужно, и арифметика с плавающей запятой должна быть достаточной, пока вы знаете о ее падениях.