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

seSquareLaplaceDirichlet(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;