如何使用隱式舍姆

嗨，我是新在這個網站來解決一些2D與時間有關的擴散功能蟒蛇迭代和我在尋找我的Python代碼的答案如何使用隱式舍姆

我的代碼只是解決地下水頭擴散（使用2d隱式分化）。換句話說，我解決依賴
與拉普拉斯方程的函數隨着時間：d^2H/DX^2 + d^2H/DY^2 = DS/dt的其中DH是頭在流路中的變化（ dx = dy） ds是一些常數值（儲存係數），dt是時間增量

所以我解決了這個例子，在邊界條件上游5米和下游3米的條件下，我的網格是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

這是我的結果與蟒蛇

Number of iterations = 132 
[[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ] 
[ 5.   4.40796148 3.88890558 3.5072111 3.24606331 3.0865587 3.  ]] 
Error = 9.17176556925e-10

與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次迭代

這裏請注意，我使用了高斯賽德爾迭代

來源

2016-06-11 Hassan Saleem

哪一個是正確的，你的python實現還是你的matlab嗎？（或不是？） –

Matlab是好的，但我有問題的Python不是。爲了解釋更多從matlab結果中獲取一行，並從python中獲得一行（同一行）（從5到3），您可以看到matlab結果更加線性，而python則是曲線 –

我想我想通了什麼是錯我的Python代碼

我應該讓

for k in range(ntm-1):

，而不是

for k in range(1, ntm-1):

由於矩陣蟒蛇第一要素是從零開始，以及我的初始條件頭值出現在哪裏。也就是說，在短暫的地下水頭將開始瀰漫達到平衡和頭部將採取直（線性）線（我一直在用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

來源

2016-06-13 07:30:47

如何使用隱式舍姆

回答

相關問題