В этом примере показано, как использовать функции fsolve
решатель, чтобы решить большие разреженные системы уравнений эффективно. Пример использует целевую функцию, заданную для системы уравнения,
Уравнения, чтобы решить , . Использование в качестве примера .
Эта целевая функция достаточно проста, что можно вычислить ее якобиан аналитически. Как объяснено в записи Векторных и Матричных Целевых функций, якобиана из системы уравнений . Обеспечьте эту производную как второй выход вашей целевой функции. nlsf1
функция помощника в конце этого примера создает целевую функцию и его якобиан .
Решите систему уравнений с помощью опций по умолчанию, которые вызывают '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)
7007
Уравнение решено точно, но 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)
6
disp(output4.cgiterations)
8
Эта опция, устанавливающая причины 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)
6
disp(output5.cgiterations)
19
Норма невязки немного ниже на этот раз, и количество вычислений функции неизменно. Количество итераций метода сопряженных градиентов увеличивается от 8 до 19, указывая, что эта установка заставляет решатель делать, больше работает.
Решите уравнения с помощью '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)
6
Этот алгоритм дает самую низкую норму невязки и использует только еще один, чем самое низкое количество вычислений функции.
Обобщите результаты.
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