category
This commit is contained in:
@@ -0,0 +1,114 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
# os.chdir('E:/data') # 设置路径
|
||||
|
||||
# 万有引力常数
|
||||
G = 1
|
||||
|
||||
# 三体的质量
|
||||
m1 = 10 # 绿色
|
||||
m2 = 1000 # 红色
|
||||
m3 = 10 # 蓝色
|
||||
|
||||
# 三体的初始位置
|
||||
x1 = 0 # 绿色
|
||||
y1 = 500
|
||||
x2 = 0 # 红色
|
||||
y2 = 0
|
||||
x3 = 0 # 蓝色
|
||||
y3 = 1000
|
||||
|
||||
# 三体的初始速度
|
||||
v1_x = 1.5 # 绿色
|
||||
v1_y = 0
|
||||
v2_x = 0 # 红色
|
||||
v2_y = 0
|
||||
v3_x = 0.8 # 蓝色
|
||||
v3_y = 0
|
||||
|
||||
# 步长
|
||||
t = 0.1
|
||||
|
||||
plt.ion() # 开启交互模式
|
||||
observation_max = 100 # 视线范围初始值
|
||||
x1_all = [x1] # 轨迹初始值
|
||||
y1_all = [y1]
|
||||
x2_all = [x2]
|
||||
y2_all = [y2]
|
||||
x3_all = [x3]
|
||||
y3_all = [y3]
|
||||
|
||||
i0 = 0
|
||||
for i in range(1000000):
|
||||
distance12 = np.sqrt((x1-x2)**2+(y1-y2)**2) # 物体1和物体2之间的距离
|
||||
distance13 = np.sqrt((x1-x3)**2+(y1-y3)**2) # 物体1和物体3之间的距离
|
||||
distance23 = np.sqrt((x2-x3)**2+(y2-y3)**2) # 物体2和物体3之间的距离
|
||||
|
||||
# 对物体1的计算
|
||||
a1_2 = G*m2/(distance12**2) # 物体2对物体1的加速度（用上万有引力公式）
|
||||
a1_3 = G*m3/(distance13**2) # 物体3对物体1的加速度
|
||||
a1_x = a1_2*(x2-x1)/distance12 + a1_3*(x3-x1)/distance13 # 物体1受到的水平加速度
|
||||
a1_y = a1_2*(y2-y1)/distance12 + a1_3*(y3-y1)/distance13 # 物体1受到的垂直加速度
|
||||
v1_x = v1_x + a1_x*t # 物体1的速度
|
||||
v1_y = v1_y + a1_y*t # 物体1的速度
|
||||
x1 = x1 + v1_x*t # 物体1的水平位置
|
||||
y1 = y1 + v1_y*t # 物体1的垂直位置
|
||||
x1_all = np.append(x1_all, x1) # 记录轨迹
|
||||
y1_all = np.append(y1_all, y1) # 记录轨迹
|
||||
|
||||
# 对物体2的计算
|
||||
a2_1 = G*m1/(distance12**2)
|
||||
a2_3 = G*m3/(distance23**2)
|
||||
a2_x = a2_1*(x1-x2)/distance12 + a2_3*(x3-x2)/distance23
|
||||
a2_y = a2_1*(y1-y2)/distance12 + a2_3*(y3-y2)/distance23
|
||||
v2_x = v2_x + a2_x*t
|
||||
v2_y = v2_y + a2_y*t
|
||||
x2 = x2 + v2_x*t
|
||||
y2 = y2 + v2_y*t
|
||||
x2_all = np.append(x2_all, x2)
|
||||
y2_all = np.append(y2_all, y2)
|
||||
|
||||
# 对物体3的计算
|
||||
a3_1 = G*m1/(distance13**2)
|
||||
a3_2 = G*m2/(distance23**2)
|
||||
a3_x = a3_1*(x1-x3)/distance13 + a3_2*(x2-x3)/distance23
|
||||
a3_y = a3_1*(y1-y3)/distance13 + a3_2*(y2-y3)/distance23
|
||||
v3_x = v3_x + a3_x*t
|
||||
v3_y = v3_y + a3_y*t
|
||||
x3 = x3 + v3_x*t
|
||||
y3 = y3 + v3_y*t
|
||||
x3_all = np.append(x3_all, x3)
|
||||
y3_all = np.append(y3_all, y3)
|
||||
|
||||
# 选择观测坐标
|
||||
axis_x = np.mean([x1, x2, x3]) # 观测坐标中心固定在平均值的地方
|
||||
axis_y = np.mean([y1, y2, y3]) # 观测坐标中心固定在平均值的地方
|
||||
while True:
|
||||
if np.abs(x1-axis_x) > observation_max or np.abs(x2-axis_x) > observation_max or np.abs(x3-axis_x) > observation_max or\
|
||||
np.abs(y1-axis_y) > observation_max or np.abs(y2-axis_y) > observation_max or np.abs(y3-axis_y) > observation_max:
|
||||
observation_max = observation_max * 2 # 有一个物体超出视线时，视线范围翻倍
|
||||
elif np.abs(x1-axis_x) < observation_max/10 and np.abs(x2-axis_x) < observation_max/10 and np.abs(x3-axis_x) < observation_max/10 and\
|
||||
np.abs(y1-axis_y) < observation_max/10 and np.abs(y2-axis_y) < observation_max/10 and np.abs(y3-axis_y) < observation_max/10:
|
||||
observation_max = observation_max / 2 # 所有物体都在的视线的10分之一内，视线范围减半
|
||||
else:
|
||||
break
|
||||
|
||||
plt.axis([axis_x-observation_max, axis_x+observation_max, axis_y-observation_max, axis_y+observation_max])
|
||||
|
||||
plt.plot(x1, y1, 'og', markersize=m1*100/observation_max) # 默认密度相同，质量越大的，球面积越大。视线范围越宽，球看起来越小。
|
||||
plt.plot(x2, y2, 'or', markersize=m2*100/observation_max)
|
||||
plt.plot(x3, y3, 'ob', markersize=m3*100/observation_max)
|
||||
plt.plot(x1_all, y1_all, '-g') # 画轨迹
|
||||
plt.plot(x2_all, y2_all, '-r')
|
||||
plt.plot(x3_all, y3_all, '-b')
|
||||
|
||||
|
||||
# plt.show() # 显示图像
|
||||
# plt.pause(0.00001) # 暂停0.00001，防止画图过快
|
||||
|
||||
if np.mod(i, 100) == 0: # 当运动明显时，把图画出来
|
||||
print(i0)
|
||||
plt.savefig(str(i0)+'.jpg') # 保存为图片可以用来做动画
|
||||
i0 += 1
|
||||
plt.clf() # 清空
|
||||
@@ -0,0 +1,110 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
# os.chdir('E:/data') # 设置路径
|
||||
|
||||
# 万有引力常数
|
||||
G = 1
|
||||
|
||||
# 三体的质量
|
||||
m1 = 3 # 绿色
|
||||
m2 = 100 # 红色
|
||||
m3 = 10 # 蓝色
|
||||
|
||||
# 三体的初始位置
|
||||
x1 = 0 # 绿色
|
||||
y1 = -100
|
||||
x2 = 0 # 红色
|
||||
y2 = 0
|
||||
x3 = 0 # 蓝色
|
||||
y3 = 50
|
||||
|
||||
# 三体的初始速度
|
||||
v1_x = 1 # 绿色
|
||||
v1_y = 0
|
||||
v2_x = 0 # 红色
|
||||
v2_y = 0
|
||||
v3_x = 2 # 蓝色
|
||||
v3_y = 0
|
||||
|
||||
# 步长
|
||||
t = 0.05
|
||||
|
||||
plt.ion() # 开启交互模式
|
||||
observation_max = 100 # 观测坐标范围初始值
|
||||
x1_all = [x1] # 轨迹初始值
|
||||
y1_all = [y1]
|
||||
x2_all = [x2]
|
||||
y2_all = [y2]
|
||||
x3_all = [x3]
|
||||
y3_all = [y3]
|
||||
|
||||
i0 = 0
|
||||
for i in range(100000000000):
|
||||
distance12 = np.sqrt((x1-x2)**2+(y1-y2)**2) # 物体1和物体2之间的距离
|
||||
distance13 = np.sqrt((x1-x3)**2+(y1-y3)**2) # 物体1和物体3之间的距离
|
||||
distance23 = np.sqrt((x2-x3)**2+(y2-y3)**2) # 物体2和物体3之间的距离
|
||||
# 对物体1的计算
|
||||
a1_2 = G*m2/(distance12**2) # 物体2对物体1的加速度（用上万有引力公式）
|
||||
a1_3 = G*m3/(distance13**2) # 物体3对物体1的加速度
|
||||
a1_x = a1_2*(x2-x1)/distance12 + a1_3*(x3-x1)/distance13 # 物体1受到的水平加速度
|
||||
a1_y = a1_2*(y2-y1)/distance12 + a1_3*(y3-y1)/distance13 # 物体1受到的垂直加速度
|
||||
v1_x = v1_x + a1_x*t # 物体1的速度
|
||||
v1_y = v1_y + a1_y*t # 物体1的速度
|
||||
x1 = x1 + v1_x*t # 物体1的水平位置
|
||||
y1 = y1 + v1_y*t # 物体1的垂直位置
|
||||
x1_all = np.append(x1_all, x1) # 记录轨迹
|
||||
y1_all = np.append(y1_all, y1) # 记录轨迹
|
||||
# 对物体2的计算
|
||||
a2_1 = G*m1/(distance12**2)
|
||||
a2_3 = G*m3/(distance23**2)
|
||||
a2_x = a2_1*(x1-x2)/distance12 + a2_3*(x3-x2)/distance23
|
||||
a2_y = a2_1*(y1-y2)/distance12 + a2_3*(y3-y2)/distance23
|
||||
v2_x = v2_x + a2_x*t
|
||||
v2_y = v2_y + a2_y*t
|
||||
x2 = x2 + v2_x*t
|
||||
y2 = y2 + v2_y*t
|
||||
x2_all = np.append(x2_all, x2)
|
||||
y2_all = np.append(y2_all, y2)
|
||||
# 对物体3的计算
|
||||
a3_1 = G*m1/(distance13**2)
|
||||
a3_2 = G*m2/(distance23**2)
|
||||
a3_x = a3_1*(x1-x3)/distance13 + a3_2*(x2-x3)/distance23
|
||||
a3_y = a3_1*(y1-y3)/distance13 + a3_2*(y2-y3)/distance23
|
||||
v3_x = v3_x + a3_x*t
|
||||
v3_y = v3_y + a3_y*t
|
||||
x3 = x3 + v3_x*t
|
||||
y3 = y3 + v3_y*t
|
||||
x3_all = np.append(x3_all, x3)
|
||||
y3_all = np.append(y3_all, y3)
|
||||
|
||||
# 选择观测坐标
|
||||
axis_x = np.mean([x1, x2, x3]) # 观测坐标中心固定在平均值的地方
|
||||
axis_y = np.mean([y1, y2, y3]) # 观测坐标中心固定在平均值的地方
|
||||
while True:
|
||||
if np.abs(x1-axis_x) > observation_max or np.abs(x2-axis_x) > observation_max or np.abs(x3-axis_x) > observation_max or\
|
||||
np.abs(y1-axis_y) > observation_max or np.abs(y2-axis_y) > observation_max or np.abs(y3-axis_y) > observation_max:
|
||||
observation_max = observation_max * 2 # 有一个物体超出视线时，视线范围翻倍
|
||||
elif np.abs(x1-axis_x) < observation_max/10 and np.abs(x2-axis_x) < observation_max/10 and np.abs(x3-axis_x) < observation_max/10 and\
|
||||
np.abs(y1-axis_y) < observation_max/10 and np.abs(y2-axis_y) < observation_max/10 and np.abs(y3-axis_y) < observation_max/10:
|
||||
observation_max = observation_max / 2 # 所有物体都在的视线的10分之一内，视线范围减半
|
||||
else:
|
||||
break
|
||||
|
||||
plt.axis([axis_x-observation_max, axis_x+observation_max, axis_y-observation_max, axis_y+observation_max])
|
||||
|
||||
plt.plot(x1, y1, 'og', markersize=m1*100/observation_max) # 默认密度相同，质量越大的，画出来的面积越大。视线范围越宽，球看起来越小。
|
||||
plt.plot(x2, y2, 'or', markersize=m2*100/observation_max)
|
||||
plt.plot(x3, y3, 'ob', markersize=m3*100/observation_max)
|
||||
plt.plot(x1_all, y1_all, '-g') # 画轨迹
|
||||
plt.plot(x2_all, y2_all, '-r')
|
||||
plt.plot(x3_all, y3_all, '-b')
|
||||
|
||||
# plt.show() # 显示图像
|
||||
# plt.pause(0.00001) # 暂停0.00001，防止画图过快
|
||||
|
||||
if np.mod(i, 200) == 0:
|
||||
print(i0)
|
||||
plt.savefig(str(i0)+'.jpg') # 保存为图片可以用来做动画
|
||||
i0 += 1
|
||||
plt.clf() # 清空
|
||||
@@ -0,0 +1,110 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
# os.chdir('E:/data') # 设置路径
|
||||
|
||||
# 万有引力常数
|
||||
G = 1
|
||||
|
||||
# 三体的质量
|
||||
m1 = 3 # 绿色
|
||||
m2 = 100 # 红色
|
||||
m3 = 10 # 蓝色
|
||||
|
||||
# 三体的初始位置
|
||||
x1 = 0 # 绿色
|
||||
y1 = -100
|
||||
x2 = 0 # 红色
|
||||
y2 = 0
|
||||
x3 = 0 # 蓝色
|
||||
y3 = 50
|
||||
|
||||
# 三体的初始速度
|
||||
v1_x = 1 # 绿色
|
||||
v1_y = 0
|
||||
v2_x = 0 # 红色
|
||||
v2_y = 0
|
||||
v3_x = 2 # 蓝色
|
||||
v3_y = 0
|
||||
|
||||
# 步长
|
||||
t = 0.1
|
||||
|
||||
plt.ion() # 开启交互模式
|
||||
observation_max = 100 # 观测坐标范围初始值
|
||||
x1_all = [x1] # 轨迹初始值
|
||||
y1_all = [y1]
|
||||
x2_all = [x2]
|
||||
y2_all = [y2]
|
||||
x3_all = [x3]
|
||||
y3_all = [y3]
|
||||
|
||||
i0 = 0
|
||||
for i in range(100000000000):
|
||||
distance12 = np.sqrt((x1-x2)**2+(y1-y2)**2) # 物体1和物体2之间的距离
|
||||
distance13 = np.sqrt((x1-x3)**2+(y1-y3)**2) # 物体1和物体3之间的距离
|
||||
distance23 = np.sqrt((x2-x3)**2+(y2-y3)**2) # 物体2和物体3之间的距离
|
||||
# 对物体1的计算
|
||||
a1_2 = G*m2/(distance12**2) # 物体2对物体1的加速度（用上万有引力公式）
|
||||
a1_3 = G*m3/(distance13**2) # 物体3对物体1的加速度
|
||||
a1_x = a1_2*(x2-x1)/distance12 + a1_3*(x3-x1)/distance13 # 物体1受到的水平加速度
|
||||
a1_y = a1_2*(y2-y1)/distance12 + a1_3*(y3-y1)/distance13 # 物体1受到的垂直加速度
|
||||
v1_x = v1_x + a1_x*t # 物体1的速度
|
||||
v1_y = v1_y + a1_y*t # 物体1的速度
|
||||
x1 = x1 + v1_x*t # 物体1的水平位置
|
||||
y1 = y1 + v1_y*t # 物体1的垂直位置
|
||||
x1_all = np.append(x1_all, x1) # 记录轨迹
|
||||
y1_all = np.append(y1_all, y1) # 记录轨迹
|
||||
# 对物体2的计算
|
||||
a2_1 = G*m1/(distance12**2)
|
||||
a2_3 = G*m3/(distance23**2)
|
||||
a2_x = a2_1*(x1-x2)/distance12 + a2_3*(x3-x2)/distance23
|
||||
a2_y = a2_1*(y1-y2)/distance12 + a2_3*(y3-y2)/distance23
|
||||
v2_x = v2_x + a2_x*t
|
||||
v2_y = v2_y + a2_y*t
|
||||
x2 = x2 + v2_x*t
|
||||
y2 = y2 + v2_y*t
|
||||
x2_all = np.append(x2_all, x2)
|
||||
y2_all = np.append(y2_all, y2)
|
||||
# 对物体3的计算
|
||||
a3_1 = G*m1/(distance13**2)
|
||||
a3_2 = G*m2/(distance23**2)
|
||||
a3_x = a3_1*(x1-x3)/distance13 + a3_2*(x2-x3)/distance23
|
||||
a3_y = a3_1*(y1-y3)/distance13 + a3_2*(y2-y3)/distance23
|
||||
v3_x = v3_x + a3_x*t
|
||||
v3_y = v3_y + a3_y*t
|
||||
x3 = x3 + v3_x*t
|
||||
y3 = y3 + v3_y*t
|
||||
x3_all = np.append(x3_all, x3)
|
||||
y3_all = np.append(y3_all, y3)
|
||||
|
||||
# 选择观测坐标
|
||||
axis_x = np.mean([x1, x2, x3]) # 观测坐标中心固定在平均值的地方
|
||||
axis_y = np.mean([y1, y2, y3]) # 观测坐标中心固定在平均值的地方
|
||||
while True:
|
||||
if np.abs(x1-axis_x) > observation_max or np.abs(x2-axis_x) > observation_max or np.abs(x3-axis_x) > observation_max or\
|
||||
np.abs(y1-axis_y) > observation_max or np.abs(y2-axis_y) > observation_max or np.abs(y3-axis_y) > observation_max:
|
||||
observation_max = observation_max * 2 # 有一个物体超出视线时，视线范围翻倍
|
||||
elif np.abs(x1-axis_x) < observation_max/10 and np.abs(x2-axis_x) < observation_max/10 and np.abs(x3-axis_x) < observation_max/10 and\
|
||||
np.abs(y1-axis_y) < observation_max/10 and np.abs(y2-axis_y) < observation_max/10 and np.abs(y3-axis_y) < observation_max/10:
|
||||
observation_max = observation_max / 2 # 所有物体都在的视线的10分之一内，视线范围减半
|
||||
else:
|
||||
break
|
||||
|
||||
plt.axis([axis_x-observation_max, axis_x+observation_max, axis_y-observation_max, axis_y+observation_max])
|
||||
|
||||
plt.plot(x1, y1, 'og', markersize=m1*100/observation_max) # 默认密度相同，质量越大的，画出来的面积越大。视线范围越宽，球看起来越小。
|
||||
plt.plot(x2, y2, 'or', markersize=m2*100/observation_max)
|
||||
plt.plot(x3, y3, 'ob', markersize=m3*100/observation_max)
|
||||
plt.plot(x1_all, y1_all, '-g') # 画轨迹
|
||||
plt.plot(x2_all, y2_all, '-r')
|
||||
plt.plot(x3_all, y3_all, '-b')
|
||||
|
||||
# plt.show() # 显示图像
|
||||
# plt.pause(0.00001) # 暂停0.00001，防止画图过快
|
||||
|
||||
if np.mod(i, 100) == 0:
|
||||
print(i0)
|
||||
plt.savefig(str(i0)+'.jpg') # 保存为图片可以用来做动画
|
||||
i0 += 1
|
||||
plt.clf() # 清空
|
||||
@@ -0,0 +1,113 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
# os.chdir('E:/data') # 设置路径
|
||||
|
||||
# 万有引力常数
|
||||
G = 1
|
||||
|
||||
# 三体的质量
|
||||
m1 = 15 # 绿色
|
||||
m2 = 12 # 红色
|
||||
m3 = 8 # 蓝色
|
||||
|
||||
# 三体的初始位置
|
||||
x1 = 300 # 绿色
|
||||
y1 = 50
|
||||
x2 = -100 # 红色
|
||||
y2 = -200
|
||||
x3 = -100 # 蓝色
|
||||
y3 = 150
|
||||
|
||||
# 三体的初始速度
|
||||
v1_x = 0 # 绿色
|
||||
v1_y = 0
|
||||
v2_x = 0 # 红色
|
||||
v2_y = 0
|
||||
v3_x = 0 # 蓝色
|
||||
v3_y = 0
|
||||
|
||||
# 步长
|
||||
t = 0.05
|
||||
|
||||
plt.ion() # 开启交互模式
|
||||
observation_max = 100 # 视线范围初始值
|
||||
x1_all = [x1] # 轨迹初始值
|
||||
y1_all = [y1]
|
||||
x2_all = [x2]
|
||||
y2_all = [y2]
|
||||
x3_all = [x3]
|
||||
y3_all = [y3]
|
||||
|
||||
i0 = 0
|
||||
for i in range(1000000):
|
||||
distance12 = np.sqrt((x1-x2)**2+(y1-y2)**2) # 物体1和物体2之间的距离
|
||||
distance13 = np.sqrt((x1-x3)**2+(y1-y3)**2) # 物体1和物体3之间的距离
|
||||
distance23 = np.sqrt((x2-x3)**2+(y2-y3)**2) # 物体2和物体3之间的距离
|
||||
|
||||
# 对物体1的计算
|
||||
a1_2 = G*m2/(distance12**2) # 物体2对物体1的加速度（用上万有引力公式）
|
||||
a1_3 = G*m3/(distance13**2) # 物体3对物体1的加速度
|
||||
a1_x = a1_2*(x2-x1)/distance12 + a1_3*(x3-x1)/distance13 # 物体1受到的水平加速度
|
||||
a1_y = a1_2*(y2-y1)/distance12 + a1_3*(y3-y1)/distance13 # 物体1受到的垂直加速度
|
||||
v1_x = v1_x + a1_x*t # 物体1的速度
|
||||
v1_y = v1_y + a1_y*t # 物体1的速度
|
||||
x1 = x1 + v1_x*t # 物体1的水平位置
|
||||
y1 = y1 + v1_y*t # 物体1的垂直位置
|
||||
x1_all = np.append(x1_all, x1) # 记录轨迹
|
||||
y1_all = np.append(y1_all, y1) # 记录轨迹
|
||||
|
||||
# 对物体2的计算
|
||||
a2_1 = G*m1/(distance12**2)
|
||||
a2_3 = G*m3/(distance23**2)
|
||||
a2_x = a2_1*(x1-x2)/distance12 + a2_3*(x3-x2)/distance23
|
||||
a2_y = a2_1*(y1-y2)/distance12 + a2_3*(y3-y2)/distance23
|
||||
v2_x = v2_x + a2_x*t
|
||||
v2_y = v2_y + a2_y*t
|
||||
x2 = x2 + v2_x*t
|
||||
y2 = y2 + v2_y*t
|
||||
x2_all = np.append(x2_all, x2)
|
||||
y2_all = np.append(y2_all, y2)
|
||||
|
||||
# 对物体3的计算
|
||||
a3_1 = G*m1/(distance13**2)
|
||||
a3_2 = G*m2/(distance23**2)
|
||||
a3_x = a3_1*(x1-x3)/distance13 + a3_2*(x2-x3)/distance23
|
||||
a3_y = a3_1*(y1-y3)/distance13 + a3_2*(y2-y3)/distance23
|
||||
v3_x = v3_x + a3_x*t
|
||||
v3_y = v3_y + a3_y*t
|
||||
x3 = x3 + v3_x*t
|
||||
y3 = y3 + v3_y*t
|
||||
x3_all = np.append(x3_all, x3)
|
||||
y3_all = np.append(y3_all, y3)
|
||||
|
||||
# 选择观测坐标
|
||||
axis_x = np.mean([x1, x2, x3]) # 观测坐标中心固定在平均值的地方
|
||||
axis_y = np.mean([y1, y2, y3]) # 观测坐标中心固定在平均值的地方
|
||||
while True:
|
||||
if np.abs(x1-axis_x) > observation_max or np.abs(x2-axis_x) > observation_max or np.abs(x3-axis_x) > observation_max or\
|
||||
np.abs(y1-axis_y) > observation_max or np.abs(y2-axis_y) > observation_max or np.abs(y3-axis_y) > observation_max:
|
||||
observation_max = observation_max * 2 # 有一个物体超出视线时，视线范围翻倍
|
||||
elif np.abs(x1-axis_x) < observation_max/10 and np.abs(x2-axis_x) < observation_max/10 and np.abs(x3-axis_x) < observation_max/10 and\
|
||||
np.abs(y1-axis_y) < observation_max/10 and np.abs(y2-axis_y) < observation_max/10 and np.abs(y3-axis_y) < observation_max/10:
|
||||
observation_max = observation_max / 2 # 所有物体都在的视线的10分之一内，视线范围减半
|
||||
else:
|
||||
break
|
||||
|
||||
plt.axis([axis_x-observation_max, axis_x+observation_max, axis_y-observation_max, axis_y+observation_max])
|
||||
|
||||
plt.plot(x1, y1, 'og', markersize=m1*100/observation_max) # 默认密度相同，质量越大的，球面积越大。视线范围越宽，球看起来越小。
|
||||
plt.plot(x2, y2, 'or', markersize=m2*100/observation_max)
|
||||
plt.plot(x3, y3, 'ob', markersize=m3*100/observation_max)
|
||||
plt.plot(x1_all, y1_all, '-g') # 画轨迹
|
||||
plt.plot(x2_all, y2_all, '-r')
|
||||
plt.plot(x3_all, y3_all, '-b')
|
||||
|
||||
# plt.show() # 显示图像
|
||||
# plt.pause(0.00001) # 暂停0.00001，防止画图过快
|
||||
|
||||
if np.mod(i, 1000) == 0: # 当运动明显时，把图画出来
|
||||
print(i0)
|
||||
plt.savefig(str(i0)+'.jpg') # 保存为图片可以用来做动画
|
||||
i0 += 1
|
||||
plt.clf() # 清空
|
||||
@@ -0,0 +1,113 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
# os.chdir('E:/data') # 设置路径
|
||||
|
||||
# 万有引力常数
|
||||
G = 1
|
||||
|
||||
# 三体的质量
|
||||
m1 = 15 # 绿色
|
||||
m2 = 12 # 红色
|
||||
m3 = 8 # 蓝色
|
||||
|
||||
# 三体的初始位置
|
||||
x1 = 300 # 绿色
|
||||
y1 = 50
|
||||
x2 = -100 # 红色
|
||||
y2 = -200
|
||||
x3 = -100 # 蓝色
|
||||
y3 = 150
|
||||
|
||||
# 三体的初始速度
|
||||
v1_x = 0 # 绿色
|
||||
v1_y = 0
|
||||
v2_x = 0 # 红色
|
||||
v2_y = 0
|
||||
v3_x = 0 # 蓝色
|
||||
v3_y = 0
|
||||
|
||||
# 步长
|
||||
t = 0.1
|
||||
|
||||
plt.ion() # 开启交互模式
|
||||
observation_max = 100 # 视线范围初始值
|
||||
x1_all = [x1] # 轨迹初始值
|
||||
y1_all = [y1]
|
||||
x2_all = [x2]
|
||||
y2_all = [y2]
|
||||
x3_all = [x3]
|
||||
y3_all = [y3]
|
||||
|
||||
i0 = 0
|
||||
for i in range(1000000):
|
||||
distance12 = np.sqrt((x1-x2)**2+(y1-y2)**2) # 物体1和物体2之间的距离
|
||||
distance13 = np.sqrt((x1-x3)**2+(y1-y3)**2) # 物体1和物体3之间的距离
|
||||
distance23 = np.sqrt((x2-x3)**2+(y2-y3)**2) # 物体2和物体3之间的距离
|
||||
|
||||
# 对物体1的计算
|
||||
a1_2 = G*m2/(distance12**2) # 物体2对物体1的加速度（用上万有引力公式）
|
||||
a1_3 = G*m3/(distance13**2) # 物体3对物体1的加速度
|
||||
a1_x = a1_2*(x2-x1)/distance12 + a1_3*(x3-x1)/distance13 # 物体1受到的水平加速度
|
||||
a1_y = a1_2*(y2-y1)/distance12 + a1_3*(y3-y1)/distance13 # 物体1受到的垂直加速度
|
||||
v1_x = v1_x + a1_x*t # 物体1的速度
|
||||
v1_y = v1_y + a1_y*t # 物体1的速度
|
||||
x1 = x1 + v1_x*t # 物体1的水平位置
|
||||
y1 = y1 + v1_y*t # 物体1的垂直位置
|
||||
x1_all = np.append(x1_all, x1) # 记录轨迹
|
||||
y1_all = np.append(y1_all, y1) # 记录轨迹
|
||||
|
||||
# 对物体2的计算
|
||||
a2_1 = G*m1/(distance12**2)
|
||||
a2_3 = G*m3/(distance23**2)
|
||||
a2_x = a2_1*(x1-x2)/distance12 + a2_3*(x3-x2)/distance23
|
||||
a2_y = a2_1*(y1-y2)/distance12 + a2_3*(y3-y2)/distance23
|
||||
v2_x = v2_x + a2_x*t
|
||||
v2_y = v2_y + a2_y*t
|
||||
x2 = x2 + v2_x*t
|
||||
y2 = y2 + v2_y*t
|
||||
x2_all = np.append(x2_all, x2)
|
||||
y2_all = np.append(y2_all, y2)
|
||||
|
||||
# 对物体3的计算
|
||||
a3_1 = G*m1/(distance13**2)
|
||||
a3_2 = G*m2/(distance23**2)
|
||||
a3_x = a3_1*(x1-x3)/distance13 + a3_2*(x2-x3)/distance23
|
||||
a3_y = a3_1*(y1-y3)/distance13 + a3_2*(y2-y3)/distance23
|
||||
v3_x = v3_x + a3_x*t
|
||||
v3_y = v3_y + a3_y*t
|
||||
x3 = x3 + v3_x*t
|
||||
y3 = y3 + v3_y*t
|
||||
x3_all = np.append(x3_all, x3)
|
||||
y3_all = np.append(y3_all, y3)
|
||||
|
||||
# 选择观测坐标
|
||||
axis_x = np.mean([x1, x2, x3]) # 观测坐标中心固定在平均值的地方
|
||||
axis_y = np.mean([y1, y2, y3]) # 观测坐标中心固定在平均值的地方
|
||||
while True:
|
||||
if np.abs(x1-axis_x) > observation_max or np.abs(x2-axis_x) > observation_max or np.abs(x3-axis_x) > observation_max or\
|
||||
np.abs(y1-axis_y) > observation_max or np.abs(y2-axis_y) > observation_max or np.abs(y3-axis_y) > observation_max:
|
||||
observation_max = observation_max * 2 # 有一个物体超出视线时，视线范围翻倍
|
||||
elif np.abs(x1-axis_x) < observation_max/10 and np.abs(x2-axis_x) < observation_max/10 and np.abs(x3-axis_x) < observation_max/10 and\
|
||||
np.abs(y1-axis_y) < observation_max/10 and np.abs(y2-axis_y) < observation_max/10 and np.abs(y3-axis_y) < observation_max/10:
|
||||
observation_max = observation_max / 2 # 所有物体都在的视线的10分之一内，视线范围减半
|
||||
else:
|
||||
break
|
||||
|
||||
plt.axis([axis_x-observation_max, axis_x+observation_max, axis_y-observation_max, axis_y+observation_max])
|
||||
|
||||
plt.plot(x1, y1, 'og', markersize=m1*100/observation_max) # 默认密度相同，质量越大的，球面积越大。视线范围越宽，球看起来越小。
|
||||
plt.plot(x2, y2, 'or', markersize=m2*100/observation_max)
|
||||
plt.plot(x3, y3, 'ob', markersize=m3*100/observation_max)
|
||||
plt.plot(x1_all, y1_all, '-g') # 画轨迹
|
||||
plt.plot(x2_all, y2_all, '-r')
|
||||
plt.plot(x3_all, y3_all, '-b')
|
||||
|
||||
# plt.show() # 显示图像
|
||||
# plt.pause(0.00001) # 暂停0.00001，防止画图过快
|
||||
|
||||
if np.mod(i, 500) == 0: # 当运动明显时，把图画出来
|
||||
print(i0)
|
||||
plt.savefig(str(i0)+'.jpg') # 保存为图片可以用来做动画
|
||||
i0 += 1
|
||||
plt.clf() # 清空
|
||||
@@ -0,0 +1,73 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1145
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import random
|
||||
import time
|
||||
|
||||
|
||||
def integral(): # 直接数值积分
|
||||
integral_value = 0
|
||||
for x in np.arange(0, 1, 1/10**7):
|
||||
integral_value = integral_value + x**2*(1/10**7) # 对x^2在0和1之间积分
|
||||
return integral_value
|
||||
|
||||
|
||||
def MC_1(): # 蒙特卡洛求定积分1：投点法
|
||||
n = 10**7
|
||||
x_min, x_max = 0.0, 1.0
|
||||
y_min, y_max = 0.0, 1.0
|
||||
count = 0
|
||||
for i in range(0, n):
|
||||
x = random.uniform(x_min, x_max)
|
||||
y = random.uniform(y_min, y_max)
|
||||
# x*x > y，表示该点位于曲线的下面。所求的积分值即为曲线下方的面积与正方形面积的比。
|
||||
if x * x > y:
|
||||
count += 1
|
||||
integral_value = count / n
|
||||
return integral_value
|
||||
|
||||
|
||||
def MC_2(): # 蒙特卡洛求定积分2：期望法
|
||||
n = 10**7
|
||||
x_min, x_max = 0.0, 1.0
|
||||
integral_value = 0
|
||||
for i in range(n):
|
||||
x = random.uniform(x_min, x_max)
|
||||
integral_value = integral_value + (1-0)*x**2
|
||||
integral_value = integral_value/n
|
||||
return integral_value
|
||||
|
||||
|
||||
print('【计算时间】')
|
||||
start_clock = time.perf_counter() # 或者用time.clock()
|
||||
a00 = 1/3 # 理论值
|
||||
end_clock = time.perf_counter()
|
||||
print('理论值（解析）：', end_clock-start_clock)
|
||||
start_clock = time.perf_counter()
|
||||
a0 = integral() # 直接数值积分
|
||||
end_clock = time.perf_counter()
|
||||
print('直接数值积分：', end_clock-start_clock)
|
||||
start_clock = time.perf_counter()
|
||||
a1 = MC_1() # 用蒙特卡洛求积分投点法
|
||||
end_clock = time.perf_counter()
|
||||
print('用蒙特卡洛求积分_投点法：', end_clock-start_clock)
|
||||
start_clock = time.perf_counter()
|
||||
a2 = MC_2()
|
||||
end_clock = time.perf_counter()
|
||||
print('用蒙特卡洛求积分_期望法：', end_clock-start_clock, '\n')
|
||||
|
||||
print('【计算结果】')
|
||||
print('理论值（解析）：', a00)
|
||||
print('直接数值积分：', a0)
|
||||
print('用蒙特卡洛求积分_投点法：', a1)
|
||||
print('用蒙特卡洛求积分_期望法：', a2, '\n')
|
||||
|
||||
print('【计算误差】')
|
||||
print('理论值（解析）：', 0)
|
||||
print('直接数值积分：', abs(a0-1/3))
|
||||
print('用蒙特卡洛求积分_投点法：', abs(a1-1/3))
|
||||
print('用蒙特卡洛求积分_期望法：', abs(a2-1/3))
|
||||
|
||||
@@ -0,0 +1,33 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1247
|
||||
"""
|
||||
|
||||
import random
|
||||
import math
|
||||
import numpy as np
|
||||
from scipy.stats import norm
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def function(x): # 期望样品的分布，target distribution function
|
||||
y = (norm.pdf(x, loc=2, scale=1)+norm.pdf(x, loc=-5, scale=1.5))/2 # loc代表了均值,scale代表标准差
|
||||
return y
|
||||
|
||||
|
||||
T = 100000 # 取T个样品
|
||||
pi = [0 for i in range(T)] # 任意选定一个马尔科夫链初始状态
|
||||
for t in np.arange(1, T):
|
||||
pi_star = np.random.uniform(-10, 10) # a proposed distribution， 例如是均匀分布的，或者是一个依赖pi[t - 1]的分布
|
||||
alpha = min(1, (function(pi_star) / function(pi[t - 1]))) # 接收率
|
||||
u = random.uniform(0, 1)
|
||||
if u < alpha:
|
||||
pi[t] = pi_star # 以alpha的概率接收转移
|
||||
else:
|
||||
pi[t] = pi[t - 1] # 不接收转移
|
||||
|
||||
pi = np.array(pi) # 转成numpy格式
|
||||
print(pi.shape) # 查看抽样样品的维度
|
||||
plt.plot(pi, function(pi), '*') # 画出抽样样品期望的分布 # 或用plt.scatter(pi, function(pi))
|
||||
plt.hist(pi, bins=100, density=1, facecolor='red', alpha=0.7) # 画出抽样样品的分布 # bins是分布柱子的个数，density是归一化，后面两个参数是管颜色的
|
||||
plt.show()
|
||||
@@ -0,0 +1,18 @@
|
||||
% This code is supported by the website: https://www.guanjihuan.com
|
||||
% The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1247
|
||||
|
||||
clc;clear all;clf;
|
||||
s=100000; % 取的样品数
|
||||
f=[1,2,3,3,3,3,6,5,4,3,2,1]; % 期望得到样品的分布函数
|
||||
d=zeros(1,s); % 初始状态
|
||||
x=1;
|
||||
for i=1:s
|
||||
y=unidrnd(12); % 1到12随机整数
|
||||
alpha=min(1,f(y)/f(x)); % 接收率
|
||||
u=rand;
|
||||
if u<alpha % 以alpha的概率接收转移
|
||||
x=y;
|
||||
end
|
||||
d(i)=x;
|
||||
end
|
||||
hist(d,1:1:12);
|
||||
@@ -0,0 +1,97 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1249
|
||||
"""
|
||||
|
||||
import random
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import copy
|
||||
import math
|
||||
import time
|
||||
|
||||
|
||||
def main():
|
||||
size = 30 # 体系大小
|
||||
T = 2 # 温度
|
||||
ising = get_one_sample(sizeOfSample=size, temperature=T) # 得到符合玻尔兹曼分布的ising模型
|
||||
plot(ising) # 画图
|
||||
energy_s = []
|
||||
for i in range(size):
|
||||
for j in range(size):
|
||||
energy_s0 = getEnergy(i=i, j=j, S=ising) # 获取格点的能量
|
||||
energy_s = np.append(energy_s, [energy_s0], axis=0)
|
||||
plt.hist(energy_s, bins=50, density=1, facecolor='red', alpha=0.7) # 画出格点能量分布 # bins是分布柱子的个数，density是归一化，后面两个参数是管颜色的
|
||||
plt.show()
|
||||
|
||||
|
||||
def get_one_sample(sizeOfSample, temperature):
|
||||
S = 2 * np.pi * np.random.random(size=(sizeOfSample, sizeOfSample)) # 随机初始状态，角度是0到2pi
|
||||
print('体系大小：', S.shape)
|
||||
initialEnergy = calculateAllEnergy(S) # 计算随机初始状态的能量
|
||||
print('系统的初始能量是:', initialEnergy)
|
||||
newS = np.array(copy.deepcopy(S))
|
||||
for i00 in range(1000): # 循环一定次数，得到平衡的抽样分布
|
||||
newS = Metropolis(newS, temperature) # Metropolis方法抽样，得到玻尔兹曼分布的样品体系
|
||||
newEnergy = calculateAllEnergy(newS)
|
||||
if np.mod(i00, 100) == 0:
|
||||
print('循环次数%i, 当前系统能量是:' % (i00+1), newEnergy)
|
||||
print('循环次数%i, 当前系统能量是:' % (i00 + 1), newEnergy)
|
||||
return newS
|
||||
|
||||
|
||||
def Metropolis(S, T): # S是输入的初始状态， T是温度
|
||||
delta_max = 0.5 * np.pi # 角度最大的变化度数，默认是90度，也可以调整为其他
|
||||
k = 1 # 玻尔兹曼常数
|
||||
for i in range(S.shape[0]):
|
||||
for j in range(S.shape[0]):
|
||||
delta = (2 * np.random.random() - 1) * delta_max # 角度变化在-90度到90度之间
|
||||
newAngle = S[i, j] + delta # 新角度
|
||||
energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 获取该格点的能量
|
||||
energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 获取格点变成新角度时的能量
|
||||
alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 这个接受率对应的是玻尔兹曼分布
|
||||
if random.uniform(0, 1) <= alpha:
|
||||
S[i, j] = newAngle # 接受新状态
|
||||
else:
|
||||
pass # 保持为上一个状态
|
||||
return S
|
||||
|
||||
|
||||
def getEnergy(i, j, S, angle=None): # 计算(i,j)位置的能量，为周围四个的相互能之和
|
||||
width = S.shape[0]
|
||||
height = S.shape[1]
|
||||
top_i = i - 1 if i > 0 else width - 1 # 用到周期性边界条件
|
||||
bottom_i = i + 1 if i < (width - 1) else 0
|
||||
left_j = j - 1 if j > 0 else height - 1
|
||||
right_j = j + 1 if j < (height - 1) else 0
|
||||
environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]]
|
||||
energy = 0
|
||||
if angle == None:
|
||||
for num_i in range(4):
|
||||
energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]])
|
||||
else:
|
||||
for num_i in range(4):
|
||||
energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]])
|
||||
return energy
|
||||
|
||||
|
||||
def calculateAllEnergy(S): # 计算整个体系的能量
|
||||
energy = 0
|
||||
for i in range(S.shape[0]):
|
||||
for j in range(S.shape[1]):
|
||||
energy += getEnergy(i, j, S)
|
||||
return energy/2 # 作用两次要减半
|
||||
|
||||
|
||||
def plot(S): # 画图
|
||||
X, Y = np.meshgrid(np.arange(0, S.shape[0]), np.arange(0, S.shape[0]))
|
||||
U = np.cos(S)
|
||||
V = np.sin(S)
|
||||
plt.figure()
|
||||
plt.quiver(X, Y, U, V, units='inches')
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
||||
@@ -0,0 +1,76 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1249
|
||||
"""
|
||||
|
||||
import random
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import copy
|
||||
import math
|
||||
import time
|
||||
|
||||
|
||||
def main():
|
||||
size = 30 # 体系大小
|
||||
for T in np.linspace(0.02, 5, 100):
|
||||
ising, magnetism = get_one_sample(sizeOfSample=size, temperature=T)
|
||||
print('温度=', T, ' 磁矩=', magnetism, ' 总能量=', calculateAllEnergy(ising))
|
||||
plt.plot(T, magnetism, 'o')
|
||||
plt.show()
|
||||
|
||||
|
||||
def get_one_sample(sizeOfSample, temperature):
|
||||
newS = np.zeros((sizeOfSample, sizeOfSample)) # 初始状态
|
||||
magnetism = 0
|
||||
for i00 in range(100):
|
||||
newS = Metropolis(newS, temperature)
|
||||
magnetism = magnetism + abs(sum(sum(np.cos(newS))))/newS.shape[0]**2
|
||||
magnetism = magnetism/100
|
||||
return newS, magnetism
|
||||
|
||||
|
||||
def Metropolis(S, T): # S是输入的初始状态， T是温度
|
||||
k = 1 # 玻尔兹曼常数
|
||||
for i in range(S.shape[0]):
|
||||
for j in range(S.shape[0]):
|
||||
newAngle = np.random.randint(-1, 1)*np.pi
|
||||
energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 获取该格点的能量
|
||||
energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 获取格点变成新角度时的能量
|
||||
alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 这个接受率对应的是玻尔兹曼分布
|
||||
if random.uniform(0, 1) <= alpha:
|
||||
S[i, j] = newAngle # 接受新状态
|
||||
else:
|
||||
pass # 保持为上一个状态
|
||||
return S
|
||||
|
||||
|
||||
def getEnergy(i, j, S, angle=None): # 计算(i,j)位置的能量，为周围四个的相互能之和
|
||||
width = S.shape[0]
|
||||
height = S.shape[1]
|
||||
top_i = i - 1 if i > 0 else width - 1 # 用到周期性边界条件
|
||||
bottom_i = i + 1 if i < (width - 1) else 0
|
||||
left_j = j - 1 if j > 0 else height - 1
|
||||
right_j = j + 1 if j < (height - 1) else 0
|
||||
environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]]
|
||||
energy = 0
|
||||
if angle == None:
|
||||
for num_i in range(4):
|
||||
energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]])
|
||||
else:
|
||||
for num_i in range(4):
|
||||
energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]])
|
||||
return energy
|
||||
|
||||
|
||||
def calculateAllEnergy(S): # 计算整个体系的能量
|
||||
energy = 0
|
||||
for i in range(S.shape[0]):
|
||||
for j in range(S.shape[1]):
|
||||
energy += getEnergy(i, j, S)
|
||||
return energy/2 # 作用两次要减半
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
||||
@@ -0,0 +1,61 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import time
|
||||
|
||||
|
||||
|
||||
time_1 = np.array([])
|
||||
time_2 = np.array([])
|
||||
time_3 = np.array([])
|
||||
n_all = np.arange(2,5000,200) # 测试的范围
|
||||
start_all = time.process_time()
|
||||
for n in n_all:
|
||||
print(n)
|
||||
matrix_1 = np.zeros((n,n))
|
||||
matrix_2 = np.zeros((n,n))
|
||||
for i0 in range(n):
|
||||
for j0 in range(n):
|
||||
matrix_1[i0,j0] = np.random.uniform(-10, 10)
|
||||
for i0 in range(n):
|
||||
for j0 in range(n):
|
||||
matrix_2[i0,j0] = np.random.uniform(-10, 10)
|
||||
|
||||
start = time.process_time()
|
||||
matrix_3 = np.dot(matrix_1, matrix_2) # 矩阵乘积
|
||||
end = time.process_time()
|
||||
time_1 = np.append(time_1, [end-start], axis=0)
|
||||
|
||||
start = time.process_time()
|
||||
matrix_4 = np.linalg.inv(matrix_1) # 矩阵求逆
|
||||
end = time.process_time()
|
||||
time_2 = np.append(time_2, [end-start], axis=0)
|
||||
|
||||
start = time.process_time()
|
||||
eigenvalue, eigenvector = np.linalg.eig(matrix_1) # 求矩阵本征值
|
||||
end = time.process_time()
|
||||
time_3 = np.append(time_3, [end-start], axis=0)
|
||||
|
||||
|
||||
|
||||
end_all = time.process_time()
|
||||
print('总共运行时间：', (end_all-start_all)/60, '分')
|
||||
|
||||
plt.subplot(131)
|
||||
plt.xlabel('n^3/10^9')
|
||||
plt.ylabel('时间（秒）')
|
||||
plt.title('矩阵乘积')
|
||||
plt.plot((n_all/10**3)*(n_all/10**3)*(n_all/10**3), time_1, 'o-')
|
||||
|
||||
plt.subplot(132)
|
||||
plt.xlabel('n^3/10^9')
|
||||
plt.title('矩阵求逆')
|
||||
plt.plot((n_all/10**3)*(n_all/10**3)*(n_all/10**3), time_2, 'o-')
|
||||
|
||||
plt.subplot(133)
|
||||
plt.xlabel('n^3/10^9')
|
||||
plt.title('求矩阵本征值')
|
||||
plt.plot((n_all/10**3)*(n_all/10**3)*(n_all/10**3), time_3, 'o-')
|
||||
|
||||
plt.rcParams['font.sans-serif'] = ['SimHei'] # 在画图中正常显示中文
|
||||
plt.rcParams['axes.unicode_minus'] = False # 中文化后，加上这个使正常显示负号
|
||||
plt.show()
|
||||
@@ -0,0 +1,73 @@
|
||||
(* Content-type: application/vnd.wolfram.mathematica *)
|
||||
|
||||
(*** Wolfram Notebook File ***)
|
||||
(* http://www.wolfram.com/nb *)
|
||||
|
||||
(* CreatedBy='Mathematica 12.0' *)
|
||||
|
||||
(*CacheID: 234*)
|
||||
(* Internal cache information:
|
||||
NotebookFileLineBreakTest
|
||||
NotebookFileLineBreakTest
|
||||
NotebookDataPosition[ 158, 7]
|
||||
NotebookDataLength[ 2028, 65]
|
||||
NotebookOptionsPosition[ 1687, 50]
|
||||
NotebookOutlinePosition[ 2075, 67]
|
||||
CellTagsIndexPosition[ 2032, 64]
|
||||
WindowFrame->Normal*)
|
||||
|
||||
(* Beginning of Notebook Content *)
|
||||
Notebook[{
|
||||
Cell[BoxData[
|
||||
RowBox[{"\[IndentingNewLine]",
|
||||
RowBox[{
|
||||
RowBox[{"Clear", "[", "\"\<`*\>\"", "]"}], "\[IndentingNewLine]",
|
||||
RowBox[{
|
||||
RowBox[{"A", "=",
|
||||
RowBox[{"{",
|
||||
RowBox[{
|
||||
RowBox[{"{",
|
||||
RowBox[{"3", ",", "2", ",",
|
||||
RowBox[{"-", "1"}]}], "}"}], ",",
|
||||
RowBox[{"{",
|
||||
RowBox[{
|
||||
RowBox[{"-", "2"}], ",",
|
||||
RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", " ",
|
||||
RowBox[{"{",
|
||||
RowBox[{"3", ",", "6", ",",
|
||||
RowBox[{"-", "1"}]}], "}"}]}], "}"}]}], ";"}], "\[IndentingNewLine]",
|
||||
RowBox[{"MatrixForm", "[", "A", "]"}], "\n",
|
||||
RowBox[{"eigenvalue", "=",
|
||||
RowBox[{"MatrixForm", "[",
|
||||
RowBox[{"Eigenvalues", "[", "A", "]"}], "]"}]}], "\n",
|
||||
RowBox[{"eigenvector", "=",
|
||||
RowBox[{"MatrixForm", "[",
|
||||
RowBox[{"Eigenvectors", "[", "A", "]"}], "]"}]}]}]}]], "Input",
|
||||
CellChangeTimes->{{3.8249089321894894`*^9, 3.8249090194456663`*^9}, {
|
||||
3.8249090545199647`*^9, 3.824909158471919*^9}, {3.8249092372038407`*^9,
|
||||
3.824909243121014*^9}},
|
||||
CellLabel->"In[24]:=",ExpressionUUID->"54dc89aa-0cf7-4c91-9e21-ea7b96cf97e8"]
|
||||
},
|
||||
WindowSize->{1128, 568},
|
||||
WindowMargins->{{Automatic, 375}, {Automatic, 189}},
|
||||
Magnification:>1.2 Inherited,
|
||||
FrontEndVersion->"12.0 for Microsoft Windows (64-bit) (2019\:5e744\:67088\
|
||||
\:65e5)",
|
||||
StyleDefinitions->"Default.nb"
|
||||
]
|
||||
(* End of Notebook Content *)
|
||||
|
||||
(* Internal cache information *)
|
||||
(*CellTagsOutline
|
||||
CellTagsIndex->{}
|
||||
*)
|
||||
(*CellTagsIndex
|
||||
CellTagsIndex->{}
|
||||
*)
|
||||
(*NotebookFileOutline
|
||||
Notebook[{
|
||||
Cell[558, 20, 1125, 28, 238, "Input",ExpressionUUID->"54dc89aa-0cf7-4c91-9e21-ea7b96cf97e8"]
|
||||
}
|
||||
]
|
||||
*)
|
||||
|
||||
@@ -0,0 +1,3 @@
|
||||
clc;clear all;
|
||||
A = [3, 2, -1; -2, -2, 2; 3, 6, -1]
|
||||
[V,D] = eig(A)
|
||||
@@ -0,0 +1,40 @@
|
||||
program main
|
||||
use lapack95
|
||||
implicit none
|
||||
integer i,j,info
|
||||
complex*16 A(3,3), eigenvalues(3), eigenvectors(3,3)
|
||||
|
||||
A(1,1)=(3.d0, 0.d0)
|
||||
A(1,2)=(2.d0, 0.d0)
|
||||
A(1,3)=(-1.d0, 0.d0)
|
||||
A(2,1)=(-2.d0, 0.d0)
|
||||
A(2,2)=(-2.d0, 0.d0)
|
||||
A(2,3)=(2.d0, 0.d0)
|
||||
A(3,1)=(3.d0, 0.d0)
|
||||
A(3,2)=(6.d0, 0.d0)
|
||||
A(3,3)=(-1.d0, 0.d0)
|
||||
|
||||
write(*,*) 'matrix:'
|
||||
do i=1,3
|
||||
do j=1,3
|
||||
write(*,"(f10.4, '+1i*',f7.4)",advance='no') A(i,j) ! <20><>ѭ<EFBFBD><D1AD>Ϊ<EFBFBD>е<EFBFBD>ָ<EFBFBD><D6B8>
|
||||
enddo
|
||||
write(*,*) ''
|
||||
enddo
|
||||
|
||||
call geev(A=A, W=eigenvalues, VR=eigenvectors, INFO=info)
|
||||
|
||||
write(*,*) 'eigenvectors:'
|
||||
do i=1,3
|
||||
do j=1,3
|
||||
write(*,"(f10.4, '+1i*',f7.4)",advance='no') eigenvectors(i,j) ! <20><>ѭ<EFBFBD><D1AD>Ϊ<EFBFBD>е<EFBFBD>ָ<EFBFBD>ꡣ<EFBFBD><EAA1A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
enddo
|
||||
write(*,*) ''
|
||||
enddo
|
||||
write(*,*) 'eigenvalues:'
|
||||
do i=1,3
|
||||
write(*,"(f10.4, '+1i*',f7.4)",advance='no') eigenvalues(i)
|
||||
enddo
|
||||
write(*,*) ''
|
||||
write(*,*) ''
|
||||
end program
|
||||
@@ -0,0 +1,21 @@
|
||||
import numpy as np
|
||||
|
||||
A = np.array([[3, 2+1j, -1], [2-1j, -2, 6], [-1, 6, 1]])
|
||||
eigenvalue, eigenvector = np.linalg.eig(A)
|
||||
print('矩阵：\n', A)
|
||||
print('特征值:\n', eigenvalue)
|
||||
print('特征向量:\n', eigenvector)
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(eigenvector.transpose().conj(), eigenvector))
|
||||
print('判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose().conj()))
|
||||
|
||||
print('特征向量矩阵的列向量模方和：')
|
||||
for i in range(3):
|
||||
print(np.abs(eigenvector[0, i])**2+np.abs(eigenvector[1, i])**2+np.abs(eigenvector[2, i])**2)
|
||||
print('特征向量矩阵的行向量模方和：')
|
||||
for i in range(3):
|
||||
print(np.abs(eigenvector[i, 0])**2+np.abs(eigenvector[i, 1])**2+np.abs(eigenvector[i, 2])**2)
|
||||
|
||||
print('\n对角化验证：')
|
||||
print(np.dot(np.dot(eigenvector.transpose().conj(), A), eigenvector))
|
||||
print(np.dot(np.dot(eigenvector, A), eigenvector.transpose().conj()))
|
||||
@@ -0,0 +1,18 @@
|
||||
import numpy as np
|
||||
|
||||
A = np.array([[3, 2, -1], [-2, -2, 2], [3, 6, -1]])
|
||||
eigenvalue, eigenvector = np.linalg.eig(A)
|
||||
print('矩阵：\n', A)
|
||||
print('特征值:\n', eigenvalue)
|
||||
print('特征向量:\n', eigenvector)
|
||||
print('特征值为-4对应的特征向量理论值:\n', np.array([1/3, -2/3, 1])/np.sqrt((1/3)**2+(-2/3)**2+1**2))
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(eigenvector.transpose(), eigenvector))
|
||||
print('判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose()))
|
||||
|
||||
print('特征向量矩阵的列向量模方和：')
|
||||
for i in range(3):
|
||||
print(eigenvector[0, i]**2+eigenvector[1, i]**2+eigenvector[2, i]**2)
|
||||
print('特征向量矩阵的行向量模方和：')
|
||||
for i in range(3):
|
||||
print(eigenvector[i, 0]**2+eigenvector[i, 1]**2+eigenvector[i, 2]**2)
|
||||
@@ -0,0 +1,21 @@
|
||||
import numpy as np
|
||||
|
||||
A = np.array([[3, 2, -1], [2, -2, 6], [-1, 6, 1]])
|
||||
eigenvalue, eigenvector = np.linalg.eig(A)
|
||||
print('矩阵：\n', A)
|
||||
print('特征值:\n', eigenvalue)
|
||||
print('特征向量:\n', eigenvector)
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(eigenvector.transpose(), eigenvector))
|
||||
print('判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose()))
|
||||
|
||||
print('特征向量矩阵的列向量模方和：')
|
||||
for i in range(3):
|
||||
print(eigenvector[0, i]**2+eigenvector[1, i]**2+eigenvector[2, i]**2)
|
||||
print('特征向量矩阵的行向量模方和：')
|
||||
for i in range(3):
|
||||
print(eigenvector[i, 0]**2+eigenvector[i, 1]**2+eigenvector[i, 2]**2)
|
||||
|
||||
print('\n对角化验证：')
|
||||
print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))
|
||||
print(np.dot(np.dot(eigenvector, A), eigenvector.transpose()))
|
||||
@@ -0,0 +1,32 @@
|
||||
import numpy as np
|
||||
|
||||
A = np.array([[0, 1, 1, -1], [1, 0, -1, 1], [1, -1, 0, 1], [-1, 1, 1, 0]])
|
||||
eigenvalue, eigenvector = np.linalg.eig(A)
|
||||
print('矩阵：\n', A)
|
||||
print('特征值:\n', eigenvalue)
|
||||
print('特征向量:\n', eigenvector)
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(eigenvector.transpose(), eigenvector))
|
||||
print('判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose()))
|
||||
|
||||
print('特征向量矩阵的列向量模方和：')
|
||||
for i in range(4):
|
||||
print(eigenvector[0, i]**2+eigenvector[1, i]**2+eigenvector[2, i]**2+eigenvector[3, i]**2)
|
||||
print('特征向量矩阵的行向量模方和：')
|
||||
for i in range(4):
|
||||
print(eigenvector[i, 0]**2+eigenvector[i, 1]**2+eigenvector[i, 2]**2+eigenvector[i, 3]**2)
|
||||
|
||||
print('\n对角化验证：')
|
||||
print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))
|
||||
print(np.dot(np.dot(eigenvector, A), eigenvector.transpose()))
|
||||
|
||||
print('\n特征向量理论值:')
|
||||
T = np.array([[1/np.sqrt(2), 1/np.sqrt(6), -1/np.sqrt(12), 1/2], [1/np.sqrt(2), -1/np.sqrt(6), 1/np.sqrt(12), -1/2], [0, 2/np.sqrt(6), 1/np.sqrt(12), -1/2], [0, 0, 3/np.sqrt(12), 1/2]])
|
||||
print(T)
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(T.transpose(), T))
|
||||
print('判断是否正交：\n', np.dot(T, T.transpose()))
|
||||
|
||||
print('\n对角化验证：')
|
||||
print(np.dot(np.dot(T.transpose(), A), T))
|
||||
print(np.dot(np.dot(T, A), T.transpose()))
|
||||
@@ -0,0 +1,44 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/10890
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def main():
|
||||
A = np.array([[0, 1, 1, -1], [1, 0, -1, 1], [1, -1, 0, 1], [-1, 1, 1, 0]])
|
||||
eigenvalue, eigenvector = np.linalg.eig(A)
|
||||
print('矩阵：\n', A)
|
||||
print('特征值:\n', eigenvalue)
|
||||
print('特征向量:\n', eigenvector)
|
||||
|
||||
print('\n判断是否正交：\n', np.dot(eigenvector.transpose(), eigenvector))
|
||||
print('判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose()))
|
||||
|
||||
print('对角化验证：')
|
||||
print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))
|
||||
|
||||
# 施密斯正交化
|
||||
eigenvector = Schmidt_orthogonalization(eigenvector)
|
||||
|
||||
print('\n施密斯正交化后，特征向量：\n', eigenvector)
|
||||
|
||||
print('施密斯正交化后，判断是否正交：\n', np.dot(eigenvector.transpose(), eigenvector))
|
||||
print('施密斯正交化后，判断是否正交：\n', np.dot(eigenvector, eigenvector.transpose()))
|
||||
|
||||
print('施密斯正交化后，对角化验证：')
|
||||
print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))
|
||||
|
||||
|
||||
def Schmidt_orthogonalization(eigenvector):
|
||||
num = eigenvector.shape[1]
|
||||
for i in range(num):
|
||||
for i0 in range(i):
|
||||
eigenvector[:, i] = eigenvector[:, i] - eigenvector[:, i0]*np.dot(eigenvector[:, i].transpose().conj(), eigenvector[:, i0])/(np.dot(eigenvector[:, i0].transpose().conj(),eigenvector[:, i0]))
|
||||
eigenvector[:, i] = eigenvector[:, i]/np.linalg.norm(eigenvector[:, i])
|
||||
return eigenvector
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -0,0 +1,45 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/15978
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from math import *
|
||||
|
||||
def main():
|
||||
a1 = [0, 1]
|
||||
a2 = [1, 0]
|
||||
b1, b2 = calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2)
|
||||
print(b1, b2)
|
||||
|
||||
|
||||
def calculate_one_dimensional_reciprocal_lattice_vector(a1):
|
||||
b1 = 2*pi/a1
|
||||
return b1
|
||||
|
||||
|
||||
def calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2):
|
||||
a1 = np.array(a1)
|
||||
a2 = np.array(a2)
|
||||
a1 = np.append(a1, 0)
|
||||
a2 = np.append(a2, 0)
|
||||
a3 = np.array([0, 0, 1])
|
||||
b1 = 2*pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3))
|
||||
b2 = 2*pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3))
|
||||
b1 = np.delete(b1, 2)
|
||||
b2 = np.delete(b2, 2)
|
||||
return b1, b2
|
||||
|
||||
|
||||
def calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3):
|
||||
a1 = np.array(a1)
|
||||
a2 = np.array(a2)
|
||||
a3 = np.array(a3)
|
||||
b1 = 2*pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3))
|
||||
b2 = 2*pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3))
|
||||
b3 = 2*pi*np.cross(a1, a2)/np.dot(a1, np.cross(a2, a3))
|
||||
return b1, b2, b3
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -0,0 +1,22 @@
|
||||
import guan
|
||||
import sympy
|
||||
|
||||
a1 = [0, 1]
|
||||
a2 = [1, 0]
|
||||
b1, b2 = guan.calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2)
|
||||
print(b1, b2, '\n')
|
||||
|
||||
print('bases in the real space')
|
||||
a = sympy.symbols('a')
|
||||
a1 = sympy.Matrix(1, 2, [3*a/2, sympy.sqrt(3)*a/2])
|
||||
a2 = sympy.Matrix(1, 2, [3*a/2, -sympy.sqrt(3)*a/2])
|
||||
print('a1:')
|
||||
sympy.pprint(a1)
|
||||
print('a2:')
|
||||
sympy.pprint(a2)
|
||||
print('\nbases in the reciprocal space')
|
||||
b1, b2 = guan.calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2)
|
||||
print('b1:')
|
||||
sympy.pprint(b1)
|
||||
print('b2:')
|
||||
sympy.pprint(b2)
|
||||
@@ -0,0 +1,55 @@
|
||||
"""
|
||||
This code is supported by the website: https://www.guanjihuan.com
|
||||
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/15978
|
||||
"""
|
||||
|
||||
import sympy
|
||||
|
||||
|
||||
def main():
|
||||
print('bases in the real space')
|
||||
a = sympy.symbols('a')
|
||||
a1 = sympy.Matrix(1, 2, [3*a/2, sympy.sqrt(3)*a/2])
|
||||
a2 = sympy.Matrix(1, 2, [3*a/2, -sympy.sqrt(3)*a/2])
|
||||
print('a1:')
|
||||
sympy.pprint(a1)
|
||||
print('a2:')
|
||||
sympy.pprint(a2)
|
||||
print('\nbases in the reciprocal space')
|
||||
b1, b2 = calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2)
|
||||
print('b1:')
|
||||
sympy.pprint(b1)
|
||||
print('b2:')
|
||||
sympy.pprint(b2)
|
||||
|
||||
|
||||
def calculate_one_dimensional_reciprocal_lattice_vector_with_sympy(a1):
|
||||
b1 = 2*sympy.pi/a1
|
||||
return b1
|
||||
|
||||
|
||||
def calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2):
|
||||
a1 = sympy.Matrix(1, 3, [a1[0], a1[1], 0])
|
||||
a2 = sympy.Matrix(1, 3, [a2[0], a2[1], 0])
|
||||
a3 = sympy.Matrix(1, 3, [0, 0, 1])
|
||||
cross_a2_a3 = a2.cross(a3)
|
||||
cross_a3_a1 = a3.cross(a1)
|
||||
b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3)
|
||||
b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3)
|
||||
b1 = sympy.Matrix(1, 2, [b1[0], b1[1]])
|
||||
b2 = sympy.Matrix(1, 2, [b2[0], b2[1]])
|
||||
return b1, b2
|
||||
|
||||
|
||||
def calculate_three_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2, a3):
|
||||
cross_a2_a3 = a2.cross(a3)
|
||||
cross_a3_a1 = a3.cross(a1)
|
||||
cross_a1_a2 = a1.cross(a2)
|
||||
b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3)
|
||||
b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3)
|
||||
b3 = 2*sympy.pi*cross_a1_a2/a1.dot(cross_a2_a3)
|
||||
return b1, b2, b3
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
Reference in New Issue
Block a user