lh2.m

二阶椭圆偏微分方程的差分算法

function lh2(p,q,f,alpha,xspan,yspan,M,N)
% 二阶椭圆型方程边值问题的差分格式
% -[d(p*du/dx)/dx+d(p*du/dy)/dy]+q*u=f;
% u(x,y)=alpha(x,y) every (x,y) belongs to d(xspan*yspan)
% p,q,f均为关于x,y的二元函数
% test: u(x,y)=x^2+y^2;f(x,y)=-4+x^2+y^2;
p=@(x,y) 1;
q=@(x,y) 0; 
f=@(x,y) 4-x.^2-y.^2;
xspan=[-1,1];yspan=[-1,1];
alpha=@(x,y) 0;
M=30;N=10;

h1=(xspan(2)-xspan(1))/M;
x=xspan(1):h1:xspan(2);
h2=(yspan(2)-yspan(1))/N;
y=yspan(1):h2:yspan(2);
u=zeros(M+1,N+1);
H=zeros((M-1)*(N-1),(M-1)*(N-1));
g=zeros(M-1,N-1);
u(1,:)=alpha(x(1),y(:));
u(M+1,:)=alpha(x(M+1),y(:));
u(:,1)=alpha(x(:),y(1));
u(:,N+1)=alpha(x(:),y(N+1));
uh=u([2:M],[2:N]);
for i=1:M-1,
    g(i,:)=f(x(i),y(2:N));
end;
[x,y]=meshgrid(x,y);
x=x';
y=y';
u=u';
uh=uh';
g=g';
for i=1:(M-1)*(N-1),
            alpha1=p(x(i)+h1/2,y(i))/h1^2;%alpha1
            alpha2=p(x(i),y(i)+h2/2)/h2^2;%alpha2
            alpha3=p(x(i)-h1/2,y(i))/h1^2;%alpha3
            alpha4=p(x(i),y(i)-h2/2)/h2^2;%alpha4
            alpha0=(alpha1+alpha2+alpha3+alpha4)+q(round(i/(M-1)),mod(i,M-1));%alpha0
            if i>1,
                H(i,i-1)=-alpha3;
            else
                g(i)=g(i)+u(i);
            end;
            if i>M-1,
                H(i,i-M+1)=-alpha4;
            else
                g(i)=g(i)+alpha3*u(i);
            end;
            H(i,i)=alpha0;
            if i<(M-1)*(N-2),
                H(i,i+M-1)=-alpha2;
            else
                g(i)=g(i)+u(i);
            end;
            if i<(M-1)*(N-1)-1,
                H(i,i+1)=-alpha1;
            else
                g(i)=g(i)+u(i);
            end;
end;
uh(:)=H^-1*g(:);
x=x';
y=y';
u=u';
uh=uh';
g=g';
%uh(:)=in_guass_seidel(H,g(:),eye(length(g(:)),1))%用高斯塞德尔迭代求解
%uh(:)=in_jacobi(H,g(:),eye(length(g(:)),1))%用雅克比迭代求解
u([2:M],[2:N])=uh;
figure;
mesh(x,y,u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%雅克比跌代
function x=in_jacobi(A,b,x1)
exps=1.0e-6;
count=0;
n=length(x1);
for row=1:n,
    A(row,[1:row-1,row+1:n])=-A(row,[1:row-1,row+1:n])/A(row,row);
    b(row,1)=b(row,1)/A(row,row);
    A(row,row)=0;
end
x=x1+2*exps;
while max(abs(x-x1))>exps & count<100,
    x=A*x1+b;
    temp=x1;
    x1=x;
    x=temp;
    count=count+1;
end;
in_jacobi=x;
disp('叠代次数是:');
count
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%高斯塞得尔迭代
function x=in_guass_seidel(A,b,x)
exps=1.0e-6;
dx=x;
count=0;
n=length(x);
while max(abs(dx))>exps & count<=100,
    for row=1:n,
        if A(row,row)==0,
            return;
        end
        dx(row)=(b(row)-A(row,[1:row-1])*x([1:row-1])-A(row,[row:n])*x([row:n]))...
            /A(row,row);
        x(row)=x(row)+dx(row);
    end
    count=count+1;
end
in_guass_seidel=x;
disp('叠代次数:');
count