@pythonnoob (noob python)posted at 2021-12-04

for文について

Q&A

解決したいこと

for文を用いて解析を行っているのですがKの値を変えていっても出てくる値がすべて一緒なので、どこがおかしいのか知りたいです。
解析したい関数は以下になります。

\dot{θ} = ω_i+\frac{K}{N}\sum_{k=1}^{N} sin(θ_j-θ_i)

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

該当するソースコード

ソースコードを入力
import numpy as np
import random
import matplotlib.pyplot as plt
import time

start=time.time()


def theta_mat(theta):
    theta=theta.reshape(1, -1)
    theta_mat=theta.copy()

    return theta_mat


N = 3#振動子数
Ks = np.arange(1.5,5,0.1) #結合強度
theta = np.array([2.0 * np.pi * random.uniform(0,1) for i in range(N)]).reshape(1, -1)  # start point

omega=np.array([np.random.standard_cauchy() for i in range(N)]).reshape(1,-1)# anglular velocity
t0=0
tf=300#測りたい時間の10倍
ts=np.arange(t0,tf,1)
h=0.1        # micro time

Rs=[]
for K in Ks:
    for t in ts:
        print(K)
        theta_now=theta[t, :].reshape(-1, 1)
        sumx=np.sum(np.cos(theta_now))
        sumy=np.sum(np.sin(theta_now))
        R=(sumx**2+sumy**2)**(1/2)/N
        Rs.append(R)

        k1=(omega+K*np.sum(np.sin(theta_mat(theta_now)-theta_now),axis=1)/N).reshape(-1, 1)
        theta4k2=theta_now+(h/2)*k1
        k2=(omega+K*np.sum(np.sin(theta_mat(theta4k2)-theta4k2),1)/N).reshape(-1, 1)
        theta4k3=theta_now+(h/2)*k2
        k3=(omega+K*np.sum(np.sin(theta_mat(theta4k3)-theta4k3),1)/N).reshape(-1, 1)
        theta4k4=theta_now+h*k3
        k4=(omega+K*np.sum(np.sin(theta_mat(theta4k4)-theta4k4),1)/N).reshape(-1, 1)

        theta_next=theta_now+(h/6)*(k1+2*k2+2*k3+k4)
        theta_next=np.mod(theta_next, 2*np.pi)
        theta=np.r_[theta, theta_next.reshape(1, N)]

   #print(sumx,sumy,R)

        t+=1
        #print(t*h)
    print(Rs)
    Rave=np.sum(Rs)/len(Rs)
    plt.scatter(K,Rave)
    #print(Rs,Rave)
    Rs.clear()
    #print(Rs)
end=time.time()
plt.ylabel('R')
plt.xlabel('K')
plt.grid()
print('処理時間',end-start)
plt.show()

自分で試したこと

printでKの値が変わっていることを確認したのですが、それがうまく代入できていない（？）ようです。
もし、お気づきの点ありましたらご教授ください。

3Answer

@bizzpaper posted at 2021-12-07

こちらで回答したものを少し変えた実装例になります。

アドバイスとしては、今回の実装のように計算と描画は分けた方が良いです。
これは、CSVなどに角度のデータを一時的に保存できるようにするためです。
CSVに一時的に保存することで、グラフの装飾を変えたい場合や解析方法を変えたい場合に、もう一度計算をやり直さなくて済むようになります。

かなり面白いテーマなので、最後まであきらめずにがんばってください！

$N=100, T=500$ の結果

$N$や$T$を大きくすれば、理論値に近づいていくと思います。

実装例

import numpy as np
import matplotlib.pyplot as plt

# 微分方程式の右辺
def func(t, theta, **params):
	K = params.get("K", 0.0)
	omega = params.get("omega", np.zeros_like(theta))
	Ti, Tj = np.meshgrid(theta, theta)
	return omega + K * np.average(np.sin(Tj - Ti), axis=0)

# 4次ルンゲクッタ
def rk4(dt, t, Xs, func, **params):
	k1 = dt * func(t, Xs, **params)
	k2 = dt * func(t + dt / 2, Xs + k1 / 2, **params)
	k3 = dt * func(t + dt / 2, Xs + k2 / 2, **params)
	k4 = dt * func(t + dt, Xs + k3, **params)
	# 更新
	new_t = t + dt
	newXs = Xs + (k1 + 2 * k2 + 2 * k3 + k4) / 6
	return new_t, newXs

# 秩序変数の計算
def calc_order(dt, T, K, omega, theta0):
	print(f"K = {K:.2f}")

	# 初期値のコピー
	theta = theta0.copy()

	# 計算の実行
	rs = []
	for t in np.arange(0, T + dt, dt):
		# 複素数で取り扱う
		r = np.average(np.cos(theta) + np.sin(theta) * 1j)
		rs.append(r)

		# 更新
		_, theta = rk4(dt, t, theta, func, K=K, omega=omega)

	return np.array(rs)

# パラメータ
N = 100
T = 500
dt = 0.1

# 結合強度
Kmin, Kmax = 1.5, 4.0
Ks = np.arange(Kmin, Kmax, 0.2)

# 初期値
omega = np.array([np.random.standard_cauchy() for _ in range(N)])
theta0 = np.array([np.random.uniform(0.0, 2 * np.pi) for _ in range(N)])

# Kの値ごとに計算
Raves = []
for K in Ks:
	r = calc_order(dt, T, K, omega, theta0)
	Raves.append(np.average(np.abs(r)))

# 描画
fig, ax = plt.subplots()
ax.scatter(Ks, Raves, facecolor="None", edgecolors="red")

ax.set_xlim(Kmin, Kmax)
ax.set_ylim(0, 1)

plt.show()

2Like

Comments

@pythonnoob
Questioner
毎度ご回答いただきありがとうございます。
自分でもコードを書いてみるのですが、解析にかかる時間が大きくなってしまうのは計算と描画の処理を同時に行わせていることも関係しているのでしょうか？
正直、自分には不相応な難しいテーマだと痛感しておりますがとても楽しいのは事実です。
引き続きがんばります！
@bizzpaper
描画が計算時間に対してクリティカルな処理かは私も検証したことがないのでわかりません。少なくとも、計算中もax.scatterなどのオブジェクトを保持し続けるのはメモリに影響がありそうです。

そもそもですが、pythonは他の言語と比べて非常に遅いです。これは、pythonが1行ずつ解釈・翻訳しながら実行するインタプリタ型の言語であるためです。大規模な数値計算で用いられるCやFortranは実行前にコードをすべて解釈・翻訳（コンパイル）し、機械語に翻訳された実行用のファイルを別に生成します。この実行用ファイルは機械にわかりやすい内容なので、インタプリタ型と比べて非常に速く処理できます。ただ、pythonは可読性や書きやすさがウリであるため、一長一短です。

私は学生時代、計算はFortranで行い、解析とグラフ作成はpythonで行って、役割分担していました(Fortranのintelコンパイラでも1日〜2日かかる計算でした...)。

それぞれの一長一短を考えてツールを適切に選ぶことが望ましいですが、まずは計算と解析・グラフ作成のコードをわけるところから始めてみてはどうでしょうか？

>正直、自分には不相応な難しいテーマだと痛感しておりますがとても楽しいのは事実です。
>引き続きがんばります！

楽しいのはよいですね。モチベーションは大事です。その気持ちを忘れずにがんばってください。
@pythonnoob
Questioner
なるほど、用途や目的に合わせて適宜使用する言語を選択する必要があるのですね。
pcで演算するとき、とても時間がかかっていたので私だけではなかったと一安心です…

自分にもpcにも読みやすいコードにできるように努力します！
ありがとうございました。

@PondVillege posted at 2021-12-05

恐らくRunge-Kutta法を用いて数値解析を試みていると思われますが，まだKに依存していない

Rs=[]
for K in Ks:
    for t in ts:
        print(K)
        theta_now=theta[t, :].reshape(-1, 1)
        sumx=np.sum(np.cos(theta_now))
        sumy=np.sum(np.sin(theta_now))
        R=(sumx**2+sumy**2)**(1/2)/N
        Rs.append(R) # ここまで

の段階でRをRsに追加，その後Runge-Kutta法を用いてKの値に依存するような解析をするも，それを何ら記録せずにplt.scatter(K,np.sum(Rs)/len(Rs))としているためだと考えられます．

0Like

Comments

@PondVillege
Qiitaでは，ソースコードの掲載のみならず，数式を書くこともできます．今回のコードで解析することを目的としている関数等の記述もあればもっとアドバイスしやすいです．ご検討をお願いします．
数式の書き方は， Markdown記法チートシート https://qiita.com/PlanetMeron/items/63ac58898541cbe81ada 等を参考にしてください．
@pythonnoob
Questioner
こんばんは。ご回答いただきありがとうございます。
ただいま関数を追加させていただきました。蔵本モデルの解析です。
Rはすべてのθから求めており、全体のまとまり具合を表しています。
もし、よろしかったらご教授いただけると幸いです。
@PondVillege
もしかして @pythonnoob さんが過去に @bizzpaper さんに回答してもらった https://qiita.com/pythonnoob/questions/ab3f2af83d81d06efe08 をKで可変して実行しようとしている感じでしょうか...
@pythonnoob
Questioner
はい。そうです！
時間発展の部分をルンゲクッタを用いて書き直そうと思ったのですがKを変えていったときにこのようなエラーが出てしまいました。
@PondVillege
ブロードキャストの演算を用いようとして，多数の初期thetaが変動する様子を観察したいものと見受けましたが，reshapeの多用もあり，かなり追いにくいコードになっています．
行列演算をしているわけではないので，reshapeはあまり用いない方がいいと思います．

実際にthetaの値を追ってみたのですが，Kが1つずつ変化するごとに，thetaのサイズが300ずつ増加しています．初期に(1x3)のサイズで3個の状態を用意したのに，各K毎に(301x3)，(601x3)，(901x3)のように列のサイズが300ずつ増加してしまっているのは意図した実装ではないと思いますが，どうなのでしょうか．

すなわち，tを用いて各K毎に0番目からアクセスしてしまっていることが原因だと思われます．初回だけKの値を用いた解析を行うものの，次のKではまた0番目からthetaにアクセスし，データは301番目以降に追加していくようなものだと思います．
@pythonnoob
Questioner
夜分にすみません、只今確認してみたところ確かにthetaのサイズ変わってしまっていました。
これでは私の意図した処理になっていないのでもう一度見直してみます...
ご指摘感謝申し上げます。

@PondVillege posted at 2021-12-05

結論としては，
各K毎にthetaを初期化する．だと思います．
したがって，ループの開始部分

for K in Ks:
    for t in ts:

にthetaの初期化文を入れて

for K in Ks:
    theta = np.array([2.0 * np.pi * random.uniform(0,1) for i in range(N)]).reshape(1, -1)
    for t in ts:

とするか，大元のthetaを全てのKで共通にするために

theta_init = np.array([2.0 * np.pi * random.uniform(0,1) for i in range(N)]).reshape(1, -1)

としておき，

for K in Ks:
    theta = theta_init.copy()
    for t in ts:

とするかのどちらかだと思います．私的には後者を用いた方が良さそうだと思います．次の画像がそのときの結果です．randomを用いているので私と同じ結果にはなかなかならないと思います．参考までに．

0Like

Comments

@PondVillege
非常に細かくて申し訳ないのですが，書いていただいた数式において，変数たるθやωが斜体ではなく立体になっています．おそらく全角で入力したからだと考えられますが，\thetaや\omegaを用いると特に気にせず斜体になってくれると思います．また，関数であるsinが斜体になってしまっています．ここは，\sinとすれば立体になってくれると思います．　
今後レポートや論文を書くのだろうとお察ししますが，このような場で慣れておくといいと思います．
@pythonnoob
Questioner
表記の面でも気にしなければいけないことがあるのですね...大変勉強になりました。

現在、thetaの初期化をどうやって表現しようか悩んでいたところだったのでとてもありがたいです。
初学者かつ無知な者の質問で色々ご不便をおかけしてしまいましたが、丁寧にご回答してくださり本当に助かりました。大変感謝申し上げます。
@PondVillege
私としても，蔵本モデルについて初めて知ることができ，とてもいい機会になりました。ありがとうございます。勉学に励んでください。

Are you sure you want to delete the question?

for文について

解決したいこと

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

該当するソースコード

自分で試したこと

3Answer

Comments

Comments

Comments

Your answer might help someone💌