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:MxUxianM(x,1)=InitialConditions((x-1)*h);
end
for t=1:NtUxianM(1,t)=BoundaryConditions(0,(t-1)*k);UxianM(Mx,t)=BoundaryConditions(1,(t-1)*k);
end
%利用显格式公式逐行递推
for t=2:Ntfor x=2:Mx-1UxianM(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:NtUxian= [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-1temp{t}=C;
end
mid = repmat({C},Nt-1,Nt-1);%对角块
for x=1:Nt-1for t=1:Nt-1if x~=tmid{x,t}=zeros(Mx-2,Mx-2); %非主对角线置0endif x==t+1mid{x,t}=D; %下次对角线上的块为Dendend
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-3f1(x)=InitialConditions(h*x);
end
fyin(1:Mx-2)=f1;%计算边值条件计算f2~fNt-1
for t=2:Nt-1fyin((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:MxU0(x)=InitialConditions((x-1)*h);
endUcnM=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-1n=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-1Ucn= [Ucn;UcnM(:,t)];
end%%
%计算精确解，用于误差分析****************************************************
U_M=zeros(Mx-2,Nt-1);
U=[];
for x=1:Mx-2for t=1:Nt-1 U_M(x,t)=ExactSolution(x*h,t*k);end
end
%将结果按Ku=f方法的u的结构重排
for t=1:Nt-1U=[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); endfunction [ Uxt ] = InitialConditions ( x )
%在此函数中定义初值条件Uxt=sin(pi*x);endfunction [ Uxt ] = ExactSolution( x,t)
%在此函数中定义待求函数的精确解u(x,t),用于比较数值解Uxt=exp(-1*pi*pi*t)*sin(pi*x);endfunction [ Uxt ] = BoundaryConditions( LeftOrRight,t )
%在此函数中定义边值条件
%LeftOrRight标识左边界条件还是右边界条件 0代表左，1代表右Uxt=-1;
if LeftOrRight==0Uxt=0;
elseif LeftOrRight==1Uxt=0;
else%不应该出现此情况，抛异常error('错误使用边界条件!');
endend
