三角関数を用いて掛け算を行う その２

三角関数を用いて掛け算を行う　その２

 

 

三角関数表と次の三角関数表を用いて、掛け算を行う方法を紹介した。

　　

 

その際、実際に、この方法を用いて５１５０×３０９０の計算を計算してみたところ、精度よく計算できなかった。

これは掛け算に使われた数の桁数が互いに４桁と非常に短いためではないかと考え、次の掛け算で再挑戦してみることにする。

　　766044443×529919264=405941707425849952

９桁同士の掛け算だにゃ。これならば、きっと正しい答えを導き出してくれるに違いない。

　　

そして、ここで三角表を見て、

　　

であるa、bを探すと、a=５０°、b=３２°、

したがって、

　　

ここで再び、三角対数表を取り出し、cos18°とcos８２°の値を調べると、0.951056516と0.139173101だから、

　　

したがって、

 

　　

左から９桁目の数字が8と7で違っているけれど（しかし、その差はわずかに１）、この計算法の有効数字がこの程度だから、この桁の数字が少し違うのはしょうがない。

科学計算で使用されるFORTRANの単精度（実数型32ビット）の計算の場合、0.01×100は1にならない。

それに比べれば、この計算法は十分すぎるほどの精度を誇っているといって過言はないと思う。

 

「ネムネコ、ちょっと待った！！」

「何だにゃ、おじゃま虫。」

「お前、計算が一致するように、三角関数表にあるsin50°とsin32°の数字を整数化し、766044443、529919264という数字を引っ張りだしただろう。だったら、最初から、三角関数を用いたこの計算が有効数字程度で一致するように仕組まれているじゃないか。」

「たまたまの一致だにゃ。これは偶然。だから、妙な言いがかりは付けないでもらおうか！！」

 

実際、おじゃま虫の指摘通りなのですが、三角関数表と三角関数の積和公式を用いることで、桁数の大きい数同士の掛け算を、簡単な足し算と引き算に還元して計算できるというお話でした。

なお、ネムネコのPCにある電卓ソフトでこの掛け算を行うと、

　　766044443×529919264=4.059417074×10¹⁷

となる。コンピュータは２進数で計算するので、32ビットだと、符号なしの整数で表せる最大の数字は2³²−1=4294967295となって、今回のような桁数の多い整数同士の掛け算の計算を整数として正確に計算できないんだケロ。64ビットでも2⁶⁴−1=1.844674407×10¹⁹だから、今回の掛け算はその限界能力に迫るものになっているのであった。たとえば、

　　7660444431×5299192645

としたら、符号なしの64ビット整数ですら、この掛け算は整数として計算できない。そして、「こんな計算はできないケロ」というエラーメッセージを表示するのであった。（あくまで、これはC言語などで整数型データとして計算した場合の話。）

対数を使うと、この掛け算は次のようになる。

　　

有効数字9桁で計算した場合、この掛け算に関しては、対数を用いた掛け算の計算よりも三角関数を用いた計算の方が精度よく計算されているにゃ。

対数を利用した計算に勝ったケロ！！

 

 

ちなみに有効数字１０桁で計算すると、4.059417074×10¹⁷になる。

 

今日のアニソン、「モンスターストライク」から『We Will Rock You』

今日のアニソンは、「モンスターストライク」から『We Will Rock You』です。

この欲は、どこかで聞いたことのある局だにゃ(^^)

［dy/dx＝0を通過する微分方程式の数値解］ その２

［dy/dx＝0を通過する微分方程式の数値解］　その２

 

 

　次にx＝0で解が止まらなかったとします。これは数値誤差によりy(0)＝ε≠0になったという事です。もしx＝0で解が止まるべきだという積極的理由があるなら、計算はx＝0で打ち切るべきです。じっさい今は(11)から、数値解はx＝0で止まっても良いとわかっています。しかし一般には、そんなこたぁ～事前にわかりません(^^;)。やってみて判断するのが普通です。だとすれば、y(0)＝ε≠0のεはルンゲクッタ法が許す精度内では妥当で、x＝0で解が止まらないのは正解かも知れないのです。もちろん疑ってかかるべきですが。

　ところがyの零点付近というのは、数値解法の生の誤差値が現れるところで、しかも相対誤差も意味をなさないので、εが非常に小さくても、それが本当に0かどうかの判定はかなり難しいのです。いずれにしろ今度もy＝x³にしか目を向けず、分岐解y(x)＝0を考えないというのは数学的に言ってまずいです。

 

　例を挙げますね。(11)は、ある力学問題を表していると仮定します。

　初期値y₀＝y(-1)＝-1，Δx＝2/80＝0.025で計算を開始したところ、y(0)＝ε≠0となりx＝0で解は止まらなかったとします。

　じっさいに計算すると、ε＝y(0)＝1.45×10^-5になります。Δx＝0.025の4次精度なので、yの精度はΔx³＝1.56×10^-5程度のはずです。これらを比べてみると、ほぼ同じです。

　これの意味は、4次のルンゲクッタ法では、ε程度の大きさは0と見分けがつかないという事です。従って正解がε＝0である可能性も非常に濃厚です。でも可能性は可能性です。

　いま我々は理論解を持っているので、初期値y(-1)＝-1を与えればy(0)＝0でなければならない事を知ってますが、普通はそうじゃありません。理論解を出せないからこそ数値計算します。そうするとy(0)＝0となるような初期値の与え方はy(-1-ε)であって、そのずれのためにy(0)＝ε≠0なのかも知れないからです。そうするとルンゲクッタ法は超正解を与えてる事になります。

 

　数学的には(x，y)＝(0，0)で解が分岐する可能性は否定できません。よって図-3のようになってしまった以上、y(0)＝0を初期値として(11)を数値積分した分岐解も探すべきです。ところがその計算は、蟻地獄です(^^;)。

 

　初期条件y(0)＝0から計算をスタートできるようにする、何らかのスターター手段が必要です。

 

　では物理的にはどうでしょう？。(11)は、ある力学問題を表しているのでした。分岐解のどちらかを選べるでしょうか？。力はyの2階微分に比例します。もしdy/dx＝0でもd²y/dx²≠0ならば、yはy(0)＝0から動くと言えます（前回は計算を間違えました(^^;)）。

　　

　(12)よりy＝0ではd²y/dx²＝0なので、速度0の時に力0なら動かないと結論できます。従ってy(0)＝0が真実なら、x＝0以後はy(x)＝0が正解です。しかし図-3は正しいのでしょうか？。その判断はできません。この段階は大学初年級レベルかな？(^^)。

 

　大学学部レベルになると、次のように言い出します。現実の物理系には減衰（摩擦）があるはずだ。よって減衰無しとして計算した結果がε＝y(0)＝1.45×10^-5と非常に小さいなら、そこで止まってしまうはずだと。何故なら(11)よりdy/dx＝0となるのは、y＝0の時のみだからy(0)＝0とみなすべきだ。よってx＝0以後はy(x)＝0が正解と。

 

 

　ところが院レベルになると、また趣向が変わります(^^)。

　摩擦があるから止まるって？。その結論は減衰定数の大きさに依存する。そんな事は、減衰項つきの微分方程式を解いてから言え。それよりも重要なことは、現実の物理系には必ず外乱がある事だ！。

　たとえx＝0以降の厳密解がy(x)＝0であったとしても(x，y)＝(0，0)は不安定停留点だ。それは位相空間で考えればわかる。

　外乱を認めるとx＝0でy＝dy/dx＝0に達しても、外乱によりy(0)は、微小に±δだけ、ほとんど必ずわけもなくジャンプすると考えられます。

　図-4に示すように、(11)のy(x)は常に単調増加です。もしジャンプ先が0＜δ＝y(0)なら、x＝0からy(x)は単調増加し続けます。よって図-3の赤点群は、定性的には正解という事になります。

　微小な外乱δは、微小という事以外は何もわかりません。従ってその理想化として、y(0)＝0を初期値とするy(x)＝0以外の数値解を計算する必要に迫られます。よってスターターが必要です。

 

　結局、数学科にいようと物理・工学系にいようと、スターターが必要になりました(^^;)。

 

　スターターとして自然なのは、件の微分方程式を高階微分方程式に持ち込むやり方です。単振動の時のように上手くいけば、解は自然にスタートしてくれます。じつは(11)に関してはdy³/dx³＝6です(^^)。3階微分まで考慮すれば、y＝0からでも数値解はスタートします。

 

　しかしこれには3つの問題があります。

　　1) 一般に高階微分であればある程、数値精度は悪くなり、計算手続きは面倒になり、計算量も増える。

　　2) 一般に何回微分すれば良いかわからない。

　　3) 最悪の事態として何回微分してもdⁿy/dxⁿが、その点で0かもしれない。

　　（※ あまり知られていない、テーラー展開不能点(^^;)）

 

　例えば(6)を(7)へ持ち込んだとたん、計算量は倍に増えました。3階微分なら3倍です。

 

　もっと泥臭いスターターを考えます。泥臭いという事は、強い方法だという意味です。単振動(6)では、y(0)＝1でなく、y(0.01)＝0.999950000416665を初期値とするだけで非常に高精度の解が得られます。

　そういうものですぐ思い付けるのは、x＝0からΔx離れたx＝Δxにおけるy(Δx)を、逆向きルンゲクッタ法に従ってy(0)＝0となるように逆算する手です。(11)の場合なら、y(Δx)＝y₁として、

 

となりますけれど、上記の非線形方程式をy₁について解くのはまぁ～、あきらめた方が無難ですな(^^;)。ところでΔxはかなり小さいのでした。という事はy₁もyの零点付近の値なので、かなり小さいはずです。Δx×y₁はすごく小さいと予想できます。Δx×y₁について線形化しましょう。

 

　上記近似は、Δxが小さいほど正確なはずです。そこでΔx＝0.001とすると、y₁＝1.259422×10^-10となり、ほぼ倍精度実数の実用精度の限界です。この結果をもとにΔx＝-0.001とした逆向きルンゲクッタを行うと、y(0)＝-1.71×10^-9になります。y(-1)＝-1を初期値とした時の結果、y(0)＝1.45×10^-5よりもだいぶ0に近いです。y₁＝y(0.001)＝1.259422×10^-10からΔx＝0.001でx＝0.025までルンゲクッタを行うと、y(0.025)＝1.459363×10^-5です。以後これを初期値とし、Δx＝0.025でx＝1までルンゲクッタで計算します。

 

 

　そして図-5をいっしょうけんめい観察するのです。そうすると、ある疑念がふつふつと沸いてきますきます(^^)。y(x)はx＝0に関して反対称ではないのか？（！）。じっさい黒点のy(1)はy(1)＝0.996619です。もしy(1)＝1が正解なら誤差は1/1000のオーダーで、これはy(0.001)を逆算した時の計算誤差の影響と考えられます。

 

 

　図-6は図-5の赤点のx≦0側の絶対値と黒点を比較したものです。両者の比はほぼ常に1です。x＝0近傍で1からはずれるのは、両者の値が小さくなりすぎて0/0状態になり比が意味をなさないからだと判断できます。y(x)は、x＝0を境に反対称である可能性濃厚です。もしもそうなら、y(0)＝0でなければなりません。もともとの条件(11)が、y＝0に関して対称だからです。

 

　こうしていちおう、数値的にも次の結論が得られます。

　(11)のy(x)は、初期値y(-1)＝-1に対してy(0)＝0で数学的には解が分岐する。物理的には、どちらを選ぶかは状況による。

 

　でもあくまで「いちおう」です(^^;)。結論を確定させるには、もっと高精度の力業が・・・(^^)。

 

ネムネコのつぶやき

なんか、微妙な表になってしまった

三角関数表を使うと、たとえば、
  3090×5150
という掛け算ができるんだよ、という話をするために、三角関数表をブログにアップしたんだけれど、何とも微妙な表になってしまった。

画像だと小さくて見づらいし、pdfファイルだと書物並みの美しさを誇るけれど、このブログはpdfファイルには対応していない。

有料版のブログにすれば、役立たずのSo-netブログに毎月お金払わないといけないし・・・。困ったもんだにゃ。

ところで、「三角関数表を使うと、掛け算ができる」ってことを知っていたケロか。
ニュートンよりも少し前の時代に、最新の計算技術として、三角関数表を利用した掛け算が実際に行われていたんだケロよ。その当時、１０桁くらいの三角関数表が出版されていたんだから、驚きだよね〜。
いったい、どうやって、計算したんだろう。
必要に迫られていたということもあるけれど、昔のヒトは、現代人と違って、智慧があるよね〜。そして、何より辛抱強い!!
だって、今と違って、電卓やコンピュータという便利なツールは何もないんだぜ〜。この１行、つまり、１°に相当する三角関数の値を求めるのに、最低でも数日、ヘタをすれば、１ヶ月、１年くらいかかったんじゃ〜ないかしら。
ただただ、感心するばかりのネムネコであった。

三角関数による計算法

三角関数表

三角関数表

角度（度）	正弦	余弦	正接	角度（度）	正弦	余弦	正接							π
0	0	1	0	45	0.7071067812	0.7071067812	1							3.1415926536
1	0.0174524064	0.9998476952	0.0174550649	46	0.7193398003	0.6946583705	1.0355303138
2	0.0348994967	0.999390827	0.0349207695	47	0.7313537016	0.6819983601	1.07236871
3	0.0523359562	0.9986295348	0.0524077793	48	0.7431448255	0.6691306064	1.1106125148
4	0.0697564737	0.9975640503	0.0699268119	49	0.7547095802	0.656059029	1.1503684072
5	0.0871557427	0.9961946981	0.0874886635	50	0.7660444431	0.6427876097	1.1917535926
6	0.1045284633	0.9945218954	0.1051042353	51	0.7771459615	0.629320391	1.2348971565
7	0.1218693434	0.9925461516	0.1227845609	52	0.7880107536	0.6156614753	1.2799416322
8	0.139173101	0.9902680687	0.1405408347	53	0.79863551	0.6018150232	1.3270448216
9	0.156434465	0.9876883406	0.1583844403	54	0.8090169944	0.5877852523	1.3763819205
10	0.1736481777	0.984807753	0.1763269807	55	0.8191520443	0.5735764364	1.4281480067
11	0.1908089954	0.9816271834	0.1943803091	56	0.8290375726	0.5591929035	1.4825609685
12	0.2079116908	0.9781476007	0.2125565617	57	0.8386705679	0.544639035	1.5398649638
13	0.2249510543	0.9743700648	0.2308681911	58	0.8480480962	0.5299192642	1.600334529
14	0.2419218956	0.9702957263	0.2493280028	59	0.8571673007	0.5150380749	1.6642794824
15	0.2588190451	0.9659258263	0.2679491924	60	0.8660254038	0.5	1.7320508076
16	0.2756373558	0.9612616959	0.2867453858	61	0.8746197071	0.4848096202	1.8040477553
17	0.2923717047	0.956304756	0.3057306815	62	0.8829475929	0.4694715628	1.8807264653
18	0.3090169944	0.9510565163	0.3249196962	63	0.8910065242	0.4539904997	1.9626105055
19	0.3255681545	0.9455185756	0.3443276133	64	0.8987940463	0.4383711468	2.0503038416
20	0.3420201433	0.9396926208	0.3639702343	65	0.906307787	0.4226182617	2.1445069205
21	0.3583679495	0.9335804265	0.383864035	66	0.9135454576	0.4067366431	2.2460367739
22	0.3746065934	0.9271838546	0.4040262258	67	0.9205048535	0.3907311285	2.3558523658
23	0.3907311285	0.9205048535	0.4244748162	68	0.9271838546	0.3746065934	2.4750868534
24	0.4067366431	0.9135454576	0.4452286853	69	0.9335804265	0.3583679495	2.6050890647
25	0.4226182617	0.906307787	0.4663076582	70	0.9396926208	0.3420201433	2.7474774195
26	0.4383711468	0.8987940463	0.4877325886	71	0.9455185756	0.3255681545	2.9042108777
27	0.4539904997	0.8910065242	0.5095254495	72	0.9510565163	0.3090169944	3.0776835372
28	0.4694715628	0.8829475929	0.5317094317	73	0.956304756	0.2923717047	3.2708526185
29	0.4848096202	0.8746197071	0.5543090515	74	0.9612616959	0.2756373558	3.4874144438
30	0.5	0.8660254038	0.5773502692	75	0.9659258263	0.2588190451	3.7320508076
31	0.5150380749	0.8571673007	0.600860619	76	0.9702957263	0.2419218956	4.0107809335
32	0.5299192642	0.8480480962	0.6248693519	77	0.9743700648	0.2249510543	4.3314758743
33	0.544639035	0.8386705679	0.6494075932	78	0.9781476007	0.2079116908	4.7046301095
34	0.5591929035	0.8290375726	0.6745085168	79	0.9816271834	0.1908089954	5.144554016
35	0.5735764364	0.8191520443	0.7002075382	80	0.984807753	0.1736481777	5.6712818196
36	0.5877852523	0.8090169944	0.726542528	81	0.9876883406	0.156434465	6.3137515147
37	0.6018150232	0.79863551	0.7535540501	82	0.9902680687	0.139173101	7.1153697224
38	0.6156614753	0.7880107536	0.7812856265	83	0.9925461516	0.1218693434	8.144346428
39	0.629320391	0.7771459615	0.8097840332	84	0.9945218954	0.1045284633	9.5143644542
40	0.6427876097	0.7660444431	0.8390996312	85	0.9961946981	0.0871557427	11.4300523028
41	0.656059029	0.7547095802	0.8692867378	86	0.9975640503	0.0697564737	14.3006662567
42	0.6691306064	0.7431448255	0.9004040443	87	0.9986295348	0.0523359562	19.0811366877
43	0.6819983601	0.7313537016	0.9325150861	88	0.999390827	0.0348994967	28.6362532829
44	0.6946583705	0.7193398003	0.9656887748	89	0.9998476952	0.0174524064	57.2899616308
45	0.7071067812	0.7071067812	1	90	1	0