FPU单链模式能量

接上一篇,FPU单链模式能量模拟得到了单链各个原子不同时刻速度位移分布,根据这些数据计算一下每个时刻,各个模式的能量大小,得到FPU经典图

模式能量含义

上一篇文章记录了FPU模型模拟时的完整过程,最后得到的是不同时刻的每个原子的位移 xst 和速度分布 vst (矩阵形状 \([(steps+2) \times n]\)

现在计算不同时刻模式能量

\(E_n=\frac{m(v_n^2 + x_n^2\times w^2)}{2}\)

注意,

  • 计算的本征值本身就是 \(\omega ^ 2\),不要把本征值在进行平方了
  • 注意投影位移和速度时,本征向量行列别反了

测试

每个模式是否正确

画出来每个模式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# 第k个模式,从1开始,j原子数从0开始到32
def plt_bz(k, bzxl, bzz):
N = len(bzxls)

x = np.linspace(1, N, N)
print(x)
js = np.linspace(0, N, N + 1)
print(js)

file = "cache/imgs/bz/%d" % k
plt.scatter(x, bzxl)

plt.plot((2 / (N + 1))**0.5 * np.sin(np.pi * k * js / (N + 1)))

title = "%d=%f" % (k, bzz)
plt.title(title)
print(title)

plt.savefig(file)
plt.close()


def plt_bzs(bzzs, bzxls):

plt.plot(bzzs)
plt.savefig("cache/imgs/bz/bz")
plt.close()

pool = newPool()
N = len(bzzs)
pool.map(plt_bz, np.linspace(1, N, N), bzxls, bzzs)
pool.close()
pool.join()

初始位移是否正确

测试本征值和本征向量行列是否正确,直接画图初始位移投影

1
2
3
4
5
6
7
8
9
def ty_xst(tys, bzxls, log=True):
plt.figure()
if (log):
plt.axes(yscale="log")
print(tys)
ys = np.dot(tys, np.transpose(bzxls))
plt.scatter(np.linspace(1, len(tys), len(tys)), ys)
plt.savefig("cache/imgs/ty_xst")
plt.close()

上一篇文章,模拟计算以后,载入相关数据,测试一下

  • 防止本征向量行列反了
  • 防止模拟时,初始位移用错了
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
file_add = "_alpha"
# file_add = "_beta"
# file_add = ""
N = 32
step_start = 0
steps = 1000000
dt = 0.01

xst = np.load("cache/data/xst%s.npy" % file_add)[step_start:step_start + steps, :]
bzxls = np.load("data/bzxls.npy")

ty_xst(xst[0, :], bzxls, False)

# 第2个模式
# ty_xst(bzxls[2], bzxls, False)

生成的图片,应该是,测试的模式再上面,其他的模式几乎为0

1

某个原子位移随时间变化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def plot_n(xst, vst, n_see=16):
legends = []
plt.ylim(-1.5, 1.5)
for i in [0]:
s = "-."
if (i == 0):
s = "-"
plt.plot(xst[:, n_see + i], s)
plt.plot(vst[:, n_see + i], s)
legends.append("n=%d" % (n_see + i))

plt.legend(legends)

# plt.plot(vst[:, 1])
plt.savefig("test/nt")
# plt.show()
plt.close()
1

模式能量

获取模式能量

1
2
3
4
5
6
7
8
9
10
11
def get_Ets(xst, vst, bzzs, bzxls):
bzxls = np.transpose(bzxls)

xst_ = np.dot(xst, bzxls)
vst_ = np.dot(vst, bzxls)

Ets = vst_ ** 2 + xst_**2 * np.abs(bzzs)

np.save("cache/data/ets.npy", Ets)

return Ets

画出某个时间,各个模式能量投影情况

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def plot_ek(Ets, t=1, onlyOdd=True, log=True):
es = Ets[t, :]
N = len(es)

xs = np.linspace(1, N, N)

if(onlyOdd):
es = es[::2]
xs = xs[::2]

print(xs)

plt.figure()
if(log):
plt.axes(yscale="log")
plt.scatter(xs, es)
for i in range(len(es)):
plt.annotate(str(int(xs[i])), xy=(xs[i], es[i]), xytext=(xs[i], es[i]))
plt.savefig("test/eks")
plt.close()

画出某些模式的能量随时间变化,注意啊,千万不能log坐标下画出来

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def plot_et(Es, ks=[0, 1, 2, 3, 4, 5], log=False):

Es = np.array(Es)

print(np.shape(Es))

plt.figure()
if (log):
plt.axes(yscale="log")

Es_ = []
legend = []

for k in ks:
k = k
print(k)
legend.append(k)
Es_.append(Es[:, k])

# 归一化
# Es_ = Es_ / Es[:, 0]

ts = np.arange(len(Es[:, 0])) * t_

for j in range(len(ks)):
y = Es_[j]
plt.plot(ts, y)

plt.legend(legend)
# plt.savefig("%s/imgs/cals/cal%d_ty.png" % (path_run, cal))
plt.savefig("test/ets")
plt.close()
# plt.show()
1

完整代码

提前创建文件夹,没有写在python,上一篇模型时已经创建了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
import numpy as np
from matplotlib import pyplot as plt
from pathos.multiprocessing import ProcessingPool as newPool


def test_bz(bzzs, bzxls, V):
print("测试")
bzz = bzzs[0]
bzxl = bzxls[0]
test_ = (bzz * bzxl) - np.dot(bzxl, V)
print(test_)

# 通过投影,也可以查看是否对
# ty_xst(xst[0, :], bzxls, False)
# ty_xst(bzxls[2], bzxls, False)

if (sum(test_) < 10**-8):
print("本征值没问题")


def plot_n(xst, vst, n_see=16):

legends = []
plt.ylim(-1.5, 1.5)
for i in [0]:
s = "-."
if (i == 0):
s = "-"
plt.plot(xst[:, n_see + i], s)
# plt.plot(vst[:, n_see + i], s)
legends.append("n=%d" % (n_see + i))

plt.legend(legends)

# plt.plot(vst[:, 1])
plt.savefig("test/nt")
# plt.show()
plt.close()


def plot_t(xst, jump=100):
for xs in xst[::jump]:
plt.plot(xs)
plt.savefig("test/xst")
plt.close()

# plt.show()


def get_Ets(xst, vst, bzzs, bzxls):
bzxls = np.transpose(bzxls)

xst_ = np.dot(xst, bzxls)
vst_ = np.dot(vst, bzxls)

Ets = vst_ ** 2 + xst_**2 * np.abs(bzzs)

np.save("cache/data/ets.npy", Ets)

return Ets


def ty_xst(tys, bzxls, log=True):
plt.figure()
if (log):
plt.axes(yscale="log")
print(tys)
ys = np.dot(tys, np.transpose(bzxls))
plt.scatter(np.linspace(1, len(tys), len(tys)), ys)
plt.savefig("cache/imgs/ty_xst")
plt.close()


def plot_ek(Ets, t=1, onlyOdd=True, log=True):
es = Ets[t, :]
N = len(es)

xs = np.linspace(1, N, N)

if(onlyOdd):
es = es[::2]
xs = xs[::2]

print(xs)

plt.figure()
if(log):
plt.axes(yscale="log")
plt.scatter(xs, es)
for i in range(len(es)):
plt.annotate(str(int(xs[i])), xy=(xs[i], es[i]), xytext=(xs[i], es[i]))
plt.savefig("test/eks")
plt.close()


def plot_et(Es, ks=[0, 1, 2, 3, 4, 5], log=False):

Es = np.array(Es)

print(np.shape(Es))

plt.figure()
if (log):
plt.axes(yscale="log")

Es_ = []
legend = []

for k in ks:
k = k
print(k)
legend.append(k)
Es_.append(Es[:, k])

# 归一化
# Es_ = Es_ / Es[:, 0]

ts = np.arange(len(Es[:, 0])) * t_

for j in range(len(ks)):
y = Es_[j]
plt.plot(ts, y)

plt.legend(legend)
# plt.savefig("%s/imgs/cals/cal%d_ty.png" % (path_run, cal))
plt.savefig("test/ets")
plt.close()
# plt.show()


# 第k个模式,从1开始,j原子数从0开始到32
def plt_bz(k, bzxl, bzz):
N = len(bzxls)

x = np.linspace(1, N, N)
print(x)
js = np.linspace(0, N, N + 1)
print(js)

file = "cache/imgs/bz/%d" % k
plt.scatter(x, bzxl)

plt.plot((2 / (N + 1))**0.5 * np.sin(np.pi * k * js / (N + 1)))

title = "%d=%f" % (k, bzz)
plt.title(title)
print(title)

plt.savefig(file)
plt.close()


def plt_bzs(bzzs, bzxls):

plt.plot(bzzs)
plt.savefig("cache/imgs/bz/bz")
plt.close()

pool = newPool()
N = len(bzzs)
pool.map(plt_bz, np.linspace(1, N, N), bzxls, bzzs)
pool.close()
pool.join()


file_add = "_alpha"
# file_add = "_beta"
# file_add = ""
N = 32
step_start = 0
steps = 1000000
t_ = 0.01

xst = np.load("cache/data/xst%s.npy" % file_add)[step_start:step_start +
steps, :]
vst = np.load("cache/data/vst%s.npy" % file_add)[step_start:step_start +
steps, :]
print(np.shape(xst))
bzzs = np.load("data/bzzs.npy")
bzxls = np.load("data/bzxls.npy")
V = np.load("data/V.npy")

test_bz(bzzs, bzxls, V)
# plt_bzs(bzzs, bzxls)

# plot_t(xst[0:20000], jump=100)
# plot_n(xst[0:200000, :], vst[0:200000, :], 15)

Ets = get_Ets(xst, vst, bzzs, bzxls)
print(np.shape(Ets))
plot_et(Ets)
# plot_ek(Ets, 100000-1, log=True)


# ty_xst(xst[0, :], bzxls, False)
# ty_xst(bzxls[2], bzxls, False)

\(\beta\) 模型代码不变,注意一下,文件名就行


本文作者:yuhldr
本文地址: https://yuhldr.github.io/posts/46949.html
版权声明:转载请注明出处!