How to change many loops to recursive function?

Hi, i would like to replace for loops to recursive function. I also would like to have my 5 variables m,al,R0,R1,e in array of (min, max, step,name) because may be in the future i might have more variables so i just add more array and change my function f.

function xm = minfricloop
    
    xm = [ Inf  Inf  Inf  Inf  Inf  Inf];   %For minimization
    
    mmin = 3;
    mmax = 12;
    almin = 0.2;
    almax = 1;
    R0min = 0.014;
    R0max = 0.020;
    R1min = 0.025;
    R1max = 0.030;
    emin = 0.5;
    emax = 0.9;
    rpm = 20000;
    ita = 0.1;
    h0 = 2.51849556599183e-05;
    n = 5;
   

    
   for v = 0 :20
       
      m = mmin + (mmax-mmin)*(v/(20));
        for w = 0 :20
             al = almin + (almax-almin)*(w/(20));
             
             for  x = 0 :20
                R0 = R0min + (R0max-R0min)*(x/(20));
                for y = 0 :20
                    R1 = R1min + (R1max-R1min)*(y/(20));
                    for z = 0 :20
                         e = emin + (emax-emin)*(z/(20));
                       
         k   =  (2*pi*e*(R0+R1)/(2*m))*tan(al*pi/180)/h0;
      

        f = (R1-R0)*(rpm*2*pi*(R0+R1)/(2*60))*ita*(2*pi*e*(R0+R1)/(2*m))*((6/(k+2))+(4*log(k+1)/k))/h0;

       
       if f < xm(6) %For minimization
           xm = [m al R0 R1 e f];
       end
                    end    
                end
             end
        end
   end
end

Answer 1

[Walter Roberson]

Code attached showing how to use ndgrid for this kind of work.

On my system execution time is less than 0.15 seconds.

[Paranee Sawad]

Thank you very much Walter Roberson. ^___^

[Walter Roberson]

Note that each array would have N^5 entries and you have N of them, each entry double precision, so with N = 21, that is 21^5 * 5 * 8 = about 160 megabytes. Plus additional temporary copies as it processes the expressions.160 megabytes is not bad, but if you were to increase N to get a finer resolution then the memory could go up quickly. By N = 31 you would be using more than 1 gigabyte.

[Paranee Sawad]

Dear Walter Roberson,

Thank you for reminding of the N I have already tried and checked the time. And i also have one more question if i would like to put a constraint about 1.2 < k_ < 1.55 . Is it possible to add this since i try to find the optimal value?

[Walter Roberson]

Modified version attached.

For situations like this, there are two basic strategies:

1. Calculate at all locations, but only take the minima over the locations that meet the criteria; or
2. Filter the coefficients down to the combinations that meet the criteria, and calculate the main expression for those locations

When most of the locations meet the criteria then calculate-anyhow-and-reject can be faster. When few of the locations meet the criteria, then reject-then-calculate is likely to be faster.

You might think that rejecting first would always be faster, but you need to apply the rejection to each of the variables individually, which has a cost that has to be balanced against the cost of the extra computation, and the cost of extra computation is not always obvious because of automatic multi-threaded calculations in some cases.

[Paranee Sawad]

Aha i got it. Thank you very much (:

[Paranee Sawad]

Dear Walter Roberson,

From previous, i set h0 to be constant but for now i also want to compute it by each value of my parameters (min, max) which will give different h0. I compute h0 by Newton Method until it makes fu-frad = 0. And at the end it will compute min f_. According to your recommendation, it seems like i should let them compute all possible value and then find the min f_. Do i understand it correctly?

[Walter Roberson]

In my original code, I used the naming convention that all of the multidimensional arrays of variables are named with trailing _ such as al_ . In the revised code that selects "good" positions, the resulting vectors are named with _g such as al_g .

As I go through your code above, I noticed that you did not use any naming convention about which variables are scalars and which are vectors (or arrays.) Because of that, you have accidentally ended up in a situation where you are calculating vectors of values but you are thinking of them as if they are scalars. In particular, when you are working with arrays, you have one error value for each array location, but your while loop is testing as if you are working with a scalar. You need to figure out what you want to do in this case: do you want to continue the loop only until the condition is false in one location (your current code), or do you want to continue the loop until the condition is false everywhere (that is, all entries fall into range) ? If you are not certain it is possible for the condition to be false everywhere, then you might need to run the calculation over all of the locations individually.

[Paranee Sawad]

Dear Walter Roberson,

In while loop, you are right since i chage the value of errF = 1e-4 to errF = 1e-9 the optimal value will be changed. So if i want h0 to be computed in all entries can i just only change the errF = 1e-20? to make sure that the err = abs(f/frad); has already run for all.

Therefore in my while loop i should also change my variables for example, fu_, kk_ right? since they are computed with variables in array.

[Walter Roberson]

The problem is that your err variable is a multidimensional array, but you are asking while err > errF . What doe you want that to mean? Do you want it to stop as soon as any member of err is <= errF ? Do you want it to continue until all members of err are <= errF? If you want it to continue until all members are <= errF, are you sure that it is possible to have that happen for all of the members of the array simultaneously ?

[Paranee Sawad]

Dear Walter Roberson,

I understand my problem now. I am really sorry because i have no background about coding before. So when i run my code it will start at err(:,:,14,18) until err(:,:,21,21) that means i lost many values of them and each of them is 21x21. When compare err > errF my errF which is just a scalar number cannot compare to err so for what i understand now i should also set the errF to be matrix of 21x21 right?

[Walter Roberson]

Comparing a matrix to a scalar is fine.

V = [1 2 3 4 5 6 7 8 9];
V >= 5

gives

[false false fasle false true true true true true]

But if you were to test

while V >= 5

then because not all of the logical results are true, the test would be considered to be false, just as if you had written

while all(V >= 5)

But sometimes what you want is

    while any(V >= 5)

But you need to decide which test is appropriate for you.

[Paranee Sawad]

Dear Walter Roberson,

since i set

   while all(err > errF) 

and when i check the value in array of

        err = abs(f/frad);

why they are not convergent to the value that will satisfy

err < errF

but if i set

       while any(err > errF)       

instead, i can notice that most of the value tend to satisfy

err < errF

Can i say because of this then i decide to use while any(err > errF)

[Walter Roberson]

Remember with the all() you will stop as soon as the condition is false for even one entry, so as soon as the very first entry gets within tolerence you would stop. Most of the entries would tend to be out of range. When you use any() you will not stop as long as even one entry is out of tolerance, so when you stop then all of them will be within tolerance.

What you have to be careful of is that it becomes possible that in order to satisfy the tolerance for one value, that another value has to be outside of tolerance. Perhaps that is not a problem for this particular set of equations, but it is a general problem when you are working on a search grid if the results at one location can affect the calculation at a different location.

[Paranee Sawad]

Dear Walter Roberson,

Aha! Thank you very much i got it ^__^

[Paranee Sawad]

Dear Walter Roberson,

since in

while any(err > errF)

i have computed to find h0 with 4 input variables

[m_, al_, R0_, R1_] = ndgrid( linspace(mmin, mmax, N), ...
                              linspace(almin, almax, N), ...
                              linspace(R0min, R0max, N), ...
                              linspace(R1min, R1max, N));

then after that i would like to find

f_ = (R1_ - R0_) .* (rpm .* 2 .* pi .* (R0_ + R1_) ./ (2 .* 60)) .* ita .* (2 .* pi  .* (R0_ + R1_) ./ (2 .* m_)) .* ((6 ./(k_ + 2)) + (4 .* log(k_+ 1) ./ k_)) ./ h0_;

but in this case [m_, al_, R0_, R1_] would not be the same value that used to compute h0 right? so in this case in

while any(err > errF)

How can i let them keep the value [m_, al_, R0_, R1_] that related to h0 so before they left the loop they can keep that value to compute f_.

[Walter Roberson]

You haven't posted your latest code.

[Paranee Sawad]

Dear Walter Roberson,

You can see the attached file. And for what i wonder is that if i change the value for

frad = 220;
rpm = 20000;
     
%     frad = [220;256;343;462;567;657;772;899;1012;1059;994;537;100];
%     rpm = [20000;22500;25000;27500;30000;32500;35000;37500;40000;42500;45000;47500;50000];

it seems like everything in

xm = [m_(idx), al_(idx), R0_(idx), R1_(idx), h0(idx) , f_(idx)];

will be the same except the value for f_(idx). because they didnt use the relative value of [m_, al_, R0_, R1_] which give h0.

[Walter Roberson]

It turns out not to be the case that the same location is always chosen. I have attached code that shows otherwise.

The location chosen can involve either the first m value, or the last m value. It is always the first al value, the first R0 value, and the last R1 value.

Five of your 13 cases involve complex-valued h0 and f_ . This happens because h0 can potentially be negative, and when it is then k can become less than -1, after which log(k+1) is log of a negative number, leading you to complex results.

[Paranee Sawad]

Dear Walter Roberson,

So i should set h0 > 0, right? since if h0 can potentially be negative it will effect k to become less than -1 and my result will be in complex number which is not allowed to happen in reallity for my case.

[Paranee Sawad]

Dear Walter Roberson,

I know it i should change the coefficient in the equation when i let them compute the h0 in the while loop

h0    = h0 - 0.01*f./jacobian;

Thank you very much (:

[Walter Roberson]

Just before you calculate error, set

h0 = max(h0, realmin);

This will prevent h0 from going negative, and so will prevent complex values from being calculated. Somehow it is able to satisfy the error conditions.

[Paranee Sawad]

Dear Walter Roberson,

Thank you very much (:

and now i try to recheck the value back, for example in attached file G.m this is the original code when i use to compute h0 if i know the value for m, al, R0, R1 but it seems like when i substitute it didnt return the value h0 same as the value that compute in file minftnormalit123.m