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=5000;%时间分割
T=0.5;%总时间
tau=T./N;
L=2;
a=1;%导热系数决定的
h=L./M;
r=tau./h./h./2.*a.^2;
disp(r);
H=diag(ones(1,M)*(2.*r+1))+diag(ones(1,M-1)*(-r),1)+diag(ones(1,M-1)*(-r),-1);
v=zeros(M,M);
%result=zeros(M,M,N);
u=zeros(M,M);%初值
% subject to
u(1,:)=0;u(M,:)=0;u(:,1)=0;u(:,M)=0;
%result(:,:,1)=u;
[X,Y] = meshgrid(h:h:L,h:h:L);
for k = 2:N
u(1,:)=0;u(M,:)=0;u(:,1)=0;u(:,M)=0;
u(M./4,M./4)=1;u(3.*M./4,M./4)=1;u(M./4,3.*M./4)=1;u(3.*M./4,3.*M./4)=1;
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';
v(M./4,M./4)=1;v(3.*M./4,M./4)=1;v(M./4,3.*M./4)=1;v(3.*M./4,3.*M./4)=1;
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;
u(M./4,M./4)=1;u(3.*M./4,M./4)=1;u(M./4,3.*M./4)=1;u(3.*M./4,3.*M./4)=1;
colormap(hot(10));
    surf(h:h:L,(h:h:L)',u,'EdgeColor','none');
    %[U,V,W] = surfnorm(X,Y,u);
    hold on;
    %quiver3(X,Y,u,U,V,W,0.01,'k');
    axis([h,L,h,L]);
    grid off
    view(0,90);
    %pause(0.01);
    colorbar
    title('热传导初边值问题');
    hold off;
     frame=getframe(gcf);
    imind=frame2im(frame); [imind,cm] = rgb2ind(imind,256);
    if k==2
       imwrite(imind,cm,'test.gif','gif', 'Loopcount',inf,'DelayTime',1e-3);
    else
        if mod(k,50)==1
       imwrite(imind,cm,'test.gif','gif','WriteMode','append','DelayTime',1e-3);
       fprintf('epoch=%d\n',k);
        end
    end
end

这里是
$\\ \varphi \left( {x,y,0} \right)=0$
下的解。