偏微分方程数值解法的MATLAB源码Word文档下载推荐.docx
《偏微分方程数值解法的MATLAB源码Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解法的MATLAB源码Word文档下载推荐.docx（16页珍藏版）》请在冰点文库上搜索。
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（'
%[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,不稳定'
）
%计算初值和边值
U=zeros（M+1,N+1）;
fori=1:
M+1
U（i,1）=phi（x（i））;
forj=1:
N+1
U（1,j）=psi1（t（j））;
U（M+1,j）=psi2（t（j））;
%逐层求解
N
fori=2:
M
U（i,j+1）=r*U（i-1,j）+r1*U（i,j）+r*U（i+1,j）;
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）
M=50;
N=50;
%[Uxt]=PDEParabolicClassicalImplicit（uX,uT,phi,psi1,psi2,M,N,C）;
Diag=zeros（1,M-1）;
%矩阵的对角线元素
Low=zeros（1,M-2）;
%矩阵的下对角线元素
Up=zeros（1,M-2）;
%矩阵的上对角线元素
M-2
Diag（i）=1+2*r;
Low（i）=-r;
Up（i）=-r;
Diag（M-1）=1+2*r;
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
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）;
古典隐式格式，一维热传导方程的解的图像'
此算法需要使用追赶法求解三对角线性方程组，这个算法在上一篇帖子中已经给出，为了方便，再给出来
追赶法解三对角线性方程组
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]'
%检查参数的输入是否正确
n=length（D）;
m=length（b）;
n1=length（L）;
n2=length（U）;
ifn-n1~=1||n-n2~=1||n~=m
输入参数有误！
'
x='
'
return;
%追的过程
fori=2:
n
L（i-1）=L（i-1）/D（i-1）;
D（i）=D（i）-L（i-1）*U（i-1）;
x=zeros（n,1）;
x
（1）=b
（1）;
x（i）=b（i）-L（i-1）*x（i-1）;
%赶的过程
x（n）=x（n）/D（n）;
fori=n-1:
-1:
1
x（i）=（x（i）-U（i）*x（i+1））/D（i）;
古典隐式格式
在以后的程序中，我们都取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<
%[Uxt]=PDEParabolicCN（uX,uT,phi,psi1,psi2,M,N）;
r=dt/dx/dx;
Diag（i）=1+r;
Low（i）=-r/2;
Up（i）=-r/2;
Diag（M-1）=1+r;
B=zeros（M-1,M-1）;
B（i,i）=1-r;
B（i,i+1）=r/2;
B（i+1,i）=r/2;
B（M-1,M-1）=1-r;
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:
Crank-Nicolson隐式格式，一维热传导方程的解的图像'
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='
%ub=4;
M=20;
%phi1=inline（'
y*（4-y）'
phi2=inline（'
sin（pi*x/4）'
%[Uxy]=PDEEllipseSquareLaplaceDirichlet（ub,phi1,phi2,psi1,psi2,M,'
ifnargin==6
type='
%步长
h=ub/M;
%横纵坐标
M）*h;
y=（0:
%差分格式的矩阵形式AU=K
%构造矩阵A
M2=（M-1）^2;
A=zeros（M2）;
M2
A（i,i）=4;
M2-1
ifmod（i,M-1）~=0
A（i,i+1）=-1;
A（i+1,i）=-1;
M2-M+1
A（i,i+M-1）=-1;
A（i+M-1,i）=-1;
U=zeros（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）;
%构造K
K=zeros（M2,1）;
M-1
K（i）=U（i+1,1）;
K（M2-i+1）=U（i+1,M+1）;
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）;
K（（M-1）*（i-1）+1）=U（1,i+1）;
K（（M-1）*i）=U（M+1,i+1）;
x0=ones（M2,1）;
switchtype
%调用Guass-Seidel迭代法求解线性方程组AU=K
case'
X=EqtsJacobi（A,K,x0）;
X=EqtsGS（A,K,x0）;
otherwise
disp（'
差分格式类型输入错误'
%把求解结果化成矩阵型式
forj=2:
U（j,i）=X（j-1+（M-1）*（i-2））;
mesh（x,y,U）;
五点差分格式Laplace方程Diriclet问题的解的图像'
x'
y'
Laplace方程Diriclet问题的解U'
正方形区域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<
=uT,0<
=uX
%初值条件:
%输出参数:
U-解矩阵，第一行表示初值，第二行表示第2个时间层……
t-纵坐标，时间
%输入参数:
uX-变量x的上界
uT-变量t的上界
M-变量x的等分区间数
N-变量t的等分区间数
C-系数
phi-初值条件函数，定义为内联函数
psi1,psi2-边值条件函数，定义为内联函数
type-差分格式，从下列值中选取
-type='
LaxFriedrichs'
，采用Lax-Friedrichs差分格式求解
CourantIsaacsonRees'
，采用Courant-Isaacson-Rees差分格式求解
LeapFrog'
，采用Leap-Frog（蛙跳）差分格式求解
LaxWendroff'
，采用Lax-Wendroff差分格式求解
CrankNicolson'
，采用Crank-Nicolson差分格式求解，此格式需调用追赶法
求解三对角线性方程组
h=uX/M;
%变量x的步长
k=uT/N;
%变量t的步长
r=k/h;
N）*k;
%初值条件
%U（1,j）=NaN;
%U（M+1,j）=NaN;
%Lax-Friedrichs差分格式
ifabs（C*r）>
|C*r|>
1，Lax-Friedrichs差分格式不稳定！
U（i,j+1）=（U（i+1,j）+U（i-1,j））/2-C*r*（U（i+1,j）-U（i-1,j））/2;
%Courant-Isaacson-Rees差分格式
ifC<
C<
0，采用前差公式'
ifC*r<
-1
Courant-Isaacson-Lees差分格式不稳定！
%逐层求解
forj=1:
U（i,j+1）=（1+C*r）*U（i,j）-C*r*U（i+1,j）;
else
C>
0，采用后差公式'
ifC*r>
U（i,j+1）=C*r*U（i-1,j）+（1-C*r）*U（i,j）;
%Leap-Frog（蛙跳）差分格式
phi2=input（'
请输入第二层初值条件函数：
psi2='
1，Leap-Frog差分格式不稳定！
升级会员 