1、原创偏微分方程数值解法的MATLAB源码【更新完毕】说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！其他的数值算法见：./Announce/Announce.asp?BoardID=209&id=82450041、古典显式格式求解抛物型偏微分方程（一维热传导方程）function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典显式格式求解抛物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1
2、,psi2,M,N,C)%方程：u_t=C*u_xx 0 = x = uX,0 = t 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
3、(空间变量 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(
4、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=in
5、line(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)=
6、-r; Up(i)=-r;endDiag(M-1)=1+2*r;%计算初值和边值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%逐层求解，需要使用追赶法（调用函数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
7、,Diag,Up,b);endU=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
8、:三对角矩阵的对角线，行向量%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)%检查参数的输入是否正确n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 = 1 | n-n2 = 1 | n = m disp(输入参数有误！) x= ; return;end%追的过程for i=2:n L(i-1)=L(i-1)/D(i-1); D(i)=D(
9、i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:n x(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i=n-1:-1:1 x(i)=(x(i)-U(i)*x(i+1)/D(i);endreturn;古典隐式格式在以后的程序中，我们都取C=1，不再作为一个输入参数处理3、Crank-Nicolson隐式格式求解抛物型偏微分方程需要调用追赶法的程序function U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%Crank-Nicolson隐式格式
10、求解抛物型偏微分方程%U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%方程：u_t=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 -边值条件，定义为内联函
11、数% M -沿x轴的等分区间数% N -沿t轴的等分区间数%应用举例：%uX=1;uT=0.2;M=50;N=50;%phi=inline(sin(pi*x);psi1=inline(0);psi2=inline(0);%U x t=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);%矩阵的上
12、对角线元素for i=1:M-2 Diag(i)=1+r; Low(i)=-r/2; Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值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);endB=zeros(M-1,M-1);for i=1:M-2 B(i,i)=1-r; B(i,i+1)=r/2; B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBac
13、kward）for j=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);endU=U;%作出图形mesh(x,t,U);title(Crank-Nicolson隐式格式，一维热传导方程的解的图像)xlabel(空间变量 x)ylabel(时间变量 t)zlabel(一维热传导方程的解 U)return;Crank-Nicolson隐式格式4、正方形
14、区域Laplace方程Diriclet问题的求解需要调用Jacobi迭代法和Guass-Seidel迭代法求解线性方程组function U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%正方形区域Laplace方程的Diriclet边值问题的差分求解%此程序需要调用Jacobi迭代法或者Guass-Seidel迭代法求解线性方程组%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%方程：u_xx+u_yy=0 0=
15、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=2
16、0;%phi1=inline(y*(4-y);phi2=inline(0);psi1=inline(sin(pi*x/4);psi2=inline(0);%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,GS);if nargin=6 type=GS;end%步长h=ub/M;%横纵坐标x=(0:M)*h;y=(0:M)*h;%差分格式的矩阵形式AU=K%构造矩阵AM2=(M-1)2;A=zeros(M2);for i=1:M2 A(i,i)=4;endfor i=1:M2-1 if mod(i,M-1)=0 A(
17、i,i+1)=-1; A(i+1,i)=-1; endendfor i=1:M2-M+1 A(i,i+M-1)=-1; A(i+M-1,i)=-1;endU=zeros(M+1);%边值条件for i=1: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);end%构造KK=zeros(M2,1);for i=1:M-1 K(i)=U(i+1,1); K(M2-i+1)=U(i+1,M+1);endK(1)=K(1)+U(1,2);K(M-1)=K(M-1)+U(M
18、+1,2);K(M2-M+2)=K(M2-M+2)+U(1,M);K(M2)=K(M2)+U(M+1,M);for i=2:M-2 K(M-1)*(i-1)+1)=U(1,i+1); K(M-1)*i)=U(M+1,i+1);endx0=ones(M2,1);switch type %调用Guass-Seidel迭代法求解线性方程组AU=K case Jacobi X=EqtsJacobi(A,K,x0); %调用Guass-Seidel迭代法求解线性方程组AU=K case GS X=EqtsGS(A,K,x0); otherwise disp(差分格式类型输入错误) return;end%
19、把求解结果化成矩阵型式for i=2:M for j=2:M U(j,i)=X(j-1+(M-1)*(i-2); endendU=U;%作出图形mesh(x,y,U);title(五点差分格式Laplace方程Diriclet问题的解的图像)xlabel(x)ylabel(y)zlabel(Laplace方程Diriclet问题的解 U)return;正方形区域Laplace方程五点差分格式5、一阶双曲型方程的差分方法function U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%一阶双曲型方程的差分格式%U x t=PDEHyperb
20、olic(uX,uT,M,N,C,phi,psi1,psi2,type)%方程:u_t+C*u_x=0 0 = t = uT, 0 = x 1 disp(|C*r|1，Lax-Friedrichs差分格式不稳定！) end %逐层求解 for j=1:N for i=2:M U(i,j+1)=(U(i+1,j)+U(i-1,j)/2-C*r*(U(i+1,j)-U(i-1,j)/2; end end %Courant-Isaacson-Rees差分格式 case CourantIsaacsonRees if C0 disp(C0，采用前差公式) if C*r0，采用后差公式) if C*r1
21、disp(Courant-Isaacson-Lees差分格式不稳定！) end %逐层求解 for j=1:N for i=2:M U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j); end end end %Leap-Frog（蛙跳）差分格式 case LeapFrog phi2=input(请输入第二层初值条件函数：psi2=); if abs(C*r)1 disp(|C*r|1，Leap-Frog差分格式不稳定！) end %第二层初值条件 for i=1:M+1 U(i,2)=phi2(x(i); end %逐层求解 for j=2:N for i=2:M U(
22、i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j); end end %Lax-Wendroff差分格式 case LaxWendroff if abs(C*r)1 disp(|C*r|1，Lax-Wendroff差分格式不稳定！) end %逐层求解 for j=1:N for i=2:M U(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j)/2+C2*r2*(U(i+1,j)-2*U(i,j)+U(i-1,j)/2; end end %Crank-Nicolson隐式差分格式，需调用追赶法求解三对角线性方程组的算法 case CrankNi
23、colson Diag=zeros(1,M-1);%矩阵的对角线元素 Low=zeros(1,M-2);%矩阵的下对角线元素 Up=zeros(1,M-2);%矩阵的上对角线元素 for i=1:M-2 Diag(i)=4; Low(i)=-r*C; Up(i)=r*C; end Diag(M-1)=4; B=zeros(M-1,M-1); for i=1:M-2 B(i,i)=4; B(i,i+1)=-r*C; B(i+1,i)=r*C; end B(M-1,M-1)=4; %逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward） for j=1:N b1=zero
24、s(M-1,1); b1(1)=r*C*(U(1,j+1)+U(1,j)/2; b1(M-1)=-r*C*(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 otherwise disp(差分格式类型输入有误！) return;endU=U;%作出图形mesh(x,t,U);title(type 格式求解一阶双曲型方程的解的图像);xlabel(空间变量 x);ylabel(时间变量 t);zlabel(一阶双曲型方程的解 U);return;附录资料：MA
25、TLAB Cell函数使用技巧谈谈MATLAB中cell函数如果p为一个数，那么h(1)=p,是没有问题的。如果p为一个向量，那么h(1,:)=p是没有问题的。如果p是一个矩阵的话，上面的两种赋值方法都是会有错误的。那么要如何处理呢？这时就用到了cell数据类型了。cell的每个单元都可以存储任何数据，比如传递函数等。当然，存储矩阵更是没有问题的了。但是用cell数据类型之前，要先初始化。a=cell(n,m)那么就把a初始化为一个n行m列的空cell类型数据。如何赋值呢？a1,1=rand(5)那么a的1行1列的单元中存储的就是一个随机的55的方阵了。那么要用第一个单元中的方阵中的某个值呢？
26、可以如下引用：a1,1(2,3)就可以了，引用cell单元时要用,再引用矩阵的某个数据就要用()了。cell单元中的每个单元都是独立的，可以分别存储不同大小的矩阵或不同类型的数据。下面举个例子：a=cell(2,2);%预分配a1,1=cellclass;a1,2=1 2 2;a2,1=a,b,c;a2,2=9 5 6;a1,1ans =cellclassa1,2ans = 1 2 2a2,:ans =abcans = 9 5 6 b=a1,1b =cellclass元胞数组：元胞数组是MATLAB的一种特殊数据类型，可以将元胞数组看做一种无所不包的通用矩阵，或者叫做广义矩阵。组成元胞数组的元
27、素可以是任何一种数据类型的常数或者常量，每一个元素也可以具有不同的尺寸和内存占用空间，每一个元素的内容也可以完全不同，所以元胞数组的元素叫做元胞（cell）。和一般的数值矩阵一样，元胞数组的内存空间也是动态分配的。（1）元胞数组的创建 a=matlab,20;ones(2,3),1:10a = matlab 20 2x3 double 1x10 double b=matlab,20;ones(2,3),1:10b = matlab 20 2x3 double 1x10 double c=10c = 10c(1,2)=2c = 10 2c(2,2)=5c = 10 2 5isequal(a,b)
28、ans = 1whosName Size Bytes Class Attributesa 2x2 388 cell ans 1x1 1 logical b 2x2 388 cell c 2x2 208 cell 用cell函数创建元胞数组，创建的数组为空元胞。cell函数创建空元胞数组的主要目的是为数组预先分配连续的存储空间，节约内存占用，提高执行效率。 a=cell(1)a = b=cell(1,2)b = c=cell(3,3)c = d=cell(2,2,2)d(:,:,1) = d(:,:,2) = whosName Size Bytes Class Attributesa 1x1 4
29、 cell ans 1x1 1 logical b 1x2 8 cell c 3x3 36 cell d 2x2x2 32 cell （2）元胞数组的数据获得从元胞数组中读取数据，可保存为一个标准的数组或一个新的单元数组，或取出数组进行计算。元胞数组中数据的访问，可通过元胞内容的下标进行，用元胞数组名加大括号。大括号中数值表示元胞的下标。如a1，2表示元胞数组中第一行第二列的元胞。 a=20,matlab;ones(2,3),1:3a = 20 matlab 2x3 double 1x3 doublestr=a(1,2)str = matlabclass(str)ans =cellstr=a1
30、,2str =matlabclass(str)ans =char（）和有着本质的区别，大括号用于表示元胞的内容，小括号表示指定的元胞。a = 20 matlab 2x3 double 1x3 doublea2,1(2,2)ans = 0.9134a2,1(2,3)ans = 0.0975a1,2(2)ans =a使用元胞的下标，可将一个元胞数组的子集赋值给另一个变量，创建新的元胞数组。 a=1,2,3;4,5,6;7,8,9a = 1 2 3 4 5 6 7 8 9 b=a(2:3,2:3)b = 5 6 8 9 c=a(1:3,2:3)c = 2 3 5 6 8 9本例使用元胞下标的方式创建了新的元胞数组b和c，通过结果看出b和c就是元胞数组a的一部分。（3）元胞数组的删除和重塑要删除单元数组中的行或列，可以用冒号表示单元数组中的行或列，然后对其赋一个空矩阵即可。a=20,matlab;ones(2,3),1:3a = 20 matlab 2x3 double 1x3 d
