solvelucgsor1219.m

自己写的求解波松方程的matlab程序 并采用两种方法求解

clear all
%---------------------------矩阵A--------------------------------------------
N=9;
n=N*N;
  A=zeros(N*N,N*N);
   A(1,1)=4;
   for i=2:N*N;
        A(i,i)=4;
        if (i>9) A(i,i-N)=-1; A(i-N,i)=A(i,i-N);end
        if mod((i-1),N)~=0
       A(i,i-1)=-1;A(i-1,i)=A(i,i-1);end
end
%--------------------------常数项B------------------------------------------
h=0.1;
F=zeros(N,N);
f=inline('-(x.^2+y.^2)*exp(x.*y)','x','y');
for y=1:N;
    for x=1:N;
        F(x,y)=h^2*f(x*h,y*h);
   end
end
 
for i=1:N;
    bottombound(i)=1;
    leftbound(i)=1;
    upbound(i)=exp(i*h);
    rightbound(i)=exp(i*h);
end
%------------------------------------上下边界 ----------------------------------
for i=1:N;
    F(i,1)=F(i,1)+bottombound(i);
    F(i,N)=F(i,N)+upbound(i);
 end
  %-----------------------------------左右边界----------------------------------
  for j=1:N;
      F(1,j)=F(1,j)+leftbound(j);
      F(N,j)=F(N,j)+rightbound(j);
  end
%   B=zeros(N*N,1);
%   for i=0:N-1
%     B(9*i+1:9*(i+1))=F(1:9,i+1);
%   end
B=[];
for i = 1:N
    B = [B; F(:,i)];   %B = [B, F(:,i)]，逗号与冒号差别
end

 %----------------------------lu分解求解---------------------------------------------
 
                         tic    %秒表开
%  tm=cputime;
                         disp('lu分解求解')
 [L,U]=lu(A);
 y=L\B;
 x=U\y;
%  cputime-tm
 uu=zeros(N,N);
   for i=0:N-1;
   uu(1:9,i+1)=x(9*i+1:9*(i+1));
   end
                            toc
 %-----------------------------------------------------------------------秒表关
                           tic %秒表开
 %----------------------------------CG迭代求解-------------------------------                       
                        disp('CG迭代求解')
   xx0=ones(81,1);
   uuu=CG(A,B,xx0);         
                          toc   
 %------------------------------------------------------------------------秒表关 
 % -----------------------------------SOR迭代求解------------------------------
                          disp('SOR迭代求解')
                          tic    %秒表开
x0=ones(81,1);
U=sor(A,B,x0);
toc    
%--------------------------------------------------------------------------秒表关 
%-----------------------可以简化------------------------------------------------------%
    BB=[1 bottombound(1:9) 1;
    leftbound(1),uu((1:9),1)',rightbound(1);
    leftbound(2),uu((1:9),2)',rightbound(2);
    leftbound(3),uu((1:9),3)',rightbound(3);
    leftbound(4),uu((1:9),4)',rightbound(4);
    leftbound(5),uu((1:9),5)',rightbound(5);
    leftbound(6),uu((1:9),6)',rightbound(6);
    leftbound(7),uu((1:9),7)',rightbound(7);
    leftbound(8),uu((1:9),8)',rightbound(8);
    leftbound(9),uu((1:9),9)',rightbound(9);
    1, upbound(1:9), 2.7818;];
%-------------------------------------------------------------------------%
  x=0.1:0.1:0.9;% 绘制9*9网格
  y=0.1:0.1:0.9;
  mesh(x,y,uu);    
  figure
  surf(uu);% 绘制11*11网格
  x=0:0.1:1;
  y=0:0.1:1;
  [X,Y]=meshgrid(x,y);
  figure
  mesh(X,Y,BB)
  figure
  surf(X,Y,BB)
  cputime
  % 计算系统运行时间的命令  tic toc cputime
  % 绘图函数   mesh surf