@pythonnoob (noob python)posted at 2021-11-05

pythonを用いた蔵本モデルのシミュレーション

Q&A

解決したいこと

現在、pythonを用いて蔵本モデルのシミュレーションを行うべくプログラミングに奮闘しております。（恥ずかしながらど素人です）
知りたいのはN個の振動子が単位円上を運動し、同期していく様子のアニメーションです。
当サイトでC言語を用いて解析している方のコードを参考にPythonに応用させたいのですが、以下のようなエラーが出てしまいます。どのようにコードを改善したらよいのかご教授いただければ幸いです。もし、可能であれば回答者様のコードも教えていただければ嬉しいです。

発生している問題・エラー

IndexError: list index out of range

該当するソースコード

import numpy as np
import random
import matplotlib.pyplot as plt
figure, axes = plt.subplots()

N = 10
K = 1
dt = 0.01
maxt=10
t = np.arange(0,maxt,dt)
theta = []
omega = []
thetanew = []

for i in range(N):
    thetanew.append(2.0*np.pi*random.uniform(0,1))
    omega.append(np.random.standard_cauchy())

for i in range(N):
    thetanew[i] += omega[i]*dt
for i in range(N):
    for j in range(N):
        thetanew[j] += K*np.sin(theta[j]-theta[i])/N
for i in range(N):
    thetanew[i] = theta[i]

試したこと

固有振動数ωはコーシー分布、i番目の振動子がスタートする位相は０～２πにランダムに存在させています。
ど素人なので質問も的を得ていないかもしれませんがよろしくお願いいたします。

2Answer

@ligun posted at 2021-11-05

ぱっと見ですがthetaに何も入れておらず空配列なのが原因な気がします。

0Like

Comments

@pythonnoob
Questioner
回答いただきありがとうございます。
自分もthetaのリストのところでエラーが出ている気がしたのですが対策できるわけでもなく...
引き続き勉強していきたいと思います！

@bizzpaper posted at 2021-11-05

こんな感じでどうでしょうか。

C言語のコードを参考に書かれているとのことですが、以下のように配列のブロードキャストやリスト内包表記などPython特有の書き方を利用するとC言語よりもスマートに書けます。

引き続きがんばってください！

calc.py

# coding: utf-8
import math
import random
import numpy as np
import pandas as pd

# <出力ファイル名>
file = "./result.csv"

# <定数の設定>
N = 40
K = 10
dt = 0.01
t0 = 0
te = 5

if __name__ == "__main__":
	# <変数と初期値の準備>
	time = np.arange(t0, te + dt, dt)
	theta = np.array([2.0 * np.pi * random.uniform(0,1) for i in range(N)])
	omega = np.array([np.random.standard_cauchy() for i in range(N) ])

	cols = ["Time"] + [f"p{i:03}" for i in range(N)]
	data = pd.DataFrame([], index=[], columns=cols)

	# <計算の実行>
	for t in time:
		print(f"{t:.2f} sec ({100 * t / (te - t0):.2f} %)\r", end="")
		# 格納
		row = dict(**{"Time": t}, **{f"p{i:03}": theta[i] for i in range(N)})
		data = data.append(row, ignore_index=True)

		# 更新
		delta = omega
		for i in range(N):
			for j in range(N):
				# 配列演算ではないので、numpyよりmathのほうが計算が高速
				delta[j] += K * math.sin(theta[j] - theta[i]) / N
		theta += delta * dt

	print(t)

	# <CSVに保存>
	data.to_csv(file, index=False)

animate.py

# coding: utf-8
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from calc import file, t0, te

# <出力ファイル名>
out = "./animation.gif"

# <アニメーションの描画間隔(msec)>
intrvl = 10

# <結果の読込>
data = pd.read_csv(file, index_col=None, header=0)
cnt, _ = data.shape

# <アニメーションの実装>
def gen():
	for _, row in data.iterrows():
		t = row["Time"]
		print(f"{t:.2f} sec ({100 * t / (te - t0):.2f} %)\r", end="")
		yield row

def update(row):
	# 図のクリア
	ax.cla()

	# 行をアンパック
	time = row["Time"]
	theta = row["p001":].values

	# デカルト座標に変換
	x = np.cos(theta)
	y = np.sin(theta)

	# 描画
	ax.scatter(x, y, facecolor="None", edgecolors="red")
	ax.text(
		0.05, 0.95, f"$t = {time:.3f}$ sec",
		verticalalignment='top', transform=ax.transAxes
	)

	rlim = 1.2
	ax.set_xlim(-rlim, rlim)
	ax.set_ylim(-rlim, rlim)
	ax.set_aspect("equal")

# <アニメーションの生成>
fig, ax = plt.subplots()
ani = FuncAnimation(fig, update, gen, interval=10, save_count=cnt)
ani.save(out)
# plt.show()

0Like

Comments

@pythonnoob
Questioner
回答いただきありがとうございます。
自分の稚拙なコーディングよりも見やすい上に、学ばせていただけることが多く感謝申し上げます。
初期値（omega,theta)はあのような感じで記述できるのですね、とても新鮮でした。
勉強に励みます。また、壁に当たったらご教授ください！
@bizzpaper
お力添えできたようで何よりです。

ちなみに、エラーの原因は @ligun さんのおっしゃられる通り、空配列であることです。
最後のfor文の代入が逆で、正しくはtheta[i] = theta_new[i]です。

今後は解析してグラフとかかなと思うので、私のコードを関数化してパラメータスタディしやすい形にできるとよいですね。
@pythonnoob
Questioner
なるほど、承知いたしました。
自分のコーディングを見直してみようと思います！

蔵本モデルにもあるように秩序定数の解析になるかと思います。
またお力添えしていただけましたら幸いです。

Are you sure you want to delete the question?

pythonを用いた蔵本モデルのシミュレーション

解決したいこと

発生している問題・エラー

該当するソースコード

試したこと

2Answer

Comments

Comments

Your answer might help someone💌