许久没接触python,又有点忘了
继续学习流沙公众号
对这个方程,非定常项,以及x和y方向的对流项
碰巧今天看安德鲁的计算流体力学偏微分方程所讲,利用克莱姆法则和特征值法来判断方程的属性
例子为一阶,在做习题碰到二阶。
此为闲话
同时,根据文章所写,稍加修改,将for循环和数组放在同一程序中进行对比,以温习python,发现遗忘甚多。
from mpl_toolkits.mplot3d import Axes3D # 新
import numpy as np
import matplotlib.pyplot as plt
import time
nx = 81
ny = 81
nt = 100 # 时间步
c = 1 # 常数
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = 0.2
dt = sigma * dx
x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)
u = np.ones((ny,nx))
un = np.ones((ny,nx))
u[int(0.5/dy):int(1/dy+1),int(0.5/dx):int(1/dx+1)] = 2
print("forloop? matrix?")
choose = input()
# 对比for循环运算和数组运算的速度
if choose == 'forloop':
# for 循环
start = time.perf_counter()
for n in range(nt + 1): ## loop across number of time steps
un = u.copy()
row, col = u.shape
for j in range(1,row):
for i in range(1, col):
u[j, i] = (un[j,i] - (c * dt / dx * (un[j, i] - un[j, i - 1])) - (c * dt / dy * (un[j, i] - un[j - 1, i])))
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
end = time.perf_counter()
else:
# 数组运算
start = time.perf_counter()
for n in range(nt + 1): ## matrix
un = u.copy()
u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, :-1])) - (c * dt / dy * (un[1:, 1:] - un[:-1, 1:])))
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
end = time.perf_counter()
time_elapsed = end - start # 普通计算耗时
print("time: %f" %(time_elapsed))
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = Axes3D(fig) # 说明此图为3D
x,y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x, y, u[:],cmap = 'viridis') # cmap 表示的名字
plt.show()
对于此程序代码的说明:
1、from mpl_toolkits.mplot3d import Axes3D ,是3维绘图的库ax = Axes3D(fig)
2、放弃time.clock(),改为time.perf_counter()计时
结果
数组计算的时间为0.0206
for循环的计算时间为3.604