Недавно я сравнивал различные нелинейные решатели от scipy и был особенно впечатлен примером Ньютона-Крылова в книге рецептов Scipy, в которой они решают уравнение дифференциального уравнения второго порядка с нелинейным членом реакции примерно в 20 строках кода.
Я изменил пример кода для решения нелинейного уравнения Пуассона ( также называемого уравнением Пуассона-Больцмана , см. Стр. 17 в этих заметках) для полупроводниковых гетероструктур, которое имеет вид,
(Это остаточная функция, которая передается решателю.)
Это проблема электростатики, где и p ( x , ϕ ) являются нелинейными функциями для вида n i ( x ) e - ( E i ( x , ϕ ) - E f ) . Детали здесь не важны, но дело в том, что нелинейная функция изменяется экспоненциально с ϕ, поэтому функция невязки может изменяться в огромном диапазоне ( 10 - 6 - 10 16 )с небольшим изменением .
Я численно решаю это уравнение с помощью Ньютона-Крылова Сципи, но оно никогда не сойдется (фактически оно всегда сообщало бы об ошибке при вычислении якобиана). Я переключился с Солвера Ньютона-Крылова на fsolve (который основан на гибриде MINPACK), и он сработал впервые!
Существуют ли общие причины, по которым Ньютон-Крылов не подходит для определенных задач? Должны ли входные уравнения быть как-то обусловлены?
Может быть, нужно больше информации, чтобы прокомментировать, но как вы думаете, почему fsolve сработал в этом случае?
источник
sol = newton_krylov(func, guess, method='gmres')
) решило проблему. Не совсем уверен, почему, но кто-то еще с этой проблемой может подумать о том же.Ответы:
Есть две проблемы, с которыми вы, вероятно, столкнетесь.
Жестокое кондиционирование
Во-первых, проблема плохо обусловлена, но если вы предоставляете только остаток, Ньютон-Крылов выбрасывает половину ваших значащих цифр путем конечного дифференцирования остатка, чтобы получить действие якобиана:
Если вы предоставите аналитический якобиан, вы получите все цифры (например, 16 с двойной точностью). Методы Крылова также полагаются на внутренние продукты, поэтому, если ваш якобиан плохо подготовлен к мелодии1016 Это действительно единственное число, и Крылов может застаиваться или возвращать ошибочные решения. Это также может предотвратить схождение прямых решателей. Иногда вы можете использовать многосеточные методы для гомогенизации в грубую сетку с гибким кондиционированием. Когда проблема не может быть сформулирована с лучшей подготовкой, стоит поработать с четверной точностью. (Это поддерживается PETSc, например.)
Обратите внимание, что те же проблемы применимы к квазиньютоновским методам, но без конечных разностей. Все масштабируемые методы для задач с некомпактными операторами (например, дифференциальные уравнения) должны использовать информацию Якоби для предобусловливания.
Вполне вероятно, что
fsolve
либо не использовали якобианскую информацию, либо использовался метод изгиба ноги или сдвиг, чтобы добиться прогресса с помощью метода «градиентного спуска», несмотря на по существу единый якобиан (т. Е. Конечная разность имела бы много «шума» от арифметика конечной точности). Это не масштабируется и,fsolve
вероятно, становится медленнее по мере увеличения размера проблемы.глобализация
Если линейные задачи решены точно, мы можем исключить проблемы, относящиеся к линейной задаче (Крылов), и сосредоточиться на ней из-за нелинейности. Локальные минимумы и негладкие признаки медленной конвергенции или вызывают застой. Пуассон-Больцман является гладкой моделью, поэтому, если вы начнете с достаточно хорошего начального предположения, Ньютон сойдется квадратично. Большинство стратегий глобализации включают в себя своего рода продолжение, чтобы произвести высококачественное начальное предположение для заключительных итераций. Примеры включают продолжение сетки (например, Full Multigrid), продолжение параметра и псевдотранзитивное продолжение. Последнее обычно применимо к установившимся задачам и предлагает некоторую теорию глобальной конвергенции, см. Coffey, Kelley, and Keyes (2003) . Поиском оказывается эта статья, которая может быть вам полезна:Шестаков, Милович, Ной (2002) Решение нелинейного уравнения Пуассона-Больцмана с использованием псевдопереходного продолжения и метода конечных элементов . Псевдотранзитивное продолжение тесно связано с алгоритмом Левенберга-Марквардта.
дальнейшее чтение
источник