How to avoid a numerical problem in a for loop
    4 vues (au cours des 30 derniers jours)
  
       Afficher commentaires plus anciens
    
    Eli Dim
 le 16 Juil 2015
  
    
    
    
    
    Réponse apportée : Star Strider
      
      
 le 16 Juil 2015
            I have two input files T and dispatch_seq. My code looks at T and sees how many rows of dispatch_seq should be taken. The result is stored in Output.
    load T.mat
    load dispatch_seq.mat
    for j = 1:length(T)
        ph1 = cumsum(cell2mat(dispatch_seq{j}(:,2)));
        DecR1 = 10*rem(ph1,1);
        for k = 1:length(ph1)
            if DecR1(k) <= 0.1           
                ph1(k) = round(ph1(k));
            end 
        end
        j_1 = find(ph1>=T(j),1);
        Output(j) = {dispatch_seq{j}(1:j_1,:)};
        if ph1(j_1) > T(j) 
            Output{j}{end,2} = (Output{j}{end,2}-ph1(j_1)+T(j));
        end 
    end
Everything works as expected with the exception of Output{28,1} and Output{77,1}. What happens there is the number in the corresponding row of the vector T equals exactly the sum of all rows of dispatch_seq, but the result I get in Output is []. How can I avoid this error?
0 commentaires
Réponse acceptée
  Star Strider
      
      
 le 16 Juil 2015
        I apologise for the delay. Life intrudes.
I discovered a problem. When I inserted this line as a trace of the arguments to the ‘j_1’ assignment, just below your ‘j_1’ assignment:
fprintf(1,'\nmax(ph1) = %23.15E, T(%d) = %23.15E, -> max(ph1)-T(%d) = %23.15E\n', max(ph1), j, T(j), j, max(ph1)-T(j))
the result was:
max(ph1) =   4.086999999999996E+03, T(28) =   4.087000000000000E+03, -> max(ph1)-T(28) =  -3.637978807091713E-12
but the difference was zero or positive for all prior values of ‘ph1’ and ‘T(j)’. So for j=28, the find call returns an empty value for ‘j_1’. That results in the empty elements in your cell array.
Since you’re testing for ph1 >= T(j), one option might be to subsequently define j_1=length(ph1) if ‘j_1’ is empty, or to introduce a tolerance, or just round ‘ph1’ such as:
j_1 = find(round(ph1)>=T(j),1);
Since I’m not certain what your code is doing (so it may not be an appropriate solution), I offer that only as a suggestion.
0 commentaires
Plus de réponses (0)
Voir également
Catégories
				En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

