[原创]偏微分方程数值解法的MATLAB源码【更新完毕】
说明:由于偏微分的程序都比较长,比其他的算法稍复杂一些,所以另开一贴,专门上传偏微分的程序 谢谢大家的支持!
其他的数值算法见:
..//Announce/Announce.asp?BoardID=209&id=8245004
1、古典显式格式求解抛物型偏微分方程(一维热传导方程)
function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典显式格式求解抛物型偏微分方程
%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %
%方程:u_t=C*u_xx 0 <= 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');
%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值 if nargin==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;
if r > 0.5
disp('r > 0.5,不稳定') end
%计算初值和边值 U=zeros(M+1,N+1); for i=1:M+1 U(i,1)=phi(x(i)); end for j=1:N+1 U(1,j)=psi1(t(j)); U(M+1,j)=psi2(t(j)); end
%逐层求解 for j=1:N for i=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 [U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典隐式格式求解抛物型偏微分方程
%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) %
%方程:u_t=C*u_xx 0 <= 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');
%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值 if nargin==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);%矩阵的上对角线元素 for i=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); for i=1:M+1 U(i,1)=phi(x(i)); end for j=1:N+1
U(1,j)=psi1(t(j)); U(M+1,j)=psi2(t(j)); end
%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackward) for j=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;
此算法需要使用追赶法求解三对角线性方程组,这个算法在上一篇帖子中已经给出,为了方便,再给出来
追赶法解三对角线性方程组
function x=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=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]'; %x=EqtsForwardAndBackward(L,D,U,b)
偏微分方程数值解法的MATLAB源码
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)