Need help solving heat equation using adi method

2024 5 12
1 个回答
更新时间：2024 5 16
3 次查看（30 天）

0 个投票

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;

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Torsten
Torsten 2024-5-12
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2024-5-12

0 个投票

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

19 个评论
显示 17更早的评论 隐藏 17更早的评论

Torsten
Torsten 2024-5-12
T = 0.00001;
dt = 0.000000001;
That's really small ...
Vard
Vard 2024-5-12
and how can I ensure stability?
Torsten
Torsten 2024-5-12
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
Vard 2024-5-12
the question is that there shouldn't be a stability problem, but I don't know how to fix the code
Vard
Vard 2024-5-12
编辑：Vard 2024-5-12
can you help me to write the solution with three-dimensional array ?
Vard
Vard 2024-5-12
or maybe you have examples?
John D'Errico
John D'Errico 2024-5-12
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
Vard 2024-5-12
but it should work for any dt value yes?
I need a ADI scheme's solution example
Torsten
Torsten 2024-5-12
编辑：Torsten 2024-5-12
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
Vard 2024-5-12
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
Vard 2024-5-12
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.
Torsten
Torsten 2024-5-12
编辑：Torsten 2024-5-12
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
Vard 2024-5-12
i don't have any document, i just have the Heat equation and should solve it using ADI scheme
Vard
Vard 2024-5-12
also using three-dimensional array
Vard
Vard 2024-5-13
@Torsten @John D'Errico don't you have a sample of solution that you can share with me?
Torsten
Torsten 2024-5-13
No. ADI for a parabolic partial differential equation in two spatial coordinates is not a standard solution method for this kind of problem.
Vard
Vard 2024-5-13
@Torsten Do you have a another solution with using three dimensional arrays?
Torsten
Torsten 2024-5-13
编辑：Torsten 2024-5-13
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

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

提问:

Vard
Vard
2024-5-12

编辑：

Vard
Vard
2024-5-16

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by