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

A(i+M-1,i)=-1; end

U=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 %构造K K=zeros(M2,1); for i=1:M-1 K(i)=U(i+1,1);

K(M2-i+1)=U(i+1,M+1); end

K(1)=K(1)+U(1,2); K(M-1)=K(M-1)+U(M+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); end

x0=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

%把求解结果化成矩阵型式 for i=2:M for j=2:M

U(j,i)=X(j-1+(M-1)*(i-2)); end end U=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]=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type) %

%方程:u_t+C*u_x=0 0 <= t <= uT, 0 <= x <= uX %初值条件:u(x,0)=phi(x) %

%输出参数:U -解矩阵，第一行表示初值，第二行表示第2个时间层……

% x -横坐标 % t -纵坐标，时间 %输入参数:uX -变量x的上界 % uT -变量t的上界 % M -变量x的等分区间数 % N -变量t的等分区间数 % C -系数

% phi -初值条件函数，定义为内联函数 % psi1,psi2 -边值条件函数，定义为内联函数 % type -差分格式，从下列值中选取

% -type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解

% -type='CourantIsaacsonRees'，采用Courant-Isaacson-Rees差分格式求解 % -type='LeapFrog'，采用Leap-Frog（蛙跳）差分格式求解 % -type='LaxWendroff'，采用Lax-Wendroff差分格式求解

% -type='CrankNicolson'，采用Crank-Nicolson差分格式求解，此格式需调用追赶法 % 求解三对角线性方程组 %

h=uX/M;%变量x的步长 k=uT/N;%变量t的步长 r=k/h;%步长比

x=(0:M)*h; t=(0:N)*k;

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)); %U(1,j)=NaN; %U(M+1,j)=NaN; end

switch type

%Lax-Friedrichs差分格式 case 'LaxFriedrichs' if abs(C*r)>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 C<0

disp('C<0，采用前差公式') if C*r<-1

disp('Courant-Isaacson-Lees差分格式不稳定！') end %逐层求解 for j=1:N for i=2:M

U(i,j+1)=(1+C*r)*U(i,j)-C*r*U(i+1,j); end end else

disp('C>0，采用后差公式') if C*r>1

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(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+C^2*r^2*(U(i+1,j)-2*U(i,j)+U(i-1,j))/2; end end

%Crank-Nicolson隐式差分格式，需调用追赶法求解三对角线性方程组的算法 case 'CrankNicolson'

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=zeros(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;