热传导方程的差分格式原理与matlab实现

function [  ] = ParabolicEquation( h,k )
%求解抛物型方程中的一种：热传导方程
%h:x轴步长
%k:t轴步长

r=k/(h*h);%网格比
Mx=floor(1.0/h)+1;%网格在x轴上的节点个数（算上0）
Nt=floor(1.0/k)+1;%网格在t轴上的节点个数（算上0）
N=(Mx-2)*(Nt-1); %U的维数

%%
%********************************古典显格式*********************************
%直接递推的方法求古典显格式

UxianM=zeros(Mx,Nt);
Uxian=[];
%先赋初值和边界值
for x=1:Mx
    UxianM(x,1)=InitialConditions((x-1)*h);
end
for t=1:Nt
    UxianM(1,t)=BoundaryConditions(0,(t-1)*k);
    UxianM(Mx,t)=BoundaryConditions(1,(t-1)*k);
end
%利用显格式公式逐行递推
for t=2:Nt
   for x=2:Mx-1
       UxianM(x,t)=r*UxianM(x-1,t-1)+(1-2*r)*UxianM(x,t-1)+r*UxianM(x+1,t-1);
   end
end
%将结果按Ku=f方法的u的结构重排
for t=2:Nt
   Uxian= [Uxian;UxianM(2:Mx-1,t)];
end


%%
%************************古典隐格式*****************************************

%求K,将K看出三对角块矩阵，形如：
% C 0           ...
% D C 0         ...
% 0 D C 0       ...
% 0 0 D C 0     ...
%       ...     ... 
% ...          D C

%先计算C（C是三对角矩阵）
C=eye(Mx-2)*(1+2*r);
C=C+diag(ones(1,Mx-3)*(-r),1);  %上次对角
C=C+diag(ones(1,Mx-3)*(-r),-1); %下次对角

%计算D
D=eye(Mx-2)*-1;

%计算K
temp={}; 
for t=1:Nt-1
   temp{t}=C;
end
mid = repmat({C},Nt-1,Nt-1);%对角块
for x=1:Nt-1
   for t=1:Nt-1
      if x~=t
          mid{x,t}=zeros(Mx-2,Mx-2); %非主对角线置0
      end
      if x==t+1
          mid{x,t}=D; %下次对角线上的块为D
      end
   end
end
Kyin=cell2mat(mid);

%求f
%将f分块
% f1
% f2
% ...
% fNt-1

%f中大多数值为0，非0值用初边值条件添加
fyin=zeros(N,1);

%先计算f1
f1=zeros(Mx-2,1);
%计算f（1，1），中心点为U11
f1=zeros(Mx-2,1);
f1(1)=r*BoundaryConditions(0,k)+InitialConditions(h);
%计算f（Mx-2，1），中心点为U（Mx-2，1）
f1(Mx-2)=r*BoundaryConditions(1,(Mx-1)*k)+InitialConditions(h*(Mx-2));
%计算f（x，1） (x!=1且x!=Mx-2)
for x=2:Mx-3
    f1(x)=InitialConditions(h*x);
end
fyin(1:Mx-2)=f1;

%计算边值条件计算f2~fNt-1
for t=2:Nt-1
   fyin((Mx-2)*(t-1)+1)=r* BoundaryConditions(0,k*(t));
   fyin((Mx-2)*t)=r* BoundaryConditions(1,k*(t));
end

%求u
Uyin=inv(Kyin)*fyin;

%%
%********************Crank-Nicolson格式*************************************
%追赶法求解
%An*Un+1=Bn*Un+en

%求An
An=eye(Mx-2)*(1+r);
An=An+diag(ones(1,Mx-3)*(-0.5*r),1);  %上次对角
An=An+diag(ones(1,Mx-3)*(-0.5*r),-1); %下次对角
InvA=inv(An);

%求Bn
Bn=eye(Mx-2)*(1-r);
Bn=Bn+diag(ones(1,Mx-3)*(0.5*r),1);  %上次对角
Bn=Bn+diag(ones(1,Mx-3)*(0.5*r),-1); %下次对角

Cn=InvA*Bn;

%追赶法
U0=[];%初值
for x=1:Mx
    U0(x)=InitialConditions((x-1)*h);
end

UcnM=zeros(Mx-2,Nt-1);
Ucn=[];

%求U1  An*Un+1=Bn*Un+en,e1恰好是0向量
UcnM(:,1)=Cn*U0(2:Mx-1)'+zeros(Mx-2,1);

%利用显格式公式逐行递推
for t=2:Nt-1
    n=t-1;
    en=zeros(Mx-2,1);
    en(1)=0.5*r*BoundaryConditions(0,n*k)+0.5*r*BoundaryConditions(0,(n+1)*k);
    en(Mx-2)=0.5*r*BoundaryConditions(1,n*k)+0.5*r*BoundaryConditions(1,(n+1)*k);
    UcnM(:,n+1)=Cn*UcnM(:,n)+en;
end

%将结果按Ku=f方法的u的结构重排
for t=1:Nt-1
   Ucn= [Ucn;UcnM(:,t)];
end



%%
%计算精确解，用于误差分析****************************************************
U_M=zeros(Mx-2,Nt-1);
U=[];
for x=1:Mx-2
   for t=1:Nt-1 
        U_M(x,t)=ExactSolution(x*h,t*k);
   end
end
%将结果按Ku=f方法的u的结构重排
for t=1:Nt-1
   U=[U;U_M(1:Mx-2,t)];
end

%%
%结果比较******************************************************************
%比较标准误差（均方根误差）

temp=(U-Uxian).^2;
temp=sum(sum(temp));
temp=temp/(Mx*Nt);
Exian=sqrt(temp);

temp=(U-Uyin).^2;
temp=sum(sum(temp));
temp=temp/(Mx*Nt);
Eyin=sqrt(temp);

temp=(U-Ucn).^2;
temp=sum(sum(temp));
temp=temp/(Mx*Nt);
Ecn=sqrt(temp);

fprintf('h: %d;k= %d\n',h,k);  
fprintf('古典显格式标准误差：%d\n',Exian); 
fprintf('古典隐格式标准误差：%d\n',Eyin); 
fprintf('CN格式标准误差：%d\n',Ecn); 

end

function [ Uxt ] = InitialConditions ( x )
%在此函数中定义初值条件

Uxt=sin(pi*x);

end

function [ Uxt ] = ExactSolution( x,t)
%在此函数中定义待求函数的精确解u(x,t),用于比较数值解

Uxt=exp(-1*pi*pi*t)*sin(pi*x);

end

function [ Uxt ] = BoundaryConditions( LeftOrRight,t )
%在此函数中定义边值条件
%LeftOrRight标识左边界条件还是右边界条件 0代表左，1代表右

Uxt=-1;
if LeftOrRight==0
    Uxt=0;
elseif LeftOrRight==1
    Uxt=0;
else
    %不应该出现此情况，抛异常
    error('错误使用边界条件!');
end

end

秒客网

热传导方程的差分格式原理与matlab实现

相关文章