热传导方程的差分格式原理与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