前面的案例大多數是一維的問題，從現在開始我們進入二維的世界。
事實上將一維問題擴展到二維甚至三維都是非常簡單的，採用相同的思路。在2D空間中，結構網格可定義為：
注意這裡所提到的結構網格，我們在後面還會詳細介紹。
因此，可定義一階差分格式：
下面來處理二維線性對流方程。
二維線性對流控制方程為：
這裡時間項採用向前差分，空間項採用向後差分，離散方程可寫成以下格式：
式中，i為x方向角標，j為y方向角標，n為時間項角標。
可得待求項：
採用初始條件：
邊界條件：
先用代碼將初始條件和邊界條件表達出來。
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
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
fig = plt.figure(figsize=(12,8))
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf = ax.plot_surface(x,y,u[:],cmap = cm.viridis)
plt.show()
初始條件如下圖所示。
下面開始迭代計算。我們可以分別採用循環和數組操作來實現，自己體會他們計算時間上的區別。
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
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x, y, u[:], cmap=cm.viridis)
plt.show()
計算結果如下圖所示。
利用for循環計算速度很慢，下面改用數組運算試試。
for n in range(nt + 1): ##loop across number of time steps
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
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x, y, u[:], cmap=cm.viridis)
plt.show()
相同的計算結果，如下圖所示。
本案例中，利用for循環所需的計算時間約為數組運算的400倍。
給大家推薦一本神書《批判性思維》