ADI交替差分隐格式求解二维热传导方程

ybutvdy

已于 2023-12-08 23:44:59 修改

阅读量6.4k

点赞数 12

分类专栏：计算数学文章标签：偏微分方程

于 2021-10-06 22:15:41 首次发布

本文链接：https://blog.csdn.net/undergraduate_/article/details/120630381

版权

计算数学专栏收录该内容

1 篇文章 1 订阅

订阅专栏

博客介绍了使用ADI(Alternating Direction Implicit)方法在MATLAB中解决二维热传导初边值问题的过程，展示了MATLAB代码实现，并给出了Neumann边界条件下使用FFT快速求解器的Python实现。内容涵盖了热传导方程、边界条件设定、MATLAB的隐式差分格式以及快速傅里叶变换(FFT)在解决边值问题的应用。

摘要由CSDN通过智能技术生成

ADI算法

在这里插入图片描述

初边值问题叙述

满足定解条件的齐次方程:

$\frac{{\partial u}}{{\partial t}} = {a^2}\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right) \tag{1}$
$\le x,y \le l;t > 0 \\ u{|_{t = 0}} = \varphi \left( {x,y} \right) \\ \forall t >0:u\left( {0,0,t} \right) = u\left( {0,l,t} \right) = u\left( {l,0,t} \right) = u\left( {l,l,t} \right) = 0$
以及为了求出比较美观的解的一个边值条件:
$\forall t>0:u\left( {\frac{l}{4},\frac{l}{4},t} \right) = u\left( {\frac{l}{4},\frac{{3l}}{4},t} \right) = u\left( {\frac{{3l}}{4},\frac{l}{4},t} \right) = u\left( {\frac{{3l}}{4},\frac{{3l}}{4},t} \right) = 1$

MATLAB实现

%% ADI
M=200;%网格数量
N=2000;%时间分割
T=0.05;%总时间
number = 50;%初始点数
tau=T./N;
L=1;
a=1;%导热系数决定的
h=L./M;
r=tau./h./h./2.*a.^2;%网比
disp(r);
H=sparse(diag(ones(1,M)*(2.*r+1)))+sparse(diag(ones(1,M-1)*(-r),1))+sparse(diag(ones(1,M-1)*(-r),-1));
v=zeros(M,M);
u=zeros(M,M);%初值
u(1,:)=0;u(M,:)=0;u(:,1)=0;u(:,M)=0;
[X,Y] = meshgrid(h:h:L,h:h:L);
RP = [randperm(M-2,number)'+1,randperm(M-2,number)'+1];%初始温度为100.0的点
for iter = 1:number
u(RP(iter,1),RP(iter,2)) = 100.0;
end
for k = 2:N
u(1,:)=0;u(M,:)=0;u(:,1)=0;u(:,M)=0;
delta_y=([zeros(M,1),u(:,1:M-1)]+[u(:,2:M),zeros(M,1)]-2*u)*r;
v=H\(u+delta_y);
v=v';
delta_x=([zeros(M,1),v(:,1:M-1)]+[v(:,2:M),zeros(M,1)]-2*v)*r;
u=H\(v+delta_x);
u=u';
u(1,:)=0;u(M,:)=0;u(:,1)=0;u(:,M)=0;
colormap(hot(10));
    contourf(h:h:L, (h:h:L)', u, 'ShowText','on');
    colorbar
    title('热传导初边值问题');
    hold off;
     frame=getframe(gcf);
    imind=frame2im(frame); [imind,cm] = rgb2ind(imind,256);
    if k==1
       imwrite(imind,cm,'test.gif','gif', 'Loopcount',inf,'DelayTime',1e-4);
    else
        if mod(k,100)==1
       imwrite(imind,cm,'test.gif','gif','WriteMode','append','DelayTime',1e-4);
       fprintf('epoch=%d\n',k);
        end
    end
end

在这里插入图片描述

这里是向单位正方形中滴入50个温度为100的点的散热建模
$l = 1, T = 0.05, a = 1$
下的解。

Neumman边界解:

对于边值问题:
$\Delta u = f,\frac{{\partial u}}{{\partial {\bf{n}}}} = 0;$ 在乘积型区域上的解,使用FFT方法可以得到快速求解器,可见方腔驱动流中压力泊松方程的做法:

https://qiita.com/eigs/items/628a2aeb3cb9eef91093

我这里写了一个helmholtz方程使用numpy求解器的零法向导数边值条件的快速求解器:
$\left( {\lambda I - \Delta } \right)u = f,\frac{{\partial u}}{{\partial {\bf{n}}}} = 0;$
因为热方程的隐式差分格式一定能化成这个格式(很显然):

import numpy as np


def fft_neumann_poisson(b, h):
    m, n = b.shape

    row = n * np.cos(np.pi * np.arange(0, n) / n)
    row[0] = row[0] + n
    col = m / 2 * np.ones(m)
    col[0] = col[0] + m / 2
    rc = np.outer(row, col)
    row = n / 2 * np.ones(n)
    row[0] = row[0] + n / 2
    col = m * np.cos(np.pi * np.arange(0, m) / m) - 2 * m
    col[0] = col[0] - m
    rc = rc + np.outer(row, col)
    
    y1 = np.fft.fft(np.concatenate([b, np.zeros_like(b, dtype=np.complex128)], axis=0),axis = 0)
#     print(np.concatenate([b, np.zeros_like(b, dtype=np.complex128)], axis=0))
#     print(y1.round(4))
    y1 = np.real(y1[:m, :] * (np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) - np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) * 1j).reshape(-1,1))
    y1 = np.fft.fft(np.concatenate([y1.T.conj(), np.zeros_like(y1.T, dtype=np.complex128)], axis=0),axis = 0)
    y1 = np.real(y1[:n, :] * (np.cos(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) - np.sin(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) * 1j).reshape(-1,1))
#     print(rc.round(4))
    y1[0, 0] = 0
    y1 = (y1 / rc).T.conj()
    y1[0, 0] = 0
#     print(y1.shape)
#     print(np.concatenate([ ((np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) ) - (np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m))) * 1j)*y1, np.zeros((m, n))], axis=0).round(2))
    y2 = np.real( np.fft.fft(np.concatenate([ ((np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) ) - (np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m))) * 1j)*y1, np.zeros((m, n))], axis=0),axis = 0))
#     print(y2.shape)
#     print(y2.round(4))
    y2 = y2[:m,:].T.conj()
#     print(y2.round(4))
#     print(y2.shape)
    y2 = np.real( np.fft.fft(np.concatenate([(np.cos(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) )*y2 -(np.sin(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) )*y2*1j, np.zeros((n, m))], axis=0),axis = 0))

    re = (h**2) * y2[:n, :].T.conj()

    return re

def fft_neumann_helmholtz(b,lamb,h):
        m, n = b.shape

        row = n * np.cos(np.pi * np.arange(0, n) / n)
        row[0] = row[0] + n
        col = m / 2 * np.ones(m)
        col[0] = col[0] + m / 2
        rc = np.outer(row, col)
        row = n / 2 * np.ones(n)
        row[0] = row[0] + n / 2
        col = m * np.cos(np.pi * np.arange(0, m) / m) - 2 * m
        col[0] = col[0] - m
        rc = rc + np.outer(row, col)
        row = n / 2 * np.ones(n)
        row[0] = n
        col = m / 2 *np.ones(m)
        col[0] = m
        rc = -1*rc + lamb*(h**2)*np.outer(row, col)

        y1 = np.fft.fft(np.concatenate([b, np.zeros_like(b, dtype=np.complex128)], axis=0),axis = 0)
    #     print(np.concatenate([b, np.zeros_like(b, dtype=np.complex128)], axis=0))
    #     print(y1.round(4))
        y1 = np.real(y1[:m, :] * (np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) - np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) * 1j).reshape(-1,1))
        y1 = np.fft.fft(np.concatenate([y1.T.conj(), np.zeros_like(y1.T, dtype=np.complex128)], axis=0),axis = 0)
        y1 = np.real(y1[:n, :] * (np.cos(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) - np.sin(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) * 1j).reshape(-1,1))
    #     print(rc.round(4))
        y1[0, 0] = 0
        y1 = (y1 / rc).T.conj()
        y1[0, 0] = 0
    #     print(y1.shape)
    #     print(np.concatenate([ ((np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) ) - (np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m))) * 1j)*y1, np.zeros((m, n))], axis=0).round(2))
        y2 = np.real( np.fft.fft(np.concatenate([ ((np.cos(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m)) ) - (np.sin(np.arange(0, m).reshape(-1, 1) * np.pi / (2 * m))) * 1j)*y1, np.zeros((m, n))], axis=0),axis = 0))
    #     print(y2.shape)
    #     print(y2.round(4))
        y2 = y2[:m,:].T.conj()
    #     print(y2.round(4))
    #     print(y2.shape)
        y2 = np.real( np.fft.fft(np.concatenate([(np.cos(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) )*y2 -(np.sin(np.arange(0, n).reshape(-1, 1) * np.pi / (2 * n)) )*y2*1j, np.zeros((n, m))], axis=0),axis = 0))
 
        return y2[:n, :].T.conj()

用来求解效果是这样的:

import matplotlib.pyplot as plt
x = np.linspace(0, 1, 224)
y = np.linspace(0, 1, 224)
h = x[1]-x[0]  # 差分步长
lamb = 1/(h**2)
x, y = np.meshgrid(x, y)

# 计算函数值
z = np.sin(np.pi * x) * np.cos(np.pi * y)
# print(x[1]-x[0])
# rep = fft_neumann_helmholtz(z,lamb,h)
rep = z.copy()
for k in range(1000):
    rep = fft_neumann_helmholtz(rep,lamb,h)
    print("{:d}".format(k+1),end = ',')
    if (k+1) % 50 == 0:
        plt.contourf(x, y, rep, cmap='viridis')  # 用contourf函数绘制填充等高线图
        plt.colorbar()  # 添加颜色条
        plt.title('time = {:.4f},initial f(x, y) = sin(pi*x)*cos(pi*y)'.format((k+1)*h**2))
        plt.xlabel('x')
        plt.ylabel('y')
        plt.show()
        print()

在这里插入图片描述