Использование fzero в Matlab или Octave, избегая цикла for и сложных решений

Я использую функцию fzero для решения нелинейного уравнения, зависящего от одного параметра, и меня не устраивает мой метод. У меня есть следующие проблемы:

1) Можно ли избежать цикла for для параметра?

2) Чтобы избежать сложных решений, я сначала должен предварительно вычислить допустимый интервал для fzero. Есть ли лучшее решение здесь?

Если я уменьшу размер шага параметра, время выполнения станет медленным. Если я не вычисляю интервал заранее, я получаю сообщение об ошибке «Значения функции в конечных точках интервала должны быть конечными и действительными». в Matlab и «fzero: неверный начальный брекетинг» в Octave.

Вот код

% 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');

matlab octave equation-solving nonlinear-functions
person miquo    schedule 01.05.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

С Matlab я получил график ниже со следующим упрощенным циклом.

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

Разница в форме уравнения: я заменил 90-y = asind(что-то) на sin(90-y) = что-то. Когда "что-то" больше 1 по абсолютной величине, предыдущая версия выдает ошибку из-за комплексного значения asind. Последний работает нормально, признавая, что это не решение (sin(90-y) не может быть равен чему-то большему, чем 1).

Предварительное вычисление скобки не требовалось, [0,90] просто работало. Еще одно изменение, которое я сделал, было в графике: plot(a0,y) вместо plot(y), чтобы получить правильную горизонтальную ось.

И вы не можете избежать цикла for здесь, и вам не следует об этом беспокоиться. Векторизация означает устранение циклов, когда содержимое является низкоуровневой операцией, которую можно выполнить в массовом порядке, работая с некоторым массивом C. Но fzero совсем не то. Если код выполняется долго, это потому, что решение множества уравнений занимает много времени, а не потому, что есть цикл for.

сюжет

person Community    schedule 02.05.2016
comment
Спасибо, теперь код выглядит проще и понятнее. Ваше решение также работает в Octave 4.0, никаких ошибок или предупреждений. - person miquo; 03.05.2016