偏微分的第4次实验，主要内容为二维抛物型方程的ADI格式求解。不足之处望读者多多指正。

使用ADI格式求解以下问题： ∂ u ∂ t − ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) = − 3 2 e 1 2 ( x + y ) − t , 0 < x , y < 1 , 0 < t ⩽ 1 \frac{\partial u}{\partial t}-\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)=-\frac{3}{2} e^{\frac{1}{2}(x+y)-t}, \quad 0C1}; %创建下对角矩阵块 A1=cell2mat(mid); %*************************************************************** %构造B2 bb1=zeros(1,M-2); for i=1:M-2 %给-1,1对角线位置赋值. bb1(i)=r./2; end bb2=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 bb2(i)=1-r; end B2=diag(bb1,-1)+diag(bb2)+diag(bb1,1); %*************************************************************** %构造A2 aa2=blkdiag(kron(eye(M-1),B2)); %创建主对角矩阵块 mid=mat2cell(aa2,(M-1)*ones(1,M-1),(M-1)*ones(1,M-1)); A2=cell2mat(mid); %*************************************************************** %构造B3 bb4=zeros(1,M-2); for i=1:M-2 %给-1,1对角线位置赋值. bb4(i)=-r./2; end bb5=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 bb5(i)=1+r; end B3=diag(bb4,-1)+diag(bb5)+diag(bb4,1); %*************************************************************** %构造A3 aa3=blkdiag(kron(eye(M-1),B3)); %创建主对角矩阵块 mid=mat2cell(aa3,(M-1)*ones(1,M-1),(M-1)*ones(1,M-1)); A3=cell2mat(mid); %*************************************************************** %构造B4 bb6=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 bb6(i)=1-r; end B4=diag(bb6); %*************************************************************** %构造C4 cc=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 cc(i)=r./2; end C4=diag(cc); %*************************************************************** %构造A4 aa4=blkdiag(kron(eye(M-1),B4)); %创建主对角矩阵块 mid=mat2cell(aa4,(M-1)*ones(1,M-1),(M-1)*ones(1,M-1)); mid(find(diag(ones(1,M-2),1)==1))={C4}; %创建上对角矩阵块 mid(find(diag(ones(1,M-2),-1)==1))={C4}; %创建下对角矩阵块 A4=cell2mat(mid); %*************************************************************** %Fk+1/2, 先构造f1 ff1=zeros(1,M-2); for i=1:M-2 %给-1,1对角线位置赋值. ff1(i)=r./4; end ff2=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 ff2(i)=(1-r)./2; end f1=diag(ff1,-1)+diag(ff2)+diag(ff1,1); %*************************************************************** %再构造f2 ff1=zeros(1,M-2); for i=1:M-2 %给-1,1对角线位置赋值. ff1(i)=-r./4; end ff2=zeros(1,M-1); %给主对角线位置赋值. for i=1:M-1 ff2(i)=(1+r)./2; end f2=diag(ff1,-1)+diag(ff2)+diag(ff1,1); %*************************************************************** U0=zeros((M-1)*(M-1),1); count1=1;%记录第0层总个数 % 初值条件如下，给第0层赋值 for i=1:M-1 %不包括边界值 for j=1:M-1 U0(count1,1)=exp(0.5*(x(i+1)+y(j+1))-t(1)); count1=count1+1; end end %*************************************************************** %先求k+1/2层 k=1;%记录行数 FF1=zeros(M-1,1); FF2=zeros(M-1,1); uk=zeros(M-1,1);%u(0,j,k)|u(M,j,k) uk1=zeros(M-1,1);%u(0,j,k+1)|u(M,j,k+1) uk_1=zeros(M-1,1);%u(0,0/M,k)|u(M,0/M,k) uk1_1=zeros(M-1,1);%u(0,0/M,k+1)|u(M,0/M,k+1) F1=zeros((M-1)*(M-1),1); while k~=N+1 %*************************************************************** %构造F(k+1/2) %FF1 for i=1:M-1 uk(i,1)=exp(0.5*(y(i+1)-t(k))); uk1(i,1)=exp(0.5*(y(i+1)-t(k+1))); end uk_1(1,1)=exp(1*t(k)); uk_1(M-1,1)=exp(0.5-t(k)); uk1_1(1,1)=exp(t(k+1)); uk1_1(M-1,1)=exp(0.5+t(k+1)); FF1=f1*uk+r./4*uk_1+f2*uk1-r./4*uk1_1; %*************************************************************** %FF2 for i=1:M-1 uk(i,1)=exp(0.5*(1+y(i+1)+t(k))); uk1(i,1)=exp(0.5*(1+y(i+1)+t(k+1))); end uk_1(1,1)=exp(0.5+t(k)); uk_1(M-1,1)=exp(1+t(k)); uk1_1(1,1)=exp(0.5+t(k+1)); uk1_1(M-1,1)=exp(1+t(k+1)); FF2=f1*uk+r./4*uk_1+f2*uk1-r./4*uk1_1; % 在下面把FF1全部赋给F1即F(k+1/2)向量的前M-1项 for j=1:M-1 F1(j,1)=FF1(j,1); end % 在下面把FF2全部赋给F1即F(k+1/2)向量的最后的M-1项 pp=1; for j=(M-1)*(M-2)+1:(M-1)*(M-1) F1(j,1)=FF2(pp,1); %Fk+1/2 pp=pp+1; end %*************************************************************** % 构造F(k) F=zeros((M-1)*(M-1),1);%Fk c=1; for j=1:M-1 F(c,1)=exp(0.5*(x(j+1)+t(k))); F(c+M-2,1)=exp(0.5*(1+x(j+1)+t(k))); c=c+M-1; end f=A2*U0+0.5.*r*F1+0.5.*r*F; U1=A1\f; %U1为第一阶段利用U(k)所求的U(k+1/2); %*************************************************************** % 继续往下求解, 利用U(k+1/2)求的U(k+1), 此为第二阶段 % U(k+1/2)已存放到U1中, F(k+1/2)已存放到F1之中, 在第二阶段需要使用两者; ff6=zeros((M-1)*(M-1),1); %Fk+1 % ff6为所需要的F(k+1); c=1; for j=1:M-1 ff6(c,1)=exp(0.5*(x(j+1)+t(k+1))); ff6(c+M-2,1)=exp(0.5*(1+x(j+1)+t(k+1))); c=c+M-1; end FF=A4*U1+0.5*r*(ff6+F1); U=A3\FF; %*************************************************************** % U为所求的U(k+1) U0=U; % 把U(k+1)作为下一次迭代的初值条件 U=reshape(U,M-1,M-1); for i=2:M for j=2:M U3(i,j,k)=U(i-1,j-1); end end for i=1:M+1 U3(i,1,k)=exp(0.5*(x(i))-t(k+1)); %左边界 U3(i,M+1,k)=exp(0.5*(1+x(i))-t(k+1)); %右边界 end for i=1:M+1 U3(1,i,k)=exp(0.5*(y(i)-t(k+1))); %上边界 U3(M+1,i,k)=exp(0.5*(1+y(i))-t(k+1)); %下边界 end k=k+1; %进入下一层 end 调用的主函数 主要实现：完成迭代矩阵U3与解析解矩阵U4的实现，主要实现图像拟合、与误差分析 求解结果（目标） 拟合图像

完成了matlab拟合与可视化的部分，误差分析并计算误差阶有待下篇文章完成。

物理老师的热传导课