偏微分方程数值解

[原创]偏微分方程数值解法的MATLAB源码【更新完毕】说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！其他的数值算法见：..//Announce/Announce.asp?BoardID=209&id=82450041、古典显式格式求解抛物型偏微分方程（一维热传导方程）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==7C=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;ifr0.5disp('r0.5,不稳定')end%计算初值和边值U=zeros(M+1,N+1);fori=1:M+1U(i,1)=phi(x(i));endforj=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解forj=1:Nfori=2:MU(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j);endendU=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==7C=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-2Diag(i)=1+2*r;Low(i)=-r;Up(i)=-r;endDiag(M-1)=1+2*r;%计算初值和边值U=zeros(M+1,N+1);fori=1:M+1U(i,1)=phi(x(i));endforj=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）forj=1:Nb1=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);endU=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~=mdisp('输入参数有误！')x='';return;end%追的过程fori=2:nL(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);fori=2:nx(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);fori=n-1:-1:1x(i)=(x(i)-U(i)*x(i+1))/D(i);endreturn;古典隐式格式在以后的程序中，我们都取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-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值U=zeros(M+1,N+1);fori=1:M+1U(i,1)=phi(x(i));endforj=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));endB=zeros(M-1,M-1);fori=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）forj=1:Nb1=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);endU=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=00=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=inli