2次のテイラー法を用いて楕円軌道を正しく描きたい

前提・実現したいこと

2物体が互いの引力によって回転する軌道を2次のテイラー法で描きたい

発生している問題・エラーメッセージ

区切りの幅を小さくしているのに、軌道が一致せずきれいな楕円を描けない

試したこと

加速度や加速度の時間微分を物体のx方向、y方向毎に記録した

補足情報（FW/ツールのバージョンなど）

python
1import numpy as np
2import math
3from matplotlib import pyplot as plt
4
5t = 0
6dt = 1/1024
7x1 = 0.0
8y1 = 1.0
9x2 = 0.0
10y2 = -1.0
11v0 = 0.25
12v1x = v0
13v1y = 0
14v2x = -v0
15v2y = 0
16tend = 30
17X1 = [0.0]
18Y1 = [1.0]
19X2 = [0.0]
20Y2 = [-1.0]
21while t < tend:
22    A1x=[]
23    A1y=[]
24    A2x=[]
25    A2y=[]
26    dA1x=[]
27    dA1y=[]
28    dA2x=[]
29    dA2y=[]
30    a1x = -(x1-x2)/(((x1-x2)**2+(y1-y2)**2)**1.5) 
31    a1y = -(y1-y2)/(((x1-x2)**2+(y1-y2)**2)**1.5)
32    a2x = -(x2-x1)/(((x1-x2)**2+(y1-y2)**2)**1.5)
33    a2y = -(y2-y1)/(((x1-x2)**2+(y1-y2)**2)**1.5)
34    da1x = -(v1x-v2x)/(((x1-x2)**2+(y1-y2)**2)**1.5)+3*np.dot(v1x-v2x,x1-x2)*(x1-x2)/(((x1-x2)**2+(y1-y2)**2)**2.5)
35    da1y = -(v1y-v2y)/(((x1-x2)**2+(y1-y2)**2)**1.5)+3*np.dot(v1y-v2y,y1-y2)*(y1-y2)/(((x1-x2)**2+(y1-y2)**2)**2.5)
36    da2x = -(v2x-v1x)/(((x1-x2)**2+(y1-y2)**2)**1.5)+3*np.dot(v2x-v1x,x2-x1)*(x2-x1)/(((x1-x2)**2+(y1-y2)**2)**2.5)
37    da2y = -(v2y-v1y)/(((x1-x2)**2+(y1-y2)**2)**1.5)+3*np.dot(v2y-v1y,y2-y1)*(y2-y1)/(((x1-x2)**2+(y1-y2)**2)**2.5)
38    A1x.append(a1x)
39    A1y.append(a1y)
40    A2x.append(a2x)
41    A2y.append(a2y)
42    dA1x.append(da1x)
43    dA1y.append(da1y)
44    dA2x.append(da2x)
45    dA2y.append(da2y)
46    x1 = x1+v1x*dt+A1x[0]*(dt**2)/2
47    y1 = y1+v1y*dt+A1y[0]*(dt**2)/2
48    x2 = x2+v2x*dt+A2x[0]*(dt**2)/2
49    y2 = y2+v2y*dt+A2y[0]*(dt**2)/2
50    v1x = v1x+A1x[0]*dt+dA1x[0]*(dt**2)/2
51    v1y = v1y+A1y[0]*dt+dA1y[0]*(dt**2)/2
52    v2x = v2x+A2x[0]*dt+dA2x[0]*(dt**2)/2
53    v2y = v2y+A2y[0]*dt+dA2y[0]*(dt**2)/2
54    X1.append(x1)
55    Y1.append(y1)
56    X2.append(x2)
57    Y2.append(y2)
58    t = t+dt
59

aは加速度、r1,r2はそれぞれ1と2の座標、v1,v2はそれぞれ速度
G=1、M=1と仮定

実際に行った結果
本来は楕円軌道になるらしいですが、どこでミスがあったのかがわかりません。
ソースコードに問題があるのではないかと考えていますが、インターネットで検索してもめぼしい情報は出てきませんでした

行動規範の内容に同意します

回答3件

自己解決

こちらのミスでした。加速度における内積の計算方法を間違えていたことによって発生していました。
これが正しいプログラミングです。

python
1import numpy as np
2import math
3from matplotlib import pyplot as plt
4
5t = 0
6dt = 1/1024
7tend = 30
8dt2 = (dt**2)/2
9
10x1 = 0.0
11y1 = 1.0
12x2 = 0.0
13y2 = -1.0
14
15v0 = 0.25
16v1x = v0
17v1y = 0
18v2x = -v0
19v2y = 0
20
21X1 = [0.0]
22Y1 = [1.0]
23X2 = [0.0]
24Y2 = [-1.0]
25
26while t < tend:
27    r3 = ((x1-x2)**2+(y1-y2)**2)**1.5 
28    r5 = ((x1-x2)**2+(y1-y2)**2)**2.5
29
30    a1x = -(x1-x2)/r3
31    a1y = -(y1-y2)/r3
32    a2x = -(x2-x1)/r3
33    a2y = -(y2-y1)/r3
34
35    da1x = -(v1x-v2x)/r3 + 3*((v1x-v2x)*(x1-x2)+(v1y-v2y)*(y1-y2))*(x1-x2)/r5
36    da1y = -(v1y-v2y)/r3 + 3*((v1x-v2x)*(x1-x2)+(v1y-v2y)*(y1-y2))*(y1-y2)/r5
37    da2x = -(v2x-v1x)/r3 + 3*((v1x-v2x)*(x1-x2)+(v1y-v2y)*(y1-y2))*(x2-x1)/r5
38    da2y = -(v2y-v1y)/r3 + 3*((v1x-v2x)*(x1-x2)+(v1y-v2y)*(y1-y2))*(y2-y1)/r5
39
40    x1 += v1x*dt+a1x*dt2
41    y1 += v1y*dt+a1y*dt2
42    x2 += v2x*dt+a2x*dt2
43    y2 += v2y*dt+a2y*dt2
44
45    v1x += a1x*dt+da1x*dt2
46    v1y += a1y*dt+da1y*dt2
47    v2x += a2x*dt+da2x*dt2
48    v2y += a2y*dt+da2y*dt2
49
50    X1.append(x1)
51    Y1.append(y1)
52    X2.append(x2)
53    Y2.append(y2)
54
55    t += dt

投稿2021/07/14 01:08

kbskyotoshinbun

総合スコア1

陽解法を使っている限り、誤差の集積は回避できません。

陰解法を使えばより正確な近似が可能ですが、計算時間は激増します。

"陽解法陰解法" でググってみればいくらでも解説が見つかります。

投稿2021/07/13 23:47

ppaul

総合スコア24670

kbskyotoshinbun

2021/07/14 00:59

検索して原理はわかったのですが、具体的なプログラム構築しようとしてもうまくいきません…

行動規範の内容に同意します

質問のコードは次のコードと同じです。

Python
1import numpy as np
2import math
3from matplotlib import pyplot as plt
4
5t = 0
6dt = 1/1024
7tend = 30
8dt2 = (dt**2)/2
9
10x1 = 0.0
11y1 = 1.0
12x2 = 0.0
13y2 = -1.0
14
15v0 = 0.25
16v1x = v0
17v1y = 0
18v2x = -v0
19v2y = 0
20
21X1 = [0.0]
22Y1 = [1.0]
23X2 = [0.0]
24Y2 = [-1.0]
25
26while t < tend:
27	r3 = ((x1-x2)**2+(y1-y2)**2)**1.5 
28	r5 = ((x1-x2)**2+(y1-y2)**2)**2.5
29
30	a1x = -(x1-x2)/r3
31	a1y = -(y1-y2)/r3
32	a2x = -(x2-x1)/r3
33	a2y = -(y2-y1)/r3
34
35	da1x = -(v1x-v2x)/r3 + 3*np.dot(v1x-v2x,x1-x2)*(x1-x2)/r5
36	da1y = -(v1y-v2y)/r3 + 3*np.dot(v1y-v2y,y1-y2)*(y1-y2)/r5
37	da2x = -(v2x-v1x)/r3 + 3*np.dot(v2x-v1x,x2-x1)*(x2-x1)/r5
38	da2y = -(v2y-v1y)/r3 + 3*np.dot(v2y-v1y,y2-y1)*(y2-y1)/r5
39
40	x1 = x1 + v1x*dt+a1x*dt2
41	y1 = y1 + v1y*dt+a1y*dt2
42	x2 = x2 + v2x*dt+a2x*dt2
43	y2 = y2 + v2y*dt+a2y*dt2
44
45	v1x = v1x + a1x*dt+da1x*dt2
46	v1y = v1y + a1y*dt+da1y*dt2
47	v2x = v2x + a2x*dt+da2x*dt2
48	v2y = v2y + a2y*dt+da2y*dt2
49
50	X1.append(x1)
51	Y1.append(y1)
52	X2.append(x2)
53	Y2.append(y2)
54
55	t += dt

ここで、x1 = x1 + v1x*dt+a1x*dt2 を
x1 += v1x*dt+a1x*dt2 に変えても同じはずなのに結果が異なります。

値を見てみると、x1 + v1x*dt+a1x*dt2 は 0.00244139664568416
そして、x1 += v1x*dt+a1x*dt2 の x1 は 0.0024413966456841595
微妙に異なります。

次のように変更すると、結果が良くなったりしませんか？

Python
1	x1 += v1x*dt+a1x*dt2
2	y1 += v1y*dt+a1y*dt2
3	x2 += v2x*dt+a2x*dt2
4	y2 += v2y*dt+a2y*dt2
5
6	v1x += a1x*dt+da1x*dt2
7	v1y += a1y*dt+da1y*dt2
8	v2x += a2x*dt+da2x*dt2
9	v2y += a2y*dt+da2y*dt2