How do I only save my solution every x iterations in a for loop?

I'm solving a PDE and iterating in time in a for loop. However, I don't want to store the solution at every single time step. Instead I want to be able to define how many times the solution should be stored.

I've tried this:

incrementer = 2
for itimestep = itimestepvector

    currentsolution = update_solution(state,currentsolution); 

    if (mod(itimestep,maxitimestep/numberoftimesIwantorecord) == 0)
        solutions(:,incrementer) = currentsolution;
        incrementer = incrementer + 1;
    end
end

The problem is this code doesn't save the solution in time "numberoftimesIwantorecord" number of times and I understand why it doesn't. How do i make the solution record "numberoftimesIwantorecord" number of times?

[Stephen23]

@J2A2B2: What exactly is the problem?

[J2A2B2]

sorry I missed that part out. updated.

Answer 1

[David Goodmanson]

Hello J,

If *numberoftimesIwantorecord* does not exactly divide *maxitimestep*, then there are two problems with this approach.

[1] Let's call the variables numof and maxitimes respectively. The time index of course runs from 1 to maxitimes. Suppose

maxitimes = 100
numof = 12
find(mod((1:maxitimes),maxitimes/numof)==0)
ans = 25    50    75   100

Only works four times, not twelve. This is happening because you are looking for a time index itime where (in this case)

mod(itime,100/12) == 0   -->   itime = m*(100/12)  % for some m
or
12*itime = m*100

3 divides the left hand side but does not divide 100, so 3 must divide m. Then m = 3*q where q is some other integer. Then

itime = q*25

and you get the result with only four answers. For equally spaced samples, many values of *numberoftimesIwantorecord* aren't even possible.

[2] There are floating point issues lurking around if maxitimes/numof is not an integer and there is a [does something ==0] test. I did get away with it here, but in general this doesn't work.

I believe you would be better off with, for example,

timespacing = 7;
find(mod(1:100,timespacing)==0) 
ans =  7    14    21    28    35    42    49    56    63    70    77    84    91   98

which is strictly an integer process.

[J2A2B2]

Thanks David

Answer 2

[Rik]

The problem is in the rounding. You can improve the robustness of your code by rounding maxitimestep/numberoftimesIwantorecord, but that still won't solve all cases.

maxitimestep=10;
for numberoftimesIwantorecord=1:10
    recorded_with_round=sum(mod(1:maxitimestep,round(maxitimestep/numberoftimesIwantorecord))==0);
    recorded_without_round=sum(mod(1:maxitimestep,(maxitimestep/numberoftimesIwantorecord))==0);
    fprintf('%2.0f times requested, %2.0f times done (%2.0f times without round)\n',...
        numberoftimesIwantorecord,recorded_with_round,recorded_without_round)
end

This results in

1 times requested,  1 times done ( 1 times without round)
2 times requested,  2 times done ( 2 times without round)
3 times requested,  3 times done ( 1 times without round)
4 times requested,  3 times done ( 2 times without round)
5 times requested,  5 times done ( 5 times without round)
6 times requested,  5 times done ( 2 times without round)
7 times requested, 10 times done ( 1 times without round)
8 times requested, 10 times done ( 2 times without round)
9 times requested, 10 times done ( 1 times without round)
10 times requested, 10 times done (10 times without round)