Fortran丸め誤差

Question

私はシンプルなコードを持っています。これは、円柱で囲まれた領域にノードをフラグします。コードを実装すると、結果は90度のシリンダ観測ケースの軽い傾きになります。

実際の問題： 上記のアルゴリズムはFortranで実装されています。コードは、円柱の内部にある場合、デカルトグリッドの点をチェックします。テストケースに続いて： シリンダは、y軸に対してyz平面内で90度の角度をなします。したがって、方位ベクトル$ \ vec {o} $は（0,1,0）です。

ケース1： オリエンテーションベクタは、$ \ vec {o} =（0.0,1.0,0.0）$で直接割り当てられます。 0.0（= 配向ベクトルが$ \ VEC {O}と倍精度精度dsinとdcosと固有Fortranの機能を指定します。これは、$の\シータ= 90 $

ケース2と完全シリンダをもたらします、\ sin（\ pi/2.0）、\ cos（\ pi/2.0））$ $ pi $の値には、20以上の重要な小数点が割り当てられています。結果として得られる円柱は軽度の傾斜を生じる。

強調表示された領域は、デカルト軸に対する円柱の傾きによる余分な材料を示します。私も試したarchitecture specific maximum precision "pi" value.これも同じ問題が発生します。

これは、シリンダによる実際の角度が90度ではないことを示しています。誰もこの問題の有効な解決策を提案することができます。私は任意の角度のためのinbuilt三角関数を使用し、正確なセルのフラグを立てる方法を探している必要があります。

注：すべての操作は倍精度の精度で行われます。

実際の機能は以下のとおりです。 rk値と8

pure logical function in_particle(p,px,x) 
    type(md_particle_type),intent(in) :: p 
    real(kind=rk),intent(in) :: px(3),x(3) 
    real(kind=rk) :: r(3),rho(3),rop(2),ro2,rdiff,u 
    rop = particle_radii(p) ! (/R_orth,R_para/) 
    ro2 = rop(1)**2 
    rdiff = rop(2) - rop(1) 

    r = x-px 

! Case 1: 
! u = dot_product((/0.0_rk,-1.0_rk,0.0_rk/),r) 
! rho = r-u*(/0.0_rk,-1.0_rk,0.0_rk/) 

! Case 2: 
    u = dot_product((/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/),r) 
    rho = r-u*(/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/) 

    if((u.le.rdiff).and.(u.ge.-rdiff)) then 
    in_particle = dot_product(rho,rho) < ro2 
    else 
    in_particle = .false. 
    end if 
end function in_particle 

注パラメータが定義される：三角の操作は、より良い問題を説明するために、コード内で行われます。しかし、元のコードは、ユーザからベクタ形式で方向を読み取ります。次に、この情報を粒子 - 粒子衝突操作のための四元数に変換する。クォータニオンを方向ベクトルに戻すと、この誤差はさらに大きくなります。衝突の開始前であっても、気筒の向きは2つの格子セルによって方向を変えられる傾向がある。

Answer 1

cos(pi/2)は必ずしも関係なく、あなたがcos計算を行う方法を正確に、あなたは正確に0を提供するつもりはなく、あなたが持っているどのように多くのpiの桁数に関係なく、理由はされていません。

• pi、無理数としてFP番号として表されるとき、最大1/2のulpの誤差を含みます。
• sinおよびcosは、IEEE-754標準によって正しく丸められ（または実装されていても）保証されていません。今

、sin(pi/2)はsin 1の周りにこのような低誘導体を有するというだけの理由にかかわらず、精度およびFPアーキテクチャの1として出てくる非常に可能性があります。単精度浮動小数点数では、正確な値がpi/2の約3e-4の範囲内にあれば、1になるはずです。問題のある呼び出しはcosです。これは、0付近で再生する精度が非常に高く、近隣では約-1の微係数です。

まだ、ここでは非常に小さな値について話しています。ここで問題を本当に増強するのは、普通のFPの丸めルールと組み合わせて、あなたがやっているイン/アウトテストです。グリッド量子の四分の一でテストポイントをバイアスすると、ボクセル化のすべての直立垂直線を見ることができます（ただし、これは副軸を中心に対称ではありません）。

ドットプロダクトを実行する前にsin/cos計算から実際にいくつかの精度を捨てて、実際に軸を量子化することもできます。

Answer 2

短い答え：一般的な角度（0、π/ 6、π/ 4、π/ 3、π/ 2、πおよびそれらの倍数）のsinとのテーブルを作成し、珍しい角度だけを計算します。一般的な角度のエラーは許容されない可能性があるのに対し、一般的でない角度のエラーはほとんどの人が許容するためです。

説明： 浮動小数点の計算が正確ではない（その性質上）ため、コードの精度と可読性の間に多少の妥協が必要になることがあります。

これを行う1つの方法は、正確にわかっているものを計算することを避けることです。これを行うには、角度の値を確認し、明らかな角度でない場合にのみ実際の計算を行います。例えば、角度0,90,180および270度は、sinおよびcosの明白な値を有する。より一般的には、一般的な角度（0、π/ 6、π/ 4、π/ 3、π/ 2、πおよびそれらの倍数）のcosおよびsinは正確に（それらが非合理的な数であっても）既知である。