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

实验内容

使用ADI格式求解以下问题：
∂u∂t−(∂2u∂x2+∂2u∂y2)=−32e12(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 0<x, y<1, \quad 0<t \leqslant 1∂t∂u−(∂x2∂2u+∂y2∂2u)=−23e21(x+y)−t,0<x,y<1,0<t⩽1
u(x,y,0)=e12(x+y),0<x,y<1u(x, y, 0)=e^{\frac{1}{2}(x+y)}, \quad 0<x, y<1u(x,y,0)=e21(x+y),0<x,y<1
u(0,y,t)=e12y−t,u(1,y,t)=e12(1+y)−t,0⩽y⩽1,0⩽t⩽1u(0, y, t)=e^{\frac{1}{2} y-t}, \quad u(1, y, t)=e^{\frac{1}{2}(1+y)-t}, \quad 0 \leqslant y \leqslant 1, \quad 0 \leqslant t \leqslant 1u(0,y,t)=e21y−t,u(1,y,t)=e21(1+y)−t,0⩽y⩽1,0⩽t⩽1
u(x,0,t)=e12x−t,u(x,1,t)=e12(1+x)−t,0<x<1,0⩽t⩽1u(x, 0, t)=e^{\frac{1}{2} x-t}, \quad u(x, 1, t)=e^{\frac{1}{2}(1+x)-t}, \quad 0<x<1, \quad 0 \leqslant t \leqslant 1u(x,0,t)=e21x−t,u(x,1,t)=e21(1+x)−t,0<x<1,0⩽t⩽1
该问题的精确解为：u(x,y,t)=e12(x+y)−tu(x, y, t)=e^{\frac{1}{2}(x+y)-t}u(x,y,t)=e21(x+y)−t

算法原理(偷懒贴图)

从第K层过渡到K+1/2层

从第K+1/2过渡到第K+1层

算法实现思想

结合上述算法思想，实现思路如下：
（1）分别给出空间区间数M，与时间间隔数N，得到空间步长h与时间步长c；
（2）构建矩阵A1、A2、A3、A4，令k=1；
（3）判断k是否等于N+1，若相等执行步骤（4）、（5），反之执行步骤（6）；
（4）构建U（k）、F（k+1/2）并求解出中间层U（k+1/2）;
（5）构建F（k+1）,求解U（K+1），将求解结果写入三维数组U3，k=k+1，并返回步骤（3）；
（6）输出U3，程序结束。

Matlab实现

ADI格式实现尝试

function U3 = ADI(M,N,h,c,x,y,t)
%***************************************************************
% 求解二维抛物型奇次方程 ut=uxx+uyy;
% 0≤x≤1 ,0≤y≤1 , 0≤t≤T;
% 该问题的定解: u(x,y,t)=exp( 1/2(x+y)-t);
% 初值条件: u(x,y,0)=exp( 1/2(x+y) );
% 边值条件: u(0,y,t)=exp( 1/2*(y)-t ); u(1,y,t)=exp( 1/2*(1+y)-t );
%  u(x,0,t)=exp( 1/2*x-t ); u(x,1,t)=exp( 1/2*(x+1)-t );
%***************************************************************
r=c./h^2;
%构造B1
bb1=zeros(1,M-1);                          %给主对角线位置赋值.
for i=1:M-1bb1(i)=1+r;
end
B1=diag(bb1);
%构造C1
cc1=zeros(1,M-1);                          %给主对角线位置赋值.
for i=1:M-1cc1(i)=-r./2;
end
C1=diag(cc1);
%***************************************************************
%构造A1
aa1=blkdiag(kron(eye(M-1),B1));             %创建主对角矩阵块
mid=mat2cell(aa1,(M-1)*ones(1,M-1),(M-1)*ones(1,M-1));
mid(find(diag(ones(1,M-2),1)==1))={C1};     %创建上对角矩阵块
mid(find(diag(ones(1,M-2),-1)==1))={C1};    %创建下对角矩阵块
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-1bb2(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-1bb5(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-1bb6(i)=1-r;
end
B4=diag(bb6);
%***************************************************************
%构造C4
cc=zeros(1,M-1);                          %给主对角线位置赋值.
for i=1:M-1cc(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-1ff2(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-1ff2(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-1U0(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-1uk(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-1uk(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/2pp=pp+1;
end
%***************************************************************
% 构造F(k)
F=zeros((M-1)*(M-1),1);%Fk
c=1;
for j=1:M-1F(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-1ff6(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:Mfor j=2:MU3(i,j,k)=U(i-1,j-1); end
end
for i=1:M+1U3(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)); %右边界
endfor i=1:M+1U3(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)); %下边界
endk=k+1;  %进入下一层
end

调用的主函数

主要实现：完成迭代矩阵U3与解析解矩阵U4的实现，主要实现图像拟合、与误差分析

求解结果（目标）

拟合图像

动态实现尝试

这里的动态实现与实验的不一样，实现思路主要来自b站的网课

实验总结

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

参考

物理老师的热传导课

偏微分方程数值解—ADI格式求解二维抛物型方程相关推荐

二维有限元方程matlab,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维 Poisson 方程的 MATLAB 实现陈莲a ,郭元辉b ,邹叶童a ( 西华师范大学 a．数学与信息学院; b．教育信息技术中心,四川南充 6437009) 摘要: ...
galerkin有限元法matlab实现,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维Poisson方程的MATLAB实现陈莲a,郭元辉b,邹叶童a [摘要]文章讨论了圆形区域上的三角形单元剖分.有限元空间,通过变分形式离散得到有限元方程. 用MATLAB编程求得数值 ...
使用matlab求解二维浅水方程的数值解（一）—浅水波
最近在读<ocean modelling for beginners>这本书,对于做海洋数值模拟工作的小白来说,这绝对是一本好书.强烈推荐给理论基础较弱的学习者,这本书循序渐进,由简入繁的 ...
有限元方法求解二维拉普拉斯方程C++实现
文章目录前言问题区域方程程序设计几何区域网格单元刚度矩阵组装数值结果问题区域网格 u 值图像(结果导出借助Matlab画图) 总结前言本文利用C++语言实现在二维任意区域(内部 ...
fortran用adi方法解二维热方程_完全耦合热弹性问题的普通态基近场动力学模拟——1 引言
Gao, Y. and Oterkus, S., 2019. Ordinary state-based peridynamic modelling for fully coupled thermoel ...
使用matlab求解二维浅水方程的数值解（二）—波浪的折射
如果大家去过海边,会有这样的感受:如果你面向大海,不管海岸是平直还是蜿蜒曲折的,你感觉到海浪总是迎着你传过来.这是波浪传播中的一种物理现象--波浪的折射.波浪由远海传入近岸的过程中,随着水深变浅,波浪 ...
[计算流体力学][Matlab] 使用 A,B,C 格式与蛙跳格式求解二维对流问题
1. 题目 2. 转述原题目要求可以简化为: 对于二维对流方程: ∂u/∂t+∂u/∂x+∂u/∂y=0 u(x,y,0)={█(1,when-4≤x≤4,-4≤y≤4@0,other)┤ 取计算范 ...
有限元方法基础-以二维拉普拉斯方程为例(附程序)
文章目录前言问题描述问题区域偏微分方程变分问题(弱形式) 有限元离散二维线性有限元 : 参考基函数 2D linear finite element : affine mapping 有限 ...
一节双曲型方程基于MATLAB的求解,二维双曲型方程的分组并行格式及其数值实验...
第 28卷第 2期 2010年 6月湖北民族学院学报 (自然科学版 ) Journal of Hubei University for Nationalities(Natural Science E ...
二维Poisson方程五点差分格式及简单求解方法Python实现
二维Poisson方程简介给出二维 PoissonPoissonPoisson 方程 DirichletDirichletDirichlet 边值问题: {−Δu=f(x,y)(x,y)∈Ωu=φ ...

偏微分方程数值解—ADI格式求解二维抛物型方程