Phystemは何をしているか?

 今回作ってもらう物理シミュレーションの「土台」としてPhystemというもの(オブジェクト)が用意されている。

 このファイルではそのPhystemの使い方と、動作原理を述べておこう。

Phystemへの描画

 コンピュータの上に物体の運動を再現したい。その為に画面に物体を「描く」のだが、ここでは2次元で直交座標で扱うことにする。

 Phystemはその中にjavascriptのcanvasを持っている。canvasの詳しい使い方を知りたい人はjavascriptに関する文献(たとえば、html5.jpのCanvasリファレンスなど、ネット上のリファレンスを読むとよいが、ここでは最低限必要なことだけを記しておく。

 コンピュータ上でよく使われる座標系は「ピクセル」(液晶やCRTの輝点一つを1とする単位)になっていて、かつ右がx軸の正、上がy軸の正と決めてある場合が多い(そして原点はたいてい、画面の左上隅である)。

 図で示すと↑のような感じになっている。しかしこれは物理に使うには使いにくいので、まず座標系の向きは右がx軸の正、上がy軸 の正とする。原点は任意の場所に置けるようにし、かつ単位系もピクセルとは関係ない自由なスケールにしておく(縦横比はピクセルの比のまま1:1である)。

 Phystemの中では座標系のいろいろな値が上のように決められている。たとえばPhystemを表す変数がpsであったなら、ps.wが横幅を、 ps.leftxが左隅x座標、ps.rightxが右隅x座標である(当然であるが、ps.w = ps.rightx - ps.leftxである)。図では原点(0,0)が中にあるが、別に中になくてもよい。ps.canvasがキャンバス、ps.ctxが描画コンテキストなので、これを 使って絵を描くことができる。

プログラムの最初で、

ps=new Phystem("canvas1",20,14,0.8);

のように書いておくと、横幅20、縦幅14の描画領域がidが「canvas1」←これは適当に変えてよいのcanvasとして作られる。

canvasの場所は、ファイルの中の

<canvas id="canvas1"></canvas>

と書いた部分である。id=の後の部分が一致していることに注意。

最後の引数の0.8は、このcanvasがブラウザの画面上の横幅の何倍になるかを指定する。

この作り方だと、中央に原点が来る。つまり、x座標は-10から10に、y座標が-7から7になる。これを変更したければ、最初を

ps=new Phystem("canvas1",20,14,0.8,0,3);

のように変更する。こうすると、左隅が0になり、下の隅が-3になる。幅は変わらないので、x座標は0から20に、y座標は-3から11ということになる。

 基本的な図形を描くためには

ps.writePoint(x,y,r,col)点(●)を(x,y)に半径r、色colで描く
ps.writeLine(x1,y1,x2,y2,col,w)(x1,y1)から(x2,y2)に色colで幅wの線を引く。
ps.writeCircle(x,y,r,col,w)(x,y)を中心として半径rの円を色colで描く。
ps.fillCircle(x,y,r,col,w)上に同じだが中身を塗りつぶす
ps.writeRect(x1,y1,x2,y2,col,w)(x1,y1)(x2,y2)を角にする長方形を描く。
ps.fillRect(x1,y1,x2,y2,col,w)上に同じだが中身を塗りつぶす
ps.writeCoordinateXY()x,y座標系を描く
ps.writeCoordinateX()x座標系のみ描く
ps.writeCoordinateY()y座標系のみ描く

などの命令があるので利用してもよい。

javascriptでcanvasに図を描く一般的な方法については、Mozilla Developer Network(MDN)の「canvasに図形を描く」のページなどを参照のこと。ただし、このページで説明されている方法で使われている「グラフィックコンテキスト」は、ps.ctx(Phystemの名前がpsでないのなら、ps.の部分を書き換えること)で定義されている。


DynamicalObjectとNonDynamicalObject

 Phystemは「系」である。物理における「系」とは互いに相互作用するいくつかの物体(このPhystemで使えるのは質点のみ)の集まりであるが、最初にnew Phystemした段階では、物体は何もない。これから物体を系に登録していかなくてはいけない。

 Phystemに「登録」することができる物体はDynamicalObjectとNonDynamicalObjectの2種類(と、これらを「継承」したオブジェクト)がある。DynamicalObjectは名前の通り「動力学的」な「もの」で、運動方程式に従って運動している物体である。NonDynamicalObjectの方は動いていないか、動いているとしても運動方程式によらない運動をしている物体である。Phystemによるシミュレーションの準備は、

  1. Phystemを初期化する。
  2. Phystemに物体を登録する。
  3. 物体が従う物理法則(働く力)を記述する。
  4. ps.start();でPhystemをスタートさせる。

と行われる。その後は物理法則にしたがって物体が運動する。

 物体の登録は、

var m=new DynamicalObject(ps,x,y,vx,vy,m,"rgba(255,255,0,0.5)");

のように行う。最初のpsが登録される「系」である。x,yが物体の初期位置、vx,vyが物体の初速度であり、mは質量である。最後の"rgba(255,255,0,0.5)"は色を指定している。

系にDynamicalObjectが登録されたことにより、以後このDynamicalObjectは一定時間ごとにそれぞれの物理法則にしたがって『次の時刻での位置』を計算され、そうやって動いた結果がcanvasに表示される。

new DynamicalObjectで起こること

 var m=new DynamicalObject(ps,x,y,vx,vy,m,"rgba(255,255,0,0.5)");を行うと、Phystemが持っているdyObjsという配列にこのmが登録される。

 dyObjsに登録されていると、

  1. Phystemが「画面を描く」時に.draw()という「自分自身を描画する関数」が呼び出される。
  2. PhystemのdrawVFlgがtrueなら、さらに.drawV()(速度を描く関数)が呼び出される。
  3. PhystemのdrawFFlgがtrueなら、さらに.drawF()(力を描く関数)が呼び出される。
  4. Phystemが時間発展を計算する時に、下で説明する方法で運動方程式に従って運動する。

という処理が行われる(1秒に何回か、決まった回数)。

 .draw().drawV()などは物体のプロパティとして定義された関数で、変更することもできる(しなければm.rで設定された半径の円がm.colで設定された色で描かれる)。

ループの中で起こっていること

 Phystemがどのように力学的計算を行っているか、その部分を抜粋しよう。

 この部分は諸君は編集する必要はない(したければしてもよいけど、その場合は課題提出の時に書き添えておいて欲しい)。どういう動きをしているのか中身を知りたい人のためにここで説明しておくだけのことである。
	for (i = 0; i < N; i++) {
		now = this.dyObjs[i];
		F = new Vector(0, 0);
		for (j = 0 ; j < N; j++) {
		  // j番目の動的オブジェクトからi番目の動的オブジェクトへの力
		  if (i !== j) {
		  f = this.dyObjs[j].interactionForce(now);
		  F.add(f);
		  }
		}
		for (j = 0 ; j < this.ndObjs.length ; j++) {
		  // j番目の【非】動的オブジェクトからi番目の動的オブジェクトへの力
		  f = this.ndObjs[j].interactionForce(now);
		  F.add(f);
		}
		now.v.add(F.prod(dt / now.mass));  // v = v+(F(x+vΔt)/m)Δt
		now.pos.copyFrom(now.npos);
	}

 この中のnow = this.dyObjs[i];nowが今着目している物体で、次のループではまずDynamicalObjectの一つ一つからこのnowに働く力をFという変数に足している。

 その計算をしている部分が f = this.dyObjs[j].interactionForce(now);の部分で、interactionForce()という、this.dyObjs[j](←j番目に登録されたDynamicalObject)のプロパティになっている関数を(nowを引数に)呼び出して力を計算させている。つまり、「this.dyObjs[j]がnowにどんな力を及ぼすか」はthis.dyObjs[j]というオブジェクトが知っている、というわけである。

 つまり、.interactionForce()という関数をオブジェクトごとに書き換えることで「いろんな力」を記述できるということである。この「自分が及ぼす力は自分が知っている」というシステムのおかげで、働く力が変わってもPhystemのコード(上に引用したループ部分)を変える必要はないのが、オブジェクト指向プログラミングの恩恵の一つ。
y

 この後次に、登録されているNonDynamicalObject(this.ndObjs[j])に対しても同様に、nowに対してどのような力を及ぼすかを計算する。これで系にある全ての物体からの力を計算したことになる。

 ところで、運動方程式は

という式であるが、これを

のように連立微分方程式として解いている。具体的には

  1. まず現在位置(pos)と現在の速度(v)から、Δt後の場所(npos)を計算する(これは上に引用したコードに入るまえに終わっている)。
  2. その場所において働く力を計算し、その力から加速度を出して速度を変化させる。now.v.add(F.prod(dt / now.mass));
  3. 現在位置(pos)に計算しておいた次の位置(npos)をコピーする。now.pos.copyFrom(now.npos);

という手順を行っている。

 先に位置を更新して、更新した位置を下に速度を更新するのは、「シンプレクティック数値積分法」という数値積分の方法で、オイラーの方法やルンゲ・クッタの方法に比べエネルギー保存則の成立において有利だと言われている。

本当のところはシンプレクティック数値積分法は解析力学のハミルトン形式を利用したもので、「速度」ではなく「運動量」を使うのが正しい方法。速度と運動量が単純なp=mvの関係なら同じことになる。Phystemでは単純な場合のみに対応している。まだハミルトン形式を知らないもしくは慣れてない人が多いと思うので、この部分については(厳密ではないけど)気にしないでやっていって欲しい。

DynamicalObjectとNonDynamicalObjectの中身

 この二つはphystem.jsの中で定義されている。phystem.jsの最初の方を見るとまず、

var PhysObject = function (ps, x, y, c, rr) {
	this.pos = new Vector(
		(x === undefined) ? 0 : x,
		(y === undefined) ? 0 : y
	);
	// 色をセット
	this.col = (c === undefined) ? "rgba(0,0,0,0.4)" : c;
	this.npos = this.pos.makeCopy();
	// npは、「次の位置」(シンプレクティック法で微分方程式を解いているのでこれが必要)
	this.r = (rr === undefined) ? 0.5 :rr;
	this.ps = ps; // 自分の属している系をセットしておく。
	// この「系」は描画用のキャンバスも含んでいる。
};
PhysObject.prototype = {
	// 描画:●を描く。
	draw: function () {
		this.ps.fillCircle(this.pos.x, this.pos.y, this.r, this.col);
	},
	// 相互作用力。「to」へと及ぼす力を返す。
	// 現在はスタブなので0ベクトルを返す。
	interactionForce: function (to) {
		return new Vector(0, 0);
	},
	// 変位べクトル。現在位置ではなく、△t後の位置で計算する。
	displacementFrom: function (from) {
		return this.npos.diff(from.npos);
	}
};

という部分がある。PhysObjectはDynamicalObjectとNonDynamicalObjectに共通する属性を持ったクラスで、このあとPhysObjectはDynamicalObjectとNonDynamicalObjectに「枝分かれ」する。

 PhysObject自体は、this.pos(現在位置のベクトル)、this.npos(しばらく後での位置のベクトル)なぜこれがいるのかは後で説明する。、this.col(色)、this.ps(自分がどの系に属しているか)、this.r(自分の半径:円だと考えている)だけを属性(プロパティ)として持つ。

NonDynamicalObjectのコンストラクタは、

// 非力学オブジェクト
// 運動しないもしくは運動方程式によらない運動をするオブジェクト
var NonDynamicalObject = function (ps, x, y, c, r) {
	PhysObject.call(this, ps, x, y, c, r);
	if (ps !== undefined) {
		ps.ndObjs.push(this);
		// 「系」に、自分自身を登録する。
	}
};

NonDynamicalObject.prototype = Object.create(PhysObject.prototype);
NonDynamicalObject.prototype.constructor = NonDynamicalObject;

のようになっている(最後の2行はjavascriptでクラスの継承を行うときのおまじないで、これでPhysObjectの持っているプロトタイプを受け継ぐことができる)。

 ps.ndOjbs.push(this);は系(ps)の持っているndOjbsという名前のリストに自分自身(this)を登録している。これで、上で説明したループの中でこの物体が他の物体に及ぼす力が計算されるようになる。

一方、DynamicalObjectの方は、

// DynamicalObject:「動く物体」のクラス
// 親クラスはPhysObject(物理的物体?)
// mass(質量)と、fList(働いている力のリスト)という属性が追加される。
// yavというオブジェクトを持っているが、外部から使う必要は多分ない。
// 初期位置、初速度、質量、色、と「どの系に属すか」が呼び出しパラメータ(初速度以降は省略可)
var DynamicalObject = function (ps, x, y, vx, vy, m, c, r) {
	PhysObject.call(this, ps, x, y, c, r);
	if (ps !== undefined) {
		// 自分自身を系の「動くものリスト」に登録する。
		ps.dyObjs.push(this);
	}
	this.v = new Vector(
		(vx === undefined) ? 0 : vx,
		(vy === undefined) ? 0 : vy
	);
	this.mass = (m === undefined) ? 1 : m;
	this.makeV(); // 動くものは速度や力を表示するので、その為の矢印を先に作っておく。
	this.fList = []; // 力のリスト。
};

のようにコンストラクタが定義されている(おまじない部分もあるが、省略した)。

ps.dyOjbs.push(this);で、今度は系の持っているdyObjsに自分自身を登録する。Dynamicalなので初速度や、働く力など、他の量も最初に用意している。

「力」のプログラミング

.interactionForce()の設定

 上で述べたように、.interactionForce()という関数が「この物体が他の物体にどんな力を及ぼすか」を設定しているので、ここを書き換えることで物理法則を変更できる。

 最初の段階では、phystem.jsの中の

      interactionForce: function (to) {
        return new Vector(0, 0);
   },

の部分でPhysObjectが出す力が「0ベクトル」に設定されている。

 「書き換える」と言ってもここを変えてしまうと、全ての物体の力がいっせいに変わってしまう。たとえばある物体mの「力」を変更するには、

m.interactionForce=function(to) {
   mからtoに働く力を計算する部分
   return 結果のVector;
}

のようなコードを付け足せばよい。

ただしこれで終わってしまうとinteraction(相互作用)にならない。to→mの方向の力を記述してないから!

 しかし、これだとたくさんの物体に力を設定するのがたいへんである。そこでここでもオブジェクト指向プログラミングの恩恵を使うことにする。

オブジェクト指向を使って力の設定を行う

 実例で示そう。

 表紙でも見せていた分子運動シミュレーションでは、

			var Molecule=function (ps,x,y,vx,vy,m,c) {
				DynamicalObject.call(this,ps,x,y,vx,vy,m,c);
			};
			Molecule.prototype=Object.create(DynamicalObject.prototype);
			Molecule.prototype.constructor=Molecule;
			Molecule.prototype.interactionForce=function (to) {
				var r=to.npos.diff(this.npos);
				var rlen=r.length();
				return r.prod(10*Math.exp(-2*rlen*rlen));
			};

というコードが書かれている。

 1行目のvar Molecule=function (ps,x,y,vx,vy,m,c) {から始まる部分は、Moleculeという新しいオブジェクトが作られた時にどのような関数が呼ばれるかが書かれている。その中身はDynamicalObject.call(this,ps,x,y,vx,vy,m,c);で、これは「DynamicalObjectが作られた時に呼ばれる関数を実行してくれ」という、言わば「丸投げ」である。

その後に書かれている、

			Molecule.prototype=Object.create(DynamicalObject.prototype);
			Molecule.prototype.constructor=Molecule;

の部分は、MoleculeのプロトタイプをDynamicalObjectのプロトタイプと一致させるための命令である(詳細は省くので呪文と思って実行すること)。

 その後のMolecule.prototype.interactionForce=function (m2) {で始まる部分がいよいよ力を定義する部分で、Molecule.protypeとあることで「Moleculeのプロトタイプ」にinteractionForceを登録する(プロトタイプに登録することで、以後new Moleculeで作られたオブジェクトにはこのinteractionForceが登録された状態になる)。

 力の中身についても説明しておこう。

var r=to.npos.diff(this.npos);では、rという変数にtoのnpos(次の段階での位置)とthis(自分自身)のnposとの差のベクトルを設定している。

次のvar rlen=r.length();はそのベクトルの長さをrlenという新しい変数に代入した。

r.prod(10*Math.exp(-2*rlen*rlen))は、差のベクトルrを(10*Math.exp(-2*rlen*rlen)倍して、これを働く力として返している。こうすることで自分から相手に向かう方向(図の緑矢印)を向いた力(大きさは最後の掛算で調整している)をtoに対して与えることになる。

この作り方では、全ての物体に同じ式の力を与えることになる。物体に応じて変えたかったら、interactionForceの中身をいじる。

新しいタイプの「物体」を作りたかったら、以下の手順を実行する。

  1. プログラムの初期化の場所に、
    			var 新クラス=function (ps,x,y,vx,vy,m,c) {
    				DynamicalObject.call(this,ps,x,y,vx,vy,m,c);
    			};
    			新クラス.prototype=Object.create(DynamicalObject.prototype);
    			新クラス.prototype.constructor=新クラス;
    を書く。新クラスの部分は自分で決めた新しいクラスの名前を入れること。
  2. その物体が他の物体に及ぼす力を、新クラス..prototype.interactionForce=function (to) {の後に、「ベクトルを返す関数」として書く。
  3. その物体を作る。たとえば、m=new 新クラス(ps,x,y,vx,vy,m,c);のように。

以下は、万有引力を及ぼし合う惑星を定義する場合のinteractionForceの例である。

Planet.prototype.interactionForce=function (m2) {
    var r=m2.npos.diff(this.npos);
    var rlen=r.length();
    var rdis=m2.r + this.r;
    if( rlen > rdis ) {
        return r.prod(-G*this.mass*m2.mass/ (rlen*rlen*rlen) );
    } else {
        return r.prod(-G*this.mass*m2.mass/ (rdis*rdis*rdis) ); 
    }
};
 万有引力は距離の自乗に反比例だが、距離が0になると発散するので、半径の和より近づいた時は定義を変えて距離に比例する力にしている。

NonDynamicalOjbectの例

重力mg

 ここではDynamicalObjectの場合を考えたが、NonDynamicalObjectも(他の物体に力を及ぼすことはできるので)同様に力を定義できる。例として、ps.makeGravity();で重力を作っている方法を説明しよう。

 makeGravityは以下のように定義されている(phystem.jsの中にある)。

// 重力を追加する
Phystem.prototype.makeGravity = function (g) {
	this.gravity = (g === undefined) ? 1 : g;
	var earth = new NonDynamicalObject(this);
	earth.interactionForce = function (to) {
		return new Vector(0, -to.mass * to.ps.gravity);
	};
	earth.draw = function () {};
};

ここでやっていることは「地球(earth)」というNonDynamicalObjectを作って(一個しか作ってないのでクラスを定義していない)、そのearthのinteractionForceを、Vector(0, -to.mass * to.ps.gravity)にしている(これは下向きmgを表している)。これでPhystemに「earth」という非力学的オブジェクトが登録され、earthからの力として重力が他の全ての物体にかかることになる。

なお最後のearth.draw = function () {};では、地球を描画する必要はないので、描画する関数を「何もしない{}」にセットしている。

 makeAir();(←空気抵抗を追加する)も同様なので、phystem.jsを見てみよう。

万有引力

 前回で見せた2個の人工衛星の追っかけっ子の中では、

  	    M=new NonDynamicalObject(gs,0,0,"rgba(0,200,0,0.8)");
		M.interactionForce=function (to) {
			var r=to.npos;
			var rlen=r.length();
			if (rlen > 0.5) {
				return r.quot(-(rlen*rlen*rlen));
			} else {
				return r.prod(-8);
			}
		}

のようにして、地球を表すオブジェクトMが人工衛星に万有引力を及ぼすようにした。

バネの力

 他に、NonDynamicalObjectの例としてはBaneもある(バネは動くのだが、その動きは運動方程式に応じたものではないので、NonDynamicalである)。bane.jsを見てみるとよい。

// m1とm2をつなぐバネ
var Bane = function (ps, m1, m2, k, l, c) {
	// バネはNonDynamicalObjectである。
	NonDynamicalObject.call(this, ps, 0, 0, c);
	this.m1 = m1;
	this.m2 = m2;
	// バネ定数(デフォルトは1)
	this.k = (k === undefined) ? 1 : k;
	// 自然長(デフォルトは0)
	this.l = (l === undefined) ? 0 : l;
};

 まずBaneはNonDynamicalObjectの一種なので、NonDynamicalObjectのコンストラクタを呼び出す(この後に、プロトタイプを一致させるためのいつもの呪文もある)。

 バネは二つのDynamicalObjectをつないでいるものなので、その二つをthis.m1とthis.m2に格納する。その後はバネ定数と自然長をセットしている。

 Baneのプロトタイプがその後につづいているが、

Bane.prototype.interactionForce = function (to) {
	if (to === this.m1) {
		return this.force1();
	} else if (to === this.m2) {
		return this.force2();
	} else {
		return new Vector(0, 0);
	}
};

の部分で.interactionForce()を定義しているが、これはto(力を呼ぼす相手)がm1かm2だった時にだけ0でない力を返すようになっている(つながっていない物体にまでバネが力を出したらたいへんだ)。力の計算自体は.force1()と.force2()という別の(やはりプロトタイプで定義された)関数に丸投げされている(どうなっているか知りたい人はソースを確認すること)。