发布时间 : 星期一 文章偏微分方程数值解法的MATLAB源码更新完毕开始阅读
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;
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个时间层……