Какие итерационные вспомогательные решатели, кроме BICG, я могу использовать для решения уравнения Ньютона в MATLAB?

Question

Я пытаюсь проанализировать результаты для различных итерационных подрешителей системы уравнений Ньютона, используя переформулировку Ньютона-Фишера LCP (линейная проблема дополнительности) . До сих пор я реализовал точный решатель - Гаусса-Зиделя и метод, модифицированный Ньютоном, который использует bicg matlab в качестве вспомогательного решателя уравнения J * h = -p (где J - это якобиан, p - это значение функции Фишера, а h - мой шаг актуализации).

Часть кода, в которой я реализовал субсольвер bicg и точный решатель:

if itt
    % use iterative solver for newton eq
    while ~all(fischer(x, A*x+b) == 0) & its < max_it
        % compute the Jacobian
        J = eval_jacobian(A, b, x);
        % compute the value of the Fischer function
        p = fischer(x, A*x + b);
        % the natural merit function for convergence measure
        residual(its) = .5*(p'*p);
        % the newton eq, solve J*h = -p
        h = bicg(J, -p, eps, s_its);
        % update the solution vector 
        x = x + h;
        % increment the iteration counter
        its = its + 1;
    end
else
    % the exact solver for newton equation
    while ~all(fischer(x, A*x+b) == 0) & its < max_it
        % compute the Jacobian
        J = eval_jacobian(A, b, x);
        % compute the value of the Fischer function
        p = fischer(x, A*x + b);
        % the natural merit function for convergence measure
        residual(its) = .5*(p'*p);
        % the newton eq, solve J*h = -p
        h = - J/p;
        % update the solution vector 
        x = x + h;
        % increment the iteration counter
        its = its + 1;
    end

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

Answer 1

В дополнение к bicg существует ряд других итерационных решателей для общих несимметричных систем, которые могут иногда обеспечивать лучшую сходимость, такие как bicgstab и gmres.Оба эти метода похожи по духу на bicg, но включают разные стратегии обновления в подпространстве Крылова.В MATLAB есть реализации для обоих этих методов, поэтому вы можете проверить их на предмет вашей проблемы.

В общем, при использовании итерационных решателей Крылова часто можно добиться лучшей производительности, включив предварительную обработку.Это довольно сложная тема - если вы еще не знакомы, здесь - это то же самое место, что и для начала.Я бы порекомендовал начать с простого предобусловливателя jacobi diag(J) и перейти оттуда.

Итеративные решатели также иногда могут получить дополнительную эффективность, основанную на так называемом подходе «без матрицы» - решающему на самом деле нужны толькочтобы вычислить матрично-векторное произведение J*p, ему не нужна полная матрица J.Может быть возможно написать эффективную подпрограмму, которая вычисляет J*p без явного формирования J, но это будет зависеть от деталей вашей конкретной проблемы.

Последнее соображение - это сама языковая ошибка программирования -MATLAB использует высокооптимизированный библиотечный код для непосредственного решения J\p, в то время как методы Крылова реализованы в виде m-кода, что обычно приводит к снижению производительности.

Надеюсь, это поможет.