旧ブログ https://blog.goo.ne.jp/yuuhae612go
https://temperament-calculation-blog2.blogspot.com/2023/05/blog-post_2.html
2023年10月22日日曜日
2023年10月7日土曜日
４点を通る球面(３次元空間内)
エクセル等の表計算ソフトの逆行列関数を使って４点を通る球面を求める。(３次元空間内)
「３点を通る円(２次元平面内)」と全く同じ方法で求める。
まず求める球面の式を
x^2+y^2+z^2+l*x+m*y+n*z+o=0
とすると
l*x+m*y+n*z+o=-(x^2+y^2+z^2)
なので４点を
p1(x1,y1,z1)、p2(x2,y2,z2)、p3(x3,y3,z3)、p4(x4,y4,z4) (表のD4:F7)
とすると
l*x1+m*y1+n*z1+o=-(x1^2+y1^2+z1^2)
l*x2+m*y2+n*z2+o=-(x2^2+y2^2+z2^2)
l*x3+m*y3+n*z3+o=-(x3^2+y3^2+z3^2)
l*x4+m*y4+n*z4+o=-(x4^2+y4^2+z4^2)
となる。行列の式で書けば(libreofficecalcでmathをオブジェクトとして使用)
となり、逆行列を使えば一発で答えが出る。
逆行列と配列定数を使用してl,m,n,oを求める。配列定数を使って前半は４行３列の行列の１～３列目の順を逆にして４列目を追加して４列目を１とした。後半は２乗して横方向の和。
D8:D11 =ArrayFormula(mmult(MINVERSE(mmult(D4:F7,{0,0,1,0;0,1,0,0;1,0,0,0})+{0,0,0,1;0,0,0,1;0,0,0,1;0,0,0,1}),-MMULT(D4:F7^2,{1;1;1})))
次に中心座標と半径を求める。
中心z D12 =-D10/2
中心y E12 =-D9/2
中心x F12 =-D8/2
半径^2 D13 =((D8^2/4+D9^2/4+D10^2/4)-D11)
半径 D14 =((D8^2/4+D9^2/4+D10^2/4)-D11)^(1/2)
検算的グラフとデータは今は省略。
2023.10.07 20:49公開 単純ミスなどは適宜修正済＆修正予定
余談１
以下テスト用
htmlの表を使った行列の式の残骸ｗ
| [x1 | y1 | z1 | 1] | [l] | [x1^2+y1^2+z1^2] | |
| [x2 | y2 | z2 | 1] | [m] | =- | [x2^2+y2^2+z2^2] |
| [x3 | y3 | z3 | 1] | [n] | [x3^2+y3^2+z3^2] | |
| [x4 | y4 | z4 | 1] | [o] | [x4^2+y4^2+z4^2] |
| [l] | [x1 | y1 | z1 | 1]^-1 | [x1^2+y1^2+z1^2] | |
| [m] | =- | [x2 | y2 | z2 | 1] | [x2^2+y2^2+z2^2] |
| [n] | [x3 | y3 | z3 | 1] | [x3^2+y3^2+z3^2] | |
| [o] | [x4 | y4 | z4 | 1] | [x4^2+y4^2+z4^2] |
2023年6月15日木曜日
銀行丸め検証結果
グーグルドライブ、リブレオフィス、エクセルでいろいろな銀行丸め(jis丸め)の式を検証してみた。結果、round系の区間交替型＆重ね合わせ型が最優秀だった。int系の区間交替型＆重ね合わせ型はシンプルで結構好きだが誤差が出てしまった。特異点型は思ったより悪くない。検証してないがround系で組めば少し改善するかも知れない。二重丸め型が最もダメだった。第二引き数を活用した式やintと併用した式は特に悪かった。
元数=$B8 (小数点以下の)桁数=$I$4 とする。
round系区間交替型 銀行丸め4c
CJ8 =if(MOD(ABS($B8)*10^$I$4,2)>=1,round($B8,$I$4),if(abs($B8)<=1/2/(10^$I$4),0,ROUNDUP($B8-sign($B8)*1/2/(10^$I$4),$I$4)))
round系重ね合わせ型 銀行丸め7
CK8 =if(abs($B8)<=1/2/10^$I$4,0,ROUNDUP($B8/2-sign($B8)*1/4/10^$I$4,$I$4)+rounddown($B8/2+sign($B8)/4/10^$I$4,$I$4))
round系重ね合わせ型 銀行丸め7b
CL8 =if(abs($B8)<=1/2/10^$I$4,0,ROUNDUP($B8/2-sign($B8)*1/4/10^$I$4,$I$4)+round($B8/2-sign($B8)/4/10^$I$4,$I$4))
検証方法
中心値xと(小数点以下の)桁数yと分割数zをセット(１件と数えた)にして、誤った値を目視でカウントしていく方法。誤差を探すための別列には「平均を使った統計的な数値」が良いと思ったが、優秀そうな式との差の絶対値の合計の10^y倍を表示した。優秀そうな式が合っていればそのまま誤差の式の個数になるから数えやすい。更に別列に誤差累計も表示したので計数途中でセルが移動しても元の位置が分かりやすい。(行の固定部分で計数した。押すと１増えるボタンを作れば効率アップだが今はそんな力はないｗ)
元の数はxより小さい方に106行、大きい方に107行。1行下に行くと1/10^y/z増える。
(x, y, z)=(0, 0, 20)の場合、元の数は-5.30～5.35。
元の数は関数ではなく値(マクロで入力するのが理想的か)の方が良さそうだが、様々な検証をするためには関数の手軽さは捨てがたい。誤差が累積しないようにrow()を元に計算した。元数由来の誤差はなるべく減らしたいが完璧は無理かも。(急ぐ時は「1行上のセル+間隔」とかよくやる)
具体的には B列(8～221行) =(row()-114)*1/10^y/z+x 114行目が中心値
最初の９件 (0, 2, 2)(0, 3, 2)(1, 3, 20)(-1, 3, 20)(1, 2, 20)(-1, 2, 20)(0, 4, 20)(1, 4, 20)(-1, 4, 20) (個数が３桁を超えた物は計数が大変なので検証に含めずスルーしてしまったものもある)
(追記 番外編 グーグルで(0, 6, 2)が175個だった。二重丸め型４つが順に0個57個59個59個。)
追加の３件 (0, 6, 20)(1, 6, 20)(-1, 6, 20)
検証十分とは言えないがまあまあの結果は出せたと思う。
列の順は検証時と同じではなく種類別に並べ替えた。CD276:CP280は誤差の個数。なので「.0」は不要だが念のため表示しておいた。
LibO7は前9件がgoogleと一致していたので省略。LibO3は範囲コピーするだけで毎回ソフトが落ちた。エクセルでもエクセルが何度か落ちた。いくつか開きっぱなしだった仕事用ファイルも一緒に落ちて「修復」させられた。検証用のファイルを別途用意するのが正解と思う。
LibO3=LibreOfficePortable3.4.5LibO7=LibreOfficePortable7.5.4
余談１
銀行丸め4bと6bは気に入ってた4と6で誤差が出てしまったから微修正して改善を期待した。結果は同数で改善なしだった。シンプルに書けるint系で誤差が出てしまったのはつくづく残念。
言い訳になるが(0, 6, 20)の中心±10行目の±0.0000005が元数の丸め結果が±0.000001の誤差だった(その値で同時に誤差だったのはマイナス側の特異点型のみ)。しかし同じ元数である(0, 6, 2)の中心±1行目の±0.0000005が元数の丸め結果はちゃんと0になった。つまり元数の設定方法由来の誤差の可能性がある。
余談２
今回最優秀な式はround系で第二引き数を使った形だった。第二引き数を使わない形の方が精度が上がるのだろうか？後で検証するかも。
余談３
最初は誤った値のつもりで「誤値」という単語を使ったが一般的ではないらしいので「誤差」に変更した(少しは電波度を減らさないとｗ)。ちょっとニュアンスが変わってしまう気がするが間違いではないから問題ない。
2023.06.15 19:26公開 単純ミスなどは適宜修正済＆修正予定
登録:
投稿 (Atom)