gaussj.m

来自「电力系统计算节点数不多的潮流程序,数度比较快」· M 代码 · 共 100 行

function x=gaussj(A,b)
format long

[n,m]=size(A);
l=size(b);

if (n==m)&(n==l(1)),

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 输出到文件的控制  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   myf=fopen('output.dat','w');
   fprintf(myf, '---------------Linear equation input data----------\n');
    	for i=1:n,
		for j=1:n,
		
			fprintf(myf, '%8.4f ', A(i,j));
        end
		fprintf(myf, '| %8.4f', b(i));
		fprintf(myf, '\n');
    end
    fclose(myf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
for k=1:n,
    temp=abs(A(k,k));
    i0=k;
    for i=k+1:n,
        if abs(A(i,k))>temp
            temp=abs(A(i,k));         %%%%%% 选主元
            i0=i;
        end 
    end
    if temp==0
        disp('A is singular');
        break
    end
    
%%%%%%    交换i0行和k行  %%%%%%%%%%
     temp=A(k,:);
     A(k,:)=A(i0,:);
     A(i0,:)=temp;
    
     temp=b(k);
     b(k)=b(i0);
     b(i0)=temp;

%%%%%%    消元过程  %%%%%%%%%%         
    for i=(k+1):n,
        A(i,k)=A(i,k)/A(k,k);
        for j=(k+1):n,
        A(i,j)=A(i,j)-A(i,k)*A(k,j);
        end
        b(i)=b(i)-A(i,k)*b(k);
        A(i,k)=0;
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  输出到文件的控制 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 	myf=fopen( 'output.dat', 'a');	
	
	fprintf(myf, '---------------Gauss Column Primary Elemination Process----------\n');
    fprintf(myf, 'The %d-th elemination\n', k);
 	for i=1:n,
		for j=1:n,
		
			fprintf(myf, '%8.4f ', A(i,j));
        end
		fprintf(myf, '| %8.4f', b(i));
		fprintf(myf, '\n');
    end
    fclose(myf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

end

%%%%%%%%%%%  回代过程  %%%%%%%%%%
x(n)=b(n)/A(n,n);

for i=n-1:-1:1,
    temp=0;
    for j=i+1:n
        temp=temp+A(i,j)*x(j);
    end
    x(i)=(b(i)-temp)/A(i,i);
end    

x=x';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 输出到文件的控制 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
myf=fopen( 'output.dat', 'a');	
fprintf(myf, '---------------Linear equation output data----------\n');
	for i=1:n,
		fprintf(myf, 'x[%d] = %8.4f\n', i, x(i));
    end
        
status=fclose(myf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

else
     disp('the matrix dimensions of A and b must agree');
end