コンテンツにスキップ

ファイル:Symplectic-method-for-keplerian-motion.svg

ページのコンテンツが他言語でサポートされていません。

元のファイル(SVG ファイル、576 × 432 ピクセル、ファイルサイズ: 71キロバイト)

概要

解説
English: Energy error in numerical solutions of the elliptic Keplerian orbit for e=0.5. The object is at the apocenter at initial time. The time step is 1/128 of the orbital period. The unit of the horizontal axis is the orbital period.
日付
原典 投稿者自身による著作物
作者 Osanshouo

ライセンス

この作品の著作権者である私は、この作品を以下のライセンスで提供します。
w:ja:クリエイティブ・コモンズ
表示
このファイルはクリエイティブ・コモンズ 表示 4.0 国際ライセンスのもとに利用を許諾されています。
あなたは以下の条件に従う場合に限り、自由に
  • 共有 – 本作品を複製、頒布、展示、実演できます。
  • 再構成 – 二次的著作物を作成できます。
あなたの従うべき条件は以下の通りです。
  • 表示 – あなたは適切なクレジットを表示し、ライセンスへのリンクを提供し、変更があったらその旨を示さなければなりません。これらは合理的であればどのような方法で行っても構いませんが、許諾者があなたやあなたの利用行為を支持していると示唆するような方法は除きます。

Source code

'''
Solve the Kepler problem $H = p^2/2 - 4 \pi / | x |$.

Code unit:
   coordinate: initial semi-major axis $a$
   time: orbital period $P$

Output file:
    "./symplectic-method-for-keplerian-motion.svg"
'''

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["text.usetex"] = True
plt.rcParams["font.family"] = "DejaVu Sans"
plt.rcParams["font.size"] = 14

fig = plt.figure()
plt.subplots_adjust(left=0.15, right=0.75)
ax = fig.add_subplot(111)

### parameters ####################
e = 0.5      # eccentricity
dt = 1./128. # time step
t_end = 4    # calculation end time
###################################

def force(x):
    r = np.linalg.norm(x)
    return -4.*np.pi**2*x/r**3

def solve(method, param):
    ts, xs, vs = [0.], [np.array([1.+e, 0])], [np.array([0., 2.*np.pi*np.sqrt((1-e)/(1+e))])]
    while ts[-1] < t_end:
        x, v = method(xs[-1], vs[-1])
        ts.append(ts[-1] + dt)
        xs.append(x)
        vs.append(v)
    ts, xs, vs = np.array(ts), np.array(xs), np.array(vs)
    energy = np.sum(vs**2, axis=1)/2. - 4.*np.pi**2/np.linalg.norm(xs, axis=1)
    ax.plot(ts, np.fabs((energy - energy[0])/energy[0]), **param)
    return 0

# Euler
def euler(x, v):
    return x + v*dt, v + force(x)*dt
solve(euler, {"label":"Euler"})

# RK4
def rk4(x, v):
    dx1, dv1 = v,           force(x)
    dx2, dv2 = v+dv1*dt/2., force(x+dx1*dt/2.)
    dx3, dv3 = v+dv2*dt/2., force(x+dx2*dt/2.)
    dx4, dv4 = v+dv3*dt,    force(x+dx3*dt   )
    return x + (dx1 + 2.*(dx2 + dx3) + dx4)*dt/6., v + (dv1 + 2.*(dv2 + dv3) + dv4)*dt/6., 
solve(rk4, {"label": "RK4"})

# symp1
def symp1(x, v):
    x = x + v*dt
    return x, v + force(x)*dt
solve(symp1, {"label":"Symp1", "ls":"-."})

# symp2
def symp2(x, v, dt=dt):
    x = x + v*dt/2.
    v = v + force(x)*dt
    return x + v*dt/2., v
solve(symp2, {"label":"Symp2", "ls":":"})

# symp4
def symp4(x, v):
    d1 = 1./(2. - np.cbrt(2.))
    d2 = 1. - 2.*d1
    x, v = symp2(x, v, dt=dt*d1)
    x, v = symp2(x, v, dt=dt*d2)
    return symp2(x, v, dt=dt*d1)
solve(symp4, {"label":"Symp4", "ls":"-."})

ax.set_title("$e={}$".format(e))
ax.set_xlabel(r"time $t/P$")
ax.set_ylabel(r"energy error $\left| \Delta E / E_0 \right|$")
ax.set_yscale("log")
ax.set_xlim([0, t_end])
ax.set_ylim([1e-7, 2])
ax.grid()
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=14)

plt.savefig("symplectic-method-for-keplerian-motion.svg")
plt.show()
plt.close()

キャプション

このファイルの内容を1行で記述してください
Energy error in numerical solutions of the elliptic Keplerian orbit for e=0.5.

このファイルに描写されている項目

題材

2 6 2020

ファイルの履歴

過去の版のファイルを表示するには、その版の日時をクリックしてください。

日付と時刻サムネイル寸法利用者コメント
現在の版2020年6月2日 (火) 13:172020年6月2日 (火) 13:17時点における版のサムネイル576 × 432 (71キロバイト)OsanshouoUploaded own work with UploadWizard

以下のページがこのファイルを使用しています:

メタデータ