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

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

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

Похожие темы