How to calculate values using for-end and if-else in Octave?

user_12918426 0
Viewed 111

I need to calculate variable Nfj in the following equation:

ea_M = (sf/E)*(2*Nfj)^b + ef*(2*Nfj)^c

so that it is equal to 4 different values of ea: 5.0900e-04, 4.3626e-04, 3.6358e-04, and 2.9084e-04. The results should be 4 values of Nfj which should be stored in results. I am rounding because I think that it is not possible to calculated exactly equal values of ea, without rounding the equality would always be false.

I wrote for it this code, but it does not work as expected: the script runs for a very long time without any results. How can I fix it to work correctly?

sf = 882.07;
ef = 0.59;
b = -0.102969;
c = -0.58;
E = 210000;
ea = [5.0900e-04; 4.3626e-04; 3.6358e-04; 2.9084e-04]

for pos = 1:length(ea)
  for Nfj = 1e3:10:1e12
    ea_M = (sf/E)*(2*Nfj)^b + ef*(2*Nfj)^c;
    if  round(ea_M * 10^5)/10^5 == round(ea(pos) * 10^5)/10^5;
      disp(ea_M)
      disp(Nfj)
      results(pos) = Nfj;
    end
  end
end
octave
1 Answers

You are trying every 10th value, from 103 to 1012: there are 1011 values here to try! This is of course a lot. You will be searching for ever, and you might skip the actual value that would make the equality true.

If you are not able to solve the equation manually, you can use a numerical solver. Let's first print the function:

sf = 882.07;
ef = 0.59;
b = -0.102969;
c = -0.58;
E = 210000;

ea_M = @(Nfj) (sf/E)*(2*Nfj).^b + ef*(2*Nfj).^c;

Nfj = logspace(3,12,1000);
plot(Nfj, ea_M(Nfj))
set(gca, 'xscale', 'log')

plot generated by code above

It looks like this function is nicely monotonic, and the four values you are looking for are in the interval of 103 to 1012 you were searching. To find where it equals one of your values, we can subtract that value and find where it equals zero. You can very quickly narrow down your search if you start at a point where the function is larger than zero, and one where it is smaller than zero. You halve the interval every time, keeping the half of the interval that contains the zero crossing. The fzero function does just this.

ea = [5.0900e-04; 4.3626e-04; 3.6358e-04; 2.9084e-04];
results = zeros(size(ea));
for pos = 1:numel(ea)
   results(pos) = fzero(@(Nfj) ea_M(Nfj) - ea(pos), [1e3,1e12]);
end
results

In MATLAB, this code runs in a small fraction of a second and outputs:

results =

   1.0e+10 *

    0.0429
    0.1849
    1.0627
    9.1919

and ea_M(results) - ea is approximately zero.

A few notes on the code I posted here:

  1. ea_M is defined as an anonymous function. This makes it easier to reuse the expression, rather than writing it over and over again. I replaced ^ with .^ to allow this function to do the computation on an array of Nfj values at once, rather than only single values. This is necessary for the fzero call.

  2. I plotted the function on a logarithmic scale, because this function calls for that. The same is not true for all functions.

  3. I preallocated the results array, you should avoid increasing the size of an array in a loop.

user_73287820
Related