Home > wafo > trgauss > mc2rfc.m

mc2rfc

PURPOSE ^

Calculates a rainflow matrix given a Markov chain with kernel f_xy;

SYNOPSIS ^

f_rfc = mc2rfc(f_xy,paramv,paramu)

DESCRIPTION ^

 MC2RFC  Calculates a rainflow matrix given a Markov chain with kernel f_xy; 
        f_rfc = f_xy + F_mc(f_xy). 
  
   CALL: f_rfc = mc2rfc(f_xy); 
  
   where 
  
         f_rfc = the rainflow matrix, 
         f_xy =  the frequency matrix of Markov chain (X0,X1) 
                 but only the triangular part for X1>X0.  
  
   Further optional input arguments; 
  
   CALL:  f_rfc = mc2rfc(f_xy,paramx,paramy); 
  
        paramx = the parameter matrix defining discretization of x-values, 
        paramy = the parameter matrix defining discretization of y-values,

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f_rfc = mc2rfc(f_xy,paramv,paramu) 
002 %MC2RFC  Calculates a rainflow matrix given a Markov chain with kernel f_xy; 
003 %       f_rfc = f_xy + F_mc(f_xy). 
004 % 
005 %  CALL: f_rfc = mc2rfc(f_xy); 
006 % 
007 %  where 
008 % 
009 %        f_rfc = the rainflow matrix, 
010 %        f_xy =  the frequency matrix of Markov chain (X0,X1) 
011 %                but only the triangular part for X1>X0.  
012 % 
013 %  Further optional input arguments; 
014 % 
015 %  CALL:  f_rfc = mc2rfc(f_xy,paramx,paramy); 
016 % 
017 %       paramx = the parameter matrix defining discretization of x-values, 
018 %       paramy = the parameter matrix defining discretization of y-values, 
019 %       
020 if nargin<2 
021 paramv=[-1, 1, length(f_xy)]; 
022 paramu=paramv; 
023 end 
024  
025 if nargin<3 
026 paramu=paramv; 
027 end 
028 dd=diag(rot90(f_xy)); 
029 N=length(f_xy); 
030 Splus=sum(f_xy'); 
031 Sminus=fliplr(sum(f_xy)); 
032 Max_rfc=zeros(N,1); 
033 Min_rfc=zeros(N,1); 
034 norm=zeros(N,1); 
035 for i=1:N 
036   Spm=Sminus(i)+Splus(i)-dd(i); 
037   if Spm>0. 
038     Max_rfc(i)=(Splus(i)-dd(i))*(Splus(i)-dd(i))/(1-dd(i)/Spm)/Spm; 
039     Min_rfc(i)=(Sminus(i)-dd(i))*(Sminus(i)-dd(i))/(1-dd(i)/Spm)/Spm; 
040     norm(i)=Spm; 
041   end 
042 end 
043  
044 %cross=zeros(N,1); 
045 %for i=2:N 
046 %   cross(N-i+1)=cross(N-i+2)+Sminus(N-i+2)-Splus(N-i+2); 
047 %end 
048  
049 f_rfc=zeros(N,N); 
050 f_rfc(N-1,1)=Max_rfc(N-1); 
051 f_rfc(1,N-1)=Min_rfc(2); 
052  
053 for k=3:N-1 
054      for i=2:k-1, 
055  
056 %       AAe= f_xy(1:N-k,1:k-i); 
057 %       SAAe=sum(sum(AAe)); 
058        AA = f_xy(N-k+1:N-k+i,k-i+1:k); 
059        RAA=f_rfc(N-k+1:N-k+i,k-i+1:k); 
060        nA=length(AA); 
061        MA= Splus(N-k+1:N-k+i); 
062        mA=Sminus(N-k+1:N-k+i); 
063        normA=norm(N-k+1:N-k+i); 
064        MA_rfc=Max_rfc(N-k+1:N-k+i); 
065        mA_rfc=Min_rfc(k-i+1:k); 
066        SA=sum(sum(AA)); 
067        SRA=sum(sum(RAA)); 
068        SMA_rfc=sum(MA_rfc); 
069        SMA=sum(MA); 
070        DRFC=SA-SMA-SRA+SMA_rfc; 
071         
072        NT=MA_rfc(1)-sum(RAA(1,:)); 
073  
074 %       if k==35 
075 %          check=[MA_rfc(1) sum(RAA(1,:))] 
076 %          pause 
077 %       end 
078  
079     NT=max(NT,0); 
080  
081     if NT>1e-6*MA_rfc(1) 
082  
083      NN=MA-sum(AA'); 
084      e=(fliplr(mA)-sum(AA))'; 
085      e=flipud(e); 
086      AA=AA+flipud(rot90(AA,-1)); 
087      AA=rot90(AA); 
088          AA=AA-0.5*diag(diag(AA)); 
089  
090  
091      for j=1:nA, 
092           if normA(j)~=0 
093             AA(j,:)=AA(j,:)/normA(j); 
094             e(j)=e(j)/normA(j); 
095           end 
096        end 
097        fx=0.;           
098  
099        if max(abs(e))>1e-7 & max(abs(NN))>1e-7*MA_rfc(1) 
100        
101   
102        I=eye(size(AA)); 
103  
104          if nA==1 
105         fx=NN/(1-AA)*e; 
106          else 
107         fx=NN*((I-AA)\e); 
108          end 
109        end 
110         
111     f_rfc(N-k+1,k-i+1)=DRFC+fx; 
112     end 
113    end 
114   m0=max(0,Min_rfc(N)-sum(f_rfc(N-k+2:N,1))); 
115   M0=max(0,Max_rfc(N-k+1)-sum(f_rfc(N-k+1,2:k))); 
116   f_rfc(N-k+1,1)=min(m0,M0); 
117 %  n_loops_left=N-k+1 
118  end 
119  
120 for k=2:N 
121   M0=max(0,Max_rfc(1)-sum(f_rfc(1,N-k+2:N))); 
122   m0=max(0,Min_rfc(k)-sum(f_rfc(2:k,N-k+1))); 
123   f_rfc(1,N-k+1)=min(m0,M0); 
124 end 
125 f_rfc=f_rfc+rot90(diag(dd),-1); 
126 clf 
127 subplot(1,2,2) 
128 pcolor(levels(paramv),levels(paramu),flipud(f_xy+flipud(rot90(f_xy,-1)))) 
129   axis([paramv(1), paramv(2), paramu(1), paramu(2)]) 
130   title('MC-kernel  f(x,y)') 
131   ylabel('y'), xlabel('x')                     
132 axis('square') 
133  
134 subplot(1,2,1) 
135 pcolor(levels(paramv),levels(paramu),flipud(f_rfc)) 
136   axis([paramv(1), paramv(2), paramu(1), paramu(2)]) 
137   title('Rainflow matrix') 
138   ylabel('max'), xlabel('rfc-min')                     
139 axis('square') 
140  
141  
142  
143  
144  
145  
146  
147  
148

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003