Need help solving heat equation using adi method

Nx = 50;

Ny = 50;

dx = 1.0 / (Nx - 1);

dy = 2.0 / (Ny - 1);

T = 0.00001;

dt = 0.000000001;

Nt = 1000; % Total number of time steps

alpha = 4; % Diffusion coefficient

x = linspace(0, 1, Nx);

y = linspace(0, 2, Ny);

[X, Y] = meshgrid(x, y);

U = Y .* X + 1; % Initial condition

function val = f(x, y, t)

val = exp(t) * cos(pi * x / 2) * sin(pi * y / 4);

end

function U = apply_boundary_conditions(U, x, y, dy, dx)

U(:, 1) = 1; % y=0

U(:, end) = U(:, end-1) + dy .* x'; % y=2

U(1, :) = U(2, :) - dx * y; % x=0

U(end, :) = y' + 1; % x=1

end

% Thomas algorithm

function x = thomas_algorithm(a, b, c, d)

n = length(d);

c_star = zeros(n-1, 1);

d_star = zeros(n, 1);

x = zeros(n, 1);

c_star(1) = c(1) / b(1);

d_star(1) = d(1) / b(1);

for i = 2:n-1

temp = b(i) - a(i-1) * c_star(i-1);

c_star(i) = c(i) / temp;

d_star(i) = (d(i) - a(i-1) * d_star(i-1)) / temp;

end

d_star(n) = (d(n) - a(n-1) * d_star(n-1)) / b(n);

x(n) = d_star(n);

for i = n-1:-1:1

x(i) = d_star(i) - c_star(i) * x(i+1);

end

end

% ADI method

for n = 1:Nt

U = apply_boundary_conditions(U, x, y, dy, dx);

% First half-step: X-direction implicit, Y-direction explicit

for j = 2:Ny-1

a = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);

b = (1 + alpha * dt / dx^2) * ones(Nx, 1);

c = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);

d = U(j, 2:end-1)' + 0.5 * alpha * dt / dy^2 * (U(j+1, 2:end-1) - 2 * U(j, 2:end-1) + U(j-1, 2:end-1))' + dt * f(x(2:end-1), y(j), n*dt);

U(j, 2:end-1) = thomas_algorithm(a, b(2:end-1), c, d);

end

% Second half-step: Y-direction implicit, X-direction explicit

for i = 2:Nx-1

a = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);

b = (1 + alpha * dt / dy^2) * ones(Ny, 1);

c = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);

d = U(2:end-1, i) + 0.5 * alpha * dt / dx^2 * (U(2:end-1, i+1) - 2 * U(2:end-1, i) + U(2:end-1, i-1)) + dt * f(x(i), y(2:end-1), n*dt);

U(2:end-1, i) = thomas_algorithm(a, b(2:end-1), c, d);

end

U = apply_boundary_conditions(U, x, y, dy, dx);

end

% Visualization

surf(X, Y, U);

title('ADI Solution at T=' + string(T));

xlabel('X');

ylabel('Y');

zlabel('U');

colorbar;

[Torsten]

Shouldn't your solution array be three-dimensional instead of two-dimensional ? The way you arranged the code, you only get one solution at time T - the complete history for U is overwritten.

Answer 1

[John D'Errico]

You say it works for sufficiently small values dt.

With that exp(t) term in there, do you seriously expect it to work well for large values of dt? I have dreams myself somedays...

[Torsten]

T = 0.00001;
dt = 0.000000001;

That's really small ...

[Vard]

and how can I ensure stability?

[Torsten]

Since you seem to use an implicit time integration method, you should be able to use larger timesteps without getting stability issues. So I'm quite sure there must be something wrong with your code.

[Vard]

the question is that there shouldn't be a stability problem, but I don't know how to fix the code

[Vard]

can you help me to write the solution with three-dimensional array ?

[Vard]

or maybe you have examples?

[John D'Errico]

There may easily be problems in your code. For example, there is no need to write a Thomas algorithm to solve a simple banded linear system. You make your code more complex, easier to fail due to a bug when you do things like that.

But debugging your code is difficult for others. We need to figure out exactly what you are doing from the few comments we see, variable names that don't make a huge amount of sense, etc.

Regardless, if your code works for smaller dt, but not larger dt, then it is by far most likely to be due to instability.

[Vard]

but it should work for any dt value yes?

I need a ADI scheme's solution example

[Torsten]

Isn't ADI iterative in each time step ? I don't see any iterations in your code.

If you write down the mathematical algorithm for the above PDE, it would help to understand your code.

[Vard]

i can also send python code

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

Nx, Ny = 50, 50

dx, dy = 1.0 / (Nx - 1), 2.0 / (Ny - 1)

T = 0.00001

dt = 0.000000001

Nt = 1000 # Total number of time steps

alpha = 4 # Diffusion coefficient

x = np.linspace(0, 1, Nx)

y = np.linspace(0, 2, Ny)

t = np.linspace(0, T, Nt+1)

X, Y = np.meshgrid(x, y)

U = np.outer(y, x) + 1 # Initial condition

def f(x, y, t):

return np.exp(t) * np.cos(np.pi * x / 2) * np.sin(np.pi * y / 4)

def apply_boundary_conditions(U):

U[:, 0] = 1 # y=0

U[:, -1] = U[:, -2] + dy * x # y=2

U[0, :] = U[1, :] - dx * y # x=0

U[-1, :] = y + 1 # x=1

return U

# Thomas algorithm

def thomas_algorithm(a, b, c, d):

n = len(d)

c_star = np.zeros(n-1)

d_star = np.zeros(n)

x = np.zeros(n)

c_star[0] = c[0] / b[0]

d_star[0] = d[0] / b[0]

for i in range(1, n-1):

temp = b[i] - a[i-1] * c_star[i-1]

c_star[i] = c[i] / temp

d_star[i] = (d[i] - a[i-1] * d_star[i-1]) / temp

d_star[n-1] = (d[n-1] - a[n-2] * d_star[n-2]) / b[n-1]

x[-1] = d_star[-1]

for i in range(n-2, -1, -1):

x[i] = d_star[i] - c_star[i] * x[i+1]

return x

# ADI method

for n in range(Nt):

U = apply_boundary_conditions(U)

# First half-step: X-direction implicit, Y-direction explicit

for j in range(1, Ny - 1):

a = np.full(Nx - 1, -alpha * dt / (2 * dx**2))

b = np.full(Nx, 1 + alpha * dt / dx**2)

c = np.full(Nx - 1, -alpha * dt / (2 * dx**2))

d = U[j, 1:-1] + 0.5 * alpha * dt / dy**2 * (U[j + 1, 1:-1] - 2 * U[j, 1:-1] + U[j - 1, 1:-1]) +dt * f(x[1:-1], y[j], n*dt)

U[j, 1:-1] = thomas_algorithm(a, b[1:-1], c, d)

# Second half-step: Y-direction implicit, X-direction explicit

for i in range(1, Nx - 1):

a = np.full(Ny - 1, -alpha * dt / (2 * dy**2))

b = np.full(Ny, 1 + alpha * dt / dy**2)

c = np.full(Ny - 1, -alpha * dt / (2 * dy**2))

d = U[1:-1, i] + 0.5 * alpha * dt / dx**2 * (U[1:-1, i + 1] - 2 * U[1:-1, i] + U[1:-1, i - 1])+dt * f(x[i], y[1:-1], n*dt)

U[1:-1, i] = thomas_algorithm(a, b[1:-1], c, d)

U = apply_boundary_conditions(U)

# Visualization

fig = plt.figure(figsize=(14, 8))

ax = fig.add_subplot(111, projection='3d')

surf = ax.plot_surface(X, Y, U, cmap='viridis', edgecolor='none')

ax.set_title('ADI Solution at T={}'.format(T))

ax.set_xlabel('X')

ax.set_ylabel('Y')

ax.set_zlabel('U')

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

[Vard]

i write Thomas algorithm for tridiagonal matrix solutions.

then loops`

For each time step n (0 to Nt - 1):

1. Apply the boundary conditions to the solution grid U.
2. First Half-Step (X-Direction Implicit):

• For each row j from 1 to Ny-2:
• Set up tridiagonal matrix coefficients for the X-direction implicit solution:
• Diagonal and off-diagonal coefficients (a, b, c).
• Right-hand side (d) based on explicit Y-direction updates and the source term.
• Solve the tridiagonal system using the Thomas algorithm.
• Update row values in the solution grid U.

1. Second Half-Step (Y-Direction Implicit):

• For each column i from 1 to Nx-2:
• Set up tridiagonal matrix coefficients for the Y-direction implicit solution:
• Diagonal and off-diagonal coefficients (a, b, c).
• Right-hand side (d) based on explicit X-direction updates and the source term.
• Solve the tridiagonal system using the Thomas algorithm.
• Update column values in the solution grid U.

[Vard]

@Torsten @John D'Errico please help me if you can

[Torsten]

I'm not sure about which implicit time-integration method you use and how the boundary conditions and the source term must be included in the course of the two steps. So I cannot help you until you give me a link to a document where these questions are answered.

[Vard]

i don't have any document, i just have the Heat equation and should solve it using ADI scheme

[Vard]

also using three-dimensional array

[Vard]

@Torsten @John D'Errico don't you have a sample of solution that you can share with me?

[Torsten]

No. ADI for a parabolic partial differential equation in two spatial coordinates is not a standard solution method for this kind of problem.

[Vard]

@Torsten Do you have a another solution with using three dimensional arrays?

[Torsten]

If your code were correct, you could simply save the solution matrix U after each time step in a three-dimensional matrix:

U_3d = zeros(Nt,Ny,Nx)
for nt = 1:Nt
    ...
   U_3d(nt,:,:) = U;     
end