Метод Зейделя для статической задачи о распространении загрязняющих веществ

a=0;b=1;T=1;N=30;M=25;mu=1; %1e-3;%ЛПЬЖЖЙГЙЕОФ ДЙЖЖХЪЙЙ РБТПЧ ЬФЙМПЧПЗП УРЙТФБ Ч ЧПЪДХИЕ
tau=T/M; h=(b-a)/N;
Matr=zeros((N-1)^2,(N-1)^2);
A=zeros(N-1,N-1);
for k=1:N-1
if k==1
    A(k,k)=h^2/tau+4;
    A(k,k+1)=-mu;
else
    if k==N-1
        A(k,k)=h^2/tau+4;
        A(k,k-1)=-mu;
    else
        A(k,k)=h^2/tau+4;
        A(k,k-1)=-mu;
        A(k,k+1)=-mu;
    end;
end;
    end;
 
for k=0:N-2
    Matr(k*(N-1)+1:k*(N-1)+N-1,k*(N-1)+1:k*(N-1)+N-1)=A;
    if k==0
        Matr(k*(N-1)+1:k*(N-1)+N-1,(k+1)*(N-1)+1:(k+1)*(N-1)+N-1)=-mu*eye(N-1);
        Matr(k*(N-1)+1:k*(N-1)+N-1,(N-1)^2+1:(N-1)^2+(N-1))=-eye(N-1);
        
    else
        if k==N-2
            Matr(k*(N-1)+1:k*(N-1)+N-1,(k-1)*(N-1)+1:(k-1)*(N-1)+N-1)=-mu*eye(N-1);
        else
            Matr(k*(N-1)+1:k*(N-1)+N-1,(k-1)*(N-1)+1:(k-1)*(N-1)+N-1)=-mu*eye(N-1);
            Matr(k*(N-1)+1:k*(N-1)+N-1,(k+1)*(N-1)+1:(k+1)*(N-1)+N-1)=-mu*eye(N-1);
            
        end;
    end;
end;
Matr((N-1)^2+1:(N-1)^2+(N-1),1:N-1)=-eye(N-1);
Matr((N-1)^2+1:(N-1)^2+(N-1),(N-1)^2+1:(N-1)^2+(N-1))=eye(N-1);
     
v0=zeros(1,(N-1)^2+N-1);
%for i=2:N-20
  %  for j=2:N-20
   %     v0(i+(j-1)*(N-1))=10;
 %   end;
%end;
x=0+h:h:1-h;
y=0+h:h:1-h;
for i=1:N-1
    for j=1:N-1 
        if ((x(i)-0.5)^2+(y(j)-0.5)^2<0.09)
            v0(i+(j-1)*(N-1))=10;
        end;
    end;
end;
 
u=zeros(N-1,N-1);
    for i=1:N-1
        for j=1:N-1
            u(i,j)=v0(i+(j-1)*(N-1));
        end;
    end;
    figure;
surf(x,y,u);colormap(jet); view(90,270); shading flat; colorbar;axis tight;title(strcat('t=',num2str(0),', max=',num2str(max(max(u))))); drawnow;
    pause(0.5);
 
f=zeros(1,(N-1)^2+N-1);
    %s='';
%for i=1:N-2
  %  s= strcat(s,'%8.6f_');
%end;
%s=strcat(s,'%8.6f\n'); 
%formatSpec=s;
%fileID = fopen('matrica.txt','w');
 
for k=1:M
    for i=1:(N-1)
        for j=1:N-1
            f(i+(j-1)*(N-1))=h^2/tau*v0(i+(j-1)*(N-1))+h^2*sin(x(i)^2)*sin(y(j));%*x(i)*y(j)%sin(x(i))*sin(y(j))%x(i)*sinh(y(j))%exp(x(i))*y(j)^8;
        end;
    end;
    v1=Matr\f';
    %v1=zeid(0.001,Matr,f);
    for i=1:N-1
        for j=1:N-1
            u(i,j)=v1(i+(j-1)*(N-1));
        end;
    end;
    set(gcf,'renderer','painters');
    surf(x,y,u);colormap(jet); view(90,270);shading interp; colorbar;axis tight;title(strcat('t=',num2str(tau*k),', max=',num2str(max(max(u))))); drawnow;
    pause(0.01);
   % fprintf(fileID, formatSpec,u);
    %fprintf(fileID,'\n\n');
 
    v0=v1';
 
end;
%fclose(fileID);

function[x2]=zeid(eps,A,b)
s=size(b);
N=s(2);
x1=zeros(1,N);
x2=zeros(1,N);
norma=100500;
while norma>eps
    x1=x2;
    for i=1:N 
        sum1=0;
        sum2=0;
        for j=1:i-1
            sum1=sum1+A(i,j)*x2(j);
        end;
        for j=i+1:N
            sum2=sum2+A(i,j)*x1(j);
        end;
        x2(i)=(b(i)-sum1-sum2)/A(i,i);
     end;
        norma=norm(x2-x1);
end;