Нелинейные уравнения с якобианом

Рассмотрите проблему нахождения решения системы нелинейных уравнений, якобиан которых разрежен. Размерность проблемы в этом примере 1000. Цель состоит в том, чтобы найти x таким образом что F (x) = 0. Принимая n = 1000, нелинейные уравнения

F(1)=3x12x122x2+1,F(i)=3xi2xi2xi12xi+1+1,F(n)=3xn2xn2xn1+1.

Чтобы решить большую нелинейную систему уравнений, F (x) = 0, можно использовать доверительную область отражающий алгоритм, доступный в fsolve, крупномасштабный алгоритм (Крупномасштабный по сравнению с Алгоритмами Средней шкалы).

Шаг 1: Запишите файл nlsf1.m, который вычисляет значения целевой функции и якобиан.

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

Шаг 2: Вызовите решить стандартную программу для системы уравнений.

xstart = -ones(1000,1);
fun = @nlsf1; 
options = optimoptions(@fsolve,'Display','iter',...
    'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'PrecondBandWidth',0);
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

Отправная точка дана, а также имя функции. Метод по умолчанию для fsolve является доверительным резким искривлением области, таким образом, необходимо задать 'Algorithm' как 'trust-region' в аргументе options в порядке выбрать алгоритм доверительной области. Установка опции Display к 'iter' заставляет fsolve отображать вывод в каждой итерации. Установка 'SpecifyObjectiveGradient' к true, fsolve причин, чтобы использовать якобиевскую информацию, доступную в nlsf1.m.

Команды отображают этот вывод:

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          1            1011                            19
     1          2         16.1942        7.91898           2.35      
     2          3       0.0228027        1.33142          0.291      
     3          4     0.000103359      0.0433329         0.0201      
     4          5      7.3792e-07      0.0022606       0.000946      
     5          6     4.02299e-10    0.000268381       4.12e-05

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.

Линейная система (приблизительно) решена в каждой главной итерации с помощью предобусловленного метода сопряженных градиентов. Установка PrecondBandWidth к 0 в options означает, что диагональный предварительный формирователь используется. (PrecondBandWidth задает пропускную способность матрицы перед созданием условий. Пропускная способность 0 средних значений там является только одной диагональю в матрице.)

От значений оптимальности первого порядка быстро происходит линейная сходимость. Количество итераций метода сопряженных градиентов (CG), требуемых на главную итерацию, является низким, самое большее пять для проблемы 1 000 размерностей, подразумевая, что линейные системы не очень трудно решить в этом случае (хотя больше, работа требуется как сходимость, прогрессирует).

Если вы хотите использовать трехдиагональный предварительный формирователь, т.е. матрицу перед созданием условий с тремя диагоналями (или пропускная способность одной), установите PrecondBandWidth на значение 1:

options = optimoptions(@fsolve,'Display','iter','SpecifyObjectiveGradient',true,...
   'Algorithm','trust-region','PrecondBandWidth',1);
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

В этом случае вывод

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          1            1011                            19
     1          2         16.0839        7.92496           1.92      
     2          3       0.0458181         1.3279          0.579      
     3          4     0.000101184      0.0631898         0.0203      
     4          5     3.16615e-07     0.00273698        0.00079      
     5          6     9.72481e-10     0.00018111       5.82e-05

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.

Обратите внимание на то, что несмотря на то, что то же количество итераций происходит, количество итераций PCG понизилось, так меньше, работа делается на итерацию. Смотрите Предобусловленный Метод сопряженных градиентов.

Установка PrecondBandWidth к Inf (это - значение по умолчанию) означает, что решатель использует факторизацию Холесского, а не PCG.