possion.txt

利用1五点差分格式解决的方程为Possion方程

function lapacezjp2(n)
f=zeros((n-1),(n-1));       %本解求的是81×81的矩阵解
f(1,1)=-(20+80);f(1,[2:end])=-20;%top
f([2:end-1],1)=-80;%left
f(n-1,1)=-(80+180);f(n-1,[2:end])=-180;%bottle
V=zeros(n+1,n+1);
V(1,[2:n])=20;V(n+1,[2:n])=180;
V([2:n],1)=80;V([2:n],n+1)=0;
D=-4*eye(n-1);
L=zeros(n-1,n-1);
for i=2:n-1
    L(i,i-1)=1;
end
U=L';
A=D+L+U;
E=eye(n-1);%上式分别求出(n-1)*(n-1)矩阵的分块矩阵A，E及f
for i=1:n-1
    c{i}=E;
    b{i}=A;
    F{i}=f(1,:);
end
  a=c;
bata{1}=c{1}*inv(b{1});%从这开始为追赶法a,b,c,f都设为元胞数组，元胞元素为矩阵
for i=2:n-2
    bata{i}=c{i}*inv(b{i}-a{i}*bata{i-1});
end
Y{1}=F{1}*inv(b{1});
for i=2:n-1
    Y{i}=(F{i}-Y{i-1}*a{i})*inv(b{i}-a{i}*bata{i-1});
end
X{n-1}=Y{n-1};
for i=n-2:-1:1
    X{i}=Y{i}-X{i+1}*bata{i};
end               %追赶法到此
x=zeros(n-1,n-1);
j=1:4:n-1;
m=length(j);
M=zeros(m,n-1);
for i=1:n-1         %取出元胞数组里面所有元素赋给矩阵
  x(i,:)=X{i}(1,:);
end
for i=1:m
     M(i,:)=X{i}(1,:);
 end
V([2:n],[2:n])=x;
M,surf(V)