偏微分方程数值解法的MATLAB源码..docx

[原创]偏微分方程数值解法的MATLAB源码【更新完毕】说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！其他的数值算法见：..//Announce/Announce.asp?BoardID=209id古典显式格式求解抛物型偏微分方程（一维热传导方程）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));endfor 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);? ? endendU=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,psi