偏微分方程数值解法的MATLAB代码.docx
《偏微分方程数值解法的MATLAB代码.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解法的MATLAB代码.docx（22页珍藏版）》请在冰豆网上搜索。
偏微分方程数值解法的MATLAB代码
[原创]偏微分方程数值解法的MATLAB源码【更新完毕】
说明：
由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序
谢谢大家的支持！
其他的数值算法见：
..//Announce/Announce.asp?
BoardID=209&id=8245004
1、古典显式格式求解抛物型偏微分方程（一维热传导方程）
function[Uxt]=PDEParabolicClassicalExplicit（uX,uT,phi,psi1,psi2,M,N,C）
%古典显式格式求解抛物型偏微分方程
%[Uxt]=PDEParabolicClassicalExplicit（uX,uT,phi,psi1,psi2,M,N,C）
%
%方程：
u_t=C*u_xx0<=x<=uX,0<=t<=uT
%初值条件：
u（x,0）=phi（x）
%边值条件：
u（0,t）=psi1（t）,u（uX,t）=psi2（t）
%
%输出参数：
U-解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……
% x-空间变量
% t-时间变量
%输入参数：
uX-空间变量x的取值上限
% uT-时间变量t的取值上限
% phi-初值条件，定义为内联函数
% psi1-边值条件，定义为内联函数
% psi2-边值条件，定义为内联函数
% M-沿x轴的等分区间数
% N-沿t轴的等分区间数
% C-系数，默认情况下C=1
%
%应用举例：
%uX=1;uT=0.2;M=15;N=100;C=1;
%phi=inline（'sin（pi*x）'）;psi1=inline（'0'）;psi2=inline（'0'）;
%[Uxt]=PDEParabolicClassicalExplicit（uX,uT,phi,psi1,psi2,M,N,C）;
%设置参数C的默认值
ifnargin==7
C=1;
end
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=（0:
M）*dx;
t=（0:
N）*dt;
r=C*dt/dx/dx;%步长比
r1=1-2*r;
ifr>0.5
disp（'r>0.5,不稳定'）
end
%计算初值和边值
U=zeros（M+1,N+1）;
fori=1:
M+1
U（i,1）=phi（x（i））;
end
forj=1:
N+1
U（1,j）=psi1（t（j））;
U（M+1,j）=psi2（t（j））;
end
%逐层求解
forj=1:
N
fori=2:
M
U（i,j+1）=r*U（i-1,j）+r1*U（i,j）+r*U（i+1,j）;
end
end
U=U';
%作出图形
mesh（x,t,U）;
title（'古典显式格式，一维热传导方程的解的图像'）
xlabel（'空间变量x'）
ylabel（'时间变量t'）
zlabel（'一维热传导方程的解U'）
return;
古典显式格式不稳定情况
古典显式格式稳定情况
2、古典隐式格式求解抛物型偏微分方程（一维热传导方程）
function[Uxt]=PDEParabolicClassicalImplicit（uX,uT,phi,psi1,psi2,M,N,C）
%古典隐式格式求解抛物型偏微分方程
%[Uxt]=PDEParabolicClassicalImplicit（uX,uT,phi,psi1,psi2,M,N,C）
%
%方程：
u_t=C*u_xx0<=x<=uX,0<=t<=uT
%初值条件：
u（x,0）=phi（x）
%边值条件：
u（0,t）=psi1（t）,u（uX,t）=psi2（t）
%
%输出参数：
U-解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……
% x-空间变量
% t-时间变量
%输入参数：
uX-空间变量x的取值上限
% uT-时间变量t的取值上限
% phi-初值条件，定义为内联函数
% psi1-边值条件，定义为内联函数
% psi2-边值条件，定义为内联函数
% M-沿x轴的等分区间数
% N-沿t轴的等分区间数
% C-系数，默认情况下C=1
%
%应用举例：
%uX=1;uT=0.2;M=50;N=50;C=1;
%phi=inline（'sin（pi*x）'）;psi1=inline（'0'）;psi2=inline（'0'）;
%[Uxt]=PDEParabolicClassicalImplicit（uX,uT,phi,psi1,psi2,M,N,C）;
%设置参数C的默认值
ifnargin==7
C=1;
end
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=（0:
M）*dx;
t=（0:
N）*dt;
r=C*dt/dx/dx;%步长比
Diag=zeros（1,M-1）;%矩阵的对角线元素
Low=zeros（1,M-2）;%矩阵的下对角线元素
Up=zeros（1,M-2）;%矩阵的上对角线元素
fori=1:
M-2
Diag（i）=1+2*r;
Low（i）=-r;
Up（i）=-r;
end
Diag（M-1）=1+2*r;
%计算初值和边值
U=zeros（M+1,N+1）;
fori=1:
M+1
U（i,1）=phi（x（i））;
end
forj=1:
N+1
U（1,j）=psi1（t（j））;
U（M+1,j）=psi2（t（j））;
end
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
forj=1:
N
b1=zeros（M-1,1）;
b1
（1）=r*U（1,j+1）;
b1（M-1）=r*U（M+1,j+1）;
b=U（2:
M,j）+b1;
U（2:
M,j+1）=EqtsForwardAndBackward（Low,Diag,Up,b）;
end
U=U';
%作出图形
mesh（x,t,U）;
title（'古典隐式格式，一维热传导方程的解的图像'）
xlabel（'空间变量x'）
ylabel（'时间变量t'）
zlabel（'一维热传导方程的解U'）
return;
此算法需要使用追赶法求解三对角线性方程组，这个算法在上一篇帖子中已经给出，为了方便，再给出来
追赶法解三对角线性方程组
functionx=EqtsForwardAndBackward（L,D,U,b）
%追赶法求解三对角线性方程组Ax=b
%x=EqtsForwardAndBackward（L,D,U,b）
%x:
三对角线性方程组的解
%L:
三对角矩阵的下对角线，行向量
%D:
三对角矩阵的对角线，行向量
%U:
三对角矩阵的上对角线，行向量
%b:
线性方程组Ax=b中的b，列向量
%
%应用举例:
%L=[-1-2-3];D=[2345];U=[-1-2-3];b=[61-21]';
%x=EqtsForwardAndBackward（L,D,U,b）
%检查参数的输入是否正确
n=length（D）;m=length（b）;
n1=length（L）;n2=length（U）;
ifn-n1~=1||n-n2~=1||n~=m
disp（'输入参数有误！
'）
x='';
return;
end
%追的过程
fori=2:
n
L（i-1）=L（i-1）/D（i-1）;
D（i）=D（i）-L（i-1）*U（i-1）;
end
x=zeros（n,1）;
x
（1）=b
（1）;
fori=2:
n
x（i）=b（i）-L（i-1）*x（i-1）;
end
%赶的过程
x（n）=x（n）/D（n）;
fori=n-1:
-1:
1
x（i）=（x（i）-U（i）*x（i+1））/D（i）;
end
return;
古典隐式格式
在以后的程序中，我们都取C=1，不再作为一个输入参数处理
3、Crank-Nicolson隐式格式求解抛物型偏微分方程
需要调用追赶法的程序
function[Uxt]=PDEParabolicCN（uX,uT,phi,psi1,psi2,M,N）
%Crank-Nicolson隐式格式求解抛物型偏微分方程
%[Uxt]=PDEParabolicCN（uX,uT,phi,psi1,psi2,M,N）
%
%方程：
u_t=u_xx0<=x<=uX,0<=t<=uT
%初值条件：
u（x,0）=phi（x）
%边值条件：
u（0,t）=psi1（t）,u（uX,t）=psi2（t）
%
%输出参数：
U-解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……
% x-空间变量
% t-时间变量
%输入参数：
uX-空间变量x的取值上限
% uT-时间变量t的取值上限
% phi-初值条件，定义为内联函数
% psi1-边值条件，定义为内联函数
% psi2-边值条件，定义为内联函数
% M-沿x轴的等分区间数
% N-沿t轴的等分区间数
%
%应用举例：
%uX=1;uT=0.2;M=50;N=50;
%phi=inline（'sin（pi*x）'）;psi1=inline（'0'）;psi2=inline（'0'）;
%[Uxt]=PDEParabolicCN（uX,uT,phi,psi1,psi2,M,N）;
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=（0:
M）*dx;
t=（0:
N）*dt;
r=dt/dx/dx;%步长比
Diag=zeros（1,M-1）;%矩阵的对角线元素
Low=zeros（1,M-2）;%矩阵的下对角线元素
Up=zeros（1,M-2）;%矩阵的上对角线元素
fori=1:
M-2
Diag（i）=1+r;
Low（i）=-r/2;
Up（i）=-r/2;
end
Diag（M-1）=1+r;
%计算初值和边值
U=zeros（M+1,N+1）;
fori=1:
M+1
U（i,1）=phi（x（i））;
end
forj=1:
N+1
U（1,j）=psi1（t（j））;
U（M+1,j）=psi2（t（j））;
end
B=zeros（M-1,M-1）;
fori=1:
M-2
B（i,i）=1-r;
B（i,i+1）=r/2;
B（i+1,i）=r/2;
end
B（M-1,M-1）=1-r;
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
forj=1:
N
b1=zeros（M-1,1）;
b1
（1）=r*（U（1,j+1）+U（1,j））/2;
b1（M-1）=r*（U（M+1,j+1）+U（M+1,j））/2;
b=B*U（2:
M,j）+b1;
U（2:
M,j+1）=EqtsForwardAndBackward（Low,Diag,Up,b）;
end
U=U';
%作出图形
mesh（x,t,U）;
title（'Crank-Nicolson隐式格式，一维热传导方程的解的图像'）
xlabel（'空间变量x'）
ylabel（'时间变量t'）
zlabel（'一维热传导方程的解U'）
return;
Crank-Nicolson隐式格式
4、正方形区域Laplace方程Diriclet问题的求解
需要调用Jacobi迭代法和Guass-Seidel迭代法求解线性方程组
function[Uxy]=PDEEllipseSquareLaplaceDirichlet（ub,phi1,phi2,psi1,psi2,M,type）
%正方形区域Laplace方程的Diriclet边值问题的差分求解
%此程序需要调用Jacobi迭代法或者Guass-Seidel迭代法求解线性方程组
%[Uxy]=PDEEllipseSquareLaplaceDirichlet（ub,phi1,phi2,psi1,psi2,M,type）
%
%方程：
u_xx+u_yy=0 0<=x,y<=ub
%边值条件：
u（0,y）=phi1（y）
% u（ub,y）=phi2（y）
% u（x,0）=psi1（x）
% u（x,ub）=psi2（x）
%
%输出参数：
U-解矩阵，第一行表示y=0时的值，第二行表示第y=h时的值……
% x-横坐标
% y-纵坐标
%输入参数：
ub-变量边界值的上限
% phi1,phi2,psi1,psi2-边界函数，定义为内联函数
% M-横纵坐标的等分区间数
% type-求解差分方程的迭代格式，若type='Jacobi'，采用Jacobi迭代格式
% 若type='GS'，采用Guass-Seidel迭代格式。
默认情况下，type='GS'
%
%应用举例：
%ub=4;M=20;
%phi1=inline（'y*（4-y）'）;phi2=inline（'0'）;psi1=inline（'sin（pi*x/4）'）;psi2=inline（'0'）;
%[Uxy]=PDEEllipseSquareLaplaceDirichlet（ub,phi1,phi2,psi1,psi2,M,'GS'）;
ifnargin==6
type='GS';
end
%步长
h=ub/M;
%横纵坐标
x=（0:
M）*h;
y=（0:
M）*h;
%差分格式的矩阵形式AU=K
%构造矩阵A
M2=（M-1）^2;
A=zeros（M2）;
fori=1:
M2
A（i,i）=4;
end
fori=1:
M2-1
ifmod（i,M-1）~=0
A（i,i+1）=-1;
A（i+1,i）=-1;
end
end
fori=1:
M2-M+1
A（i,i+M-1）=-1;
A（i+M-1,i）=-1;
end
U=zeros（M+1）;
%边值条件
fori=1:
M+1
U（i,1）=psi1（（i-1）*h）;
U（i,M+1）=psi2（（i-1）*h）;
U（1,i）=phi1（（i-1）*h）;
U（M+1,i）=phi2（（i-1）*h）;
end
%构造K
K=zeros（M2,1）;
fori=1:
M-1
K（i）=U（i+1,1）;
K（M2-i+1）=U（i+1,M+1）;
end
K
（1）=K
（1）+U（1,2）;
K（M-1）=K（M-1）+U（M+1,2）;
K（M2-M+2）=K（M2-M+2）+U（1,M）;
K（M2）=K（M2）+U（M+1,M）;
fori=2:
M-2
K（（M-1）*（i-1）+1）=U（1,i+1）;
K（（M-1）*i）=U（M+1,i+1）;
end
x0=ones（M2,1）;
switchtype
%调用Guass-Seidel迭代法求解线性方程组AU=K
case'Jacobi'
X=EqtsJacobi（A,K,x0）;
%调用Guass-Seidel迭代法求解线性方程组AU=K
case'GS'
X=EqtsGS（A,K,x0）;
otherwise
disp（'差分格式类型输入错误'）
return;
end
%把求解结果化成矩阵型式
fori=2:
M
forj=2:
M
U（j,i）=X（j-1+（M-1）*（i-2））;
end
end
U=U';
%作出图形
mesh（x,y,U）;
title（'五点差分格式Laplace方程Diriclet问题的解的图像'）
xlabel（'x'）
ylabel（'y'）
zlabel（'Laplace方程Diriclet问题的解U'）
return;
正方形区域Laplace方程五点差分格式
5、一阶双曲型方程的差分方法
function[Uxt]=PDEHyperbolic（uX,uT,M,N,C,phi,psi1,psi2,type）
%一阶双曲型方程的差分格式
%[Uxt]=PDEHyperbolic（uX,uT,M,N,C,phi,psi1,psi2,type）
%
%方程:
u_t+C*u_x=0 0<=t<=uT,0<=x<=uX
%初值条件:
u（x,0）=phi（x）
%
%输出参数:
U-解矩阵，第一行表示初值，第二行表示第2个时间层……
% x-横坐标
% t-纵坐标，时间
%输入参数:
uX-变量x的上界
% uT-变量t的上界
% M-变量x的等分区间数
% N-变量t的等分区间数
% C-系数
% phi-初值条件函数，定义为内联函数
% psi1,psi2-边值条件函数，定义为内联函数
% type-差分格式，从下列值中选取
% -type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解
% -type='CourantIsaacsonRees'，采用Courant-Isaacson-Rees差分格式求解
% -type='LeapFrog'，采用Leap-Frog（蛙跳）差分格式求解
% -type='LaxWendroff'，采用Lax-Wendroff差分格式求解
% -type='CrankNicolson'，采用Crank-Nicolson差分格式求解，此格式需调用追赶法
% 求解三对角线性方程组
%
h=uX/M;%变量x的步长
k=uT/N;%变量t的步长
r=k/h;%步长比
x=（0:
M）*h;
t=（0:
N）*k;
U=zeros（M+1,N+1）;
%初值条件
fori=1:
M+1
U（i,1）=phi（x（i））;
end
%边值条件
forj=1:
N+1
U（1,j）=psi1（t（j））;
U（M+1,j）=psi2（t（j））;
%U（1,j）=NaN;
%U（M+1,j）=NaN;
end
switchtype
%Lax-Friedrichs差分格式
case'LaxFriedrichs'
ifabs（C*r）>1
disp（'|C*r|>1，Lax-Friedrichs差分格式不稳定！
'）
end
%逐层求解
forj=1:
N
fori=2:
M
U（i,j+1）=（U（i+1,j）+U（i-1,j））/2-C*r*（U（i+1,j）-U（i-1,j））/2;
end
end
%Courant-Isaacson-Rees差分格式
case'CourantIsaacsonRees'
ifC<0
disp（'C<0，采用前差公式'）
ifC*r<-1
disp（'Courant-Isaacson-Lees差分格式不稳定！
'）
end
%逐层求解
forj=1:
N
fori=2:
M
U（i,j+1）=（1+C*r）*U（i,j）-C*r*U（i+1,j）;
end
end
else
disp（'C>0，采用后差公式'）
ifC*r>1
disp（'Courant-Isaacson-Lees差分格式不稳定！
'）
end
%逐层求解
forj=1:
N
fori=2:
M
U（i,j+1）=C*r*U（i-1,j）+（1-C*r）*U（i,j）;
end
end
end
%Leap-Frog（蛙跳）差分格式
case'LeapFrog'
phi2=input（'请输入第二层初值条件函数：
psi2='）;
ifabs（C*r）>1
disp（'|C*r|>1，Leap-Frog差分格式不稳定！
'）
end
%第二层初值条件
fori=1:
M+1
U（i,2）=phi2（x（i））;
end
%逐层求解
升级会员 