How to avoid a numerical problem in a for loop

Question

Eli Dim 2015-7-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/230097-how-to-avoid-a-numerical-problem-in-a-for-loop

回答： Star Strider 2015-7-16

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2015-7-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/230097-how-to-avoid-a-numerical-problem-in-a-for-loop#answer_186295

在 MATLAB Online 中打开

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.