Как перебрать функции?

Question

Я хотел бы применить цикл над функцией. У меня есть следующий «мать» код:

v = 1; 
fun = @root; 
x0 = [0,0] 
options = optimset('MaxFunEvals',100000,'MaxIter', 10000); 
x = fsolve(fun,x0, options) 

Кроме того, у меня есть следующие функции в отдельном файле:

function D = root(x) 
v = 1; 
D(1) = x(1) + x(2) + v - 2; 
D(2) = x(1) - x(2) + v - 1.8; 
end 

Теперь я хотел бы найти корни, когда v = sort(rand(1,1000)). Другими словами, я хотел бы перебрать функцию для каждого значения v.

Answer 1

Вам нужно будет изменить root принять дополнительную переменную (v), а затем изменить функцию ручки укоренить в анонимную функцию, которая питает в v, что вы хотите

function D = root(x, v) 
    D(1) = x(1) + x(2) + v - 2; 
    D(2) = x(1) - x(2) + v - 1.8; 
end 

% Creates a function handle to root using a specific value of v 
fun = @(x)root(x, v(k)) 

Answer 2

vvec = sort(rand(1,2)); 

x0 = [0,0]; 
for v = vvec, 
    fun = @(x) root(v, x); 
    options = optimset('MaxFunEvals',100000,'MaxIter', 10000); 
    x = fsolve(fun, x0, options); 
end 

с функцией определение:

function D = root(v, x) 
D(1) = x(1) + x(2) + v - 2; 
D(2) = x(1) - x(2) + v - 1.8; 
end 

Answer 3

Только в том случае, уравнение ваш фактическое уравнения (а не фиктивный пример): что уравнение линейного, смысл, вы можете решить для всех v с простой mldivide:

v = sort(rand(1,1000)); 
x = [1 1; 1 -1] \ bsxfun(@plus, -v, [2; 1.8]) 

И, в случае, если те не ваших реальных уравнения , вам не нужно в цикле, вы можете векторизации все это:

function x = solver() 

    options = optimset('Display' , 'off',... 
         'MaxFunEvals', 1e5,... 
         'MaxIter' , 1e4); 

    v = sort(rand(1, 1000)); 
    x0 = repmat([0 0], numel(v), 1);  
    x = fsolve(@(x)root(x,v'), x0, options); 

end 

function D = root(x,v) 

    D = [x(:,1) + x(:,2) + v - 2 
     x(:,1) - x(:,2) + v - 1.8]; 

end 

Это может или не может быть быстрее, чем зацикливание, это зависит от ваших действительных уравнений.

Это может быть медленнее, потому что fsolve необходимо будет вычислить якобиан 2000 × 2000 (элементы 4M) вместо 2 × 2, 1000 раз (4k элементов).

Но это может быть быстрее, потому что стоимость запуска fsolve может быть большой, то есть накладные расходы многих вызовов могут фактически перевесить затраты на вычисление более крупного якобиана.

В любом случае, при условии, что якобиан в качестве второго выхода будет достаточно огромной скорости все вверх:

function solver() 

    options = optimset('Display' , 'off',... 
         'MaxFunEvals', 1e5,... 
         'MaxIter' , 1e4,...      
         'Jacobian' , 'on'); 

    v = sort(rand(1, 1000)); 
    x0 = repmat([1 1], numel(v), 1);  
    x = fsolve(@(x)root(x,v'), x0, options); 

end 

function [D, J] = root(x,v) 

    % Jacobian is constant: 
    persistent J_out 
    if isempty(J_out) 
     one = speye(numel(v));  
     J_out = [+one +one 
       +one -one]; 
    end 

    % Function values at x 
    D = [x(:,1) + x(:,2) + v - 2 
     x(:,1) - x(:,2) + v - 1.8]; 

    % Jacobian at x: 
    J = J_out; 

end