情報教育システムにはv3.3.5のProcessingが既にインストールされてい る。この実習ではこれを利用して頂く。
もし自宅のＰＣ等で使用したい場合は、 オリジナルサイトからダ ウンロードして使う事が可能である。但し、オリジナルサイトでは最新 のv3.5.3がある。バージョンが異なると動作や使い方が変わることがあ るので、注意が必要である。皆さんが各自でより進んだ内容を勉強する 事は大変よいことであるが、試験、又はレポートで使用するバージョン は公平性を保つため、情報教育システムと同じバージョンとする。
アルゴリズムを改良する事で計算精度を上げる方法がある。
前回の対面授業の課題４では、以下のプログラムを用いて地球の公転をシュミュレーションしてもらった。
def setup():
size(600, 600)
frameRate(1000)
colorMode(RGB)
background(0,0,0)
stroke(255,0,0)
fill(255,0,0)
ellipse(sx(0),sy(0),5,5)
stroke(0,255,0)
fill(0,255,0)
x1=-1.6e11
y1=-1.6e11
x2=+1.6e11
y2=+1.6e11
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
t=0.0
dt=40000.0
x=1.496e11
y=0.0
vx=0.0
vy=1.5e4
G=6.67384e-11
Ms=1.989e30
i=0
def draw():
global x
global y
global t
global i
global vx
global vy
global G
global Ms
print("vx="+str(vx)+" vy="+str(vy)+" x="+str(x)+" y="+str(y)+" t="+str(t))
point(sx(x), sy(y))
t=t+dt
r3=(x**2+y**2)**1.5
vx=vx-G*Ms*dt*x/r3
vy=vy-G*Ms*dt*y/r3
x=x+vx*dt
y=y+vy*dt
i=i+1
if i>10000:
noLoop()
前回からの変更として、ｄｔを4００００.０、ｖｙの初期値を、１．5ｅ４としている。復習のため、このプログラムを実行してみよ。本日の授業ではコピー＆ペーストを使用してもよい。楕円軌道がだんだんずれていく事が分かる。厳密な解ではこれは起きない。
次に、同じ初期条件と時間ステップで、４次のルンゲクッタ法を用いたプログラムを実行してみる。
以下のプログラムを入力して実行してみよ。
def setup():
size(600, 600)
frameRate(1000)
colorMode(RGB)
background(0,0,0)
stroke(255,0,0)
fill(255,0,0)
ellipse(sx(0),sy(0),5,5)
stroke(0,255,0)
fill(0,255,0)
x1=-1.6e11
y1=-1.6e11
x2=+1.6e11
y2=+1.6e11
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
t=0.0
dt=40000.0
x=1.496e11
y=0.0
vx=0.0
vy=1.5e4
G=6.67384e-11
Ms=1.989e30
i=0
def draw():
global t
global dt
global x
global y
global vx
global vy
global G
global Ms
global i
point(sx(x), sy(y))
print("i="+str(i)+" t="+str(t)+" Vx="+str(vx)+" Vy="+str(vy)+" X="+str(x)+" Y="+str(y))
# for j in range (10000):
t=t+dt
#step1
r3=(x**2+y**2)**1.5
ax=-G*Ms*x/r3
ay=-G*Ms*y/r3
k00=vx*dt
k10=vy*dt
k20=ax*dt
k30=ay*dt
#step2
r3=((x+0.5*k00)**2+(y+0.5*k10)**2)**1.5
ax=-G*Ms*(x+0.5*k00)/r3
ay=-G*Ms*(y+0.5*k10)/r3
k01=(vx+0.5*k20)*dt
k11=(vy+0.5*k30)*dt
k21=ax*dt
k31=ay*dt
#step3
r3=((x+0.5*k01)**2+(y+0.5*k11)**2)**1.5
ax=-G*Ms*(x+0.5*k01)/r3
ay=-G*Ms*(y+0.5*k11)/r3
k02=(vx+0.5*k21)*dt
k12=(vy+0.5*k31)*dt
k22=ax*dt
k32=ay*dt
#step4
r3=((x+k02)**2+(y+k12)**2)**1.5
ax=-G*Ms*(x+k02)/r3
ay=-G*Ms*(y+k12)/r3
k03=(vx+k22)*dt
k13=(vy+k32)*dt
k23=ax*dt
k33=ay*dt
#last
vx=vx+(k20+2*k21+2*k22+k23)/6
vy=vy+(k30+2*k31+2*k32+k33)/6
x=x+(k00+2*k01+2*k02+k03)/6
y=y+(k10+2*k11+2*k12+k13)/6
i=i+1
if i>10000:
noLoop()
ルンゲクッタ法の場合は、見える範囲内では楕円はずれず、計算精度が
上がっている事が分かる。
一方で計算量は大幅に増えているが、画面表示にかかる時間の方 が大きいので、実際の待ち時間は変わらない。
オイラー法でもｄｔを２０００とすれば楕円がずれない精度が出 せる。しかし同じ時間をシミュレートするための計算量は５倍となる。
オイラー法は、計算量が少なく、精度が低い。ルンゲクッタ法は、 計算量が多く、精度が高い。今回の例では非常に簡単な例であったため、 両者の違いは大きくなかったが、ほとんどの科学技術計算において、ル ンゲクッタ法（更にはその改良版）が使われている。
以下のプログラムを入力して実行してみよ。コピー＆ペーストを使用してもよい。
def setup():
size(600, 600)
frameRate(100)
background(0,0,0)
stroke(255)
colorMode(HSB,360,100,100)
x1=+6.784e11
y1=-1.e11
x2=+8.784e11
y2=+1.e11
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
G=6.674e-11
pi=3.14159265358979323846
######## Parameters to chanege an initial condition of the probe ########
theta=pi/3 #Direction of momentum
disp=1.8e10 #Radial diplacement between Jupiter abd the probe
#########################################################################
def accel(m,r):
return G*m/(r**2)
t=0.0
dt=10000.0
x_jup=7.784e11*cos(-pi/30)
y_jup=7.784e11*sin(-pi/30)
x_prb=7.784e11-disp
y_prb=-7.0e10
vx_jup=1.306e4*(-sin(-pi/30))
vy_jup=1.306e4*cos(-pi/30)
vx_prb=1.0e4*cos(theta)
vy_prb=1.0e4*sin(theta)
m_sol=1.989e30
m_jup=1.898e27
m_prb=721.9 #voyger 1 & 2
ax_jup=0.0
ay_jup=0.0
ax_prb=0.0
ay_prb=0.0
## the sun stand still ##
i=0
def draw():
global t
global dt
global x_jup
global y_jup
global x_prb
global y_prb
global vx_jup
global vy_jup
global vx_prb
global vy_prb
global m_sol
global m_jup
global m_prb
global ax_jup
global ay_jup
global ax_prb
global ay_prb
global i
print("vx_prb="+str(vx_prb)+" vy_prb="+str(vy_prb)+" x_prb="+str(x_prb)+" y_prb="+str(y_prb)+" t="+str(t))
V_prb=vx_prb**2+vy_prb**2
stroke(0,0,99)
point(sx(x_jup),sy(y_jup))
stroke(180,99,99)
point(sx(x_prb),sy(y_prb))
r_sol_jup=sqrt(x_jup**2+y_jup**2)
r_sol_prb=sqrt(x_prb**2+y_prb**2)
r_jup_prb=sqrt((x_jup-x_prb)**2+(y_jup-y_prb)**2)
ax_jup=-x_jup/r_sol_jup*accel(m_sol,r_sol_jup)
ay_jup=-y_jup/r_sol_jup*accel(m_sol,r_sol_jup)
ax_prb=-x_prb/r_sol_prb*accel(m_sol,r_sol_prb)+(x_jup-x_prb)/r_jup_prb*accel(m_jup,r_jup_prb)
ay_prb=-y_prb/r_sol_prb*accel(m_sol,r_sol_prb)+(y_jup-y_prb)/r_jup_prb*accel(m_jup,r_jup_prb)
vx_jup=vx_jup+ax_jup*dt
vy_jup=vy_jup+ay_jup*dt
vx_prb=vx_prb+ax_prb*dt
vy_prb=vy_prb+ay_prb*dt
x_jup=x_jup+vx_jup*dt
y_jup=y_jup+vy_jup*dt
x_prb=x_prb+vx_prb*dt
y_prb=y_prb+vy_prb*dt
i=i+1
if i>100000:
noLoop()
このプログラムは、スイングバイをシミュレーションしている。スイン
グバイについての資料は、CLE上の「第１３回 追加講義資料」の「Processingでスイングバイ」の章を参照せよ。
前回の対面授業の課題４は ２体問題（２つの質点の運動を扱う）であったが、これは、太陽、木星、探査機の３体問題である。 尚、数学的には３体以上の厳密解はなく、高精度の結果を得るには数値的に解くしか方法がない。また、 ３体ではあるが、太陽は座標原点に固定してある事に留意せよ。
「第１３回 追加講義資料」の「Processing でスイングバイ」を 参照して、初期条件を変え、スイングバイが色々なパターンを取る事を 確認せよ。地球に帰還する軌道に載せたり、惑星軌道に載せたり、次の 惑星に向かったり、太陽系を脱出したりする事ができる。
以下のプログラムを入力して実行してみよ。コピー＆ペーストを使用してもよい。
i=0
t=0.0
dt=15.0
x=[0,0,0]
y=[0,0,0]
vx=[0,0,0]
vy=[0,0,0]
ax=[0,0,0]
ay=[0,0,0]
r=[[0,0,0],[0,0,0],[0,0,0]]
m=[5.0e30, 5.0e30, 5.0e30]
pi=3.14159265358979323846
def setup():
x[0]=0.0
y[0]=5.0e10
x[1]=x[0]*cos(2*pi/3)-y[0]*sin(2*pi/3)
y[1]=x[0]*sin(2*pi/3)+y[0]*cos(2*pi/3)
x[2]=x[0]*cos(-2*pi/3)-y[0]*sin(-2*pi/3)
y[2]=x[0]*sin(-2*pi/3)+y[0]*cos(-2*pi/3)
vx[0]=5.0e4
vy[0]=2.0e4
vx[1]=vx[0]*cos(2*pi/3)-vy[0]*sin(2*pi/3)
vy[1]=vx[0]*sin(2*pi/3)+vy[0]*cos(2*pi/3)
vx[2]=vx[0]*cos(-2*pi/3)-vy[0]*sin(-2*pi/3)
vy[2]=vx[0]*sin(-2*pi/3)+vy[0]*cos(-2*pi/3)
size(600, 600)
print("width="+str(width)+" height="+str(height))
background(0)
stroke(255)
colorMode(HSB ,360,100,100)
frameRate(500)
G=6.674e-11
x1=-1.0e11
y1=-1.0e11
x2=+1.0e11
y2=+1.0e11
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
def accel(m, r):
return G*m/sq(r)
def draw():
global i
global t
global dt
global x
global y
global vx
global vy
global ax
global ay
global r
global m
for j in range(500):
t=t+dt
r[0][1]=sqrt((x[0]-x[1])**2+(y[0]-y[1])**2)
r[1][2]=sqrt((x[1]-x[2])**2+(y[1]-y[2])**2)
r[2][0]=sqrt((x[2]-x[0])**2+(y[2]-y[0])**2)
ax[0]=(x[1]-x[0])/r[0][1]*accel(m[1],r[0][1])+(x[2]-x[0])/r[2][0]*accel(m[2],r[2][0])
ay[0]=(y[1]-y[0])/r[0][1]*accel(m[1],r[0][1])+(y[2]-y[0])/r[2][0]*accel(m[2],r[2][0])
ax[1]=(x[2]-x[1])/r[1][2]*accel(m[2],r[1][2])+(x[0]-x[1])/r[0][1]*accel(m[0],r[0][1])
ay[1]=(y[2]-y[1])/r[1][2]*accel(m[2],r[1][2])+(y[0]-y[1])/r[0][1]*accel(m[0],r[0][1])
ax[2]=(x[0]-x[2])/r[2][0]*accel(m[0],r[2][0])+(x[1]-x[2])/r[1][2]*accel(m[1],r[1][2])
ay[2]=(y[0]-y[2])/r[2][0]*accel(m[0],r[2][0])+(y[1]-y[2])/r[1][2]*accel(m[1],r[1][2])
vx[0]=vx[0]+dt*ax[0]
vy[0]=vy[0]+dt*ay[0]
vx[1]=vx[1]+dt*ax[1]
vy[1]=vy[1]+dt*ay[1]
vx[2]=vx[2]+dt*ax[2]
vy[2]=vy[2]+dt*ay[2]
x[0]=x[0]+dt*vx[0]
y[0]=y[0]+dt*vy[0]
x[1]=x[1]+dt*vx[1]
y[1]=y[1]+dt*vy[1]
x[2]=x[2]+dt*vx[2]
y[2]=y[2]+dt*vy[2]
stroke(0,99,99)
point(sx(x[0]), sy(y[0]))
stroke(120,99,99)
point(sx(x[1]), sy(y[1]))
stroke(240,99,99)
point(sx(x[2]), sy(y[2]))
print("X0="+str(x[0])+" Y0="+str(y[0])+" X1="+str(x[1])+" Y1="+str(y[1])+" X2="+str(x[2])+" Y2="+str(y[2]))
i=i+1
if i >10000:
noLoop()
このプログラムは、連星系のシミュレーションをしている。連星系につ いての資料は、CLE上の「第１３回 追加講義資料」の「連星系のシミュレーショ ン」の章を参照せよ。
課題２では太陽の座標を原点に固定していたが、この計算では３ 体とも固定されておらず、全て計算している。
配列が使用されている。
高精度の計算が必要になってくるので、ｄｔは非常に小さくなっ ている。そのため、５００回の繰り返しを計算して、１回の描画を行う ようになっている。
ニュートンの万有引力の法則下では、２体問題（解くべき物体が ２個）の場合は、数学的な厳密解があり、それが楕円軌道（と彗星のよ うな双曲線軌道）になる事が分かっている。しかし、３体以上になると 数学的な厳密解がない事も分かっている。（解がないという意味ではな い。式に書けないという意味である。）このような場合、数学的に近似 する（つまりモデルを変更する）方法もあるが、それを行ってしまうと 得られた数値の精度がモデルの精度以上にはできないという限界がある。 より高精度な結果を得るには、数学的な近似をせず、計算上の近似を高 めていく手法しかない。それがシミュレーションが必要な理由の１つで ある。
CLE上の「第１３回追加講義資料」の「円周率ってホンマに3.14なん？」 の章を読んでから、以下のプログラムを入力して実行してみよ。コピー ＆ペーストを使用してもよいが、プログラムの中身をしっかりと 読んで、理解する事。
import random
def setup():
size(600,600)
print("width="+str(width)+" height="+str(height))
frameRate(500)
colorMode(RGB)
x1=-1.0
y1=-1.0
x2=+1.0
y2=+1.0
def sx(x):
return width*(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
random.seed()
n=50000
i=0
s=0
def f(x,y):
return x**2+y**2
def draw():
global s
global n
global i
x=random.uniform(-1,1)
y=random.uniform(-1,1)
if f(x,y)<=1:
s=s+1
stroke(0,255,0)
ellipse(sx(x),sy(y),2,2)
else:
stroke(0,0,255)
ellipse(sx(x),sy(y),2,2)
i=i+1
if i>=n:
noLoop()
calc_pi=4.0*s/i
print("pi="+str(calc_pi))
このプログラムは、規則性のない数(乱数)を[-1,1]の区間で２つ(x,y)生成しその数の二乗和、
つまり原点からの距離が円の半径(1)よりも大きいかを判定を行なっている。
この事象判定を複数回(50000)程行う事で、50000回中何回円の内部に点が生成されたか(確率)を計算している。 今回の例では、まず上の確率を計算し、乱数(x,y)の定義域である正方形2×2の面積で規格化する事で円の面積を求めている。 円の半径は1であるので、円の面積が"π"になることを利用することで円周率を見積もっている。
→(乱数の定義域で表される正方形の面積 : 4)/(乱数の生成回数 : 50000)・(円内側に入った回数 : s)
このように、乱数を複数回生成する事で面積を見積もったり、自然界の崩壊現象などを擬似的に再現する事がある。 例えば、放射線検出器に放射線が入射する確率は、放射線検出器と放射線線源の距離や放射線検出器の形などに依存し、手計算で求める時が難しい。 その時、放射線検出器の形を再現し、放射線線源の崩壊を乱数を用いて再現する事で擬似的に放射線が放射線検出器に入射することを再現する事が可能である。
これで今学期の対面授業は全て終了です。お疲れ様でした。来週の期末テストに向けて頑張ってください。第１４回の E-learning は来週の期末テストまでに終わらせてください。