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-1A2(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);%数值解图