计算科学 - 使用连续过度松弛技术计算泊松方程的误差，Python - 吾爱随笔录

我正在尝试解决泊松方程

$\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = 32(x(x-1) + y(y-1))$

对于使用 Python3 的 61x61 网格，边界条件为 $T=0$ 在所有四个边界上并初步猜测 $T=50$ . 我正在考虑一个松弛参数 $w=1.6$ .

为了使用 SOR 技术求解泊松方程，我使用有限差分法进行离散化并执行以下操作：

$T^{(n+1)}_{i,j} = 0.25w{(T^{n}_{i+1, j} + T^{n}_{i-1, j} + T^{n}_{i, j+1} +T^{n}_{i, j-1} - 32(x_i(x_i - 1) + y_i(y_i-1))dx^2)} + T^n_{i,j}(1-w)$

我想通过简单地找到旧迭代和新迭代之间的差异（通过找到数组中的最大差异）来找到迭代之间的错误，并检查该错误是否小于我的允许错误（ $10^{-7}$ ）。如果它较小，则迭代停止，否则它们继续。

但是，即使经过 2000 次迭代，这种方法似乎也不会收敛。我想知道我在哪里犯了错误（我怀疑它在错误计算中）。为什么这个简单的误差计算不起作用？我可以实施哪些更改以使其正常工作？

这是我的代码：

import numpy as np 
import matplotlib.pyplot as plt 
 
# Set Dimension and delta
nx = 61 #grid size
my = 61
x = 1.0 #total x length
y = 1.0 #total y length
dx = x/(nx-1)
dy = y/(my-1)
xarr = np.linspace(0,x,nx)
yarr = np.linspace(0,y,my)
print(xarr)

"""
Output:
[0.         0.01666667 0.03333333 0.05       0.06666667 0.08333333
 0.1        0.11666667 0.13333333 0.15       0.16666667 0.18333333
 0.2        0.21666667 0.23333333 0.25       0.26666667 0.28333333
 0.3        0.31666667 0.33333333 0.35       0.36666667 0.38333333
 0.4        0.41666667 0.43333333 0.45       0.46666667 0.48333333
 0.5        0.51666667 0.53333333 0.55       0.56666667 0.58333333
 0.6        0.61666667 0.63333333 0.65       0.66666667 0.68333333
 0.7        0.71666667 0.73333333 0.75       0.76666667 0.78333333
 0.8        0.81666667 0.83333333 0.85       0.86666667 0.88333333
 0.9        0.91666667 0.93333333 0.95       0.96666667 0.98333333
 1.        ]

 """
# Boundary condition
Ttop = 0
Tbottom = 0
Tleft = 0
Tright = 0

# Initial guess of interior grid
Tguess = 50

# Set colour interpolation and colour map
colorinterpolation = 50
colourMap = plt.cm.jet

# Set meshgrid
X, Y = np.meshgrid(np.arange(0, nx), np.arange(0, my))

# Set array size and set the interior value with Tguess
T = np.empty((nx, my))
T.fill(Tguess)

#Boundary conditions
T[(my-1):, :] = Ttop
T[:1, :] = Tbottom
T[:, (nx-1):] = Tright
T[:, :1] = Tleft
T_init=T
print("The initial matrix is: \n", T_init)

"""
Output:
The initial matrix is: 
 [[ 0.  0.  0. ...  0.  0.  0.]
 [ 0. 50. 50. ... 50. 50.  0.]
 [ 0. 50. 50. ... 50. 50.  0.]
 ...
 [ 0. 50. 50. ... 50. 50.  0.]
 [ 0. 50. 50. ... 50. 50.  0.]
 [ 0.  0.  0. ...  0.  0.  0.]]

"""

#SOR Technique
def SORAlgo(error, w, T, MaxIter):
    for n in range(MaxIter):
        Tn=T.copy()
        #Solving the Poisson Equation using array operations
        T[1:-1, 1:-1] = w*0.25*((Tn[2:, 1:-1] + Tn[:-2, 1:-1])*dy**2 + (Tn[1:-1, 2:] + Tn[1:-1, :-2])*dx**2 - 32*((xarr[1:-1]*(xarr[1:-1]-1) + yarr[1:-1]*(yarr[1:-1]-1)))*dx**2*dy**2)/(2*(dx**2 + dy**2)) + (1-w)*Tn[1:-1, 1:-1]

        #Tn will be the older value, T will be the newer value. Finding the max difference in corresponding values of both arrays
        max_error = (abs(T-Tn)).max()

        if max_error<error or n==MaxIter-2:
            print("The relaxation parameter is: ", w)
            print("The number of iterations taken is: ", n)
            print("The error is: ", max_error)
            break

    return Tn 

#Taking error = 10^-7, relaxation parameter=1.6 and maximum iterations=2000
T_final = SORAlgo(0.0000001, 1.6, T_init, 2000)
print(T_final)

#print("The final matrix is: ", T_final)

cp = plt.contourf(X, Y, T_final, colorinterpolation, cmap=colourMap)
plt.colorbar()
plt.show()

这是我的输出：

The relaxation parameter is:  1.6
The number of iterations taken is:  1998
The error is:  2.35862095815448e-05
[[0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 2.25244388e-05 4.41640095e-05 ... 4.41640095e-05
  2.25244388e-05 0.00000000e+00]
 [0.00000000e+00 2.40057069e-05 4.75951601e-05 ... 4.75951601e-05
  2.40057069e-05 0.00000000e+00]
 ...
 [0.00000000e+00 2.40057069e-05 4.75951601e-05 ... 4.75951601e-05
  2.40057069e-05 0.00000000e+00]
 [0.00000000e+00 2.25244388e-05 4.41640095e-05 ... 4.41640095e-05
  2.25244388e-05 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
  0.00000000e+00 0.00000000e+00]]

请注意它是如何进行到近 2000 次迭代并停止的，只是因为我明确要求循环在 1999 次迭代时中断并且错误高于指定的（ $10^{-7}$ ）。这是情节的链接，因为我无法将其直接粘贴到此处：

https://drive.google.com/file/d/1xLbVSp9XA92saZr26b0gMXYZZ_YsR2GQ/view?usp=sharing

谢谢你。

编辑：我接受了评论中的建议。

我调整了源项并简化了方程，所以不是

Tn[1:-1, 1:-1] = w*0.25*(Tn[2:, 1:-1] + Tn[:-2, 1:-1] + Tn[1:-1, 2:] + Tn[1:-1, :-2])/(dx**2) - 32*((xarr[1:-1]*(xarr[1:-1]-1) + yarr[1:-1]*(yarr[1:-1]-1)))*dx**2*dy**2 + (1-w)*Tn[1:-1, 1:-1]

我做了：

b  = np.zeros((my, nx)) #The new source term

for i,j in zip(range(nx), range(my)):
    b[i, j] = 32*(xarr[i]*(xarr[i]-1) + yarr[j]*(yarr[j]-1))
...
T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1] - dx**2*b[1:-1, 1:-1])

当我将源项设为 0 并将松弛参数设为 1 时，即

T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1])       #w=1

我得到了正确的情节，解决方案在 349 次迭代中收敛： https ://drive.google.com/file/d/1IFqzcjkCxmlQjJwGdfwp7ACZVsP9Dd-Q/view?usp=sharing

当我使用源项 b 并保持松弛参数 w=1 时，即

T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1] - dx**2*b[1:-1, 1:-1])

我得到这个情节，需要很长时间才能收敛：https ://drive.google.com/file/d/1diOWXXY1lnP9IEgJsz7Ih7mIBAhxdtBk/view?usp=sharing

最后，当我将松弛参数 w 作为 (1.1, 1.2, 1.3....1,9) 中的任何值引入时，这些值和误差会非常高，并且温度矩阵值显示为“inf”。

当我使用松弛参数时，似乎出现了问题，但我无法弄清楚为什么会这样。