我想由已知的点应力矩阵，计算出三个主应力和主方向，理论很成熟就是代公式即可，但是结果上出现了问题，怀疑为计算精度问题，具体如下：
计算主应力：
for i=1:1:r      
T{i}=[xx(i),xy(i),zx(i);xy(i),yy(i),yz(i);zx(i),yz(i),zz(i)];
[~,principalstress{i}]=eig(T{i});
p1(i)=principalstress{i}(1,1);
p2(i)=principalstress{i}(2,2);
p3(i)=principalstress{i}(3,3);
end
计算主方向：（以主应力1为例）
l1=1;
m1=l1*(zx*xy-yz*(xx-p1))/(xy*yz-zx*(yy-p1));
n1=l1*((xx-p1)*(yy-p1)-xy*xy)/(xy*yz-zx*(yy-p1));
abcd1=sqrt(l1^2+m1^2+n1^2);%归一化
l1=l1/abcd1;
m1=m1/abcd1;
n1=n1/abcd1;
这套方法计算文献算例也好，自己随便输入的点应力矩阵（类似[1，2，3；2，4，5；3，5，6]这种）也好，都没问题；但是计算我的实验数据时得到的主应力还算合理，得到的主方向问题很大，即三个方向矢量点乘的结果完全不能认为是0，这不符合物理规律。
怀疑是程序计算精度不足导致的，因为实验测得的点应力，六个分量之间数量级从e-4到e-12不等，差距很大，公式又有乘除运算，可能导致浮点计算误差很大。
所以想请教各位大佬，有什么办法可以提高运算精度，或者修改算式回避这个问题，再或者有没有擅长应力理论的大佬，发现是其他问题导致的结果错误。

我想由已知的点应力矩阵，计算出三个主应力和主方向，理论很成熟就是代公式即可，但是结果上出现了问题，怀疑为计算精度问题，具体如下：
计算主应力：
for i=1:1:r      
T{i}=[xx(i),xy(i),zx(i);xy(i),yy(i),yz(i);zx(i),yz(i),zz(i)];
[~,principalstress{i}]=eig(T{i});
p1(i)=principalstress{i}(1,1);
p2(i)=principalstress{i}(2,2);
p3(i)=principalstress{i}(3,3);
end
计算主方向：（以主应力1为例）
l1=1;
m1=l1*(zx*xy-yz*(xx-p1))/(xy*yz-zx*(yy-p1));
n1=l1*((xx-p1)*(yy-p1)-xy*xy)/(xy*yz-zx*(yy-p1));
abcd1=sqrt(l1^2+m1^2+n1^2);%归一化
l1=l1/abcd1;
m1=m1/abcd1;
n1=n1/abcd1;
这套方法计算文献算例也好，自己随便输入的点应力矩阵（类似[1，2，3；2，4，5；3，5，6]这种）也好，都没问题；但是计算我的实验数据时得到的主应力还算合理，得到的主方向问题很大，即三个方向矢量点乘的结果完全不能认为是0，这不符合物理规律。
怀疑是程序计算精度不足导致的，因为实验测得的点应力，六个分量之间数量级从e-4到e-12不等，差距很大，公式又有乘除运算，可能导致浮点计算误差很大。
所以想请教各位大佬，有什么办法可以提高运算精度，或者修改算式回避这个问题，再或者有没有擅长应力理论的大佬，发现是其他问题导致的结果错误。