Need help solving heat equation using adi method

11 次查看（过去 30 天）

显示更早的评论

Vard 2024-5-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2117816-need-help-solving-heat-equation-using-adi-method

编辑： Vard 2024-5-16

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

% 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 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 2024-5-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2117816-need-help-solving-heat-equation-using-adi-method#answer_1456131

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更早的评论

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 2024-5-13

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

Torsten 2024-5-13

编辑：Torsten 2024-5-13

在 MATLAB Online 中打开

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

请先登录，再进行评论。

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

类别

Mathematics and Optimization Partial Differential Equation Toolbox Domain-Specific Modeling Structural Mechanics

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by