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

Этот пример показывает, как обеспечить якобианский шаблон разреженности, чтобы эффективно решить большую разреженную систему уравнений. В примере используется целевая функция, заданная для системы 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, 1in. Пример использует n=1000. The 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

См. также

Похожие темы