Using fzero in Matlab or Octave, avoiding for loop and complex solutions

Question

I'm using fzero function to solve a non-linear equation depending on one parameter and I'm not satisfied with my method. I have these issues:

1) Can for loop for the parameter be avoided ?

2) In order to avoid complex solutions I first have to pre-compute valid interval for fzero. Is there is a better solution here ?

If I reduce the parameter step size the execution time becomes slow. If I don’t pre-compute the interval I get an error "Function values at interval endpoints must be finite and real." in Matlab and "fzero: not a valid initial bracketing" in Octave.

Here is the code

% solve y = 90-asind(n*(sind(90-asind(sind(a0)/n)-y)))

% set the equation paramaters
n=1.48; a0=0:0.1:60;

% loop over a0
for i = 1:size(a0,2)

  % for each a0 find where the argument of outer asind() 
  % will not give complex solutions, i.e. argument is between 1 and -1 

  fun1 = @(y) n*(sind(90-asind(sind(a0(i))/n)-y))-0.999;
  y1 = fzero(fun1,[0 90]);  
  fun2 = @(y) n*(sind(90-asind(sind(a0(i))/n)-y))+0.999;
  y2 = fzero(fun2,[0 90]);


  % use y1, y2 as limits in fzero interval

  fun3 = @(y) 90-asind(n*(sind(90-asind(sind(a0(i))/n)-y)))-y;
  y(i) = fzero(fun3, [y1 y2]);

end

% plot the result
figure; plot(y); grid minor;
xlabel('Incident ray angle [Deg]');
ylabel('Lens surface tangent angle');

Answer 1

With Matlab, I obtained the plot below with the following simplified loop.

for i = 1:size(a0,2)
  fun3 = @(y) sind(90-y) - n*(sind(90-asind(sind(a0(i))/n)-y));
  y(i) = fzero(fun3, [0,90]);
end

The difference is in the form of equation: I replaced 90-y = asind(something) with sin(90-y) = something. When "something" is greater than 1 in absolute value, the former version throws an error due to complex value of asind. The latter proceeds normally, recognizing that this is not a solution (sin(90-y) can't be equal to something that is greater than 1).

No precomputing of the bracket was necessary, [0,90] simply worked. Another change I made was in the plot: plot(a0,y) instead of plot(y), to get the correct horizontal axis.

And you can't avoid for loop here, nor should you worry about it. Vectorization means eliminating loops where the content is a low-level operation that can be done en masse by operating on some C array. But fzero is totally not that. If the code takes long to run, it's because solving a bunch of equations takes long, not because there's a for loop.