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

Question

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

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

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

Answer 1

С 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-у = asind (то) с грехом (90-у) = что-то. Когда «что-то» больше 1 по абсолютной величине, первая версия выдает ошибку из-за сложного значения asind. Последнее протекает нормально, признавая, что это не решение (sin (90-y) не может быть равно тому, что больше 1).

Не требуется предварительный расчет кронштейна, [0,90] просто работал. Другое изменение, которое я сделал, было в сюжете: plot(a0,y) вместо plot(y), чтобы получить правильную горизонтальную ось.

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