3
$\begingroup$

I am trying to solve the Poisson Equation

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

for a 61x61 grid using Python3 with boundary conditions being $T=0$ on all four boundaries and taking an initial guess of $T=50$. I am considering a relaxation parameter $w=1.6$.

To solve the Poisson Equation using SOR technique, I discretize using the Finite Difference method and do:

$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)$

I would like to find the error between the iterations by simply finding the difference between the older iteration and the newer one (by finding the maximum difference in the array), and checking if that error is lesser than my permissible error ($10^{-7}$). If it is lesser, the iterations stop, otherwise they continue.

However, this method doesn't seem to converge even after 2000 iterations. I would like to know where I am making the mistake (I suspect it is in the error calculation). Why doesn't this simple error calculation work? What changes can I implement to make it work?

Here is my code:

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()

And this is my output:

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]]

Notice how it goes until nearly 2000 iterations and stops only because I explicitly ask the loop to break at 1999 iterations and the error is higher than specified ($10^{-7}$). This is a link to the plot since I am unable to directly paste it here:

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

Thank you.

Edit: I took the suggestions in the comments.

I tweaked the source term and simplified the equation, so instead of

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]

I did:

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])

When I make the source term 0 and take the relaxation parameter as 1, i.e.

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

I get the correct plot and the solution converges at 349 iterations: https://drive.google.com/file/d/1IFqzcjkCxmlQjJwGdfwp7ACZVsP9Dd-Q/view?usp=sharing

When I use the source term b and keep the relaxation parameter w=1, i.e.

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])

I get this plot, which takes a long time to converge: https://drive.google.com/file/d/1diOWXXY1lnP9IEgJsz7Ih7mIBAhxdtBk/view?usp=sharing

Finally, when I introduce the relaxation parameter w as any value in (1.1, 1.2, 1.3....1,9) the values and the errors shoot up really high and the temperature matrix values shows 'inf'.

It seems like the problem occurs when I use the relaxation parameter, but I cannot figure out why it's happening.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ You should use some simple right-hand side, e.g., rhs=0 or 1 (so that an exact analytic solution is possible for checking), and see if your numerical solution is going in the right direction. $\endgroup$
    – Maxim Umansky
    Oct 20, 2020 at 4:14
  • $\begingroup$ @MaximUmansky I will do that, thank you. But does my error term seem right? $\endgroup$
    – user37306
    Oct 20, 2020 at 4:20
  • $\begingroup$ You are using one possible metric for the error term. There are several other metrics; but if a solution (assuming it is reasonable - smooth and finite) converges by one metric then it converges by any other. So yes, your measure of convergence looks right, I am not seeing anything obviously wrong there. But rather than staring at the code, a much more efficient techniques for debugging problems of this kind is to apply it an exactly solvable problem; then you'll see immediately where the problem is. $\endgroup$
    – Maxim Umansky
    Oct 20, 2020 at 4:43
  • 2
    $\begingroup$ Actually, something seems wrong with the SOR formula used, with RHS=0 it should be just $T_{i,j}^{n+1}$ equals the average of the four surrounding values at the previous step, but the formula used here seems to have some extra factors. Just try it for dx=dy and $\omega$=1. $\endgroup$
    – Maxim Umansky
    Oct 20, 2020 at 4:55
  • $\begingroup$ @MaximUmansky thanks for the suggestions! I implemented them and updated my question, have a look. So the initial SOR formula I used was a general one, but I have updated it now to the simpler average of 4 surrounding values formula that you suggested. $\endgroup$
    – user37306
    Oct 20, 2020 at 13:10

1 Answer 1

Reset to default
2
$\begingroup$

This formula is correct on uniform grid (even for $\omega \neq $1 although convergence depends on $\omega$):

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

This formula is not correct (even for $\omega = $1), the average over old values has a factor $1/dx**2$ which should not be there:

$ 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] $

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for the answer, but what about the source term in the first equation of your answer? And I changed the formula to (I have added this in the question as well): 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*(Tn[1:-1, 2:] + Tn[1:-1, :-2] + Tn[2:, 1:-1] + Tn[:-2, 1:-1] - dx**2*b[1:-1, 1:-1]) And it didn't work after that either. $\endgroup$
    – user37306
    Oct 20, 2020 at 16:31
  • $\begingroup$ I used the equation in your question, after the words "so instead of", and I showed there is a mistake in the equation, presumably that resolves the problem. There is no source term in that equation that you wrote, am I wrong? In the first equation which is correct on uniform grid, if there was a source term S then it would be $T_{i,j}^{new} = (1/4) \sum_{i,j} T_{i\pm1,j\pm1} - S_{i,j} dx^2$. $\endgroup$
    – Maxim Umansky
    Oct 20, 2020 at 17:03
  • $\begingroup$ Yes, but in the question edit added after your initial comments on the question, I tried keeping source=0 and w=1 and the equation worked correctly. Then, introduced source term and w=1, and still got the solution. The problem is when there is a source and w is not 1. In the edit, the equation I used is the same as the first equation in your answer (or am I missing something?) $\endgroup$
    – user37306
    Oct 20, 2020 at 17:42
  • $\begingroup$ So, the reason I am not able to figure out where I am making a mistake is that I get the correct answers (I think) for S=0 and w=1 case, and $S \neq 0$ and w=1 case as well. If you'd like to look at these answers, the plots are linked in the question edit. So if these answers are correct, why is it that simply introducing a factor of w and (1-w) where $w \neq 1$ makes the values and the errors go to infinity? Can't seem to figure this out. But I appreciate your continued help! $\endgroup$
    – user37306
    Oct 20, 2020 at 18:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.