2014-11-07 4 views
0

Предположим, что я решаю систему нелинейных уравнений. Простым примером может служить:Определение скорости сходимости fsolve - Matlab

function example 
    x0 = [15; -2]; 
    options = optimoptions('fsolve','Display','iter','TolFun',eps,'TolX',eps); 
    [x,fval,exitflag,output] = fsolve(@P1a,x0,options); 
end 

function f1 = P1a(x) 
    f1 = [x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29]; 
end 

Как определить скорость конвергенции? 'Display', 'iter' показывает мне норму на каждом шагу, но я не могу найти способ извлечь эти значения. (Для этого конкретного примера, я считаю, что fsolve не сходится к правильному решению, а скорее к локальному минимуму. Однако это не проблема. Я просто хочу найти способ оценить скорость конвергенции.)

+0

Может 'output.iterations' поможет? – David

+0

Я использовал это, чтобы сообщить мне количество итераций. Я хочу знать, могу ли я найти порядок конвергенции? –

+0

Не могли бы вы объединить количество итераций и ошибку на каждой итерации, чтобы дать вам представление о конвергенции? Я не уверен, что еще больше вы сможете выйти из 'fsolve'. – David

ответ

0

Вы можете получить много из fsolve. Однако вам нужно будет выполнить определенную работу. Прочитайте опцию 'OutputFcn' и напишите output functions для методов оптимизации Matlab. Это отличается от варианта с тем же именем, что и для решателей ODE Matlab. Вот пример, который воспроизводит 'Display','iter' вариант стоимости для fsolve (для алгоритма по умолчанию 'trust-region-dogleg' специально):

function stop = outfun(x,optimValues,state) 
% See private/trustnleqn 
stop = false; 

switch state 
    case 'init' 
     header = sprintf(['\n           Norm of  First-order Trust-region\n',... 
          ' Iteration Func-count  f(x)   step   optimality radius']); 
     disp(header); 
    case 'iter' 
     iter = optimValues.iteration;    % Iteration 
     numFevals = optimValues.funccount;   % Func-count 
     F = optimValues.fval;      % f(x) 
     normd = optimValues.stepsize;    % Norm of step 
     normgradinf = optimValues.firstorderopt; % First-order optimality 
     Delta = optimValues.trustregionradius;  % Trust-region radius 

     if iter > 0 
      formatstr = ' %5.0f  %5.0f %13.6g %13.6g %12.3g %12.3g'; 
      iterOutput = sprintf(formatstr,iter,numFevals,F'*F,normd,normgradinf,Delta); 
     else 
      formatstr0 = ' %5.0f  %5.0f %13.6g     %12.3g %12.3g'; 
      iterOutput = sprintf(formatstr0,iter,numFevals,F'*F,normgradinf,Delta); 
     end 
     disp(iterOutput); 
    case 'done' 

    otherwise 

end 

Вы можете назвать это с помощью:

function example 
[email protected](x)[x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29]; 
x0 = [15; -2]; 
opts = optimoptions('fsolve','Display','off','OutputFcn',@outfun,'TolFun',eps,'TolX',eps); 
[x,fval,exitflag,output] = fsolve(P1a,x0,opts); 

Это еще раз печатает в окне команд , Отсюда необходимо создать функцию вывода, которая может записывать данные в массив, файл или другую структуру данных. Вот как вы можете сделать это с global variable (в общем, не очень хорошая идея):

function stop = outfun2(x,optimValues,state) 
stop = false; 
global out; % Global variable, define in main function too 

switch state 
    case 'init' 
     out = []; 
    case 'iter' 
     iter = optimValues.iteration;    % Iteration 
     numFevals = optimValues.funccount;   % Func-count 
     F = optimValues.fval;      % f(x) 
     normd = optimValues.stepsize;    % Norm of step 
     normgradinf = optimValues.firstorderopt; % First-order optimality 
     Delta = optimValues.trustregionradius;  % Trust-region radius 

     out = [out;iter numFevals F'*F normd normgradinf Delta]; 
    case 'done' 

    otherwise 

end 

Тогда просто объявить global out; в своей основной функции перед вызовом fsolve. Вы также можете выполнить это, сделав свою выходную функцию nested function, и в этом случае массив out будет совместно использоваться внешней основной функцией.

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

Я уеду «определение скорости сходимости» к вам теперь, что у вас есть инструменты в руках ...

+0

Большое вам спасибо за это! Это действительно обширно, и я ценю вашу помощь. –

Смежные вопросы