雷诺方程的数值解法

				function [A,b]=add_ps()
				A=evalin('base','A1');
				b=evalin('base','b1');
				M=evalin('base','M');
				N=evalin('base','N');
				pa=evalin('base','pa');
				p0=evalin('base','p0');
				ps_location=evalin('base','ps_location');
				beta=evalin('base','beta');
				psbeta=[ps_location,beta*p0/pa];    %供气点位置，以及压力比
				num=size(psbeta);   %num(1)供气点个数
				for k=1:num(1)     %列供气点处的方程
				    row=[];
				    row(M*(N+1))=0;
				    i=psbeta(k,1);          %得到供气点i坐标
				    j=psbeta(k,2);          %得到供气点j坐标
				    ps_beta=psbeta(k,3);    %得到供气点的供气压力比
				    row((i-1)*(N+1)+j)=1;
				    A((i-1)*(N+1)+j,:)=row;
				    b((i-1)*(N+1)+j)=ps_beta^2;
				end
				end