NO.1748 アイゼンスタイン三角形 2008.10.26. 水の流れ
第215回数学的な応募問題
皆さん、高校数学で習う余弦定理をご存知ですか。三角形の3辺の長さをa、b、cとしたときに成り立つ定理です。
c2=a2+b2−2abcosθのことです。
辺cに対する角Cが120°のとき、 c2=a2+b2+abとなる。
例えば、a=3、b=5、c=7がそうです。別名七五三の三角形と私は呼んでいます。
最近、「整数とあそぼう」(日本評論社、一松信著)を読んでいたら、このような3辺の長さが整数で、角が120°である三角形を「アイゼンスタインの三角形」と知りました。
ここで、x=a+b,y=a(またはb),z=cである三角形はz2=x2+y2―xyを満たし、角が60°になります。これがx=8、y=5、z=7となり、別名ナゴヤ(名古屋)の三角形と私は呼んでいます。
さて、c2=a2+b2+ab整数解ですが、
a=m2−n2,b=2mn+n2,c=m2+mn+n2となります。
ただし、mとnは互いに素で、m−nは3の倍数ない。
ここから、問題です。
問題1:m、nに自然数を小さい順に代入して、(a、b、c、a+b)組を幾つか求めてください。
問題2:a=m2−n2,b=2mn+n2,c=m2+mn+n2のとき、c2=a2+b2+abを満たしていることを確かめてください。
問題3:積abc(a+b)は必ず840の倍数になることを確かめてください。
問題4:可能なら、c2=a2+b2+ab整数解が、
a=m2−n2,b=2mn+n2,
c=m2+mn+n2となることを求めてください。
注:アイゼンスタイン(1823〜1852)はドイツの天才数学者です。
注:この記事に関する投稿の掲載は、2008年11月17日以降とします。
NO.1747 楕円積分について 2008.10.26. 三角定規
NO.1741において,本多欣亮さんはおそらくすべてをご存じの上で書いておられるのではないか,と推察致しますが…
各点及び座標,角θを右 <図1> のよう に定めると,
ですから,楕円上の弧 PR (青線部分) の長さ s は
で与えられます。これの積分部分が,第2種の(不完全)楕円積分 E (k,θ) です。
よく知られているように,この積分は解析的に実行することができません。私が学生の頃は,関数電卓が高嶺(高値)の花,世にパソコンなどありませんでしたので,E (k,θ) を母数 k について級数展開し,
としたり,収束が遅い(3)に「ランデン変換」とよばれる母数の変換を施して収束を速めたりといろいろな工夫をしたものでした。ですが,近頃では 高精度計算サイト のような数値計算のためのサイトが手軽に利用できるようになり,「昔の苦労は何だったのか!」と感慨ひとしおではあります。
ところで,本多さんがお書きの,楕円上の点における法線が x 軸となす角φは,下 <図2> の楕円を地球の断面とするとき,天文緯度 または 測地緯度 とよばれるものです。図の θ は球面の 地心緯度 とよ
ばれ,φとの間には,
の関係があります。ですから,赤道Sから緯度φの地点Pまでの距離を求めるには,
<1> (4)から でθを求める。
<2> 計算サイトを利用し を求める。
<3>
とするのが現実的だと思います。(2)をφについての積分に変換しても, 最終的には級数展開に頼らざるを得ないのですから。
話は変わりますが,かのニュートンは,力学的な考察から「地球は赤道方向にふくらんだ回転楕円体」と考え,著書『プリンキピア』の中で述べています。それを実証するには,(天文)緯度差1°対応する子午線の長さが,高緯度地方ほど長くなることを測量で示せばよいのですが,これは容易なことではなかったようです。
フランスは,1735年から大がかりな観測隊を北極に近いラプランドと赤道に近い南アメリカのペルーに派遣し,精密観測を敢行しました。そして,苦節8年の後,ようやく緯度差1°に対応する子午線の長さが前者の方が長いことを確認し,ニュートンの考えが実証されました。この測量は,18世紀科学史上の一大壮挙と言われています。
NO.1746 極限値は?(3) 2008.10.26. DDT
漸化式、
の極限値は?、という問題で、「上に有界な単調増加列は収束する」というヒントが付いています。
まず{an}が単調増加になるan の範囲を求めます。an>0 に注意して、
とすれば、
x2 - x - 1 ≦ 0 (1)
なので、
α2 - α - 1 ≦ 0 (2)
より、
で、an が単調増加になる範囲は、
とする.ここで2<51/2 ,(1−51/2)/2<0なので、
から、a1=1 はこの範囲に含まれます。そうすると、
1 ≦an≦α
であれば、an は単調増加です。そこで、n≦mではan が単調増加で、かつ、
となるようなmが存在する、とします。このようなam の範囲は、
でなければならないので、
α2 - 1 < x
より、
α2 - 1 < x ≦α
が必要です。ところがαは(2)を満たすので、α2−1=α であり、
α < x ≦α
となるので、このようなam の範囲はありません。従ってan はαを越える事はなく、
αを一つの上界として持つので、an は上に有界な単調増加列で収束します。
とすれば、
(√の双連続性を使っています)
なので
すなわち、
β2 - β - 1 = 0
ですが、またまた(2)が現れて、
が得られます。
ところで何故、{an} が単調増加になる範囲を定める(2)が、極限値を定める条件にも出てくるのか?、
というのは気になるところです。ところですが、図を描けば一目瞭然です。
この図を見ていると、{an} が単調増加になる範囲を定める(2)が、極限値を定める条件にもなる事で本質的なのは、数列を定義する関数、
と、数列の再帰処理を定義する関数、
g (x) = x
のグラフの形だとわかります。よって今回の問題は、次のように一般化できると思います。
[定理?]
f:R→Rが、
α2 = f (α2) (3)
を満たして、α1≦x≦α2 の範囲(α1<α2,α1=−∞も許す)でf'(x)>0 とする。
さらに、
x ≦ f (x) (4)
も成り立つとする(この条件でf(x)は、α1≦x≦α2 の範囲でf(x)≧x の上に凸な単調増加関数になる)。
数列{an}n=1,2,・・・ の初期値a1を、α1≦a1≦α2 の範囲にとり、
an+1 = f (an) (5)
で{an} を定義するなら、
(6)
が成り立つ。
[証明]
α1≦a1≦α2 なので(4),(5)より、
α1≦an≦α2 の範囲で、an は単調増加。
そこで、n≦mではan が単調増加でかつ、
am≦α2 , α2≦am+1
となるようなmが存在すると仮定する。このようなx=am の範囲は、
α2<f (x)
でなければならない。x≠α2 とすれば、f'(x)>0なのでf はx で真に単調増加で逆関数を持ち、
f-1(α2) < x
が得られる。従って、
f-1(α2) < x ≦α2
が必要。ところがα2 は(3)を満たし、x=α2 でもf'(α2)>0 なのでやはり逆関数があり、(3)より、
f-1(α2) = α2 (3)
も成り立つ。従って、
α2 < x ≦α2
となるが、このようなam の範囲はない。よってan はα2 を越える事はなく、
α2 を一つの上界として持つので、
an は上に有界な単調増加列で収束する。
x=am=α2 の場合は(3)より、
am+1=f (am)
となって、明らかに、
が成り立つ。よってam≠α2 の場合に、上記を示す。
am≠α2 の場合、任意のan でα1≦an<α2。
よって任意のx=an でf(x)はf'(x)>0 であり、x=an で連続な単調増加関数なので逆関数を持ち、逆関数も連続。従ってf はx=an で双連続(同相)なので、
limn→∞ an が収束するなら、上記も成り立つ。limn→∞ an の収束は、
上で示した。
とすれば、
β=f (β)
となる。an はα1≦an<α2 であるから、an は、α1≦β≦α2 以外のβには収束できない。従って(3)より、
β=α2
が必要で、
が得られる。
[証明終]
上記証明が正しいとしてですが、この図式って、弾塑性解析で標準的に用いられる、修正ニュートン・ラプソン法の根拠だよな、って思いました。勉強になりました・・・。
図-1は弾塑性解析におけるニュートン法の概略図で、図-2は、修正ニュートン・ラプソン法の概略図です。前の図と良く似ています。
上記証明において、f(x)が単調増加である事は本質的です。この条件を外すと「どえらい目」に会います。図-3は、歪み軟化を起こしてf(x)に(応力−歪み曲線に)単調性がなくなった場合です。この場合は、解は2つになるは、収束経路がぐるぐる回って安定しないはで、自分でやろうとすると「こっぴどい目」に会います。
そこで解析ソフトを使うわけですが、増分ステップや計算方法、収束パラメータの与え方によって、なかなか言う事を聞いてくれません。この「電子計算機全盛の時代」に、「勘と経験だけが頼り」という、なかなかトホホな事態に巻き込まれます。でも答えは出さないと、給料出ないし・・・。
NO.1745 極限値は?(2) 2008.10.26. 夜ふかしのつらいおじさん
● 数列{an}の各項の値は正であること
漸化式の成り立ちから、この数列の第2項以降の値は正です。
(正の平方根をとっている)
● 階差について
右のグラフより、 が大小の境目になります。
(点AのX座標、この値をαとする)つまり、
● 数列{an}が単調で有界であること
これを次のように略記します。
この式について、分子のつなぎの符号が「+」であるものを分子・分母に次々とかけ漸化式を用いて変形していけば、
この式の分母は、正です。
k=1,2,3,4,・・・,n-1 とします。
すると、この結果は、an−αとak−αとが同符号であることを示しています。
つまり、数列{an}の各項の値は、αよりすべて大きいか、小さいかとなります。
だから、この結果と階差の結果をあわせると、数列{an}が単調であることが分かります。
つまり、a1>αなら単調減少、a1<αなら単調増加となります。
a1>αのとき、この数列は単調減少です。
このとき、a1−α>0です。
単調減少なのに、an−αの符号が正のままであることは、
{an}が下に有界であることを示しています。
a1<αのとき、この数列は単調増加です。
このとき、a1−α<0です。
単調増加なのに、an−αの符号が負のままであることは、
{an}が上に有界であることを示しています。
● 収束する値のこと
次の関係を考慮します。
すると、
この式で、n→∞とすると、 なので、
anは、αに収束します。
問題は、a1=1ですが、初項の値にかかわらず、 この数列{an}は、αに収束します。
◎ 結果だけ出す場合
収束することを前提にすれば、nが大きなときは、anもan-1も差がなくなり、
次の方程式の解が収束値となります。
(ただし、正の解)
◎ 収束の様子を視覚的に見る場合
NO.1744 ベクトル解析Minimum(1) 2008.10.26. DDT
最近ようやっと、電磁気学を少し本気で勉強する気になりました。・・・って、今まで知ったかぶりして色々書いてきたくせに、本気じゃなかったのかよ!・・・という自己突っ込みは、いちおう入れておきます。現状はと言えば、まぁ、大学の教養時代に取った物理と、電磁気にも関係する計算もしていたうちに、なんとな〜く、簡単な問題くらいなら計算できるようになった、という程度です。
本格的な電磁気学の本って、なかなか読む気になれないんですよね。代表的なところでは、吉岡書店から出ているジャクソンの電磁気学(上下)がありますが、9800円×2(700ページ×2)なんてのを見せられると、かなりやる気を削がれます。とは言うものの、もっとあっさりした本を選んで、それなりに本気で読んでみると、電磁気学って基本原理に限ってなら、ものすごくコンパクトにまとまるのがわかります。例えば静電現象なんて、本当に基本原理だけなら数ページで済みます。そこで自習ノートの一環として、ベクトル解析Minimum,電磁気学Minimumというのを、勝手に企画しました。
電磁気学は、力学以上に初心者に対する敷居が高いと言われますが、それはベクトル解析の知識がないと、もう計算が鬱陶しくてしょうがないからだと思います。また理論形式をすっきりさせるためには、デルタ関数の扱いも必要です。これらは慣れてしまえば単なる計算道具として納得もできますが、最初はやっぱり、全然イメージが沸かなくて、見通しが立たなくて難議します。というわけで、最初にベクトル解析とデルタ関数のイメージから入ります。
1.∇(勾配).
2.ポテンシャル.
3.流量保存則.
4.div(発散)とGaussの発散定理,Δ(ラプラシアン).
5.δ関数(デルタ関数)とラプラス方程式の基本解.
6.偏微分の順序交換とrot,保存場,積分記号下の微分.
その後で、電磁気学Minimum(1)
という順序です。
1.∇(勾配)
3次現空間の点(x,y,z) を独立変数とする関数φ(x,y,z) を考えます。φ(x,y,z) 自体はある数値(スカラー関数)とします。∇はナブラと読みますが、gradと書かれる時もあります。日本語訳では勾配です(英語でgradiant)。
∇とは、φ(x,y,z) の微分をカッコ良くベクトル形式にまとめただけの話です。とりえず1変数φ(x)で考えます。y=φ(x)の微分はよく、
と書かれますが、厳密に言うと、これは微分ではなく微分係数です。その点での接線の傾きを表します。微分係数の定義に戻れば、
というイメージです。
(1)
が関数φ(x)の微分です。(dx,dy)とは、(h,φ(x0+h)−φ(x0))の極限で、h=(x0+h)−x0 である事から、それが、関数の一点(x0,φ(x0))からの微小な関数の増分を表すのは明らかと思います。(1)は、微小な増分dxに関しては、φ(x)は一次関数になっていると言っています。その一次関数の傾きが、
です。このイメージは膨らます事が出来ます。(1)の(dx,dy)を微小とせず、無制限に延長してやった場合、y0=φ(x0)として、増分(dx,dy)は、(x−x0,y−y0)と書けます。その時(1)は、
(2)
となり、関数の一点(x0,y0)における接線になります。接線とは関数の一点を目に見えるように拡大したものであり、逆に接線があるなら(微分可能なら)、関数の一点には微小な接線(2)が埋め込まれていると考えても良いと、(1)は言っています。このイメージをもろに書けば、
(3)
という事になります。一次関数なら(3)は、直線上の任意の点で無制限に成り立ちます。
(3)のイメージでφ(x,y,z) の微分を扱います。φ(x,y,z) の微分の目的も一変数の場合と同じです。局所的に関数を一次関数化(線形化)する事です。(dx,dy,dz) が微小だとして、増分、
(4)
を考えます。一変数の場合と同じように、(4)をまずxについて線形化します。
となるのは明らかと思います(一変数の場合と同じです)。ここで何故、偏微分∂φ/∂xかと言うと、
と、xだけが変化したと考えているからです。次にφ(x,y+dy,z+dz)のyについても同じ事ができます。
今度は、
と考えている事になります。zについても同じです。
となり、
(5)
が得られます。(dx,dy,dz) が微小な事を考慮して、(dx,dy,dz)→(0,0,0) の極限で考えると(ちょっとずるいですけど)、
なので、これを省略して、
(6)
と書き、全微分と呼びます。意味は微分と同じです。ここで(5)において、(dx,dy,dz)→(0,0,0) の極限を取る時、(x,y+dy,z+dz) のdyやdzは0にするくせに、なんで∂φ/∂xにかかるdxなんかは0にしないんだ?、という疑問は至極当然だと思います(ここがずるい所です)。
じつはここには、「dx,dy,dzなどの微小要素は残したい!」という目的意識が働いています。というのは、dx,dy,dzを残しておくと、直感的な計算を行えるからです。
明らかな例として、ある体積中Vの質量mと、その密度ρとの関係をあげます。
(7)
(7)は、意味を考えれば明らかな式ですが、「意味を考えなければならない」というところがポイントです。一方、微小体積dv の微小質量dm は、
(8-1)
なんだから(密度×体積)、全質量mは、それらの合計で、
(8-2)
となるに決まってる、という考えもできます。この例では僅かな差かも知れませんが、(7)の意味をいきなり考えるより、(8-1)→(8-2) と機械的に計算する方が、概念として遥かに易しいとは思いませんか?。もっと複雑なケースでは、微小要素をあからさまに残しておかないと、意味をつかむのに非常に難儀する場合もありえます。このために(6)でも、dx,dy,dzを残します。ただし(7)から明らかなように、最後の結論mにまで、微小要素dv が残るのは不可です。
そこで微小要素部分と、最後の結論にまで持ち越して良い部分とを、最初から分けておきます。(6)は、
(9)
と書けます。ここでは、横ベクトルと縦ベクトルをこの順に並べたものが、ベクトルの内積を表すものだと、約束します(ただし直後から、内積記号・も使います)。
(9)の横ベクトルと縦ベクトルを、
と書く事にすれば、
(10)
となって、微小要素部分とそうでない部分を明確に分離できます。このためには、ベクトル記法の導入が必要でした。これがベクトル解析を使う理由の一つです。
1°ベクトル記法によって、微小要素部分とそうでない部分とを分離する.
ところが(10)には、もっと意味を付与する余地があります。∇自体を形式的に、
というベクトル(ベクトル作用素と呼びます)だと定義すると、∇φとは、ベクトルかけるスカラーで、
と書けるので、(10)はベクトルかけるスカラーのベクトル∇φと、ベクトルdrとの内積なので、結果はあるスカラー dφ になると、一瞬で読めるようになります(慣れてくればですが)。ベクトル解析を使う二つ目の理由が、このカッコ良さです。
2°ベクトル記法によって、複雑な式の判読が容易になる.
ところで
には、φ(x,y,z) がその各点で、どれだけ急峻に変化するか(どれだけ傾いているか)の情報が全て含まれています。何せこれは、微分係数ですから。そういうわけで、∇は勾配と呼ばれます(1次元なら傾き)。
2.ポテンシャル
例としては、やはり流体力学(水の流れ)がわかりやすいと思います。流体力学においておよそ知りたい事とは、各時刻tにおける空間の各点(x,y,z) での流速(流れの速さ)、u=(vx,vy,vz) です。これさえわかれば、水の圧力から波の形にいたるまで、原理的には全て計算可能です。で、その運動方程式は・・・、
(1)
なんて事になります。ρは流体の密度で、pは圧力です。これにもベクトル解析が使われていますが、いずれにしても見たくもない式である事には、変わりありません。実際(1)を成分で書き下せば、(vx,vy,vz,p) の4変数の連立偏微分方程式になります(見たくない、見たくない・・・)。
流体力学で最初に発展した分野は、非圧縮性渦なし完全流体と呼ばれる領域でした。それには理由があります。速度ポテンシャルφの存在が明らかになったからです。速度ポテンシャルφとは、
(2)
というものです。(1)では、4変数連立偏微分方程式であったものが、速度ポテンシャルφの導入により、たった一個のスカラー関数を求める問題に帰着できます。しかもその支配方程式たるや、(1)に比べれば、
という、あっけないくらい見やすいものでした(Δについては後述します)。
一般に(2)の形で、求めたいベクトル値を与えてくれるスカラー関数を、ポテンシャルと言います(かなり数学よりの定義です)。求めたいベクトル値に対してポテンシャルが存在すると、計算が格段に楽になります。
またポテンシャルはふつう、何らかの保存則と関係があります。物理では、保存則によって理論の大半以上が決まったりします。なので、保存則に関係したポテンシャルはどしどし導入し、計算の大半を楽に済まそう、となるのは必定です。ポテンシャルは大変重宝がられます。
3.流量保存則
再び流体力学から例を取ります。
2次元で一方向だけに一様に水が流れ、その方向にx軸を取ったとします。x軸に平行に長方形V(体積V)をとります。また一方向だけの一様な流速はu(スカラー)とします。
図-1に示したような流れを、一様等速流と言いますが、これを「ところてん流れ」と称す人もいます。ところてん流れによって、体積Vから1秒間に流れ出す流量を考えます。縦辺の長さをL(m)だとすると、それは点線で囲った面積になりますが、その一部の幅wを取り出してみると、
です。最後の ×1(m) は、奥行き1(m) で考えてるという事です。一様等速流なので縦辺L上で、そこからの流れ出しの強度(密度)は、L上の場所によらず一定のはずです。
従って、体積Vから1秒間に流れ出す流量速度は、
と計算できます。次に体積Vに対して、流れがθ傾いている場合はどうでしょう?。
図-1と図-2のハッチ部が横切る辺の長さを比べればわかるように、同じ流出面積(体積)wuに対して、流出長さ(断面積)は、1/cosθ倍に増えている事がわかります。これは単位長さ(単位断面積)当たりの流出強度(密度)がcosθ倍に減った事を意味します。1秒間の流出効率もcosθ倍に減るはずです。実際、図-2の点線部の縦辺の長さは、Lcosθです。よってこの場合に、縦辺Lを通じて1秒間にVから流れ出す流量は、
となります。ところでLは任意に小さくとれます。Lが微小な場合を考え、それをdLで書くことにします。
さらに、
(1)
と変形してみると、u cosθは左辺の意味から、dLから流れ出す流量の、流出密度(強度)である事がわかります。ところでθとは、一様流速uと縦辺dLの成す角です。uとdLをベクトルと考え、dLに(Lに)右回りに直行する単位ベクトルn を導入すれば、(1)は、
と書けます。再びdLを右辺に移項すると、
となります。
|dL|nを線素ベクトルと呼び、dsで表わすのが慣わしです。dUは、微小な直線dLから流れ出す、単位時間当たりの流量だと考えられます。3次元の場合は、2次元のとき|dL|が、微小な直線(曲線)の長さを表わすのに対応して、微小な平面(曲面)の面積を表わすので、dsは面素ベクトルと呼ばれます。dsは、2次元でも3次元でも、考えている面積や体積の外法線方向に取られないと、以上の議論は成り立たない事は注意します。内法線方向だと、流出が流入になってしまうからです。
ここで任意の曲線や曲面は、微分の考えを用いて線形化すれば、局所的には微小直線や微小平面と考えられる事に注意します。従って、
(2)
は、任意の曲線や曲面の各点で成り立っている状況を表わした式だと言えます。
以上より(2)があれば、任意の閉曲線で囲まれた面積や、任意の閉曲面で囲まれた体積から流れ出す、任意の流れに対する流量速度を表わす式は、容易に導けます。それには、1.の(8-1)→(8-2)の過程で使った質量と密度の関係に関する、機械的な計算手順を思い出すだけで十分です。明らかに、
(3)
となります。ここでSは、流量速度Uを考えている任意の閉曲線や閉曲面です。
(3)は、ある体積Vを囲む閉曲面Sから流れ出す流量速度を表しています。今まで流れ出すと言ってきましたが、u・dsは符号を持っているので、場所によって流出だったり流入だったりします。
(3)は、それらを閉じた曲面上で合計したものですから、Uは体積Vからの流出,体積Vへの流入の収支を表しています。
体積Vは一定です。従って、水はVに入った分だけ、出て行かなければなりません(水は圧縮されない)。よって通常は、
(4)
です。もし、
なら、Vのどこかで水は吸い込まれている事になります。また、
なら、Vのどこかに蛇口が付いていて、そこから水は湧いてる事になります。湧き出し(蛇口)の符号を変えれば吸い込みにもなるので、以後は湧き出しのみ考えます。
図-3の状態を考えます。蛇口から湧いた水は、もとからVにあった水を押し出すので、必ず湧いた分だけ流出が生じます。
一方、外部からの流入も同じ理由で、流入した水がV内でどんな経路を流れようと、入った分だけ出て行きます(要するに(4)です。逆に言えば、出て行った分だけ入って来る)。
よって蛇口からの湧き出し速度をΡとした場合、
になると結論できます。
また、湧き出しが何個あっても、結局は体積の収支で、それは一定(V)でなければならないので、
です。
さらにΡが連続的に分布している場合も考えられます。この場合は、Ρでなく、湧き出し密度ρを使うべきです、微小体積からの湧き出し量は秒当たり、
(5)
に違いないから(1.(8-1)と同じ)、全体としては、
(6)
に違いないとなります(1.(8-2)と同じ)。結局(6)は、湧き出しもある場合の流量の収支を表しているので、流量保存則と呼ばれます。
4.div(発散)とGaussの発散定理,Δ(ラプラシアン)
3.(6)で、V(S)は任意でした。Vを(Sを)任意に小さく取った極限を考えれば、3.(5)を再現できるはずです。3.(6)の最右辺は、明らかに3.(5)の右辺になります。では最左辺は?。そのために図-1を考えます。
図-1では、点(x,y) における流速を(vx,vy) とし、各辺の長さdvx,dy は、最初から微小としておきます。図-1で、3.(6)の最左辺の計算を行います。図-1の長方形Vは、3.(6)の最右辺に出てくる積分領域Vに相当します。従って、3.(6)の最左辺に出てくる積分領域Sとは、図の辺@〜Cです。
3.(6)の u・ds のuとは、S上の流速です。(x,y) における流速は(vx,vy) ですが、dxとdyが微小な事から、V全体で流速は一定だとみなしても構わないくらいです(どうせ極限をとる)。そうすると@上での流速uも、u=(vx,vy) とおけます。
3.(6)の u・ds のds とは、Sの線素ベクトルです。線素ベクトルの定義は、ds=|dL|nで、|dL|は考えている曲線の微小要素(線素)でした。いま辺@は微小なので、|dL|=dy とおけます。一方nは、考えている曲線の外法線方向の単位ベクトルでした。辺@のそれは、(−1,0) です。よって@上では、
となりますが、実際には、書いていないz方向のdz もあるので、
が、3次元での結果になります。
辺Aについて、同じ事をすれば、
となります。ここでvx が、
に変わったのは、1.で述べた微分の考えを使っています。@+Aから、
が得られます。B+Cからも、
が得られ、実際にはz方向に、D,Eもあるので(書いてないけど)、
も得られます。3.(6)の最左辺の積分とは、これらの合計です。何故なら、2次元の長方形は、3次元では直方体になり、直方体は6面体なので、@+A+B+C+D+Eです。
一方、3.(6)の最右辺の積分は、Vを任意に小さくした場合、ρdvになりますが、
図-1においてdv とは、dv=dxdydz の事なので、
すなわち、
(1)
が結論です(ここで極限を取った事になります)。
(1)の左辺をカッコ良くしましょう。
これは、ベクトル作用素div を、
で定義した、という事です。よって(1)より、
(2)
と書けます。(2)と3.(6)の物理的意味は、明らかに同じなので、これも流量保存則と呼ばれます。
(2)は本質的に、3.(6)、
と同じです。そうすると(2)においては、ρが一点の湧き出し密度であるように、div・u は、一点の流出速度密度である事になります。一点に水が、強度ρで供給されるなら、それはその一点から、div・u の勢いで散逸していく(吹き出していく)というのが(2)です。散逸(吹き出し)を発散と言ってもいいでしょう。divはdivergenceの略で、日本語訳では発散です。(2)と3.(6)と結び付けます。(2)を体積Vで積分すると、
すなわち、
(3)
が得られます。(3)を Gaussの発散定理と言います。ところで、どうでしょうか?。微小要素を0とおかないで残す理由は、わかって頂けたでしょうか?。
div とラプラシアンΔとの関係は、次のようになります。2.で述べたように、速度ポテンシャルφとは、
(4)
となるスカラー関数です。(4)を(2)に代入します。
(5)
が得られます。中辺の、
を、
と書き、Δをラプラシアンと呼ぶ事にすると、(5)から、非圧縮性渦なし完全流体の支配方程式、
(6)
と、ベクトル解析の公式、
(7)
が得られます。
(7)からは自然に、ベクトル作用素間の内積、
(8)
が導かれ、この結果として、Δをスカラー作用素とみなせる事、
も出てきて、スカラーかけるスカラーのΔφは、スカラーρになると、読めるようになります。
えっ?、div=∇ ではないかって?。・・・その通りです((8)を再度見てみて下さい)。よって、
(9)
もベクトル解析の公式です。ただし、(9)の最後の∇はスカラーに作用する事を想定し、最初の∇はベクトルに作用するのを前提としています。内積を通じてベクトルに作用する∇は、divと書くのが慣わしです。それは、流量保存則(2)の導出過程において現れました。
(2)は流量保存則という、流体力学の根幹をなす法則を表わし、しかも今までやって来たように、その導出は非常に明確な物理的イメージに支えられたもので、その過程で div は開発されたのでした。ここにベクトル解析を用いる第三の理由があります。
B 式を見ただけで、物理的イメージが沸く.
がその理由です。これも慣れれば、ではありますが。
以上をまとめます。
ポテンシャルはふつう、何らかの保存則と関係すると、2.で述べました。そして今それは流量保存則で、Gaussの発散定理という数学的定理を導くものでした。しかもその過程で、div という新たなベクトル作用素も開発されました。よって、
・ポテンシャルφの勾配、∇φ
・ 保存則
または、
・ Gaussの発散定理
・ ベクトル解析の公式
は一般的に、組になって頻繁に現れます。
NO.1743 極限値は? 2008.10.8. 水の流れ
第214回数学的な応募問題
注:この記事に関する投稿の掲載は、2008年10月27日以降とします。
NO.1742 最小公倍数(2) 2008.10.8. 夜ふかしのつらいおじさん
先ず、易しい場合を考えます。
(1)
pを素数とします。
A=p の形のとき、求める整数の組は、(1,p) のみです。
N(p)=1
(2)
A=pk の形のとき、
求める整数の組は、次の通りです。
(1,pk), (p,pk), (p2,pk), ・・・・・・, (pk-1,pk)
N(pk)=k
(3)
一般の場合が見通せるように、Aが素数7個の場合を見てみます。
a,b,c,・・・・・・を素数とします。
A=abcdefg とします。
1°Aの約数は、27個あります。
求める整数の組で、(Aの約数, A) の形のものは、
(1,A),(a,A),(b,A),・・・,(ab,A),(ac,A),・・・,(bcdefg,A) [(abcdefg,A) 以外]
のように、27−1 個あります。
2°
・第1成分が[1]個の素数、第2成分が[6]個の素数の積の場合は、
(a,bcdefg), (b,acdefg), (c,abdefg), ・・・, (g,abcdef) のように7C1個あります。
・第1成分が[2]個の素数、第2成分が[5]個の素数の積の場合は、
(ab,cdefg), (ac,bdefg), ・・・, (bc,adefg), (bd,acefg), ・・・, (fg,abcde) のように7C2個あります。
・第1成分が[3]個の素数、第2成分が[4]個の素数の積の場合は、
(abc,defg), (abd,cefg), ・・・, (acd,befg), (ace,bdfg), ・・・, (efg,abcd) のように
7C3個あります。
※ 積の値が第1成分の方が大きくなれば、第2成分と入れ替えると考えておきます。
さて7C0+7C1+7C2+7C3+7C4+7C5+7C6+7C7=27 です。
二項係数の値は左右対称で等しく、両端の値は1です。
だから、[2]の3つの値の合計は、
7C1+7C2+7C3=27/2−1=26−1 です。
3°
2つの成分に共通因数がある場合を考えます。
a)
・共通素因数が1個で、その部分の素因数の個数が、第1成分[1]個と第2成分[5]個の場合は、
(a・b, a・cdefg), (a・c, a・bdefg), ・・・, (a・g, a・bcdef)
(b・a, b・cdefg), (b・c, b・adefg), ・・・, (b・g, b・acdef)
・・・・・・
(g・a, g・bcdef), (g・b, g・acdef), ・・・, (g・f, g・abcde)
たては共通の因数の選び方なので7C1個、横は残りの第1成分の因数の選び方なので6C1個あります。
だから、7C1×6C1個あります。
・共通素因数が1個で、素の部分の素因数の個数が、第1成分[2]個と第2成分[4]個の場合は、
(a・bc, a・defg), (a・bd, a・cefg), ・・・, (a・fg, a・bcde)
(b・ac, b・defg), (b・ad, b・cefg), ・・・, (b・fg, b・acde)
・・・・・・
(g・ab, g・cdef), (g・ac, g・bdef), ・・・, (g・ef, g・abcd)
たては共通の因数の選び方なので7C1個、横は残りの第1成分の因数の選び方なので6C2個あります。
だから、7C1×6C2個あります。
・共通素因数が1個で、素の部分の素因数の個数が、第1成分[3]個と第2成分[3]個の場合は、
(a・bcd, a・efg), (a・bce, a・dfg), ・・・, (a・bfg, a・cde)
(b・acd, b・efg), (b・ace, b・dfg), ・・・, (b・afg, b・cde)
・・・・・・
(g・abc, g・def), (g・abd, g・cef), ・・・, (g・aef, g・bcd)
たては共通の因数の選び方なので7C1個、横は残りの第1成分の因数の選び方なので6C3/2個あります。
※ 横は注意が必要です。同じ個数に分けるときは、総入れ替えしたものも同じなので2で割ります。
だから、7C1×6C3/2個あります。
さて6C0+6C1+6C2+6C3+6C4+6C5+6C6=26 です。
二項係数の値は左右対称で等しく、両端の値は1です。
だから6C1+6C2+6C3/2=26/2−1です。
よって、[3]a)の3つの値の合計は、
7C1×6C1+7C1×6C2+7C1×6C3/27=7C1×(6C1+6C2+6C3/2)=7C1×(25−1)
b)
・共通素因数が2個で、素の部分の素因数の個数が、第1成分[1]個と第2成分[4]個の場合は、
(ab・c, ab・defg), (ab・d, ab・defg), ・・・, (ab・g, ab・cdef)
(ac・b, ac・defg), (ac・d, ac・defg), ・・・, (ac・g, ac・bdef)
・・・・・・
(fg・a, fg・bcde), (fg・b, fg・acde), ・・・, (fg・e, fg・abcd)
たては共通の因数の選び方なので7C2個、
横は残りの第1成分の因数の選び方なので5C1個あります。
だから、7C2×5C1個あります。
・共通素因数が2個で、素の部分の素因数の個数が、第1成分[2]個と第2成分[3]個の場合は、
(ab・cd, ab・efg), (ab・ce, ab・cfg), ・・・, (ab・fg, ab・cde)
(ac・bd, ac・efg), (ac・be, ac・bfg), ・・・, (ac・fg, ac・bde)
・・・・・・
(f g・ab, fg・cde), (fg・ac, fg・bde), ・・・, (fg・de, fg・abc)
たては共通の因数の選び方なので7C2個、
横は残りの第1成分の因数の選び方なので5C2個あります。
だから、7C2×5C2個あります。
よって、[3]b)の2つの値の合計は、
7C2×5C1+7C2×5C2=7C2×(5C4+5C2)=7C2×(24−1)
c)
・共通素因数が3個で、素の部分の素因数の個数が、第1成分[1]個と第2成分[3]個の場合は、
(abc・d, abc・efg), (abc・e, abc・dfg), ・・・, (abc・g, abc・def)
(abd・c, abd・efg), (abd・e, abd・cfg), ・・・, (abd・g, abd・cef)
・・・・・・
(efg・a, efg・bcd), (efg・b, efg・acd), ・・・, (efg・d, efg・abc)
たては共通の因数の選び方なので7C3個、
横は残りの第1成分の因数の選び方なので4C1個あります。
だから、7C3×4C1個あります。
・共通素因数が3個で、素の部分の素因数の個数が、第1成分[2]個と第2成分[2]個の場合は、
(abc・de, abc・fg), (abc・df, abc・eg), ・・・, (abc・fg, abc・de)
(abd・ce, abd・fg), (abd・cf, abd・eg), ・・・, (abd・fg, abd・ce)
・・・・・・
(efg・ab, efg・cd), (efg・ac, efg・bd), ・・・, (efg・cd, efg・ab)
たては共通の因数の選び方なので7C3個、
横は残りの第1成分の因数の選び方なので4C2/2個あります。
だから、7C3×4C2/2個あります。
よって、[3]c)の2つの値の合計は、
7C3×4C1+7C3×4C2/2=7C3×(4C1+4C2/2)=7C3×(23−1)=7×7C3
d)
・共通素因数が4個で、素の部分の素因数の個数が、第1成分[1]個と第2成分[2]個の場合は、
(abcd・e, abcd・fg), (abcd・f, abcd・eg), (abcd・g, abcd・ef)
(abce・d, abce・fg), (abce・f, abce・dg), (abce・g, abce・df)
・・・・・・
(defg・a, defg・bc), (defg・b, defg・ac), (defg・c, defg・ab)
たては共通の因数の選び方なので7C4個、
横は残りの第1成分の因数の選び方なので3C1個あります。
だから、7C4×3C1個あります。
7C4×3C1=7C4×(22−1)=3×7C4
e)
・共通素因数が5個で、素の部分の素因数の個数が、第1成分[1]個と第2成分[1]個の場合は、
(abcde・f, abcde・g)
たては共通の因数の選び方なので7C5個、
横は残りの第1成分の因数の選び方なので2C1/2個あります。
だから、7C5×/2個あります。
7C5×2C1/2=7C5×(21−1)=1×7C5
よって、(3)[1]から[3]e)までを加えると、
(27−1)+7C0(26−1)+7C1(25−1)+7C2(24−1)+7C3(23−1)+7C4(22−1)+7C5(21−1)
=27+(7C026+7C125+7C224+7C323+7C422+7C521)−(7C0+7C1+7C2+7C3+7C4+7C5)
=27+(7C026+7C125+7C224+7C323+7C422+7C521)−(7C0+7C1+7C2+7C3+7C4+7C5+7C6+7C7)+7C6+7C7
=27+(7C027+7C126+7C225+7C324+7C423+7C522+7C621+7C720)/2−27−1/2
=(7C027+7C126+7C225+7C324+7C423+7C522+7C621+7C720)/2−1/2
=(2+1)7/2−1/2
=(37−1)/2
です。
つまり、[n] 個の素因数からなるAについて、正の整数の組をN([n])で表すと、
N([n])=(3n−1)/2
また、
(3n+1−1)/2
=3×(3n−3/3)/2
=3×(3n−1/3)/2
=3×(3n−1+2/3)/2
=3×(3n−1)/2+3×(2/3)/2
=3×(3n−1)/2+1
より
N([n+1])=3N([n])+1
です。
以上から
問題1 N(6)=N(2×3)=N([2])=(32−1)/2=4
問題2 N(15)=N(3×5)=N([2])=(32−1)/2=4
問題3 N(35)=N(5×7)=N([2])=(32−1)/2=4
問題4 N(105)=N(3×5×7)=N([3])=3N([2])+1=13
問題5 N(210)=N(2×3×5×7)=N([4])=3N([3])+1=40
問題6 N(2310)=N(2×3×5×7×11)=N([5])=3N([4])+1=121
となります。
次に、一般の場合を考えます。
一般の場合が見通せるように、Aが素数5種類の場合を見てみます。
a,b,c,・・・・・・を素数とします。
ア,イ,ウ,・・・・・・を正の整数とします。
A=aアbイcウdエeオ とします。
1° Aの約数は、(ア+1)(イ+1)(ウ+1)(エ+1)(オ+1) 個あります。
求める整数の組で、(Aの約数, A) の形のものは、
(ア+1)(イ+1)(ウ+1)(エ+1)(オ+1)−1 個あります。
2° 次に整数の組の第2成分がAではなく、第1成分と第2成分の各素数の次数が異なるものを数えます。
ここで、i は0〜ア-1、j は0〜イ-1、k は0〜ウ-1、l は0〜エ-1、m は0〜オ-1のある整数とします。
・(aア・bjckdlem, bイcウdエeオ・ai)の形のものは、
(aア・bjckdlem, bイcウdエeオ・ai)
(bイ・aickdlem, aアcウdエeオ・bj)
(cウ・aibjdem, aアbイdエeオ・ck)
・・・・・・・
合計 5C1×アイウエオ個
・(aアbイ・ckdlem, cウdエeオ・aibj)の形のものは、
(aアbイ・ckdlem, cウdエeオ・aibj)
(aアcウ・bjdlem, bイdエeオ・aick)
(aアdエ・bjckem, bイcウeオ・aidl)
・・・・・・・
合計5C2×アイウエオ 個あります。
2°の合計は、
5C1×アイウエオ+5C2×アイウエオ
=(5C1+5C2)×アイウエオ
=(25/2−1)×アイウエオ
=(24−1)×アイウエオ
3° 次に整数の組の第2成分がAではなく、第1成分と第2成分にaア,bイ,・・・の形の共通の因数があるものを数えます。
ここで、i は0〜ア-1、j は0〜イ-1、k は0〜ウ-1、l は0〜エ-1、m は0〜オ-1のある整数とします。
a) 共通因数が1種類のもの
・(aア・bイ・ckdlem, aア・cウdエeオ・bj)の形のものは、
(aア・bイ・ckdlem, aア・cウdエeオ・bj)
(aア・cイ・bjdlem, aア・bイdエeオ・ck)
・・・・・・・
(bイ・aイ・ckdlem, bイ・cウdエeオ・ai)
(bイ・cウ・aidlem, bイ・aアdエeオ・ck)
・・・・・・・
・・・・・・・
2つめの因数の選び方は4C1なので、
4C1×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
・(aア・bイcウ・dlem, aア・dエeオ・bjck)の形のものは、
(aア・bイcウ・dlem, aア・dエeオ・bjck)
(aア・bイdエ・ckem, aア・cウeオ・bjdl)
・・・・・・・
・・・・・・・
2つめの因数の選び方は4C2/2なので、
4C2/2×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
3°a)の合計は、
(4C1+4C2/2)×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
=(24/2−1)×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
=(23−1)×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
b) 共通因数が2種類のもの
・(aアbイ・cウ・dlem, aアbイ・dエeオ・ck)の形のものは、
(aアbイ・cウ・dlem, aアbイ・dエeオ・ck)
(aアbイ・dエ・ckem, aアbイ・cウeオ・dl)
・・・・・・・
(aアcウ・bイ・dlem, aアcウ・dエeオ・bj)
(aアcウ・dエ・bjem, aアcウ・bイeオ・dl)
・・・・・・・
・・・・・・・
2つめの因数の選び方は3C1なので、
3C1×(ウエオ+イエオ+アエオ+イウオ+アウオ+アイオ+イウエ+アウエ+アイエ+アイウ)
=(23/2−1)×(ウエオ+イエオ+アエオ+イウオ+アウオ+アイオ+イウエ+アウエ+アイエ+アイウ)
=(22−1)×(ウエオ+イエオ+アエオ+イウオ+アウオ+アイオ+イウエ+アウエ+アイエ+アイウ)
c) 共通因数が3種類のもの
・(aアbイcウ・dエ・em, aアbイcウ・eオ・dl)の形のものは、
(aアbイcウ・dエ・em, aアbイcウ・eオ・dl)
(aアbイdエ・cウ・em, aアbイdエ・eオ・ck)
・・・・・・・
・・・・・・・
2つめの因数の選び方は2C1/2なので、
2C1/2×(エオ+ウオ+イオ+アオ+ウエ+イエ+アエ+イウ+アウ+アイ)
=(21−1)×(エオ+ウオ+イオ+アオ+ウエ+イエ+アエ+イウ+アウ+アイ)
よって1°、2°、3°を合わせると、
N(A) | ={(ア+1)(イ+1)(ウ+1)(エ+1)(オ+1)-1}
|
| +(24−1)×アイウエオ
|
| +(23−1)×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
|
| +(22−1)×(ウエオ+イエオ+アエオ+イウオ+アウオ+アイオ+イウエ+アウエ+アイエ+アイウ)
|
| +(21−1)×(エオ+ウオ+イオ+アオ+ウエ+イエ+アエ+イウ+アウ+アイ)
|
| =24×アイウエオ
|
| +23×(イウエオ+アウエオ+アイエオ+アイウオ+アイウエ)
|
| +22×(ウエオ+イエオ+アエオ+イウオ+アウオ+アイオ+イウエ+アウエ+アイエ+アイウ)
|
| +21×(エオ+ウオ+イオ+アオ+ウエ+イエ+アエ+イウ+アウ+アイ)
|
| +20×(オ+エ+ウ+イ+ア)
|
よって、一般には、p1,p2,p3,・・・を素数、m1,m2,m3,・・・を正の整数として、
A=p1m1p2m2p3m3・・・・・・pnmn のとき
N(A) | =((各指数+1)のn次の基本対称式−1)
|
| +(2n-1−1)×(各指数のn次の基本対称式)
|
| +(2n-2−1)×(各指数のn-1次の基本対称式)
|
| +(2n-3−1)×(各指数のn-2次の基本対称式)
|
| +・・・・・・
|
| +(21−1)×(各指数の2次の基本対称式)
|
| =2n-1×(各指数のn次の基本対称式)
|
| +2n-2×(各指数のn-1次の基本対称式)
|
| +2n-3×(各指数のn-2次の基本対称式)
|
| +・・・・・・
|
| +21×(各指数の2次の基本対称式)
|
| +20×(各指数の1次の基本対称式)
|
NO.1741 楕円積分に関する” 3種の勘違い” 標準形 2008.10.1. 本多欣亮
地図投影のプログラムを作成中に楕円積分の計算に接する機会があり,多くの教科書やサイトに(私を含
む)初心者がつまずきやすい(不親切・不明瞭な)記載があると感じましたので投稿しました。私はこの分
野に関しても素人同然ですので,専門の方のフォローがあれば幸いです。
楕円積分の” 第1種勘違い” 標準形 → 【ねぇ,どこに3〜4次式が入ってるの?】
楕円積分の定義を正しく読まないと「う〜む,標準形って根号のなかのどこが3〜4次式なの?」と困惑
してしまいます。楕円積分の定義に合った積分⇒ Legendre-Jacobi の3つの標準形に直せる(これらが
三角関数を含んだ形で表現できる)という理屈です。見かけが楕円積分の定義に合っていても初等的に求め
られる積分もありますし(擬楕円積分,この場合上記標準形に直さずとも初等的に解ける),また,楕円積
分の標準形はLegendre-Jacobi 流の1セットだけではありません。
楕円積分の” 第2種勘違い” 標準形 → 【楕円の弧長が求まるんですよね?】
楕円4分の1周(または全体)の弧長が,第2種完全楕円積分の値そのもの(またはその定数倍)として
求まる解説がなされていることが多いですが,それより短い部分の弧長を求める式をあまり見かけません。
実用上・教育上どちらの観点からも,後者の提示を希望したいものです。
x2/a2 + y2/b2 = 1(a>b)で,点(a,0)
からスタートし,周上の各点における曲率半径を描いた線がx 軸となす角がφとなる点までの周長s( φ ) は:
ここで,k が楕円の離心率,E( φ ,k) が第2種楕円積分にあたる。確かに,φ = π /2 まで行くと[ ]内
の第2項が0 となり,見かけ上は(完全)楕円積分だけになります。
楕円積分の” 第3種勘違い” 標準形 → 【第3種超難しい!仲間はずれじゃん?】
第1種と第2種は級数展開すると互いにとても似た形を持っていますが,第3種は古典的な教科書では
楕円関数を用いた大変面倒な計算が提示されています。ところが,Carlson という研究者が1970 年代
に開発した手法により,第1種〜第3種を見通し良く簡潔に計算できるようになりました。現在の実用的
な数値計算プログラムでは,このやりかたの方が従来のものよりも主流になっています。興味ある方は,
"Numerical Recipes" のサイト等をご覧ください。
参考C プログラム
古典的(級数展開)な方法による第2種楕円積分elle( φ ,k2)。第2引数にはk でなくk2 を入れて使う。
計算機や引数の値にもよると思うが,有効数字10 〜14 桁程度は正しく出る。elle() 内のa[] の各値およ
びreturn の行の符号1か所を換えるだけで第1種楕円積分ellf( φ ,k2) となる。epoly(x,c[],d) は,実係
数c[] を持つ1変数d 次多項式の値をx において評価する。c[0] に定数項,c[d] に最高次の係数を格納(つ
まり,呼び出し側の配列の大きさはd+1 以上必要)。なお,引数のエラーチェックは一切ないので注意願い
ます。
double elle(double phi,double ks) {
double a[16],b[16],c[16] ;
double sp,cp,pw = 1. ;
int m = 7,i ;
c[ 0] = 1. ;
c[ 1] = 6.66666666666666666666666666666666666666666666667e-1 ;
c[ 2] = 5.33333333333333333333333333333333333333333333333e-1 ;
c[ 3] = 4.571428571428571428571428571428571428571428571429e-1 ;
c[ 4] = 4.063492063492063492063492063492063492063492063492e-1 ;
c[ 5] = 3.694083694083694083694083694083694083694083694084e-1 ;
c[ 6] = 3.409923409923409923409923409923409923409923409923e-1 ;
c[ 7] = 3.182595182595182595182595182595182595182595182595e-1 ;
c[ 8] = 2.995383701266054207230677818913113030760089583619e-1 ;
c[ 9] = 2.837731927515209248955378986338738660720084868692e-1 ;
c[10] = 2.702601835728770713290837129846417772114366541611e-1 ;
c[11] = 2.585097408088389377930365950287877868978959300672e-1 ;
c[12] = 2.4816935117648538028131513122763627542198009286446e-1 ;
c[13] = 2.3897789372550444027089605229327937633227712646208e-1 ;
c[14] = 2.307372767004870457787961884210973288725434324461e-1 ;
c[15] = 2.232941387424068184956092146010619311669775152705e-1 ;
a[ 0] = 0. ;
a[ 1] = 2.5e-1 ;
a[ 2] = 4.6875e-2 ;
a[ 3] = 1.953125e-2 ;
a[ 4] = 1.068115234375e-2 ;
a[ 5] = 6.7291259765625e-3 ;
a[ 6] = 4.62627410888671875e-3 ;
a[ 7] = 3.3752918243408203125e-3 ;
a[ 8] = 2.571023069322109222412109375e-3 ;
a[ 9] = 2.02349037863314151763916015625e-3 ;
a[10] = 1.633968480746261775493621826171875e-3 ;
a[11] = 1.34701120623503811657428741455078125e-3 ;
a[12] = 1.12952502189500592066906392574310302734375e-3 ;
a[13] = 9.607646266118763378472067415714263916015625e-4 ;
a[14] = 8.2718893235078638781487825326621532440185546875e-4 ;
a[15] = 7.19654371145184157398944080341607332229614257812e-4 ;
b[0] = epoly(ks,a,m-1) ;
for(i=1; i<m; i++)
b[i] = b[i-1] - a[i]*(pw*=ks) ;
for(i=1; i<m; i++)
b[i] *= c[i] ;
sp = sin(phi) ;
cp = cos(phi) ;
return(phi*(1. - b[0]) + sp*cp*epoly(sp*sp,b,m-1)) ;
}
double epoly(double x,double *c,int d) {
double e = 0. ;
while(d > 0)
e = (e + c[d--])*x ;
return(e+c[0]) ;
}