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

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 

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)

Initial values After one iteration

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?


1 Answer 1


The Peaceman-Rachford ADI scheme is unconditionally stable in 2D. See e.g. [1].


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.