sor.m

计算方法中的解线性方程组的迭代法

clear;
close all;
h=10;k=30;md=7300;cp=800;t0=1000;tw=0;tsp=600;le=0.1;%tw为环境温度
w=0.5;%此W的值不是最优,有待改善~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n_space = input('n_space?   ');
n_space=n_space+1;
n_time = input('n_time?   ');
alfa=tsp/n_time/(le*le/n_space/n_space)*k/md/cp;
alfa1=1+2*alfa;

beta=(le/n_space)*h/k;

x(1,1)=tw*beta;
x(n_space,1)=tw*beta;
x(2:(n_space-1),1)=t0;
b=x;
%矩阵计算:
D(n_space,n_space)=1+beta;
D(1,1)=1+beta;
L(n_space,n_space)=0;
L(n_space,n_space-1)=1;
U(n_space,n_space)=0;
U(1,2)=1;

for ii=2:(n_space-1)
    D(ii,ii)=alfa1;
    L(ii,ii-1)=alfa;
    U(ii,ii+1)=alfa;   
end
MM=pinv(D-w*L);
loopmatrix=MM*((1-w)*D+w*U);     %%%%%%%%%%%loopmatrix
MM=w*MM;
    
    for tt=1:n_time%对时分的循环
        
    
    
          for    ii=1:100 %某时间片一次SOR循环 循环控制条件 有待改进~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
              x=loopmatrix*x+MM*b;
              %x=x1;    
          end
          
          if tt<n_time
              b=x;
              b(1,1)=0;
              b(n_space,1)=0;
              
          end
    
   end