偏微分方程数值解法的MATLAB代码.docx
- 文档编号:1323156
- 上传时间:2022-10-20
- 格式:DOCX
- 页数:19
- 大小:323.90KB
偏微分方程数值解法的MATLAB代码.docx
《偏微分方程数值解法的MATLAB代码.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值解法的MATLAB代码.docx(19页珍藏版)》请在冰豆网上搜索。
偏微分方程数值解法的MATLAB代码
偏微分方程数值解法的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;
%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackwardforj=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]=PDEEllipseSquareLaplaceDi
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法 MATLAB 代码
![提示](https://static.bdocx.com/images/bang_tan.gif)