前面的案例大多数是一维的问题，从现在开始我们进入二维的世界。

事实上将一维问题扩展到二维甚至三维都是非常简单的，采用相同的思路。在2D空间中，结构网格可定义为：

注：

注意这里所提到的结构网格，我们在后面还会详细介绍。

”

因此，可定义一阶差分格式：

下面来处理二维线性对流方程。

1 二维线性对流模型

二维线性对流控制方程为：

这里时间项采用向前差分，空间项采用向后差分，离散方程可写成以下格式：

式中，i为x方向角标，j为y方向角标，n为时间项角标。

可得待求项：

采用初始条件：

边界条件：

2 Python代码

先用代码将初始条件和边界条件表达出来。

using PyPlot
 
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = 0.2
dt = sigma * dx
x = range(0, 2, length=nx)
y = range(0, 2, length=ny)
 
u = ones(ny, nx)
un = ones(ny, nx)
 
# 指定初始值
u[floor(Int, 0.5 / dy):floor(Int, 1 / dy),floor(Int, 0.5 / dx):floor(Int, 1 / dx)] .= 2
PyPlot.surf(x,y,u)
PyPlot.show()

初始条件如下图所示。

下面开始迭代计算。我们可以分别采用循环和数组操作来实现，自己体会他们计算时间上的区别。

 
for n in 1:nt
    un = copy(u)
    row, col = size(u)
    for j in 2:row - 1
        for i in 2:col - 1
            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[1,:] .= 1
            u[end,:] .= 1
            u[:,1] .= 1
            u[:,end] .= 1
        end
    end
end
fig = figure(figsize=(10, 6), dpi=100)
PyPlot.surf(x,y,u)
PyPlot.show()

计算结果如下图所示。

利用for循环计算速度很慢，下面改用数组运算试试。

for n in 1:nt
    un = copy(u)
    u[2:end,2:end] = un[2:end,2:end] - (c * dt / dx * (un[2:end,2:end] - un[2:end,1:end - 1])) - (c * dt / dy * (un[2:end,2:end] - un[1:end - 1,2:end]))
    u[1,:] .= 1
    u[end,:] .= 1
    u[:,1] .= 1
    u[:,end] .= 1
end
 
fig = figure(figsize=(10, 6), dpi=100)
PyPlot.surf(x, y, u)
PyPlot.show()

相同的计算结果，如下图所示。

完整程序代码如下所示。

using PyPlot
 
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = 0.2
dt = sigma * dx
x = range(0, 2, length=nx)
y = range(0, 2, length=ny)
 
u = ones(ny, nx)
un = ones(ny, nx)
 
# 指定初始值
u[floor(Int, 0.5 / dy):floor(Int, 1 / dy),floor(Int, 0.5 / dx):floor(Int, 1 / dx)] .= 2
# PyPlot.surf(x,y,u)
# PyPlot.show()
 
##循环方式进行计算
# for n in 1:nt
#     un = copy(u)
#     row, col = size(u)
#     for j in 2:row - 1
#         for i in 2:col - 1
#             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[1,:] .= 1
#             u[end,:] .= 1
#             u[:,1] .= 1
#             u[:,end] .= 1
#         end
#     end
# end
 
# 数组方式运行
for n in 1:nt
    un = copy(u)
    u[2:end,2:end] = un[2:end,2:end] - (c * dt / dx * (un[2:end,2:end] - un[2:end,1:end - 1])) - (c * dt / dy * (un[2:end,2:end] - un[1:end - 1,2:end]))
    u[1,:] .= 1
    u[end,:] .= 1
    u[:,1] .= 1
    u[:,end] .= 1
end
 
fig = figure(figsize=(10, 6), dpi=100)
PyPlot.surf(x, y, u)
PyPlot.show()