I am writing an Alternating-Directions Implicit Method for simple 2D diffusion ( \begin{equation*} \frac{df(x,y,t)}{dt}=D\Delta u \end{equation*}). Tridiagonal matrices are solved via Thomas algorithm (LU decomposing).
def calculateADI(u,delta_t, delta_x, size, D):
a1 = np.ones((size), dtype=np.float64)
b1 = np.zeros((size-1), dtype=np.float64)
c1 = np.zeros((size-1), dtype=np.float64)
q1 = np.zeros((size), dtype=np.float64)
a2 = np.ones((size), dtype=np.float64)
b2 = np.zeros((size-1), dtype=np.float64)
c2 = np.zeros((size-1), dtype=np.float64)
q2 = np.zeros((size), dtype=np.float64)
y0 = np.copy(u)
y = np.copy(y0)
gamma = D
# first direction
for j in range(1,size-1):
for i in range(0,size):
a1[i] = 2./ (delta_t) + 2.*gamma/(delta_x)**2
if i != (size-1):
b1[i] = -1*gamma/(delta_x)**2
c1[i] = -1*gamma/(delta_x)**2
if (j==0):
q2[i] = 2.*y0[i,j]/(delta_t) +gamma*(y0[i,j+1]-2.*y0[i,j])/delta_x**2
elif (j==(size-1)):
q2[i] = 2.*y0[i,j]/(delta_t) +gamma*(y0[i,j-1]-2.*y0[i,j])/delta_x**2
else:
q1[i] = 2.*y[i,j]/(delta_t) + gamma*(y[i,j+1]-2.*y[i,j]+y[i,j-1])/delta_x**2
y0[:,j] = solveLU3(a1,b1,c1,q1)
# second direction
for i in range(1,size-1):
for j in range(0,size):
a2[j] = 2./ (delta_t) + 2.*gamma/delta_x**2
if j != (size-1):
b2[j] =-1.*gamma/delta_x**2
c2[j] = -1.*gamma/delta_x**2
if (i==0):
q2[j] = 2.*y0[i,j]/(delta_t) +gamma*(y0[i+1,j]-2.*y0[i,j])/delta_x**2
elif (i==(size-1)):
q2[j] = 2.*y0[i,j]/(delta_t) +gamma*(y0[i-1,j]-2.*y0[i,j])/delta_x**2
else:
q2[j] = 2.*y0[i,j]/(delta_t) +gamma*(y0[i-1,j]-2.*y0[i,j]+y0[i+1,j])/delta_x**2
y[i,:] = solveLU3(a2,b2,c2,q2)
return y, y0 # returns both half-iterations for more information
The decomposition and solution are as follows:
@numba.njit
def decLU3(a,b,c):
n = len(a)
l = np.zeros(len(a))
l[0] = a[0]
u = np.zeros(len(b))
u[0] = b[0]/a[0]
for i in range(1,n-1):
l[i] = a[i]-c[i-1]*u[i-1]
u[i] = b[i]/l[i]
l[n-1] = a[n-1]-c[n-2]*u[n-2]
return u,l
@numba.njit
def solveLU3(a,b,c,f):
n = len(a)
u, l = decLU3(a, b, c)
z = np.copy(f)
z[0] = f[0]/l[0]
for i in range(1,n):
z[i] = (f[i]-c[i-1]*z[i-1])/l[i]
w = np.copy(z)
w[n-1] = z[n-1]
for i in range(n-2,-1,-1):
w[i] = (z[i] - u[i]*w[i+1])
return w
As an example, this is what one iteration does to a plate that has a 50x50 area with values 1 and outsides are of value 0.
uk = u = np.zeros((100,100))
u[25:75,25:75] = 1
uk, uk0 = calculateADI(uk, 1, 0.2, len(u),1)
As you can see, white lines appear that, from my understanding, are negatives, and therefore are not in range of plot values. In addition, it was found out by me, that with a $\Delta x = 1$ the program seems to do well.
So does ADI method actually have some kind of condition for stability?
And if not, what could be the problem with algorithm?
