有限差分法解第二类边界条件的椭圆形偏微分方程

有限差分法解第二类边界条件的椭圆形偏微分方程

clear
clc
ax=[0,1];
ay=[0,1];
x1=0.1;
xh=ax(1):x1:ax(2);%保持步长相同
yh=ay(1):x1:ay(2);
n1=length(xh);
%下面输入第二类边界条件，即neumann边界条件,为了保证唯一解，下面赋予函数在原点的初值为0；
syms x y
un=x+y;
A1=zeros(n1,n1);
A2=zeros(n1,n1);
I=zeros(n1,n1);
B=zeros(n1,n1);
for i=1:n1
    A1(i,i)=1;
    I(i,i)=1;
    B(i,i)=4;
end
for i=1:n1-1
    A2(i,i+1)=-1;
    B(i,i+1)=-1;
end
A2(1,2)=-2;
B(1,2)=-2;
for i=2:n1
    A2(i,i-1)=-1;
     B(i,i-1)=-1;
end
B(n1,n1-1)=-2;
A2(n1,n1-1)=-2;
A=kron(A1,B)+kron(A2,I);
%%下面构造向量；
g=zeros(n1,n1);
for i=1:n1
  for j=1:n1
      g(i,j)=subs(un,{x,y},{xh(i),yh(j)}); 
  end
end
gg=zeros(n1^2,1);
p=0;
for i=n1+1:n1:n1^2-n1
    p=p+1;
    gg(i)=g(1,p+1);
    gg(i+n1-1)=g(n1,p+1);
end
for i=1:n1
    gg(i)=g(i,1); 
end
gg(1)=2*gg(1);
gg(n1)=2*gg(n1);
p=0;
for i=n1^2-n1+1:n1^2
    p=p+1;
    gg(i)=g(p,n1); 
end
gg(n1^2-n1+1)=2*gg(n1^2-n1+1);
gg(n1^2)=2*gg(n1^2);
gg1=gg(2:end);
AA1=A(2:end,2:end);
U=AA1\gg1;
U=[0;U];
UU=zeros(n1,n1);
for i=1:n1
    for j=1:n1
        UU(i,j)=U((i-1)*n1+j);
    end
end
UU=UU';
mesh(xh,yh,UU);%数值解图