e471.m

matlab算法集 matlab算法集

%----------------------------------------------------------------------
% Example 4.7.1: Pressure-Temperature Curves  
%----------------------------------------------------------------------

% Initialize
 
   clc                  % clear screen
   clear                % clear variables
   n = 11;              % number of data points 
   q = 10*n;            % number of plot points 	
   m = 4;               % number of coefficients 
   f = 1.e5;            % conversion factor  
   vg = 1.673;          % specific volume (gas) 
   vf = 0.00104;        % specific volume (liquid) 
   b = zeros (m,1);
   d = zeros (m-1,1);
   C = zeros (m,m);
   X = zeros (q,2);
   Y = zeros (q,2);
   T = [50 60 70 80 90 100 110 120 130 140 150]';
   p = [.1235 .1994 .3119 .4739 0.7014 1.014 ...
        1.433 1.958 2.701 3.613 4.758]';
   
% Find least-squares polynomial using equations (4.5.9)-(4.5.10) */

   fprintf ('Example 4.7.1: Pressure-Temperature Curves\n');
   for k = 1 : m
      for i = 1 : n
         b(k) = b(k) +  p(i)*T(i)^(k-1);
      end
      for j = 1 : m
         for i = 1 : n
            C(k,j) = C(k,j) + T(i)^(k+j-2);
         end
      end
   end 
   a = gauss (C,b);

% Find slope at T = 100    
   
   show ('Polynomial coefficients',a)
   for i = 1 : m-1
      d(i) = i*a(i+1);
   end
   show ('Derivative coefficients',d)
   s = polynom (d,100);
   show ('Slope (bars/K)',s);

% Find change in enthalpy 

   dh = 373.15*(vg - vf)*f*s;
   show ('Change in enthalpy (kJ/kg)',dh/1000);

% Plot curve and data 
   
   for i = 1 : q
      j = min ([i,n]);
      X(i,1) = T(1) + (T(n) - T(1))*(i-1)/(q-1);
      Y(i,1) = polynom (a,X(i,1));
      X(i,2) = T(j);
      Y(i,2) = p(j);
   end
   graphxy (X,Y,'Pressure-Temperature Curves','T (^oC)','p (bars)','s')
%----------------------------------------------------------------------