exponenta event banner

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

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

F (1) = 3x1-2x12-2x2 + 1, F (i) = 3xi-2xi2-xi-1-2xi + 1 + 1, F (n) = 3xn-2xn2-xn-1 + 1.

Решаемые уравнения - Fi (x) = 0, 1≤i≤n. В примере используется n = 1000 . nlsf1a вспомогательная функция в конце этого примера реализует целевую функцию F (x).

В примере Большая разреженная система нелинейных уравнений с якобианом, которая решает ту же систему, целевая функция имеет явный якобиан. Однако иногда нельзя вычислить якобиан явно, но можно определить, какие элементы якобиана ненулевые. В этом примере якобиан тридиагональен, потому что единственными переменными, которые появляются в определении F (i), являются xi-1, xi и xi + 1. Так что единственные ненулевые части якобиана находятся на главной диагонали и двух соседних с ней диагоналях. Узор якобийской разреженности представляет собой матрицу, ненулевые элементы которой соответствуют (потенциально) ненулевым элементам в якобиане.

Создайте разреженную n-на-n тридиагональную матрицу из единиц, представляющих узор якобийской разреженности. (Пояснение к этому коду см. в разделе spdiags.)

n = 1000;
e = ones(n,1);
Jstr = spdiags([e e e],-1:1,n,n);

Просмотр структуры Jstr.

spy(Jstr)

Figure contains an axes. The axes contains an object of type line.

Набор fsolve параметры для использования 'trust-region' алгоритм, который является единственным fsolve алгоритм, который может использовать шаблон якобианской разреженности.

options = optimoptions('fsolve','Algorithm','trust-region','JacobPattern',Jstr);

Установка начальной точки x0 к вектору из -1 записей.

x0 = -1*ones(n,1);

Решите проблему.

[x,fval,exitflag,output] = fsolve(@nlsf1a,x0,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.

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

disp(norm(fval))
   1.0522e-09
disp(output.funcCount)
    25

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

options = optimoptions('fsolve','Algorithm','trust-region');
[x,fval,exitflag,output] = fsolve(@nlsf1a,x0,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(fval))
   1.0659e-09
disp(output.funcCount)
        5005

Решатель достигает по существу того же решения, что и раньше, но для этого требуются тысячи оценок функций.

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

function F = nlsf1a(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;
end

См. также

Связанные темы