３体問題の実験がうまくいきません

　　３体問題の実験がうまくいきません 会社員 2005/12/18 22:11:10 
　　├３体問題の実験がうまくいきません 青木太一 2005/12/19 07:49:24 
　　│├接近してしまうと、ものすごく大きい力にな... 青木太一 2005/12/19 08:18:52 
　　│└会社員さんのコードをちょっと書き換えまし... 青木太一 2005/12/19 08:31:56 
　　│　├このコードでも86万ぐらいになったところで... 青木太一 2005/12/19 08:57:01 
　　│　└ありがとうございました 会社員 2005/12/19 21:15:58 
　　└エネルギー保存の法則を使って 会社員 2005/12/24 01:56:49 
　　　├前半コード(1/2) 会社員 2005/12/24 02:00:15 
　　　├後半コード(2/2) 会社員 2005/12/24 02:01:31 
　　　└計算量が減りません 青木太一 2005/12/26 21:42:31 
　　　　└三回目以降のdtの値が小さすぎます。 会社員 2005/12/29 15:50:50

　　３体問題の実験がうまくいきません 会社員 2005/12/18 22:11:10  ツリーへ

３体問題の実験がうまくいきません 返事を書く

会社員 2005/12/18 22:11:10

係数がだめなんでしょうか？

!*************３体問題の実験 *****************************

!******画面の準備**********************
LET mado=1000
SET window -mado,mado,-mado,mado
DRAW grid
SET TEXT COLOR 4
PLOT TEXT ,AT -0.8*mado,0.8*mado:"TEST"
!************************************

!■■■初期設定■■■■■■■■■■
LET K=0.1 !場の係数
DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3)
!*****Aの初期位置と速度**********
LET X(1)=0.3*mado
LET Y(1)=-0.1*mado
!*****Bの初期位置と速度**********
LET X(2)=-0.35*mado
LET Y(2)=0.25*mado
!*****Cの初期位置と速度**********
LET X(3)=0.2*mado
LET Y(3)=0.35*mado
!********************************
FOR I=1 TO 3
LET XV(I)=0
LET YV(I)=0
NEXT I
!********************************

!******メインルーチン***************
FOR T=0 TO 1000000
gosub 800 !表示
NEXT T
!******************************************

!****サブルーチン**********************
800
!■場の計算距離の２乗に反比例
FOR I=0 TO 2
LET A=MOD(I,3)+1
LET B=MOD(I+1,3)+1
LET C=MOD(I+2,3)+1
LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)*(X(A)-X(B))/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^0.5 !Ｘ成分
LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)*(X(A)-X(C))/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^0.5+XG(A) !Ｘ成分
LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)*(Y(A)-Y(B))/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^0.5 !Ｙ成分
LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)*(Y(A)-Y(C))/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^0.5+YG(A) !Ｙ成分
NEXT I

!■軌跡をプロット
SET POINT COLOR 3 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 5 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I

!■単位時間後の位置
FOR I=1 TO 3
LET X(I)=XV(I)+X(I)
LET Y(I)=YV(I)+Y(I)
NEXT I

!■単位時間後の速度
FOR I=1 TO 3
LET XV(I)=XV(I)+XG(I)
LET YV(I)=YV(I)+YG(I)
NEXT I

!■移動後の位置をプロット
SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 5 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I
! WAIT DELAY 0.001
RETURN

END

　　├３体問題の実験がうまくいきません 青木太一 2005/12/19 07:49:24  ツリーへ

Re: ３体問題の実験がうまくいきません 返事を書く

青木太一 2005/12/19 07:49:24

>３体問題の実験がうまくいきません
なぜうまくいっていないと思われるのでしょうか？

もしTが6万になったころ1番と3番が接近して、お互いを吹き飛ばすような動きをするのがおかしいと思われるのでしたら、それは正常な動きだと思います。
重力は距離の二乗に反比例するので、接近した二体間ではものすごい力が働くからです。

今回の話とはちょっとズレますが、実際に太陽系ができるとき、小さい星は木星などによって重力的にボンボン吹き飛ばされたいたという話もあります。
http://www.moriyama.com/netscience/Ida_Shigeru/Ida-1.html#03

　　│├接近してしまうと、ものすごく大きい力にな... 青木太一 2005/12/19 08:18:52  ツリーへ

Re: ３体問題の実験がうまくいきません 返事を書く

青木太一 2005/12/19 08:18:52

接近してしまうと、ものすごく大きい力になってしまうため、単位時間あたりの移動量が大きくなって精度が落ちるという問題もあります。
http://jun.artcompsci.org/~makino/papers/ode/node5.html
上記URLによると解決方法として
・重力の式を変えてしまう
・時間刻み幅を変える
などがあるようです。

　　│└会社員さんのコードをちょっと書き換えまし... 青木太一 2005/12/19 08:31:56  ツリーへ

Re: ３体問題の実験がうまくいきません 返事を書く

青木太一 2005/12/19 08:31:56

会社員さんのコードをちょっと書き換えました。
変更点は
・gotoを使わない
・3体以上の問題も解けるように、力の計算のforループを変更
・時間刻み(単位時間)dtを変更できるようにした。
・位置の更新より速度の更新を先にした

最後のは
http://jun.artcompsci.org/~makino/papers/ode/ode-e.html
の1章の最後と2章の最初で言っている一次のシンプレクティック法というやつみたいです。オイラー法とどう違うのか、私にはよくわかりません(笑)。

で、シンプレクティック法にしたら、Tが6万のころ吹き飛ぶ現象はなくなりました。これはシンプレクティック法にして精度があがったおかげなのか、たまたまなのかはわかりません。(シンプレクティック法だって、吹き飛ぶときは吹き飛ぶだろうし)

!******画面の準備**********************
LET mado=1000
SET window -mado,mado,-mado,mado
DRAW grid
SET TEXT COLOR 4
PLOT TEXT ,AT -0.8*mado,0.8*mado:"TEST"
!************************************

!■■■初期設定■■■■■■■■■■
LET K=0.1 !場の係数
LET dt=1!単位時間
DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3)
!*****Aの初期位置と速度**********
LET X(1)=0.3*mado
LET Y(1)=-0.1*mado
!*****Bの初期位置と速度**********
LET X(2)=-0.35*mado
LET Y(2)=0.25*mado
!*****Cの初期位置と速度**********
LET X(3)=0.2*mado
LET Y(3)=0.35*mado
!********************************
FOR I=1 TO 3
LET XV(I)=0
LET YV(I)=0
NEXT I
!********************************

!******メインルーチン***************
FOR T=0 TO 1000000
call main_loop
NEXT T
!******************************************

!****サブルーチン**********************
sub main_loop

!■場の計算距離の２乗に反比例
for i=1 to 3
LET xg(i)=0
LET yg(i)=0
next i
for i=1 to 3
for j=i+1 to 3
LET rr=(x(i)-x(j))^2 + (y(i)-y(j))^2
LET r=sqr(rr)
LET f_ij=-K/rr ! iがjに引かれる引力
LET fx_ij=f_ij*(x(i)-x(j))/r !x成分
LET xg(i)=xg(i)+fx_ij
LET xg(j)=xg(j)-fx_ij
LET fy_ij=f_ij*(y(i)-y(j))/r !y成分
LET yg(i)=yg(i)+fy_ij
LET yg(j)=yg(j)-fy_ij
next j
next i

!■軌跡をプロット
SET POINT COLOR 3 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 5 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I

!■単位時間後の速度
FOR I=1 TO 3
LET XV(I)=XV(I)+XG(I)*dt
LET YV(I)=YV(I)+YG(I)*dt
NEXT I

!■単位時間後の位置
FOR I=1 TO 3
LET X(I)=X(I)+XV(I)*dt
LET Y(I)=Y(I)+YV(I)*dt
NEXT I

!■移動後の位置をプロット
SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 5 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I
! WAIT DELAY 0.001

end sub

END

　　│　├このコードでも86万ぐらいになったところで... 青木太一 2005/12/19 08:57:01  ツリーへ

Re: 会社員さんのコードをちょっと書き換えまし... 返事を書く

青木太一 2005/12/19 08:57:01

このコードでも86万ぐらいになったところで、吹き飛んでしまいました。軌跡の描き方をちょっと変更しましたが、こんな感じです。
http://inet-lab.naist.jp/~taichi-a/3n.png
それでは。

　　│　└ありがとうございました 会社員 2005/12/19 21:15:58  ツリーへ

Re: 会社員さんのコードをちょっと書き換えまし... 返事を書く

会社員 2005/12/19 21:15:58

ありがとうございました

時分割が少なすぎたのかもしれません。
青木太一の案の位置の更新より速度の更新を先にしました。

併せて
・色をつけ
・○表示で物体らしくしました

!******３体問題の実験 **********
!******画面の準備********************
LET mado=1000
SET window -mado,mado,-mado,mado
SET TEXT COLOR 4
!************************************

!■■■初期設定■■■■■■■■■■
LET K=1 !場の係数
LET dt=1/5 !時間の分割数
DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3)
!*****Aの初期位置と速度**********
LET X(1)=0.3*mado
LET Y(1)=-0.1*mado
!*****Bの初期位置と速度**********
LET X(2)=-0.35*mado
LET Y(2)=0.25*mado
!*****Cの初期位置と速度**********
LET X(3)=0.2*mado
LET Y(3)=0.35*mado
!********************************
!*****初期速度*************
FOR I=1 TO 3
LET XV(I)=0
LET YV(I)=0
NEXT I
!**************************

!******メインルーチン***************
FOR T=0 TO 10000000
DRAW grid
PLOT TEXT ,AT -0.8*mado,0.8*mado:"TEST"
GOSUB 800 !表示
NEXT T
!**********************************

stop

!****サブルーチン**********************
800
!■場の計算距離の２乗に反比例
FOR I=0 TO 2
LET A=MOD(I,3)+1
LET B=MOD(I+1,3)+1
LET C=MOD(I+2,3)+1
LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5*(X(A)-X(B)) !Ｘ成分
LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5*(X(A)-X(C))+XG(A) !Ｘ成分
LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5*(Y(A)-Y(B)) !Ｙ成分
LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5*(Y(A)-Y(C))+YG(A) !Ｙ成分
NEXT I

!■軌跡をプロット
FOR I=1 TO 3
SET POINT COLOR I+1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 4 ! 1・2＋3*4○5×
PLOT POINTS:X(I),Y(I)
NEXT I

!■単位時間後の速度
FOR I=1 TO 3
LET XV(I)=XV(I)+XG(I)*dt
LET YV(I)=YV(I)+YG(I)*dt
NEXT I

!■単位時間後の位置
FOR I=1 TO 3
LET X(I)=X(I)+XV(I)*dt
LET Y(I)=Y(I)+YV(I)*dt
NEXT I

!■移動後の位置をプロット
SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 4 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I
RETURN

END

　　└エネルギー保存の法則を使って 会社員 2005/12/24 01:56:49  ツリーへ

Re: ３体問題の実験がうまくいきません 返事を書く

会社員 2005/12/24 01:56:49

エネルギー保存の法則を使って
計算量を減らす目的で
エネルギー変化の大きい場合だけ
精度を10倍100倍1000倍・・に上げて再計算させてみました。

計算量が減りません。

何かよいアイデアありませんか？

　　　├前半コード(1/2) 会社員 2005/12/24 02:00:15  ツリーへ

Re: エネルギー保存の法則を使って 返事を書く

会社員 2005/12/24 02:00:15

前半コード(1/2)

!３体問題の実験
!画面の準備
LET mado=1000
SET window -mado,mado,-mado,mado
SET TEXT COLOR 4
!

!■■■初期設定■■■■■■■■■■
LET K=1 !場の係数
LET dt0=1/5 !時間の分割数
LET err_hi=10^(-5)

DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3)
DIM TMP_X(3),TMP_Y(3),TMP_XV(3),TMP_YV(3)
!*Aの初期位置と速度
LET X(1)=0.3*mado
LET Y(1)=-0.1*mado
!*Bの初期位置と速度
LET X(2)=-0.35*mado
LET Y(2)=0.25*mado
!*Cの初期位置と速度
LET X(3)=0.2*mado
LET Y(3)=0.35*mado
!
!初期速度
FOR I=1 TO 3
LET XV(I)=0
LET YV(I)=0
NEXT I
!
CALL E_SUM !エネルギー総和
PRINT "E_SUM=";E_SUM

!メインルーチン*
FOR T=0 TO 1000000000
DRAW grid
PLOT TEXT ,AT -0.8*mado,0.8*mado:"TEST"
LET dt=dt0
CALL main_loop !表示
NEXT T
!

!サブルーチン
SUB main_loop
CALL ST_SAVE !位置と速度とエネルギーのセーブ TMP_X(I)<=X(I)
CALL PLT_CLOR !軌跡をプロット

LET N=1
DO
CALL DLOAD !データをロード TMP_X(I)=>X(I)
! PRINT "N=";N
FOR H=1 TO N
CALL GRAVTY !場の計算
CALL v_POS !単位時間後の速度・位置
NEXT H
CALL E_SUM !エネルギーの計算
! PRINT "ABS((E_SUM-TMP_E_SUM)/E_SUM)=";ABS((E_SUM-TMP_E_SUM)/E_SUM)
IF ABS((E_SUM-TMP_E_SUM)/E_SUM)<err_hi THEN
! PRINT "○ ";
ELSE
PRINT "×";E_SUM
LET N=10*N
LET dt=dt/N
END IF
LOOP UNTIL ABS((E_SUM-TMP_E_SUM)/E_SUM)<err_hi !誤差以内

CALL PLT_BLK !先端の位置をプロット

!置き換え処理
LET TMP_E_SUM=E_SUM
FOR I=1 TO 3
LET TMP_X(I)=X(I)
LET TMP_Y(I)=Y(I)
LET TMP_XV(I)=XV(I)
LET TMP_YV(I)=YV(I)
NEXT I

END SUB

　　　├後半コード(2/2) 会社員 2005/12/24 02:01:31  ツリーへ

Re: エネルギー保存の法則を使って 返事を書く

会社員 2005/12/24 02:01:31

後半コード(2/2)

SUB DLOAD
!置き換え処理
LET E_SUM=TMP_E_SUM
FOR I=1 TO 3
LET X(I)=TMP_X(I)
LET Y(I)=TMP_Y(I)
LET XV(I)=TMP_XV(I)
LET YV(I)=TMP_YV(I)
NEXT I
END SUB

SUB GRAVTY !距離の２乗に反比例
FOR I=0 TO 2
LET A=MOD(I,3)+1
LET B=MOD(I+1,3)+1
LET C=MOD(I+2,3)+1
LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5*(X(A)-X(B)) !Ｘ成分
LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5*(X(A)-X(C))+XG(A) !Ｘ成分
LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5*(Y(A)-Y(B)) !Ｙ成分
LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5*(Y(A)-Y(C))+YG(A) !Ｙ成分
NEXT I
END SUB

SUB v_POS
!■単位時間後の速度・位置
FOR I=1 TO 3
LET XV(I)=XV(I)+XG(I)*dt
LET YV(I)=YV(I)+YG(I)*dt
NEXT I
FOR I=1 TO 3
LET X(I)=X(I)+XV(I)*dt
LET Y(I)=Y(I)+YV(I)*dt
NEXT I
END SUB

SUB E_SUM !■総運動エネルギーの計算
LET SUM_mvv=0
FOR I=1 TO 3
LET SUM_mvv=SUM_mvv+1/2*1*1*(SQR(XV(I)^2+YV(I)^2))^2
NEXT I
!総位置エネルギー
LET SUM_ichi=0
LET SUM_ichi=-K*1*1/SQR((X(1)-X(2))^2+(Y(1)-Y(2))^2)
LET SUM_ichi=-K*1*1/SQR((X(2)-X(3))^2+(Y(2)-Y(3))^2)+SUM_ichi
LET SUM_ichi=-K*1*1/SQR((X(3)-X(1))^2+(Y(3)-Y(1))^2)+SUM_ichi
LET E_SUM=SUM_mvv+SUM_ichi
! PRINT "合計=";E_SUM;"運動E=";SUM_mvv;"位置E=";SUM_ichi
END SUB

SUB ST_SAVE !計算前の状態保存
LET TMP_E_SUM=E_SUM
FOR I=1 TO 3
LET TMP_X(I)=X(I)
LET TMP_Y(I)=Y(I)
LET TMP_XV(I)=XV(I)
LET TMP_YV(I)=YV(I)
NEXT I
END SUB

SUB PLT_BLK !■先端の位置をプロット
SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 4 ! 1・2＋3*4○5×
FOR I=1 TO 3
PLOT POINTS:X(I),Y(I)
NEXT I
END SUB

SUB PLT_CLOR
!■軌跡をプロット
FOR I=1 TO 3
SET POINT COLOR I+1 ! 白0黒1青2緑3赤4緑5黄6
SET POINT STYLE 4 ! 1・2＋3*4○5×
PLOT POINTS:X(I),Y(I)
NEXT I
END SUB

END

　　　└計算量が減りません 青木太一 2005/12/26 21:42:31  ツリーへ

Re: エネルギー保存の法則を使って 返事を書く

青木太一 2005/12/26 21:42:31

>計算量が減りません
「時間刻みdt0の間の、力学的エネルギーの変化率をerr_hiに抑える」という条件を満たした状態で計算量を抑えたいということでしょうか？

あまりいいアイディアはないのですが、
・細かい最適化を行う
たとえば重力の相互計算ですが、二体間に働く力は同じなのに、会社員さんのコードではわざわざ二体の立場を入れ違えた場合も再計算しています。
三体問題のとき、会社員さんのコードでは6回重力計算がありますが、私のコードなら3回ですみます。
http://freebbs.around.ne.jprdeggo/ufspkc.html#ufspkc
ここを私の用に二重のfor文にするとか。
また、速度エネルギーの計算でSQRかけたあと二乗していあすが、これは無駄な処理ではないでしょうか。
どちらも、全体の計算量のうちどのていどを占めているの調べていないので、効果的かどうかはわかりませんが。
・もっとよい数値計算アルゴリズムに変更する
たとえばルンゲ・クッタ法とか。

ところでエネルギー変化率が大きいときのNとdtの更新に疑問があります。

会社員さんのコードでは、変化率が大きいとき
> LET N=10*N
> LET dt=dt/N
を行いますが、これですと

一回目
N=1
dt=dt0
失敗すると二回目は
N=10
dt=dt0/10
さらに失敗すると三回目が
N=100
dt=dt0/10/100
となってしまいます。
もし、全体的にdt0時間がたつようにFOR H=1 TO Nを回しているのでしたら、これでは三回目以降のdtの値が小さすぎます。

> LET N=10*N
> LET dt=dt/N
を
> LET N=10*N
> LET dt=dt0/N
か
> LET N=10*N
> LET dt=dt/10
のように計算するべきではないでしょうか。

会社員さんのコードの意図を読み違えていたらすいません。

　　　　└三回目以降のdtの値が小さすぎます。 会社員 2005/12/29 15:50:50  ツリーへ

Re: 計算量が減りません 返事を書く

会社員 2005/12/29 15:50:50

>三回目以降のdtの値が小さすぎます。
気づいてませんでした。
LET N=10*N
LET dt=dt0/N
が正解です。
ありがとうございました。

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

３体問題の実験がうまくいきません	返事を書く
会社員 2005/12/18 22:11:10
係数がだめなんでしょうか？ !***********３体問題の実験 ************************* !**画面の準備******************** LET mado=1000 SET window -mado,mado,-mado,mado DRAW grid SET TEXT COLOR 4 PLOT TEXT ,AT -0.8mado,0.8mado:"TEST" !********************************** !■■■初期設定■■■■■■■■■■ LET K=0.1 !場の係数 DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3) !*Aの初期位置と速度******** LET X(1)=0.3mado LET Y(1)=-0.1mado !***Bの初期位置と速度******** LET X(2)=-0.35mado LET Y(2)=0.25mado !***Cの初期位置と速度******** LET X(3)=0.2mado LET Y(3)=0.35mado !****************************** FOR I=1 TO 3 LET XV(I)=0 LET YV(I)=0 NEXT I !**************************** !**メインルーチン*********** FOR T=0 TO 1000000 gosub 800 !表示 NEXT T !************************************** !サブルーチン******************** 800 !■場の計算距離の２乗に反比例 FOR I=0 TO 2 LET A=MOD(I,3)+1 LET B=MOD(I+1,3)+1 LET C=MOD(I+2,3)+1 LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)(X(A)-X(B))/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^0.5 !Ｘ成分 LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)(X(A)-X(C))/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^0.5+XG(A) !Ｘ成分 LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)(Y(A)-Y(B))/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^0.5 !Ｙ成分 LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)(Y(A)-Y(C))/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^0.5+YG(A) !Ｙ成分 NEXT I !■軌跡をプロット SET POINT COLOR 3 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 5 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I !■単位時間後の位置 FOR I=1 TO 3 LET X(I)=XV(I)+X(I) LET Y(I)=YV(I)+Y(I) NEXT I !■単位時間後の速度 FOR I=1 TO 3 LET XV(I)=XV(I)+XG(I) LET YV(I)=YV(I)+YG(I) NEXT I !■移動後の位置をプロット SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 5 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I ! WAIT DELAY 0.001 RETURN END

Re: ３体問題の実験がうまくいきません	返事を書く
青木太一 2005/12/19 07:49:24
>３体問題の実験がうまくいきませんなぜうまくいっていないと思われるのでしょうか？もしTが6万になったころ1番と3番が接近して、お互いを吹き飛ばすような動きをするのがおかしいと思われるのでしたら、それは正常な動きだと思います。重力は距離の二乗に反比例するので、接近した二体間ではものすごい力が働くからです。今回の話とはちょっとズレますが、実際に太陽系ができるとき、小さい星は木星などによって重力的にボンボン吹き飛ばされたいたという話もあります。 http://www.moriyama.com/netscience/Ida_Shigeru/Ida-1.html#03

Re: ３体問題の実験がうまくいきません	返事を書く
青木太一 2005/12/19 08:18:52
接近してしまうと、ものすごく大きい力になってしまうため、単位時間あたりの移動量が大きくなって精度が落ちるという問題もあります。 http://jun.artcompsci.org/~makino/papers/ode/node5.html 上記URLによると解決方法として・重力の式を変えてしまう・時間刻み幅を変えるなどがあるようです。

Re: ３体問題の実験がうまくいきません	返事を書く
青木太一 2005/12/19 08:31:56
会社員さんのコードをちょっと書き換えました。変更点は・gotoを使わない・3体以上の問題も解けるように、力の計算のforループを変更・時間刻み(単位時間)dtを変更できるようにした。・位置の更新より速度の更新を先にした最後のは http://jun.artcompsci.org/~makino/papers/ode/ode-e.html の1章の最後と2章の最初で言っている一次のシンプレクティック法というやつみたいです。オイラー法とどう違うのか、私にはよくわかりません(笑)。で、シンプレクティック法にしたら、Tが6万のころ吹き飛ぶ現象はなくなりました。これはシンプレクティック法にして精度があがったおかげなのか、たまたまなのかはわかりません。(シンプレクティック法だって、吹き飛ぶときは吹き飛ぶだろうし) !****画面の準備******************** LET mado=1000 SET window -mado,mado,-mado,mado DRAW grid SET TEXT COLOR 4 PLOT TEXT ,AT -0.8mado,0.8mado:"TEST" !********************************** !■■■初期設定■■■■■■■■■■ LET K=0.1 !場の係数 LET dt=1!単位時間 DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3) !*Aの初期位置と速度******** LET X(1)=0.3mado LET Y(1)=-0.1mado !***Bの初期位置と速度******** LET X(2)=-0.35mado LET Y(2)=0.25mado !***Cの初期位置と速度******** LET X(3)=0.2mado LET Y(3)=0.35mado !****************************** FOR I=1 TO 3 LET XV(I)=0 LET YV(I)=0 NEXT I !**************************** !**メインルーチン*********** FOR T=0 TO 1000000 call main_loop NEXT T !************************************** !サブルーチン******************** sub main_loop !■場の計算距離の２乗に反比例 for i=1 to 3 LET xg(i)=0 LET yg(i)=0 next i for i=1 to 3 for j=i+1 to 3 LET rr=(x(i)-x(j))^2 + (y(i)-y(j))^2 LET r=sqr(rr) LET f_ij=-K/rr ! iがjに引かれる引力 LET fx_ij=f_ij(x(i)-x(j))/r !x成分 LET xg(i)=xg(i)+fx_ij LET xg(j)=xg(j)-fx_ij LET fy_ij=f_ij(y(i)-y(j))/r !y成分 LET yg(i)=yg(i)+fy_ij LET yg(j)=yg(j)-fy_ij next j next i !■軌跡をプロット SET POINT COLOR 3 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 5 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I !■単位時間後の速度 FOR I=1 TO 3 LET XV(I)=XV(I)+XG(I)dt LET YV(I)=YV(I)+YG(I)dt NEXT I !■単位時間後の位置 FOR I=1 TO 3 LET X(I)=X(I)+XV(I)dt LET Y(I)=Y(I)+YV(I)dt NEXT I !■移動後の位置をプロット SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 5 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I ! WAIT DELAY 0.001 end sub END

Re: 会社員さんのコードをちょっと書き換えまし...	返事を書く
青木太一 2005/12/19 08:57:01
このコードでも86万ぐらいになったところで、吹き飛んでしまいました。軌跡の描き方をちょっと変更しましたが、こんな感じです。 http://inet-lab.naist.jp/~taichi-a/3n.png それでは。

Re: 会社員さんのコードをちょっと書き換えまし...	返事を書く
会社員 2005/12/19 21:15:58
ありがとうございました時分割が少なすぎたのかもしれません。青木太一の案の位置の更新より速度の更新を先にしました。併せて・色をつけ・○表示で物体らしくしました !****３体問題の実験 ****** !**画面の準備**************** LET mado=1000 SET window -mado,mado,-mado,mado SET TEXT COLOR 4 !******************************** !■■■初期設定■■■■■■■■■■ LET K=1 !場の係数 LET dt=1/5 !時間の分割数 DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3) !*Aの初期位置と速度******** LET X(1)=0.3mado LET Y(1)=-0.1mado !***Bの初期位置と速度******** LET X(2)=-0.35mado LET Y(2)=0.25mado !***Cの初期位置と速度******** LET X(3)=0.2mado LET Y(3)=0.35mado !****************************** !*初期速度********* FOR I=1 TO 3 LET XV(I)=0 LET YV(I)=0 NEXT I !********************** !**メインルーチン************* FOR T=0 TO 10000000 DRAW grid PLOT TEXT ,AT -0.8mado,0.8mado:"TEST" GOSUB 800 !表示 NEXT T !******************************** stop !サブルーチン******************** 800 !■場の計算距離の２乗に反比例 FOR I=0 TO 2 LET A=MOD(I,3)+1 LET B=MOD(I+1,3)+1 LET C=MOD(I+2,3)+1 LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5(X(A)-X(B)) !Ｘ成分 LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5(X(A)-X(C))+XG(A) !Ｘ成分 LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5(Y(A)-Y(B)) !Ｙ成分 LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5(Y(A)-Y(C))+YG(A) !Ｙ成分 NEXT I !■軌跡をプロット FOR I=1 TO 3 SET POINT COLOR I+1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 4 ! 1・2＋34○5× PLOT POINTS:X(I),Y(I) NEXT I !■単位時間後の速度 FOR I=1 TO 3 LET XV(I)=XV(I)+XG(I)dt LET YV(I)=YV(I)+YG(I)dt NEXT I !■単位時間後の位置 FOR I=1 TO 3 LET X(I)=X(I)+XV(I)dt LET Y(I)=Y(I)+YV(I)dt NEXT I !■移動後の位置をプロット SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 4 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I RETURN END

Re: ３体問題の実験がうまくいきません	返事を書く
会社員 2005/12/24 01:56:49
エネルギー保存の法則を使って計算量を減らす目的でエネルギー変化の大きい場合だけ精度を10倍100倍1000倍・・に上げて再計算させてみました。計算量が減りません。何かよいアイデアありませんか？

Re: エネルギー保存の法則を使って	返事を書く
会社員 2005/12/24 02:00:15
前半コード(1/2) !３体問題の実験 !画面の準備 LET mado=1000 SET window -mado,mado,-mado,mado SET TEXT COLOR 4 ! !■■■初期設定■■■■■■■■■■ LET K=1 !場の係数 LET dt0=1/5 !時間の分割数 LET err_hi=10^(-5) DIM X(3),Y(3),XV(3),YV(3),XG(3),YG(3) DIM TMP_X(3),TMP_Y(3),TMP_XV(3),TMP_YV(3) !Aの初期位置と速度 LET X(1)=0.3mado LET Y(1)=-0.1mado !Bの初期位置と速度 LET X(2)=-0.35mado LET Y(2)=0.25mado !Cの初期位置と速度 LET X(3)=0.2mado LET Y(3)=0.35mado ! !初期速度 FOR I=1 TO 3 LET XV(I)=0 LET YV(I)=0 NEXT I ! CALL E_SUM !エネルギー総和 PRINT "E_SUM=";E_SUM !メインルーチン FOR T=0 TO 1000000000 DRAW grid PLOT TEXT ,AT -0.8mado,0.8mado:"TEST" LET dt=dt0 CALL main_loop !表示 NEXT T ! !サブルーチン SUB main_loop CALL ST_SAVE !位置と速度とエネルギーのセーブ TMP_X(I)<=X(I) CALL PLT_CLOR !軌跡をプロット LET N=1 DO CALL DLOAD !データをロード TMP_X(I)=>X(I) ! PRINT "N=";N FOR H=1 TO N CALL GRAVTY !場の計算 CALL v_POS !単位時間後の速度・位置 NEXT H CALL E_SUM !エネルギーの計算 ! PRINT "ABS((E_SUM-TMP_E_SUM)/E_SUM)=";ABS((E_SUM-TMP_E_SUM)/E_SUM) IF ABS((E_SUM-TMP_E_SUM)/E_SUM)<err_hi THEN ! PRINT "○ "; ELSE PRINT "×";E_SUM LET N=10*N LET dt=dt/N END IF LOOP UNTIL ABS((E_SUM-TMP_E_SUM)/E_SUM)<err_hi !誤差以内 CALL PLT_BLK !先端の位置をプロット !置き換え処理 LET TMP_E_SUM=E_SUM FOR I=1 TO 3 LET TMP_X(I)=X(I) LET TMP_Y(I)=Y(I) LET TMP_XV(I)=XV(I) LET TMP_YV(I)=YV(I) NEXT I END SUB

Re: エネルギー保存の法則を使って	返事を書く
会社員 2005/12/24 02:01:31
後半コード(2/2) SUB DLOAD !置き換え処理 LET E_SUM=TMP_E_SUM FOR I=1 TO 3 LET X(I)=TMP_X(I) LET Y(I)=TMP_Y(I) LET XV(I)=TMP_XV(I) LET YV(I)=TMP_YV(I) NEXT I END SUB SUB GRAVTY !距離の２乗に反比例 FOR I=0 TO 2 LET A=MOD(I,3)+1 LET B=MOD(I+1,3)+1 LET C=MOD(I+2,3)+1 LET XG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5(X(A)-X(B)) !Ｘ成分 LET XG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5(X(A)-X(C))+XG(A) !Ｘ成分 LET YG(A)=-K/((X(A)-X(B))^2+(Y(A)-Y(B))^2)^1.5(Y(A)-Y(B)) !Ｙ成分 LET YG(A)=-K/((X(A)-X(C))^2+(Y(A)-Y(C))^2)^1.5(Y(A)-Y(C))+YG(A) !Ｙ成分 NEXT I END SUB SUB v_POS !■単位時間後の速度・位置 FOR I=1 TO 3 LET XV(I)=XV(I)+XG(I)dt LET YV(I)=YV(I)+YG(I)dt NEXT I FOR I=1 TO 3 LET X(I)=X(I)+XV(I)dt LET Y(I)=Y(I)+YV(I)dt NEXT I END SUB SUB E_SUM !■総運動エネルギーの計算 LET SUM_mvv=0 FOR I=1 TO 3 LET SUM_mvv=SUM_mvv+1/211(SQR(XV(I)^2+YV(I)^2))^2 NEXT I !総位置エネルギー LET SUM_ichi=0 LET SUM_ichi=-K11/SQR((X(1)-X(2))^2+(Y(1)-Y(2))^2) LET SUM_ichi=-K11/SQR((X(2)-X(3))^2+(Y(2)-Y(3))^2)+SUM_ichi LET SUM_ichi=-K11/SQR((X(3)-X(1))^2+(Y(3)-Y(1))^2)+SUM_ichi LET E_SUM=SUM_mvv+SUM_ichi ! PRINT "合計=";E_SUM;"運動E=";SUM_mvv;"位置E=";SUM_ichi END SUB SUB ST_SAVE !計算前の状態保存 LET TMP_E_SUM=E_SUM FOR I=1 TO 3 LET TMP_X(I)=X(I) LET TMP_Y(I)=Y(I) LET TMP_XV(I)=XV(I) LET TMP_YV(I)=YV(I) NEXT I END SUB SUB PLT_BLK !■先端の位置をプロット SET POINT COLOR 1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 4 ! 1・2＋34○5× FOR I=1 TO 3 PLOT POINTS:X(I),Y(I) NEXT I END SUB SUB PLT_CLOR !■軌跡をプロット FOR I=1 TO 3 SET POINT COLOR I+1 ! 白0黒1青2緑3赤4緑5黄6 SET POINT STYLE 4 ! 1・2＋3*4○5× PLOT POINTS:X(I),Y(I) NEXT I END SUB END

Re: エネルギー保存の法則を使って	返事を書く
青木太一 2005/12/26 21:42:31
>計算量が減りません「時間刻みdt0の間の、力学的エネルギーの変化率をerr_hiに抑える」という条件を満たした状態で計算量を抑えたいということでしょうか？あまりいいアイディアはないのですが、・細かい最適化を行うたとえば重力の相互計算ですが、二体間に働く力は同じなのに、会社員さんのコードではわざわざ二体の立場を入れ違えた場合も再計算しています。三体問題のとき、会社員さんのコードでは6回重力計算がありますが、私のコードなら3回ですみます。 http://freebbs.around.ne.jprdeggo/ufspkc.html#ufspkc ここを私の用に二重のfor文にするとか。また、速度エネルギーの計算でSQRかけたあと二乗していあすが、これは無駄な処理ではないでしょうか。どちらも、全体の計算量のうちどのていどを占めているの調べていないので、効果的かどうかはわかりませんが。・もっとよい数値計算アルゴリズムに変更するたとえばルンゲ・クッタ法とか。ところでエネルギー変化率が大きいときのNとdtの更新に疑問があります。会社員さんのコードでは、変化率が大きいとき > LET N=10N > LET dt=dt/N を行いますが、これですと一回目 N=1 dt=dt0 失敗すると二回目は N=10 dt=dt0/10 さらに失敗すると三回目が N=100 dt=dt0/10/100 となってしまいます。もし、全体的にdt0時間がたつようにFOR H=1 TO Nを回しているのでしたら、これでは三回目以降のdtの値が小さすぎます。 > LET N=10N > LET dt=dt/N を > LET N=10N > LET dt=dt0/N か > LET N=10N > LET dt=dt/10 のように計算するべきではないでしょうか。会社員さんのコードの意図を読み違えていたらすいません。

Re: 計算量が減りません	返事を書く
会社員 2005/12/29 15:50:50
>三回目以降のdtの値が小さすぎます。気づいてませんでした。 LET N=10*N LET dt=dt0/N が正解です。ありがとうございました。