PythonとProcessingの学習のために、これらを用いて数学アニメーションを作成した。Pythonは数値計算ライブラリScipyを用いて手軽に数値計算ができるし、Processingは計算結果を可視化することに長けているから、両者の良い所取りをしようというわけだ。内容としては、①PythonでLorenz方程式を数値的に解き、結果をCSVファイルに保存する、②ProcessingでCSVファイルの内容を描画する、の2段階になっている。備忘録を兼ねてコード及び内容の説明を残しておく。ちなみに当方、プログラミングは初心者なのでもっと効率の良い方法があるかもしれないことは先に断っておく。
最初に実行結果を示すと以下のようになった。やはり視覚的にわかるアウトプットがあるとやりがいがあるものだ。
Processingについて
Processingは、Javaをベースとしたグラフィックデザイン機能に特化したプログラミング言語及びその総合開発環境のことだ。インストール方法などは検索すると多くの情報が得られるのでここでは割愛する。Processingを使用すると簡単なコードでいろいろな図形の描画や視覚的効果を実現することができるため、Pythonの数値計算のデータを読み込ませてグラフィックを生成して動画を作成してみた。
Lorenz方程式について
ローレンツ方程式は、カオス理論の分野で有名な常微分方程式だ。決定論的な方程式であるが初期値鋭敏性を持つため、無限の精度で初期値を知ることができない限り未来を予測できないことになる。あるパラメータにおけるこの方程式の解軌道は、ストレンジアトラクタと呼ばれ、その名の通り奇妙な軌道を描く。単純な規則から生成されたとは思えないほど複雑な形状が得られるのでよく数値計算の題材になる。題材は何でも良かったが、ビジュアル的に面白いのでこの方程式を採用した。方程式は以下のような形をしていて、パラメータはp=10, r=28, b=8/3として計算している。
\begin{cases}
\displaystyle \frac {dx}{dt}=-px+py\\
\displaystyle \frac {dy}{dt}=-xz+rx-y\\
\displaystyle \frac {dz}{dt}=xy-bz\\
\end{cases}
Python側コード
コード全体
import numpy as np
from scipy.integrate import odeint
import csv
p = 10.0
r = 28.0
b = 8/3
def func(state, t):
dydt = np.zeros_like(state)
dydt[0] = -p*state[0] + p*state[1]
dydt[1] = -state[0]*state[2] + r*state[0] -state[1]
dydt[2] = state[0]*state[1] - b*state[2]
return dydt
# 初期状態
x0 = 0.0
y0 = 0.1
z0 = 50
state = [x0, y0, z0]
#計算ステップ
dt = 0.01
t = np.arange(0.0, 1000, dt)
sol = odeint(func, state, t)
with open('sol', 'wt') as fout:
csvout=csv.writer(fout)
csvout.writerows(sol)
#plt.show()
説明
以下の部分ではローレンツ方程式を定義している。
p = 10.0
r = 28.0
b = 8/3
def func(state, t):
dydt = np.zeros_like(state)
dydt[0] = -p*state[0] + p*state[1]
dydt[1] = -state[0]*state[2] + r*state[0] -state[1]
dydt[2] = state[0]*state[1] - b*state[2]
return dydt
初期条件および計算のステップを設定して、odeintという常微分方程式のソルバーを使用してローレンツ方程式を数値的に解く。
# 初期状態 x0 = 0.0 y0 = 0.1 z0 = 50 state = [x0, y0, z0] #計算ステップ dt = 0.01 t = np.arange(0.0, 1000, dt) sol = odeint(func, state, t)
最後にCSV形式で計算結果を保存している。
with open('sol', 'wt') as fout:
csvout=csv.writer(fout)
csvout.writerows(sol)
#plt.show()
Processing側コード
コード全体
ArrayList<Float> x = new ArrayList<Float>();
ArrayList<Float> y = new ArrayList<Float>();
ArrayList<Float> z = new ArrayList<Float>();
int i = 0;
float rot = 0;
void setup(){
background(0); //背景黒
size(500, 500, P3D); //3Dモード
colorMode(HSB); //カラーモードHSB
String csvDataLine[] = loadStrings("sol"); //CSVファイルを読み込む
//項目ごとにバラす
for(int i = 0; i < csvDataLine.length; i=i+2){
String [] csvDataRow = split(csvDataLine[i],','); //行毎にバラす
for(int j = 0; j < csvDataRow.length; j=j+3){
x.add(float(csvDataRow[j])); //1列目がx
y.add(float(csvDataRow[j+1])); //2列目がy
z.add(float(csvDataRow[j+2])); //1列目がz
}
}
}
void draw(){
scale(5);//スケール5倍
translate(width/10, height/10, 0);//原点を画面中心にする
camera(75*cos(rot), 75*sin(rot), 0, 0, 0, 25, 0, 0, -1);//カメラの回転
rot += 0.01;//カメラ回転のパラメータ
background(0);//残像の消去用
noFill(); //必要？
float hu = 0;
beginShape(POINTS); //点で構成
for ( int n = 0; n <= i ; n++){
vertex(x.get(n), y.get(n), z.get(n)); //頂点を描画
stroke(hu, 255, 255); //頂点の色
strokeWeight(5); //点の大きさ
hu = hu + 0.1;
if (hu > 255){
hu = 0;
}
}
endShape();
//動画として保存するときはコメント外す
saveFrame("frames/####.png");
//drawのループ回数をカウント
i = i + 1;
if(i > 100000){
i = 0;
}
}
説明
setup関数の中身を説明する。CSVファイルを読み込み、各行をcsvDataRowに格納していく。さらにcsvDataRowの1列目をx,、2列目をy、3列目をzに格納する。
String csvDataLine[] = loadStrings("sol"); //CSVファイルを読み込む
//項目ごとにバラす
for(int i = 0; i < csvDataLine.length; i=i+2){
String [] csvDataRow = split(csvDataLine[i],','); //行毎にバラす
for(int j = 0; j < csvDataRow.length; j=j+3){
x.add(float(csvDataRow[j])); //1列目がx
y.add(float(csvDataRow[j+1])); //2列目がy
z.add(float(csvDataRow[j+2])); //1列目がz
draw関数の中身を説明する。ここではカメラの回転と解曲線の描画を主に行っている。cameraの位置座標に三角関数を入れることで回転を行っている。回転を行った後、backgraoundを再描画しないと残像が残ってしまう点は注意が必要。
beginShapeからendShapeの間が解曲線の描画部分となる。setup関数で定義したx, y, z座標を用いて頂点を描画する。strokeとstrokeWeightは頂点の色と大きさをそれぞれ設定している。
saveFrame(“frames/####.png”);
の記載により0001.pngから連番のpngファイルがdrawが実行される毎に作成される。drawの実行回数はiでカウントを行っている。
void draw(){
scale(5);//スケール5倍
translate(width/10, height/10, 0);//原点を画面中心にする
camera(75*cos(rot), 75*sin(rot), 0, 0, 0, 25, 0, 0, -1);//カメラの回転
rot += 0.01;//カメラ回転のパラメータ
background(0);//残像の消去用
noFill(); //必要？
float hu = 0;
beginShape(POINTS); //点で構成
for ( int n = 0; n <= i ; n++){
vertex(x.get(n), y.get(n), z.get(n)); //頂点を描画
stroke(hu, 255, 255); //頂点の色
strokeWeight(5); //点の大きさ
hu = hu + 0.1;
if (hu > 255){
hu = 0;
}
}
endShape();
//動画として保存するときはコメント外す
saveFrame("frames/####.png");
//drawのループ回数をカウント
i = i + 1;
if(i > 100000){
i = 0;
}
}
まとめ
pythonによる計算結果をprocessingによって可視化した。今回は、Lorenz方程式に関して可視化したが、どのような方程式に関しても応用が利くのでいろいろ作ってみると面白いかもしれない。
