如何用隐式模式修正二维时变扩散函数的python迭代

Question

嗨，我是这个网站的新人，我正在寻找我的python代码的答案

我的代码是简单地求解地下水水头扩散（使用二维隐式微分）。换言之，我将函数解为时间相关的
拉普拉斯方程：d^2h/dx^2+d^2h/dy^2=ds/dt 式中，dh为空间中水头过流路径的变化（dx=dy） ds是某个定值（存储系数），dt是时间增量

因此，在上游5m和下游3m的边界条件下，我的网格是7x7个单元

我用C#和matlab解决了这个问题，在转换后两个答案都很接近（对于某些精确的解析解达到平衡），但是当我使用Python时，我的解没有达到平衡

下面是我的python代码

 import numpy as np

h = np.zeros((7,7,10))
h0 = 5.
K = 2.
DT = 1.
DX = 100.
S = 0.0005
B = 8.
D = (K*B)/S
r = (D*DT)/(DX*DX)

dummy = h.shape
nrow = dummy[0]
ncol = dummy[1]
ntm = dummy[2]
print 'Groundwater Head matrix is a ', nrow, ' by ', ncol, ' matrix.','with',ntm,'increments'
for i in range(nrow):
    for j in  range(ncol):
        for k in range(1):
         h[i,j,k] = h0
         continue

        h[0:,0,1:] = 5.
        h[0:,6,1:] = 3.

        ni = 1
conv_crit = 1e-9
converged = False
while (not converged):
  max_err = 0
  for i in range(1, nrow - 1):
      for j in range(1,ncol - 1):
          for k in range(1,ntm-1):

            h_old = h[i,j,k+1]
            h[i,j,k+1] = (r*h[i-1,j,k+1]+r*h[i+1,j,K+1]+r*h[i,j-1,K+1]+
                          r*h[i,j+1,k+1]+h[i,j,k])/(1+4*r)

            h[0,0:6,k+1] = h[2,0:6,k+1]
            h[6,0:6,k+1] = h[4,0:6,k+1]

            diff = h[i,j,k+1] - h_old
            if (diff > max_err):
             max_err = diff
  if (max_err <= conv_crit):
      converged = True
  ni = ni + 1
  print 'Number of iterations = ', ni - 1
  print h[0:,0:,9]
  print ' Error =', max_err

这是我对python的结果

^{pr2}$

这里有matlab

5.0000    4.6599    4.3229    3.9892    3.6583    3.3291    3.0000
    5.0000    4.6592    4.3217    3.9879    3.6573    3.3285    3.0000
    5.0000    4.6599    4.3229    3.9892    3.6583    3.3291    3.0000
    5.0000    4.6606    4.3240    3.9904    3.6593    3.3296    3.0000
    5.0000    4.6613    4.3250    3.9915    3.6602    3.3301    3.0000
    5.0000    4.6619    4.3259    3.9925    3.6610    3.3305    3.0000
    5.0000    4.6613    4.3250    3.9915    3.6602    3.3301    3.0000

matlab只使用41次迭代

注意，我使用了高斯-赛德尔迭代法

Answer 1

我想我已经知道我的python代码出了什么问题

我应该让

for k in range(ntm-1):

而不是

^{pr2}$

因为python矩阵中的第一个元素是从零开始的，而我的初始条件head值就是在这个位置出现的。也就是说，在瞬变过程中，地下水位开始扩散，达到平衡，水头呈直线（直线）形状 python与一个保持混合索引（从matlab开始） 这是我的最终结果

   5.          4.67229917  4.34411362  4.01349084  3.67944236  3.34142973  3.        ]
 [ 5.          4.67207003  4.34374374  4.01309249  3.6791206   3.3412565   3.        ]
 [ 5.          4.67229917  4.34411362  4.01349084  3.67944236  3.34142973  3.        ]
 [ 5.          4.67251271  4.34445832  4.01386208  3.67974223  3.34159116  3.        ]
 [ 5.          4.67271172  4.34477956  4.01420805  3.68002168  3.34174161  3.        ]
 [ 5.          4.67289718  4.34507894  4.01453048  3.68028211  3.34188181  3.        ]
 [ 5.          4.67271172  4.34477956  4.01420805  3.68002168  3.34174161  3.        ]]
 Error = 0.00088610885882

即使是迭代次数也有所提高，只有49次迭代

谢谢大家

***请注意，我在之前的代码中将公差值从1e-9更改为1e-3