情報教育システムにはv3.3.5のProcessingが既にインストールされてい る。この実習ではこれを利用して頂く。
もし自宅のＰＣ等で使用したい場合は、 オリジナルサイトからダ ウンロードして使う事が可能である。但し、オリジナルサイトでは最新 のv3.5.3がある。バージョンが異なると動作や使い方が変わることがあ るので、注意が必要である。皆さんが各自でより進んだ内容を勉強する 事は大変よいことであるが、試験、又はレポートで使用するバージョン は公平性を保つため、情報教育システムと同じバージョンとする。
スタートボタンから「Devel」→「processing」を選択して Processing を起動する。
次に以下のプログラムを入力して実行してみよ。 コピー＆ペーストは使用しないこと。
size(800,600)
for i in range(350):
print("X="+str(i*2)+" Y="+str(i))
point(i*2,i)
size 文で、描画する範囲を横８００ピクセル、縦６００ピクセルに指定している。
i を０から３４９まで、１つづ増やしながら、 Ｘ座標はｉの２倍、Ｙ座標はｉの位置に、点を描画している。 Ｘ，Ｙの値は数値でも出力している。
点が小さくて見にくい場合は、point の代わりに、ellipse を使うと見やすくなる。
ellipse(i*2, i, 1, 1)
ellipse は楕円を描画する関数である。半径の小さな楕円を描く事で、大きな点を描画している。
以下の事をやってみよ。
尚、Processingで使用できる関数等の説明が
以下のプログラムを入力して実行してみよ。コピー＆ペーストは使用しないこと。
def setup():
size(600,600)
print("width="+str(width)+" height="+str(height))
frameRate(30)
line(sx(-100),sy(0),sx(+100),sy(0))
line(sx(0),sy(-100),sx(0),sy(100))
x1=float(-100.0)
y1=float(-100.0)
x2=float(+100.0)
y2=float(+100.0)
def sx(x):
return width*(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
x=float(0.0)
def draw():
global x
global i
y=0.01*x*x
print("X="+str(x)+" Y="+str(y))
point(sx(x),sy(y)) # elipse(sx(x),sy(y),1,1)
x+=1.0
if x>100.0:
noLoop()
９～１３行は、座標変換のお膳立てをしている。
size 文で宣言された６００×６００ピクセルのスクリーンは、いわばＸが１～６００、Ｙが１～６００の
座標を持っていると言える。但し、Ｙは下が大きい事に留意せよ。
これだと我々が自由に座標系を使えない。そこで、計算で使用する座標系からスクリーンに表示する座標系への
変換を考える。このプログラムでは計算で使用するＸは、－１００．０から＋１００．０まで、
Ｙも同じ範囲とし、Ｙ軸は上が大きくなるように方向を逆転させる、一次変換を行う関数を定義している。
この関数はＸ，Ｙ毎に１つづつ定義されている。 計算に使用する座標系は９～１０行目で定義されている。 変換に必要なスクリーンのピクセル数は、size 文の実行後に、自動的に、width, height という変数に代入されている。
この関数は、下から４行目の point 関数で使用されている。 この様にすると、point関数が計算座標系のＸ、Ｙを引数に持っているかのごとく使える。
プログラムでは、２次関数を計算、描画している。
Ｘを０から１づつ増やし、１００を超えたらループを終了している。 Ｘの初期値、終了条件、増分、関数形等を変更して、色々試してみよ。
setup関数はプログラムの最初に１度だけ実行される。
draw関数は、何もしなければ、自動的に無限回、実行される。 繰り返し計算し、その結果を描画する場合には、draw関数の中に、１回分の描画するプログラムを書けば、 それを繰り返してくれる。但し、パラメータの更新や、終了条件の判断は自分で書く必要がある。
関数形を変えて、異なるグラフを描け。
初期値、終わり値を変えてみよ。
スクリーンの大きさを変えてみよ。また、計算に使う座標系も変えてみよ。
draw関数はsetup関数の中でframeRate関数で設定された頻度で繰り返される。 このプログラムでは３回／秒としてある。この数値を変えてみよ。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
import random
def setup():
size(600,600)
print("width="+str(width)+" height="+str(height))
frameRate(1000)
colorMode(RGB)
x1=float(-1.0)
y1=float(-1.0)
x2=float(+1.0)
y2=float(+1.0)
def sx(x):
return width*(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
random.seed()
s=float(0)
n=int(50000)
i=int(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+=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+=1
if i>=n:
noLoop()
pi=float(4*s/n)
print("pi="+str(pi))
このプログラムは、規則性のない数(乱数)を[-1,1]の区間で２つ(x,y)生成しその数の二乗和、
つまり原点からの距離が円の半径(1)よりも大きいかを判定を行なっている。
この事象判定を複数回(50000)程行う事で、50000回中何回円の内部に点が生成されたか(確率)を計算している。 今回の例では、まず上の確率を計算し、乱数(x,y)の定義域である正方形2×2の面積で規格化する事で円の面積を求めている。 円の半径は1であるので、円の面積が"π"になることを利用することで円周率を見積もっている。
→(乱数の定義域で表される正方形の面積 : 4)/(乱数の生成回数 : 50000)・(円内側に入った回数 : s)
このように、乱数を複数回生成する事で面積を見積もったり、自然界の崩壊現象などを擬似的に再現する事がある。 例えば、放射線検出器に放射線が入射する確率は、放射線検出器と放射線線源の距離や放射線検出器の形などに依存し、手計算で求める時が難しい。 その時、放射線検出器の形を再現し、放射線線源の崩壊を乱数を用いて再現する事で擬似的に放射線が放射線検出器に入射することを再現する事が可能である。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
def setup():
size(600, 600)
#frameRate(3)
line(sx(-1000),sy(0),sx(+1000),sy(0))
line(sx(0),sy(-1000),sx(0),sy(+1000))
x1=float(-1.0)
y1=float(-100.0)
x2=float(+10.0)
y2=float(+10.0)
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
m=float(1) #weight
x0=float(0.0) #initial position of x
y0=float(0.0) #initial position of y
vx0=float(10.0) #initial velocity of x
vy0=float(10.0) #initial velocity of y
dt=float(0.0050)
x=x0
y=y0
vx=vx0
vy=vy0
t=0.0
i=int(0)
g=float(-9.8)
def draw():
strokeWeight(3)
global g
global x
global y
global vx
global vy
global t
global i
# point(sx(t),sy(vx))
# point(sx(t),sy(vy))
# point(sx(t),sy(x))
point(sx(t),sy(y))
# point(sx(x),sy(y))
print("vx="+str(vx)+" vy="+str(vy)+" x="+str(x)+" y="+str(y)+" t="+str(t))
t=t+dt
vx=vx
vy=vy+g*dt
x=x+vx*dt
y=y+vy*dt
i=i+1
if i>2000:
noLoop()
このプログラムは、質点の放物運動を数値計算で解いている。
物理学で言う、数値計算で問題が「解けた」とは、
重力加速度は、－９．８ｍ／ｓ／ｓ（下方へなので符号は負）とし、 質点を座標（０ｍ，０ｍ）から初速度（１０ｍ／ｓ、１０ｍ／ｓ）で飛ばした時の運動を追跡している。
時間の刻みは０．０１秒とし、ある時刻から次の時刻の速度と座標は以下のように求めている。
Vx ← Vx
Vy ← Vy + g・Δt
X ← X + Vx・Δt
Y ← Y + Vy・Δt
これらの式は、速度と加速度の定義から分かると思う。
Vx(t) = ( X(t+Δt) - X(t) ) / Δt （Δt→０）
Ax(t) = ( Vx(t+Δt) - Vx(t) ) / Δt （Δt→０）
この様に１次の極限での近似を用いて数値積分する方法を、オイラー法という。
プログラムの最後の方に、Ｘ軸とＹ軸に何を描画するかが、何行かに渡り書いてあり、 １つを除きコメントアウトされている。これらをそれぞれ描画し、計算の結果を確認せよ。 初速度を変えて試してみよ。
このように、物理現象をシミュレーションする場合には、物理法則に関わる部分には一切修正を加えず、 初期値等の条件を変更する事で、様々な物理現象を再現できる。
Δtは０．０１秒で計算したが、これを０．０２秒に粗くしてみよ。 グラフでは分からない程度の差しかでないので、出力された数値を比較して、計算誤差が変化する事を確認せよ。
#（シャープ）は、コメント機能である。この記号から行末までは、コメントとなり、プログラムの実行に何の影響も与えない。 また、''' コメント ''' という構造を使用すれば、複数行にまたがったコメントも作成できる。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
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=float(-1.6e11)
y1=float(-1.6e11)
x2=float(+1.6e11)
y2=float(+1.6e11)
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
t=float(0.0)
dt=float(80000.0)
x=float(1.496e11)
y=float(0.0)
vx=float(0.0)
vy=float(2.979e4)
G=float(6.67384e-11)
Ms=float(1.989e30)
i=int(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+=1
if i>10000:
noLoop()
このプログラムは、太陽の周りを回る地球の動きをシミュレーションしている。
最初の方で、色を指定している。
初期条件として、太陽と地球の距離、地球の速度は、現実的な値を入れてあるが、 この場合、綺麗に円軌道が再現できる。
初速度を変えると楕円軌道になる筈なので、やってみよ。
初速を半分（1.5e4）にすると楕円軌道がだんだんずれていってしまう。 これは計算誤差が蓄積するからである。ｄｔの値を小さくすると計算誤差が減り、 より正しい結果を得られる事を確かめよ。
ｄｔの値を小さくした場合は、frameRate で大きな値を指定すると速く計算する。 指定しなかった時の値は６０である。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
def setup():
size(600, 600)
frameRate(100)
background(0,0,0)
stroke(255)
colorMode(HSB,360,100,100)
x1=float(6.784e11)
y1=float(-1.e11)
x2=float(+8.784e11)
y2=float(+1.e11)
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
G=float(6.674e-11)
pi=float(math.pi)
######## Parameters to chanege an initial condition of the probe ########
theta=float(pi/3) #Direction if momentum
disp=float(1.8e10) #Radial diplacement between Jupiter abd the probe
#########################################################################
def accel(m,r):
return G*m/(r**2)
t=float(0.0)
dt=float(10000.0)
x_jup=float(7.784e11*math.cos(-pi/30))
y_jup=float(7.784e11*math.sin(-pi/30))
x_prb=float(7.784e11-disp)
y_prb=float(-7.0e10)
vx_jup=float(1.306e4*(-math.sin(-pi/30)))
vy_jup=float(1.306e4*math.cos(-pi/30))
vx_prb=float(1.0e4*math.cos(theta))
vy_prb=float(1.0e4*math.sin(theta))
m_sol=float(1.989e30)
m_jup=float(1.898e27)
m_prb=float(721.9) #voyger 1 & 2
ax_jup=float(0.0)
ay_jup=float(0.0)
ax_prb=float(0.0)
ay_prb=float(0.0)
## the sun stand still ##
i=int(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=float(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=float(math.sqrt(x_jup**2+y_jup**2))
r_sol_prb=float(math.sqrt(x_prb**2+y_prb**2))
r_jup_prb=float(math.sqrt((x_jup-x_prb)**2+(y_jup-y_prb)**2))
ax_jup=float(-x_jup/r_sol_jup*accel(m_sol,r_sol_jup))
ay_jup=float(-y_jup/r_sol_jup*accel(m_sol,r_sol_jup))
ax_prb=float(-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=float(-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=float(vx_jup+ax_jup*dt)
vy_jup=float(vy_jup+ay_jup*dt)
vx_prb=float(vx_prb+ax_prb*dt)
vy_prb=float(vy_prb+ay_prb*dt)
x_jup=float(x_jup+vx_jup*dt)
y_jup=float(y_jup+vy_jup*dt)
x_prb=float(x_prb+vx_prb*dt)
y_prb=float(y_prb+vy_prb*dt)
i+=1
if i>100000:
noLoop()
このプログラムは、スイングバイをシミュレーションしている。
スイングバイについては、「講義資料（物理）追加」を参照せよ。
課題５は２体問題（２つの質点の運動を扱う）であったが、これは、太陽、木星、探査機の ３体問題である。尚、数学的には３体以上の厳密解はなく、高精度の結果を得るには数値的に 解くしか方法がない。また、３体ではあるが、太陽は座標原点に固定してある事に留意せよ。
「講義資料（物理）追加」を参照して、初期条件を変え、スイングバイが色々なパターンを 取る事を確認せよ。地球に帰還する軌道に載せたり、惑星軌道に載せたり、 次の惑星に向かったり、太陽系を脱出したりする事ができる。
アルゴリズムを改良する事で計算精度を上げる方法がある。
まず、課題５のプログラムを修正して、ｄｔを１００００、ｖｙの初期値を、１．０ｅ４にして、 実行してみよ。実行に時間がかかるようになるので、コメントアウトしてある frameRate(1000) を有効にして、 速く計算させよ。 楕円軌道がだんだんずれていく事が分かる。厳密な解ではこれは起きない。
次に、同じ初期条件と時間ステップで、４次のルンゲクッタ法を用いたプログラムを実行してみる。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
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=float(-1.6e11)
y1=float(-1.6e11)
x2=float(+1.6e11)
y2=float(+1.6e11)
def sx(x):
return width *(x-x1)/(x2-x1)
def sy(y):
return height*(y2-y)/(y2-y1)
t=float(0.0)
dt=float(10000.0)
x=float(1.496e11)
y=float(0.0)
vx=float(0.0)
vy=float(1.0e4)
G=float(6.67384e-11)
Ms=float(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+=1
if i>10000:
noLoop()
ルンゲクッタ法の場合は、見える範囲内では楕円はずれず、計算精度が上がっている事が分かる。
一方で計算量は大幅に増えているが、画面表示にかかる時間の方が大きいので、実際の待ち時間は 変わらない。
オイラー法でもｄｔを２０００とすれば楕円がずれない精度が出せる。 しかし同じ時間をシミュレートするための計算量は５倍となる。
オイラー法は、計算量が少なく、精度が低い。ルンゲクッタ法は、計算量が多く、精度が高い。 今回の例では非常に簡単な例であったため、両者の違いは大きくなかったが、 ほとんどの科学技術計算において、ルンゲクッタ法（更にはその改良版）が使われている。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
import math
i=int(0)
t=float(0.0)
dt=float(10.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.e30,5.e30,5.e30]
PI=math.pi
def setup():
x[0]=0.
y[0]=5.e10
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.e4
vy[0]=2.e4
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(1000)
G=float(6.674e-11)
x1=float(-1.e11)
y1=float(-1.e11)
x2=float(+1.e11)
y2=float(+1.e11)
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
global PI
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+=1
if i >10000:
noLoop()
このような場合、数学的に近似する（つまりモデルを変更する）方法もあるが、
それを行ってしまうと得られた数値の精度がモデルの精度以上にはできないという限界がある。
より高精度な結果を得るには、数学的な近似をせず、計算上の近似を高めていく手法しかない。
それがシミュレーションが必要な理由の１つである。