程式碼：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

nx = 31 # 定義x方向節點數
ny = 31 # 定義y方向節點數
nt = 17 # 定義時間步長
nu = 0.05 # 擴散係數
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = 0.25 # 庫朗數
dt = sigma * dx * dy / nu   # 非常不理解，nu是擴散係數，不是傳播速度？還是類比擴散這個傳播速度

x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)
# 邊界條件
u = np.ones((nx,ny))
un = np.ones((nx,ny))
# 初始條件
u[int(0.5/dx):int(1/dx+1),int(0.5/dy):int(1/dy+1)] = 2

nt= int(input("Input:"))

# 定義函式進行求解
for n in range(nt + 1):
	un = u.copy()
	u[1:-1,1:-1] = (un[1:-1,1:-1] +
		             nu * dt / dx**2 * (un[2: ,1:-1] - 2 * un[1:-1,1:-1] + un[0:-2,1:-1]) +
		             nu * dt / dx**2 * (un[1:-1,2:] - 2 * un[1:-1,1:-1] + un[1:-1, 0:-2]) )
	u[0, :] = 1
	u[-1, :] = 1
	u[:, 0] = 1
	u[:, -1] = 1

fig = plt.figure()
ax = Axes3D(fig)
x, y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x,y,u[:],rstride=1,cstride=1,cmap='viridis',linewidth=0,antialiased=False)
ax.set_xlim(0,2)
ax.set_ylim(0,2)
ax.set_zlim(1,2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.show()
 

不理解：
1、dt = sigma * dx * dy / nu # 非常不理解，nu是擴散係數，不是傳播速度？還是類比擴散這個傳播速度

結果



時間分別為10， 15， 50