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

%检查参数的输入是否正确 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(i)-L(i-1)*U(i-1); end x=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); end return;

古典隐式格式

在以后的程序中，我们都取C=1，不再作为一个输入参数处理

3、Crank-Nicolson隐式格式求解抛物型偏微分方程 需要调用追赶法的程序

function [U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N) %Crank-Nicolson隐式格式求解抛物型偏微分方程 %[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 -边值条件，定义为内联函数 % 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);%矩阵的上对角线元素 for i=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); 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

B=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; end

B(M-1,M-1)=1-r;

%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward） 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); 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 [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<=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=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 %构造矩阵A M2=(M-1)^2; A=zeros(M2); for i=1:M2 A(i,i)=4; end

for i=1:M2-1 if mod(i,M-1)~=0 A(i,i+1)=-1; A(i+1,i)=-1; end end

for i=1:M2-M+1 A(i,i+M-1)=-1;