風の迷路

移転しました(でも風の迷路は閉鎖しないよ):http://fermiumbay13.hatenablog.com/

全体表示

[ リスト ]

arctan近似式

ルジャンドル多項式を利用して、定積分を求めるガウス求積という方法を使えるようになりました。
これを利用して、アークタンジェントの近似式を作ります。
アークタンジェントは、1/(1+x^2)の積分で求まり、
3次のルジャンドル多項式を利用したガウス求積では零点が0,±√(3/5)なので、
二乗とルートで打ち消されて簡単な式になるだろうという推測であります。

∫[0→a]1/(1+x^2) dx=arctan(a)

となりますから、左辺の積分を求めることを考えます。
ガウス求積を行うには積分範囲を−1〜1にしなければならないのですが、
0〜aの積分範囲のものを−1〜1に改めるにはx=a/2 (t+1)と置換しなければならないので、
a/2 ∫[-1→1]1/{1+(a/2 (t+1))^2} dt
という積分になってしまいます。
こうなると、零点の√(3/5)をtに代入して二乗を展開したとき、
ルートの項が残ってしまうことになります。結局面倒になります。

x=At+Bと置換したとき、このがあるからこんなことが起こるのです。
よって、B=0になるための条件を求めてみましょう。積分範囲をα〜βとします。
B=0を代入して、αとβの関係式を導きます。
x=αのとき、t=−1より、α=−A
x=βのとき、t=1より、β=A
よって、α=−βであればBの項がなくなります。
積分範囲は絶対値が等しくなるようにしなければならないということのようです。
それでは試しに
∫[-a→a]1/(1+x^2) dx
とするとこれはどうなるかというと、偶関数ですから、2arctan(a)になります。
即ち、
arctan(a)=1/2 ∫[-a→a]1/(1+x^2) dx
として求めることが可能なわけです。
これなら計算がラクになることでしょう。
x=atと置けば、積分範囲は−1〜1に改まります。
dx=adtなので、
arctan(a)=a/2 ∫[-1→1]1/(1+(at)^2) dt
となるわけです。それでは3次のガウス求積を行ってみましょうね。

3次のガウス求積は、重みAi、零点xiとして、
∫[-1→1]f(x)dx≒A1f(x1)+A2f(x2)+A3f(x3)
となるものでした。ここで、
零点はx1=0,x2=√(3/5),x3=−√(3/5)で、
それに対する重みはA1=8/9,A2=5/9,A3=5/9でした。
よって後は代入するだけです。
a/2 ∫[-1→1]1/(1+(at)^2) dt
≒a/2 〔8/9・1/(1+0^2)+5/9・1/{1+(a√(3/5))^2}+5/9・1/{1+(−a√(3/5))^2}〕
=a/9 {4+5/(1+3/5 a^2)}
=a/9 (20+12a^2+25)/(5+3a^2)
(12a^2+45)a/(27a^2+45)
となりました。つまり、

arctan(a)≒(12a^2+45)/(27a^2+45) a[rad]

として求まることになります。二次式÷二次式×aという形です。
実際に何か求めてみると、結構精度よく値が求まることが分かります。
例えばa=1/2のときarctan(1/2)≒0.46365になりますが、この近似式に代入すると
(12/4+45)/(27/4+45) 1/2
=24/51.75≒0.46377 となります。かなり近いじゃないですか。
がんばれば十分手計算で求まる程度の計算量で、精度よくアークタンジェントの値が求められるわけです。
ねえ、きっと使えるでしょう。

aの値が大きくなればなるほど当然誤差が大きくなっていきますが、
a>1の場合は、arctan(a)=π/2−arctan(1/a)の関係式を用いれば良いでしょう。
例えばarctan(50)を求めたいときは、
arctan(50)=π/2−arctan(1/50)ですから、
arctan(1/50)=arctan(0.02)を上の近似式で求めればいいことになります。
すると、arctan(0.02)≒0.0199973と求まりますので、
arctan(50)≒3.141592/2−0.0199973=1.5507987と求まります。
真値はarctan(50)≒1.5507990ですので、これまた近いことが分かります。


ここで与えられる値はラジアンですが、ふだんラジアンなど使うことないでしょうから、
六十分法に改めた版も用意しておいた方が何かと便利だと思われます。
180掛けてπで割れば六十分法になるわけですが、πで割るのがめんどくさいものですよ。
そこで、更に近似していくことにします。

(12a^2+45)/(27a^2+45) a ×180/π
=(2160a^2+8100)/{π(27a^2+45)} a
=(80a^2+300)/{π(a^2+5/3)} a
分子の値をπで割ります。
≒(25.46a^2+95.49)/(a^2+5/3) a
分母分子2倍してみると、
=(50.92a^2+190.98)/(2a^2+10/3) a
≒(51a^2+191)/(2a^2+10/3) a
分母分子を3倍して、
=(153a^2+573)/(6a^2+10) a
≒(150a^2+570)/(6a^2+10) a
(75a^2+285)/(3a^2+5) a

後半は結構強引ですが、つまりこういうことになります。

arctan(a)≒(75a^2+285)/(3a^2+5) a[°]

ラジアン版と似てる形になりました。
やたら近似したけれど、大丈夫なのでしょうか。
例えばarctan(0.8)を求めてみると、
arctan(0.8)≒(75×0.8^2+285)/(3×0.8^2+5)×0.8≒38.497°
真値はarctan(0.8)≒38.660°
小数点以下がちょっとずれてるっぽいですが……でも、思ったより近いです。

グラフを描いてみるとこうなります。

イメージ 1

青い実線が近似式(y=(75x^2+285)/(3x^2+5) x)、
赤い点線が真のアークタンジェント(y=arctan(x))です。
私にはもう、ぴったりくっついているようにしか見えないわけですよ。
最大誤差は約0.18°。普通の人間が0.18°のズレなんて分かりますかね。

イメージ 2

おそらくアークタンジェントの近似式の中では計算しやすい方でしょう。
何度か眺めていれば覚えてしまいますね^^!
これを使って日常の色々なものの角度を測ってみましょう。


本当に、何に使えるんだろう……

この記事に

閉じる コメント(2)

顔アイコン

arctanというとSTGには必須というぐらいの認識でしかないです(^_^;)
ルートは面倒なので、消えてくれるのは嬉しいです……と言うと可哀想ですね。

ガウス求積の前に大変な計算になっておられます@@ 邪魔者を消すのは簡単ではないのですね
むしろ前々回の記事で重みが求めていた、ガウス求積の方が簡単なのでしょうか……
おおっ、完成した式はとてもわかりやすい形をしています^^
これなら無理なく暗記もできそうなので、日常的に使えちゃいますね(?)
誤差も人の感覚では分からないほど小さく、グラフで見ても同じ形にしか見えませぬ@@

日常的ではないですが、シンプルな式なのでツクールで作りやすそうな気はします。

2013/1/5(土) 午前 2:46 [ Dik ] 返信する

顔アイコン

まさにそんな感じです^^; 縦横の長さ比から、成す角度を求めるものです。
ルートは消えてくれた方がやはり嬉しいですよ。可哀想ですが……;;

下準備がやたら面倒なのがネックですね……まあ、結果のためです。
この計算はほとんど近似なので、難しいも何もありません。すべてノリです。
ラジアンのは分母と分子の定数項がそれぞれ45になってくれてきれいになりました;
思い出したときに、定規と計算紙があれば分度器なしで角度が求まります^^;
何かやるときにはこの誤差は気になりますが、角度求めるぐらいなら……

ををを、良いかもしれません。かなり軽く計算できそうです^^;

2013/1/5(土) 午後 10:30 フェルミウム湾 返信する

コメント投稿

顔アイコン

顔アイコン・表示画像の選択

名前パスワードブログ
絵文字
×
  • オリジナル
  • SoftBank1
  • SoftBank2
  • SoftBank3
  • SoftBank4
  • docomo1
  • docomo2
  • au1
  • au2
  • au3
  • au4
投稿

.


プライバシー -  利用規約 -  メディアステートメント -  ガイドライン -  順守事項 -  ご意見・ご要望 -  ヘルプ・お問い合わせ

Copyright (C) 2019 Yahoo Japan Corporation. All Rights Reserved.

みんなの更新記事