ケプラー運動

　　ケプラー運動 Frog 2005/09/20 20:09:50 
　　├Frogさんの入手したプログラムがわからない... 青木太一 2005/09/21 06:44:17 
　　└青木さんの仰る通りプログラムを見ないこと... 落合稔彦 2005/09/21 09:54:51 
　　　└すみません。こんなプログラムです。彗星の... Frog 2005/09/21 20:43:16 
　　　　├ケプラー運動について詳しくないのですが 青木太一 2005/09/22 06:54:08 
　　　　└プログラムを拝見しました。御免なさい、元... 落合稔彦 2005/09/22 08:58:38 
　　　　　├落合さんのプログラムにバグがあると思いま... 青木太一 2005/09/22 16:16:43 
　　　　　│├二人とも有り難うございます。色々試してみ... Frog 2005/09/23 06:50:41 
　　　　　│└天文計算では普通、度を単位として計算しま... 落合稔彦 2005/09/23 09:12:53 
　　　　　└収束判定の問題です。 河川屋 2005/09/23 09:14:50 
　　　　　　├河川屋さんから２つ指摘を受けました。どち... 落合稔彦 2005/09/24 08:25:33 
　　　　　　│├初期値を180度にして小数を七桁にしたら、離... Frog 2005/09/24 20:09:02 
　　　　　　│└何か、勘違いしていません？ 河川屋 2005/09/25 11:38:02 
　　　　　　└河川屋さんに間違ったことを書いてしまいま... 落合稔彦 2005/09/24 20:47:13

　　ケプラー運動 Frog 2005/09/20 20:09:50  ツリーへ

ケプラー運動 返事を書く

Frog 2005/09/20 20:09:50

先日ケプラー運動のプログラムを入手したんですが、離心率が大きいと上手くいきません。大きな離心率用の計算法はありますか。

　　├Frogさんの入手したプログラムがわからない... 青木太一 2005/09/21 06:44:17  ツリーへ

Re: ケプラー運動 返事を書く

青木太一 2005/09/21 06:44:17

Frogさんの入手したプログラムがわからないことには、それよりよいプログラムを紹介できないと思うのですが……

一般的に、用いられている数値計算法をより精度のよいものに変更すればいいのではないでしょうか。

　　└青木さんの仰る通りプログラムを見ないこと... 落合稔彦 2005/09/21 09:54:51  ツリーへ

Re: ケプラー運動 返事を書く

落合稔彦 2005/09/21 09:54:51

　青木さんの仰る通りプログラムを見ないことには確かなことが言えないのですが、私の思うには計算式自体が e≒0 の時の為の略算式なのではないでしょうか。もし略算式でないとすれば、
Kepler 方程式を解く際に第一近似からの近似の上げ方が不充分なのではないかと思われます(Newton 近似が収束する前にやめてしまう、等)。
　自作の日月惑星の位置計算プログラムがありますので差し上げても宜しいのですが、さて方法はどうしましょうか。お持ちのウェブサイトがメイルを受け附けるようになっているのでしたら、URL を(ここに)公開して下さったらお送り出来ます。

　　　└すみません。こんなプログラムです。彗星の... Frog 2005/09/21 20:43:16  ツリーへ

Re: 青木さんの仰る通りプログラムを見ないこと... 返事を書く

Frog 2005/09/21 20:43:16

すみません。こんなプログラムです。彗星の軌道を計算しようとしたら途中で止まります。

OPTION ANGLE DEGREES
PRINT "平均近点離隔(M)"
INPUT M
PRINT "離心率(e)"
INPUT EC
PRINT "機動半径(AU)"
INPUT A
2 LET EI=M+EC*SIN(M)*180/PI
1 LET EII=M+EC*SIN(EI)*180/PI
IF ABS(EII-EI)>0.000000000001 THEN
LET EI=EII
GOTO 1
END IF
LET V=ACOS((COS(EII)-EC)/(1-EC*COS(EII)))
LET R=A*(1-EC*COS(EII))
PRINT "離心近点離隔(v)";V
PRINT "距離(r)";R
LET M=M+1
GOTO 2
END

　　　　├ケプラー運動について詳しくないのですが 青木太一 2005/09/22 06:54:08  ツリーへ

Re: すみません。こんなプログラムです。彗星の... 返事を書く

青木太一 2005/09/22 06:54:08

ケプラー運動について詳しくないのですが
・平均近点離隔→平均近点離角
・離心近点離隔→離心近点離角
・機動半径→軌道半径
ですよね。プログラムについては考え中、調べ中です。

　　　　└プログラムを拝見しました。御免なさい、元... 落合稔彦 2005/09/22 08:58:38  ツリーへ

Re: すみません。こんなプログラムです。彗星の... 返事を書く

落合稔彦 2005/09/22 08:58:38

　プログラムを拝見しました。御免なさい、元になったプログラムから大分写し違いがあったように思えます。パラグラフ２の LET EI の式と末尾の LET M=M+1, GOTO 2 の式が意味を為しま
せん。
　近点離角θ(あなたのプログラムでの V)を計算するには、Kepler 方程式 M=u-e sin u を解いて離心離角 u を出し、それ
からθを計算する手順になります。
　Kepler 方程式を解くやり方は色々ありますが、最も簡便には Newton 近似で解きます。即ち f(u)=M-u+e sin u の零点の第一
近似を u1 としますと u2=u1-f(u1)/f′(u1) がより良い第二近
似になる、というやり方です。実際に計算すれば
u2=(M+e sin u1-eu cos u1)/(1-e cos u1)
となるので、収束するまでこれを繰り返し計算すれば良いのです。
　あなたのプログラムの８行目 "2 LET EI=" 以下を次のように
してみて下さい。きっと瞬時に答えが出ます。

LET EI=180
1 LET EII=(M+EC*SIN(EI)*180/PI-EI*EC*COS(EI))/(1-EC*COS(EI))
IF ABS(EII-EI)>1E-7 THEN
LET EI=EII
GOTO 1
END IF
LET V=SGN(M)*ACOS((COS(EII)-EC)/(1-EC*COS(EII)))
LET R=A*(1-EC*COS(EII))
PRINT "離心近点離隔(v)";V
PRINT "距離(r)";R
END

　なお、Kepler 方程式の意味その他に就て(正確にでも、遠くから漠然とでも)知っておきたいとお思いでしたら、図書館に行って大学一二年生を対象とした物理または力学の本をお開き下さい。大体のことは書いてあります。

　　　　　├落合さんのプログラムにバグがあると思いま... 青木太一 2005/09/22 16:16:43  ツリーへ

Re: プログラムを拝見しました。御免なさい、元... 返事を書く

青木太一 2005/09/22 16:16:43

落合さんのプログラムにバグがあると思います。
Newton近似の式は
「u2=(M+e sin u1-eu cos u1)/(1-e cos u1)」
であっていると思いますが、コード上の
「1 LET EII=(M+EC*SIN(EI)*180/PI-EI*EC*COS(EI))/(1-EC*COS(EI))」
は、180/PIの係数をつけないといけない項が足りません。(詳細は下記コードに)

また、もともと時間変化率が一定な平均近点離角を変えてみることを
意図したプログラムだったようなので、そのようにし、
グラフィカルに面積速度一定を示すように改造したみました。
(平均近点離角mmを-720度から+720度まで回しています)

この際、落合さんのコードでは真近点離角Vを求める際の
sgnの使い方だけでは、十分な範囲の平均近点離角を受け入れられないので
「 LET m=mod(mm+180,360)-180 」
として平均近点離角を調整しました。

ケプラー運動については昨日からにわかに勉強しただけなので、
(一応力学で習った記憶もかすかにありますが)
なにか大間違いしている可能性も十分あります。
間違っていたらごめんなさい。

---コード

set window -20,20,-20,20
OPTION ANGLE DEGREES
!PRINT "離心率(e)"
!INPUT EC
LET ec=0.8
!PRINT "機動半径(AU)"
!INPUT A
LET a=10

def f(q)=q-ec*sin(q)
for mm=-360*2 to 360*2 step 0.01 !平均近点離角

LET m=mod(mm+180,360)-180 !これを行わないとVを求めるときのsgnが有効に働かない範囲がある
LET EI=180

!1 LET EII=(M+EC*SIN(EI)*180/PI-EI*EC*COS(EI))/(1-EC*COS(EI))
!上の式は「180/PI」の係数を分子のMと分母の1にもつけないといけない。
!実際得られたEIでMを逆算してみると、間違っているのがわかる。
!下が正しい式
!1 LET EII=( M*180/PI + EC*SIN(EI)*180/PI - EI*EC*COS(EI) ) / (180/PI-EC*COS(EI)) !(180/PI)を分子のMと分母の1にもかけてやらないといけない
!しかし、見難いのでニュートン近似の式を通分しないでおいた式
1 LET EII=ei - ( M -(ei + ec*sin(ei))) / (EC*COS(EI)*pi/180-1)

IF ABS(EII-EI)>1E-7 THEN
LET EI=EII
GOTO 1
END IF

!検算
if abs(m-(ei+ec*sin(ei)))>0.1 then
print "TOO BIG ERROR!!"
stop
end if

LET V=SGN(M)*ACOS((COS(EII)-EC)/(1-EC*COS(EII)))
LET R=A*(1-EC*COS(EII))
! PRINT "真近点離隔(v)";V
! PRINT "距離(r)";R

!面積速度一定が見えるはず
set line color mm/(180/3)+.5
plot lines:0,0;r*cos(V),r*sin(V)

next mm
END

　　　　　│├二人とも有り難うございます。色々試してみ... Frog 2005/09/23 06:50:41  ツリーへ

Re: 落合さんのプログラムにバグがあると思いま... 返事を書く

Frog 2005/09/23 06:50:41

二人とも有り難うございます。色々試してみます。

　　　　　│└天文計算では普通、度を単位として計算しま... 落合稔彦 2005/09/23 09:12:53  ツリーへ

Re: 落合さんのプログラムにバグがあると思いま... 返事を書く

落合稔彦 2005/09/23 09:12:53

　天文計算では普通、度を単位として計算します。それなのに式の微分計算ではラジアンを使うのですから紛らわしいですね。
　私の式がその紛れを冒したようだという指摘を受けました。でも、あれで良いのだと思います。
　Kepler の方程式を度を単位として書けば(変数 u が度を単位として動く)、
　M=u-180/π・sin(π/180・u)
となります。ここの sin はラジアンに応ずる函数ですから、度に応ずる函数を区別して仮に Sin と書くことにすれば、
　M=u-180/π・Sin u
となります。そして d/du Sin u=π/180・Cos u ですので、これを代入すれば前回示した式になると思います。

　青木さんは力学で習った記憶が不確かだそうですが、同様こちらも余り心丈夫ではありません。十年程前までポケコンで頻りに天文計算をやっていました。この春パソコンを始めたとき白石先生の十進 BASIC を知って、
早速昔のプログラムを書き換えたのは良いのですが、パソコンは習得することが多くてしかもみんな面白いものですからプログラムはそれきり放置してありました。昔のノートを引っ張り出して眺めて見ても、式計算の感覚が戻って来ないのは困ったものです。
　ですから念の為に比較検証してみました。全部をラジアンでやれば紛れが有りませんので、それで書いたプログラムと度で書いたプログラムとで出て来た結果を較べます。すると下記したようになって、矢張り私の式で良さそうに思われます。
　最後に符号の件ですが、あれは SGN(SIN(M)) と書く積りでいたのを書き損なってしまったのです。同じサインだものですから一つ書いた所で安心したのでしょう。駄目ですね。

! プログラム／ラジアン
1 OPTION ANGLE RADIANS
PRINT "平均離角(M)"
INPUT M
LET M=M*PI/180
PRINT "離心率(e)"
INPUT E
LET P=PI
2 LET Q=(M+E*SIN(P)-E*P*COS(P))/(1-E*COS(P))
IF ABS(Q-P)<1E-7 THEN 3
LET P=Q
GOTO 2
3 LET T=SGN(SIN(M))*180/PI*ACOS((COS(Q)-E)/(1-E*COS(Q)))
PRINT "離心離角=";T
END

! プログラム／度(O)
1 OPTION ANGLE DEGREES
PRINT "平均離角(M)"
INPUT M
PRINT "離心率(e)"
INPUT E
LET P=180
2 LET Q=(M+180/PI*E*SIN(P)-E*P*COS(P))/(1-E*COS(P))
IF ABS(Q-P)<1E-7 THEN 3
LET P=Q
GOTO 2
3 LET T=SGN(SIN(M))*ACOS((COS(Q)-E)/(1-E*COS(Q)))
PRINT "離心離角=";T
END

! プログラム／度(A)
! 次の行以外は (O) と同じ。
2 LET Q=(M*180/PI+E*SIN(P)*180/PI-E*P*COS(P))/(180/PI-E*COS(P))
　

　結果は次のようになります。

プログラム E=0.4, M=28°
／ラジアン T=63.2185525800302
／度(O) 63.2185525800303
／度(A) 41.9667607435834

E=0.4, M=111°
／ラジアン 145.206954356869
／度(O) 145.206954356869
／度(A) 131.849360814699

E=0.9, M=7°
／ラジアン 116.202775492282
／度(O) 116.202775492282
／度(A) 30.3105582449822

E=0.9, M=155°
／ラジアン 176.955888137840
／度(O) 176.955888137840
／度(A) 174.266830192197

　　　　　└収束判定の問題です。 河川屋 2005/09/23 09:14:50  ツリーへ

Re: プログラムを拝見しました。御免なさい、元... 返事を書く

河川屋 2005/09/23 09:14:50

収束判定の問題です。

>　元になったプログラムから大分写し違いがあったように思えます。
> パラグラフ２の LET EI の式と
> 末尾の LET M=M+1, GOTO 2 の式が意味を為しません
いえ、そこは問題無いのです。
M=M+1, GOTO 2は、平均近点離角を変えながらループを回しているため。
（しかし、－１８０＜Ｍ≦１８０以外の範囲を計算する意味は無いと思うが。）
パラグラフ２の LET EI の式は、方程式M=u-e sin uを逐次代入法で解いており、
初期値としてu=Mと置いているため。
反復法は、方程式の性質により発散してしまうこともありますが、一応こういう解き方もアリです。
で、元プログラムが解けなかった理由はココ。

IF ABS(EII-EI)>0.000000000001 THEN

において、0.000000000001という値が小さすぎ（＝計算誤差のほうが大きい）たために、
いつまでも収束判定していないと判断され、永久ループに陥っています。
ですから、0.000000000001の値をもっと大きくすれば答は出ます。
落合さん自身、こっそり、
IF ABS(EII-EI)>1E-7 THEN
と、収束判定値を変えています。そればズルイんでは？
（今回の質問は、「プログラムのどこを直せば良いか」であるため。）

>　プログラムの以下を次のように
>してみて下さい。きっと瞬時に答えが出ます。

LET EI=180
1 LET EII=(M+EC*SIN(EI)*180/PI-EI*EC*COS(EI))/(1-EC*COS(EI))
IF ABS(EII-EI)>1E-7 THEN
LET EI=EII
GOTO 1

それって、たぶん遅くなっています。（収束精度を揃えて比較したとして。）
理由は２つ。
１つめ。
ニュートン法の初期値を、EI=180としている点。
EI=M（元プログラム）としたほうが、解に近い初期値となっているから、
収束回数が少ない。実際、EC≒0であれば、解はEI≒Mであるので、
元プログラムは収束ループをほとんど回らないで解にたどり着く。
２つめ。
ニュートン法は、１回のループの中で、方程式とその微分を計算をしないとならない。
一方、元プログラムは、方程式の計算だけ。
したがって、１回のループ内の計算量は、元プログラムの２倍となっている。
（２倍というのは、方程式が複雑な場合の一般論としての値。）
よって、ループの回数が半分になって、計算時間はほぼ等しいため、
見かけ（収束ループの回数）ほど計算時間は減らない。

　　　　　　├河川屋さんから２つ指摘を受けました。どち... 落合稔彦 2005/09/24 08:25:33  ツリーへ

Re: 収束判定の問題です。 返事を書く

落合稔彦 2005/09/24 08:25:33

　河川屋さんから２つ指摘を受けました。どちらの点に就ても私は別の考えをもっています。
　(1) 収束性の問題。発端となった方の質問は、Kepler 運動の計算をしてみた所離心率の大きい場合に上手く行かない、それは何故か、という内容のものでした。使用された式が e≒0 の場合の略算式だったので、それが原因だろうとお報せしました。e が大きい場合に略算式を使えば、例えそれが収束しても意味のある結果になりません。ですから収束性が問題の要点ではなかったと考えます。従って収束の判定条件を変えたのも、単に有効数字は
12～13 桁程度が適当だろうと思ってそうしたまでで、５桁であろうと 50 桁であろうと私は構わなかったのです。
　(2) 初期条件の問題。昔惑星運動のプログラムを作ったとき、初め初期値を０として計算していました。所が、何彗星でしたか、離心率の大きい天体の軌道を計算したとき不可解な計算結果が何度も出て来たのです。大分調べて理由がやっと分りました。Newton 近似で不適当な初期値から出発した為に、遠くへ飛んでから澄して収束してしまうことがあります。あれが起っていたのです。対策を考えた末、初期値をθ＝180°とすれば大丈夫であることを突き止めました。必ず収束するという数学的証明も与えました。従って私が初期値をθ＝180°としたのは漫然と決めたのではありません。反復回数の増えることを承知の上で、安全の為にそうしたのです。もし、以上のことに疑問をお持ちでしたらHalley 彗星か何かで軌道の追跡をなさることを勧めます。或る時突然懸け離れた値が飛び出して来て首を傾げることと思います。

　　　　　　│├初期値を180度にして小数を七桁にしたら、離... Frog 2005/09/24 20:09:02  ツリーへ

Re: 河川屋さんから２つ指摘を受けました。どち... 返事を書く

Frog 2005/09/24 20:09:02

初期値を180度にして小数を七桁にしたら、離心率が0.99999でちょっと引っ掛かるみたいに遅くなりますが、うまくいきました。
ありがとうございます。

　　　　　　│└何か、勘違いしていません？ 河川屋 2005/09/25 11:38:02  ツリーへ

Re: 河川屋さんから２つ指摘を受けました。どち... 返事を書く

河川屋 2005/09/25 11:38:02

何か、勘違いしていません？

>私が初期値をθ＝180°としたのは漫然と決めたのではありません。
>反復回数の増えることを承知の上で、安全の為にそうしたのです。
私はそういうことは言っていません。
初期値をどう与えても収束することを前提に（かなり怪しい仮定ですが。）、
初期値はなるべくそれらしい値（１８０に決め打ちするより、初期値＝Ｍ。）とするほうが
収束回数が減ると言っています。
　※まあ、逐次代入法は、NEWTON法やセカント法に比べ、
　　収束回数で不利なことが多いばかりか、発散する可能性も大きいのも確かだけど。

>有効数字は12～13 桁程度が適当だろうと思ってそうしたまでで、
>５桁であろうと 50 桁であろうと私は構わなかったのです。
これは誤りです。（ただし、１０進1000桁モードなら誤りにあらず。）
理由：コンピュータ演算における実数は、数学が意味する実数でなく、「有限の有効数字」を伴う概念であるため、
離散化現象が起きてしまうため。
具体的には、10進モードにける有効桁は14桁(EPS関数参照)であるので、
DX=EI-EIIを計算すると、10^(-12)の下は、突然ゼロになってしまう(-180≦EI≦180として。)という事態が発生するため。
したがって、収束誤差を10^-13以下に設定するということと、収束誤差＝０と設定することは等しい。
また、収束過程において、このような「解の１ビットずれ」が
実現できるか、というと、関数次第のところがあるものの、一般的には実現できないと見るしかないのです。
（NEWTON法や逐次代入法などのような反復法の場合であり、２分法(かつ、２進モード)はこの限りにあらず。）
よって、収束判定の閾値は、一般的に、計算機の有効数字より大きい値とする必要があります。
　※あくまで一般的な話で、閾値がゼロでも良い場合もありうる。

オマケ。

>使用された式が e≒0 の場合の略算式だったので、それが原因だろうとお報せしました。
NEWTON法をすすめているのはわかったけど、どこをどう読めば略算式が原因というふうに読めるのかわかりません。
(9/24の記事で初めてわかった。)
>e が大きい場合に略算式を使えば、例えそれが収束しても意味のある結果になりません。
そもそも、略算式など使っていませんが。
第一近似から第二近似を求める際に略算しているという意味？
それなら、程度の違いだけで、NEWTON法だって略算（＝近似）ですが。
もう１つ。元々の質問は、「妙な値に収束する」ではなく、「途中で止まってしまう」です。
これは、元プログラムに適当な値を入れてを動かしてみればわかりますが、
・EIとEIIはほとんど同じ値になっている、でもループを抜けない
という事態が起きていることを指しているんです。
ですから、
・収束性が問題の要点ではなかったと考えます。
と書かれても困ります。
「あなたはそう考えたかもしれないけど、そうじゃないよ」としか言い様がないです。
考えただけで終わりにしないで、それが正しいかどうか確かめてみないと。
※確かに、元プログラムで、変な所に収束するという事態も生じるようです。（かなり特殊条件で。）
　でも、それは頻発するわけでないから、２次的な問題です。

　　　　　　└河川屋さんに間違ったことを書いてしまいま... 落合稔彦 2005/09/24 20:47:13  ツリーへ

Re: 収束判定の問題です。 返事を書く

落合稔彦 2005/09/24 20:47:13

　河川屋さんに間違ったことを書いてしまいました。ちょっと慌てています。
　初期値を180°にすれば収束が保証されると言うのは、別の話と混同してしまったので、全然間違いです。Halley彗星の計算なんかしないで下さい。御免なさい。
　面積速度の式を逆に解いて時間を出せば Kepler 方程式を使わずに位置が出せる。しかしやってみると離心率の大きい場合に奇妙な値に収束することがある。それを無くすには或る値のアークコサインを初期値にすれば良い、必ず収束する。そういう話だったのです。
　前に書いた私の式で初期値を180°にしたのは、単に分母を大きくする為だけだったろうと思います。
　以上、訂正します。

　

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

ケプラー運動	返事を書く
Frog 2005/09/20 20:09:50
先日ケプラー運動のプログラムを入手したんですが、離心率が大きいと上手くいきません。大きな離心率用の計算法はありますか。

Re: ケプラー運動	返事を書く
青木太一 2005/09/21 06:44:17
Frogさんの入手したプログラムがわからないことには、それよりよいプログラムを紹介できないと思うのですが…… 一般的に、用いられている数値計算法をより精度のよいものに変更すればいいのではないでしょうか。

Re: ケプラー運動	返事を書く
落合稔彦 2005/09/21 09:54:51
青木さんの仰る通りプログラムを見ないことには確かなことが言えないのですが、私の思うには計算式自体が e≒0 の時の為の略算式なのではないでしょうか。もし略算式でないとすれば、 Kepler 方程式を解く際に第一近似からの近似の上げ方が不充分なのではないかと思われます(Newton 近似が収束する前にやめてしまう、等)。　自作の日月惑星の位置計算プログラムがありますので差し上げても宜しいのですが、さて方法はどうしましょうか。お持ちのウェブサイトがメイルを受け附けるようになっているのでしたら、URL を(ここに)公開して下さったらお送り出来ます。

Re: 青木さんの仰る通りプログラムを見ないこと...	返事を書く
Frog 2005/09/21 20:43:16
すみません。こんなプログラムです。彗星の軌道を計算しようとしたら途中で止まります。 OPTION ANGLE DEGREES PRINT "平均近点離隔(M)" INPUT M PRINT "離心率(e)" INPUT EC PRINT "機動半径(AU)" INPUT A 2 LET EI=M+ECSIN(M)180/PI 1 LET EII=M+ECSIN(EI)180/PI IF ABS(EII-EI)>0.000000000001 THEN LET EI=EII GOTO 1 END IF LET V=ACOS((COS(EII)-EC)/(1-ECCOS(EII))) LET R=A(1-EC*COS(EII)) PRINT "離心近点離隔(v)";V PRINT "距離(r)";R LET M=M+1 GOTO 2 END

Re: すみません。こんなプログラムです。彗星の...	返事を書く
青木太一 2005/09/22 06:54:08
ケプラー運動について詳しくないのですが・平均近点離隔→平均近点離角・離心近点離隔→離心近点離角・機動半径→軌道半径ですよね。プログラムについては考え中、調べ中です。

Re: すみません。こんなプログラムです。彗星の...	返事を書く
落合稔彦 2005/09/22 08:58:38
プログラムを拝見しました。御免なさい、元になったプログラムから大分写し違いがあったように思えます。パラグラフ２の LET EI の式と末尾の LET M=M+1, GOTO 2 の式が意味を為しません。　近点離角θ(あなたのプログラムでの V)を計算するには、Kepler 方程式 M=u-e sin u を解いて離心離角 u を出し、それからθを計算する手順になります。　Kepler 方程式を解くやり方は色々ありますが、最も簡便には Newton 近似で解きます。即ち f(u)=M-u+e sin u の零点の第一近似を u1 としますと u2=u1-f(u1)/f′(u1) がより良い第二近似になる、というやり方です。実際に計算すれば u2=(M+e sin u1-eu cos u1)/(1-e cos u1) となるので、収束するまでこれを繰り返し計算すれば良いのです。　あなたのプログラムの８行目 "2 LET EI=" 以下を次のようにしてみて下さい。きっと瞬時に答えが出ます。 LET EI=180 1 LET EII=(M+ECSIN(EI)180/PI-EIECCOS(EI))/(1-ECCOS(EI)) IF ABS(EII-EI)>1E-7 THEN LET EI=EII GOTO 1 END IF LET V=SGN(M)ACOS((COS(EII)-EC)/(1-ECCOS(EII))) LET R=A(1-EC*COS(EII)) PRINT "離心近点離隔(v)";V PRINT "距離(r)";R END 　なお、Kepler 方程式の意味その他に就て(正確にでも、遠くから漠然とでも)知っておきたいとお思いでしたら、図書館に行って大学一二年生を対象とした物理または力学の本をお開き下さい。大体のことは書いてあります。

Re: プログラムを拝見しました。御免なさい、元...	返事を書く
青木太一 2005/09/22 16:16:43
落合さんのプログラムにバグがあると思います。 Newton近似の式は「u2=(M+e sin u1-eu cos u1)/(1-e cos u1)」であっていると思いますが、コード上の「1 LET EII=(M+ECSIN(EI)180/PI-EIECCOS(EI))/(1-ECCOS(EI))」は、180/PIの係数をつけないといけない項が足りません。(詳細は下記コードに) また、もともと時間変化率が一定な平均近点離角を変えてみることを意図したプログラムだったようなので、そのようにし、グラフィカルに面積速度一定を示すように改造したみました。 (平均近点離角mmを-720度から+720度まで回しています) この際、落合さんのコードでは真近点離角Vを求める際の sgnの使い方だけでは、十分な範囲の平均近点離角を受け入れられないので「 LET m=mod(mm+180,360)-180 」として平均近点離角を調整しました。ケプラー運動については昨日からにわかに勉強しただけなので、 (一応力学で習った記憶もかすかにありますが) なにか大間違いしている可能性も十分あります。間違っていたらごめんなさい。 ---コード set window -20,20,-20,20 OPTION ANGLE DEGREES !PRINT "離心率(e)" !INPUT EC LET ec=0.8 !PRINT "機動半径(AU)" !INPUT A LET a=10 def f(q)=q-ecsin(q) for mm=-3602 to 3602 step 0.01 !平均近点離角 LET m=mod(mm+180,360)-180 !これを行わないとVを求めるときのsgnが有効に働かない範囲がある LET EI=180 !1 LET EII=(M+ECSIN(EI)180/PI-EIECCOS(EI))/(1-ECCOS(EI)) !上の式は「180/PI」の係数を分子のMと分母の1にもつけないといけない。 !実際得られたEIでMを逆算してみると、間違っているのがわかる。 !下が正しい式 !1 LET EII=( M180/PI + ECSIN(EI)180/PI - EIECCOS(EI) ) / (180/PI-ECCOS(EI)) !(180/PI)を分子のMと分母の1にもかけてやらないといけない !しかし、見難いのでニュートン近似の式を通分しないでおいた式 1 LET EII=ei - ( M -(ei + ecsin(ei))) / (ECCOS(EI)pi/180-1) IF ABS(EII-EI)>1E-7 THEN LET EI=EII GOTO 1 END IF !検算 if abs(m-(ei+ecsin(ei)))>0.1 then print "TOO BIG ERROR!!" stop end if LET V=SGN(M)ACOS((COS(EII)-EC)/(1-ECCOS(EII))) LET R=A(1-ECCOS(EII)) ! PRINT "真近点離隔(v)";V ! PRINT "距離(r)";R !面積速度一定が見えるはず set line color mm/(180/3)+.5 plot lines:0,0;rcos(V),r*sin(V) next mm END

Re: 落合さんのプログラムにバグがあると思いま...	返事を書く
Frog 2005/09/23 06:50:41
二人とも有り難うございます。色々試してみます。

Re: 落合さんのプログラムにバグがあると思いま...	返事を書く
落合稔彦 2005/09/23 09:12:53
天文計算では普通、度を単位として計算します。それなのに式の微分計算ではラジアンを使うのですから紛らわしいですね。　私の式がその紛れを冒したようだという指摘を受けました。でも、あれで良いのだと思います。　Kepler の方程式を度を単位として書けば(変数 u が度を単位として動く)、　M=u-180/π・sin(π/180・u) となります。ここの sin はラジアンに応ずる函数ですから、度に応ずる函数を区別して仮に Sin と書くことにすれば、　M=u-180/π・Sin u となります。そして d/du Sin u=π/180・Cos u ですので、これを代入すれば前回示した式になると思います。　青木さんは力学で習った記憶が不確かだそうですが、同様こちらも余り心丈夫ではありません。十年程前までポケコンで頻りに天文計算をやっていました。この春パソコンを始めたとき白石先生の十進 BASIC を知って、早速昔のプログラムを書き換えたのは良いのですが、パソコンは習得することが多くてしかもみんな面白いものですからプログラムはそれきり放置してありました。昔のノートを引っ張り出して眺めて見ても、式計算の感覚が戻って来ないのは困ったものです。　ですから念の為に比較検証してみました。全部をラジアンでやれば紛れが有りませんので、それで書いたプログラムと度で書いたプログラムとで出て来た結果を較べます。すると下記したようになって、矢張り私の式で良さそうに思われます。　最後に符号の件ですが、あれは SGN(SIN(M)) と書く積りでいたのを書き損なってしまったのです。同じサインだものですから一つ書いた所で安心したのでしょう。駄目ですね。 ! プログラム／ラジアン 1 OPTION ANGLE RADIANS PRINT "平均離角(M)" INPUT M LET M=MPI/180 PRINT "離心率(e)" INPUT E LET P=PI 2 LET Q=(M+ESIN(P)-EPCOS(P))/(1-ECOS(P)) IF ABS(Q-P)<1E-7 THEN 3 LET P=Q GOTO 2 3 LET T=SGN(SIN(M))180/PIACOS((COS(Q)-E)/(1-ECOS(Q))) PRINT "離心離角=";T END ! プログラム／度(O) 1 OPTION ANGLE DEGREES PRINT "平均離角(M)" INPUT M PRINT "離心率(e)" INPUT E LET P=180 2 LET Q=(M+180/PIESIN(P)-EPCOS(P))/(1-ECOS(P)) IF ABS(Q-P)<1E-7 THEN 3 LET P=Q GOTO 2 3 LET T=SGN(SIN(M))ACOS((COS(Q)-E)/(1-ECOS(Q))) PRINT "離心離角=";T END ! プログラム／度(A) ! 次の行以外は (O) と同じ。 2 LET Q=(M180/PI+ESIN(P)180/PI-EPCOS(P))/(180/PI-E*COS(P)) 　　結果は次のようになります。プログラム E=0.4, M=28° ／ラジアン T=63.2185525800302 ／度(O) 63.2185525800303 ／度(A) 41.9667607435834 E=0.4, M=111° ／ラジアン 145.206954356869 ／度(O) 145.206954356869 ／度(A) 131.849360814699 E=0.9, M=7° ／ラジアン 116.202775492282 ／度(O) 116.202775492282 ／度(A) 30.3105582449822 E=0.9, M=155° ／ラジアン 176.955888137840 ／度(O) 176.955888137840 ／度(A) 174.266830192197

Re: プログラムを拝見しました。御免なさい、元...	返事を書く
河川屋 2005/09/23 09:14:50
収束判定の問題です。 >　元になったプログラムから大分写し違いがあったように思えます。 > パラグラフ２の LET EI の式と > 末尾の LET M=M+1, GOTO 2 の式が意味を為しませんいえ、そこは問題無いのです。 M=M+1, GOTO 2は、平均近点離角を変えながらループを回しているため。（しかし、－１８０＜Ｍ≦１８０以外の範囲を計算する意味は無いと思うが。）パラグラフ２の LET EI の式は、方程式M=u-e sin uを逐次代入法で解いており、初期値としてu=Mと置いているため。反復法は、方程式の性質により発散してしまうこともありますが、一応こういう解き方もアリです。で、元プログラムが解けなかった理由はココ。 IF ABS(EII-EI)>0.000000000001 THEN において、0.000000000001という値が小さすぎ（＝計算誤差のほうが大きい）たために、いつまでも収束判定していないと判断され、永久ループに陥っています。ですから、0.000000000001の値をもっと大きくすれば答は出ます。落合さん自身、こっそり、 IF ABS(EII-EI)>1E-7 THEN と、収束判定値を変えています。そればズルイんでは？（今回の質問は、「プログラムのどこを直せば良いか」であるため。） >　プログラムの以下を次のように >してみて下さい。きっと瞬時に答えが出ます。 LET EI=180 1 LET EII=(M+ECSIN(EI)180/PI-EIECCOS(EI))/(1-EC*COS(EI)) IF ABS(EII-EI)>1E-7 THEN LET EI=EII GOTO 1 それって、たぶん遅くなっています。（収束精度を揃えて比較したとして。）理由は２つ。１つめ。ニュートン法の初期値を、EI=180としている点。 EI=M（元プログラム）としたほうが、解に近い初期値となっているから、収束回数が少ない。実際、EC≒0であれば、解はEI≒Mであるので、元プログラムは収束ループをほとんど回らないで解にたどり着く。２つめ。ニュートン法は、１回のループの中で、方程式とその微分を計算をしないとならない。一方、元プログラムは、方程式の計算だけ。したがって、１回のループ内の計算量は、元プログラムの２倍となっている。（２倍というのは、方程式が複雑な場合の一般論としての値。）よって、ループの回数が半分になって、計算時間はほぼ等しいため、見かけ（収束ループの回数）ほど計算時間は減らない。

Re: 収束判定の問題です。	返事を書く
落合稔彦 2005/09/24 08:25:33
河川屋さんから２つ指摘を受けました。どちらの点に就ても私は別の考えをもっています。　(1) 収束性の問題。発端となった方の質問は、Kepler 運動の計算をしてみた所離心率の大きい場合に上手く行かない、それは何故か、という内容のものでした。使用された式が e≒0 の場合の略算式だったので、それが原因だろうとお報せしました。e が大きい場合に略算式を使えば、例えそれが収束しても意味のある結果になりません。ですから収束性が問題の要点ではなかったと考えます。従って収束の判定条件を変えたのも、単に有効数字は 12～13 桁程度が適当だろうと思ってそうしたまでで、５桁であろうと 50 桁であろうと私は構わなかったのです。　(2) 初期条件の問題。昔惑星運動のプログラムを作ったとき、初め初期値を０として計算していました。所が、何彗星でしたか、離心率の大きい天体の軌道を計算したとき不可解な計算結果が何度も出て来たのです。大分調べて理由がやっと分りました。Newton 近似で不適当な初期値から出発した為に、遠くへ飛んでから澄して収束してしまうことがあります。あれが起っていたのです。対策を考えた末、初期値をθ＝180°とすれば大丈夫であることを突き止めました。必ず収束するという数学的証明も与えました。従って私が初期値をθ＝180°としたのは漫然と決めたのではありません。反復回数の増えることを承知の上で、安全の為にそうしたのです。もし、以上のことに疑問をお持ちでしたらHalley 彗星か何かで軌道の追跡をなさることを勧めます。或る時突然懸け離れた値が飛び出して来て首を傾げることと思います。

Re: 河川屋さんから２つ指摘を受けました。どち...	返事を書く
Frog 2005/09/24 20:09:02
初期値を180度にして小数を七桁にしたら、離心率が0.99999でちょっと引っ掛かるみたいに遅くなりますが、うまくいきました。ありがとうございます。

Re: 河川屋さんから２つ指摘を受けました。どち...	返事を書く
河川屋 2005/09/25 11:38:02
何か、勘違いしていません？ >私が初期値をθ＝180°としたのは漫然と決めたのではありません。 >反復回数の増えることを承知の上で、安全の為にそうしたのです。私はそういうことは言っていません。初期値をどう与えても収束することを前提に（かなり怪しい仮定ですが。）、初期値はなるべくそれらしい値（１８０に決め打ちするより、初期値＝Ｍ。）とするほうが収束回数が減ると言っています。　※まあ、逐次代入法は、NEWTON法やセカント法に比べ、　　収束回数で不利なことが多いばかりか、発散する可能性も大きいのも確かだけど。 >有効数字は12～13 桁程度が適当だろうと思ってそうしたまでで、 >５桁であろうと 50 桁であろうと私は構わなかったのです。これは誤りです。（ただし、１０進1000桁モードなら誤りにあらず。）理由：コンピュータ演算における実数は、数学が意味する実数でなく、「有限の有効数字」を伴う概念であるため、離散化現象が起きてしまうため。具体的には、10進モードにける有効桁は14桁(EPS関数参照)であるので、 DX=EI-EIIを計算すると、10^(-12)の下は、突然ゼロになってしまう(-180≦EI≦180として。)という事態が発生するため。したがって、収束誤差を10^-13以下に設定するということと、収束誤差＝０と設定することは等しい。また、収束過程において、このような「解の１ビットずれ」が実現できるか、というと、関数次第のところがあるものの、一般的には実現できないと見るしかないのです。（NEWTON法や逐次代入法などのような反復法の場合であり、２分法(かつ、２進モード)はこの限りにあらず。）よって、収束判定の閾値は、一般的に、計算機の有効数字より大きい値とする必要があります。　※あくまで一般的な話で、閾値がゼロでも良い場合もありうる。オマケ。 >使用された式が e≒0 の場合の略算式だったので、それが原因だろうとお報せしました。 NEWTON法をすすめているのはわかったけど、どこをどう読めば略算式が原因というふうに読めるのかわかりません。 (9/24の記事で初めてわかった。) >e が大きい場合に略算式を使えば、例えそれが収束しても意味のある結果になりません。そもそも、略算式など使っていませんが。第一近似から第二近似を求める際に略算しているという意味？それなら、程度の違いだけで、NEWTON法だって略算（＝近似）ですが。もう１つ。元々の質問は、「妙な値に収束する」ではなく、「途中で止まってしまう」です。これは、元プログラムに適当な値を入れてを動かしてみればわかりますが、・EIとEIIはほとんど同じ値になっている、でもループを抜けないという事態が起きていることを指しているんです。ですから、・収束性が問題の要点ではなかったと考えます。と書かれても困ります。「あなたはそう考えたかもしれないけど、そうじゃないよ」としか言い様がないです。考えただけで終わりにしないで、それが正しいかどうか確かめてみないと。 ※確かに、元プログラムで、変な所に収束するという事態も生じるようです。（かなり特殊条件で。）　でも、それは頻発するわけでないから、２次的な問題です。

Re: 収束判定の問題です。	返事を書く
落合稔彦 2005/09/24 20:47:13
河川屋さんに間違ったことを書いてしまいました。ちょっと慌てています。　初期値を180°にすれば収束が保証されると言うのは、別の話と混同してしまったので、全然間違いです。Halley彗星の計算なんかしないで下さい。御免なさい。　面積速度の式を逆に解いて時間を出せば Kepler 方程式を使わずに位置が出せる。しかしやってみると離心率の大きい場合に奇妙な値に収束することがある。それを無くすには或る値のアークコサインを初期値にすれば良い、必ず収束する。そういう話だったのです。　前に書いた私の式で初期値を180°にしたのは、単に分母を大きくする為だけだったろうと思います。　以上、訂正します。