Большая разреженная система нелинейных уравнений с якобианом

Скрипт Open Live Script

В этом примере показано, как использовать функции fsolve решатель, чтобы решить большие разреженные системы уравнений эффективно. Пример использует целевую функцию, заданную для системы $n$ уравнения,

$\begin{array}{l} F (1) = 3 x_{1} - 2 x_{1}^{2} - 2 x_{2} + 1, \\ F (i) = 3 x_{i} - 2 x_{i}^{2} - x_{i - 1} - 2 x_{i + 1} + 1, \\ F (n) = 3 x_{n} - 2 x_{n}^{2} - x_{n - 1} + 1 . \end{array}$

Уравнения, чтобы решить $F_{i} (x) = 0$ , $1 \leq i \leq n$ . Использование в качестве примера $n = 1000$ .

Эта целевая функция достаточно проста, что можно вычислить ее якобиан аналитически. Как объяснено в записи Векторных и Матричных Целевых функций, якобиана $J (x)$ из системы уравнений $F (x)$ $J_{i j} (x) = \frac{\partial F_{i} (x)}{\partial x_{j}}$ . Обеспечьте эту производную как второй выход вашей целевой функции. nlsf1 функция помощника в конце этого примера создает целевую функцию $F (x)$ и его якобиан $J (x)$ .

Решите систему уравнений с помощью опций по умолчанию, которые вызывают 'trust-region-dogleg' алгоритм. Запустите с точки xstart(i) = -1.

n = 1000;
xstart = -ones(n,1);
fun = @nlsf1; 
[x,fval,exitflag,output] = fsolve(fun,xstart);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Отобразите качество решения и количество взятых вычислений функции.

disp(norm(fval))

   2.8577e-13

disp(output.funcCount)

fsolve решает уравнение точно, но берет тысячи вычислений функции, чтобы сделать так.

Решите уравнение с помощью якобиана и в значении по умолчанию и в 'trust-region' алгоритмы.

options = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,fval2,exitflag2,output2] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

options.Algorithm = 'trust-region';
[x3,fval3,exitflag3,output3] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

disp([norm(fval2),norm(fval3)])

   1.0e-08 *

    0.0000    0.1065

disp([output2.funcCount,output3.funcCount])

     7     5

Оба алгоритма берут крошечную часть количества вычислений функции по сравнению с номером, не используя якобиан. Алгоритм по умолчанию берет еще несколько вычислений функции, чем 'trust-region' алгоритм, но алгоритм по умолчанию достигает более точного ответа.

Смотрите ли установка 'PrecondBandWidth' опция к 1 изменению 'trust-region' ответ или КПД. Эта установка вызывает fsolve использовать трехдиагональный предварительный формирователь, который должен быть эффективным для этой трехдиагональной системы уравнений.

options.PrecondBandWidth = 1;
[x4,fval4,exitflag4,output4] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

disp(norm(fval4))

   3.1185e-05

disp(output4.funcCount)

disp(output4.cgiterations)

'PrecondBandWidth' опция, устанавливающая причины fsolve дать немного менее точный ответ, как измерено нормой невязки. Количество вычислений функции увеличивается немного, от 5 до 6. Решатель имеет меньше чем десять итераций метода сопряженных градиентов как часть его процесса решения.

Смотрите как хорошо fsolve выполняет с диагональным предварительным формирователем.

options.PrecondBandWidth = 0;
[x5,fval5,exitflag5,output5] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

disp(norm(fval5))

   2.0057e-05

disp(output5.funcCount)

disp(output5.cgiterations)

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

Решите уравнения с помощью 'levenberg-marquardt' алгоритм.

options = optimoptions('fsolve','SpecifyObjectiveGradient',true,'Algorithm','levenberg-marquardt');
[x6,fval6,exitflag6,output6] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

disp(norm(fval6))

   7.4905e-15

disp(output6.funcCount)

Этот алгоритм дает самую низкую норму невязки и использует только еще один, чем самое низкое количество вычислений функции.

Обобщите результаты.

t = table([norm(fval);norm(fval2);norm(fval3);norm(fval4);norm(fval5);norm(fval6)],...
    [output.funcCount;output2.funcCount;output3.funcCount;output4.funcCount;output5.funcCount;output6.funcCount],...
    'VariableNames',["Residual" "Fevals"],...
    'RowNames',["Default" "Default+Jacobian" "Trust-Region+Jacobian" "Trust-Region+Jacobian,BW=1" "Trust-Region+Jacobian,BW=0" "Levenberg-Marquardt+Jacobian"])

t=6×2 table
                                     Residual     Fevals
                                    __________    ______

    Default                         2.8577e-13     7007 
    Default+Jacobian                2.5886e-13        7 
    Trust-Region+Jacobian           1.0646e-09        5 
    Trust-Region+Jacobian,BW=1      3.1185e-05        6 
    Trust-Region+Jacobian,BW=0      2.0057e-05        6 
    Levenberg-Marquardt+Jacobian    7.4905e-15        6

Этот код создает nlsf1 функция.

function [F,J] = nlsf1(x)
% Evaluate the vector function
n = length(x);
F = zeros(n,1);
i = 2:(n-1);
F(i) = (3-2*x(i)).*x(i)-x(i-1)-2*x(i+1) + 1;
F(n) = (3-2*x(n)).*x(n)-x(n-1) + 1;
F(1) = (3-2*x(1)).*x(1)-2*x(2) + 1;
% Evaluate the Jacobian if nargout > 1
if nargout > 1
   d = -4*x + 3*ones(n,1); D = sparse(1:n,1:n,d,n,n);
   c = -2*ones(n-1,1); C = sparse(1:n-1,2:n,c,n,n);
   e = -ones(n-1,1); E = sparse(2:n,1:n-1,e,n,n);
   J = C + D + E;
end
end

Смотрите также

fsolve

Документация

Большая разреженная система нелинейных уравнений с якобианом

Смотрите также

Похожие темы

Документация Optimization Toolbox

Поддержка