Основная идея состоит в том, чтобы использовать итерации Ньютона Гаусса, чтобы решить нелинейные уравнения. Скажите, что вы пытаетесь решить уравнение

В установке FEM вы решаете слабую форму r (u) = 0. Установите, как обычно,

где x представляет 2D или 3-D точку. Затем умножьте уравнение на произвольную тестовую функцию _ϕi, объединяйтесь на области Ω и используйте формулу Грина и граничные условия, чтобы получить

который должен содержать для всех индексов i.

Вектор невязок ρ (U) может быть легко вычислен как

где матрицы K, M, Q и векторы F и G производятся путем сборки проблемы

Примите, что у вас есть предположение U ⁽ⁿ⁾ решения. Если U ⁽ⁿ⁾ достаточно близок к точному решению, улучшенное приближение, U ^{(n +1)} получен путем решения линеаризовавшей задачи

где $$\alpha$$ положительное число. (Не необходимо, чтобы ρ (U) = 0 имел решение, даже если ρ (u) = 0 имеет.) В этом случае итерация Ньютона Гаусса имеет тенденцию быть минимизатором невязки, т.е. решением min_U $$\Vert \rho (U)\Vert$$.

Это известно это за достаточно маленький $$\alpha$$

и

называется направлением спуска для $$\Vert \rho (U)\Vert$$, где $$\Vert \cdot \Vert$$ _2-норма L. Итерация

где $$\alpha$$ ≤ 1 выбран как можно больше таким образом, что шаг имеет разумный спуск.

Метод Ньютона Гаусса локален, и сходимость гарантируют только, когда U ⁽⁰⁾ достаточно близок к решению. В общем случае первое предположение может найтись вне области сходимости. Чтобы улучшить сходимость от плохих исходных предположений, стратегия затухания реализована для выбора α, поиск линии Армиджо-Голдстайна. Это выбирает самый большой коэффициент демпфирования α из последовательности 1, 1/2, 1/4... таким образом, что следующее неравенство содержит:

который гарантирует сокращение нормы невязки по крайней мере 1 – $$\alpha$$/2. Каждый шаг алгоритма поиска линии требует оценки невязки $$\rho \left(U^{(n)}+\alpha p_n\right)$$.

Важный момент этой стратегии - это, когда U ⁽ⁿ⁾ приближается к решению, затем $$\alpha$$→1 и таким образом повышения ставки сходимости. Если существует решение ρ (U) = 0, схема в конечном счете восстанавливает квадратичный уровень сходимости стандарта итерация Ньютона.

Тесно связанный с предыдущей проблемой выбор исходного предположения U ⁽⁰⁾. По умолчанию решатель устанавливает U ⁽⁰⁾ и затем собирает матрицы FEM K и F и вычисляет

Ослабленная итерация Ньютона Гаусса затем начата с U ⁽¹⁾, который должен быть лучшим предположением, чем U ⁽⁰⁾. Если граничные условия не зависят от решения u, то U ⁽¹⁾ удовлетворяет им, даже если U ⁽⁰⁾ не делает. Кроме того, если уравнение линейно, то U ⁽¹⁾ является точным решением FEM, и решатель не вводит цикл Ньютона Гаусса.

Существуют ситуации, где U ⁽⁰⁾ = 0 не имеет никакого смысла, или сходимость невозможна.

В некоторых ситуациях у вас может уже быть хорошее приближение, и нелинейный решатель может быть начат с него, избежав медленного режима сходимости. Эта идея используется в адаптивном генераторе mesh. Это вычисляет решение $$\tilde{U}$$ на mesh, оценивает ошибку и может совершенствовать определенные треугольники. interpolant $$\tilde{U}$$ очень хорошее начальное приближение для решения на усовершенствованной mesh.

В целом точный якобиан

не доступно. Приближение _Jn конечными разностями следующим образом является дорогим, но выполнимым. i th столбец _Jn может быть аппроксимирован

который подразумевает сборку матриц FEM для треугольников, содержащих узел решетки i. Очень простое приближение к _Jn, который дает итерацию фиксированной точки, также возможно следующие. По существу, для данного U ⁽ⁿ⁾, вычислите матрицы FEM K и F и установите

Это эквивалентно аппроксимации якобиана с матрицей жесткости. Действительно, начиная с ρ (U ⁽ⁿ⁾) = KU ⁽ⁿ⁾ – F, помещая _Jn = выражения K

Во многих случаях уровень сходимости является медленным, но стоимость каждой итерации является дешевой.

Partial Differential Equation Toolbox™ нелинейный решатель также предусматривает компромисс между этими двумя экстремальными значениями. Чтобы вычислить производную отображения U →KU, продолжите можно следующим образом. Термин a был не использован для ясности, но появляется снова в конечном результате.

Первый интегральный срок является не чем иным как _Ki,j.

Второй срок “смешан”, т.е. заменен диагональной матрицей, которая содержит суммы строки. С тех пор _Σjϕj = 1, второй срок аппроксимирован

который является i th компонент K ^(c') U, где K ^(c') является матрицей жесткости, сопоставленной с коэффициентом ∂c / ∂ u, а не c. То же обоснование может быть применено к производной отображения U →MU. Производная отображения U → –F точно

который является большой матрицей, сопоставленной с коэффициентом ∂f / ∂ u. Таким образом якобиан остаточного ρ (U) аппроксимирован

где дифференцирование относительно u, K и M определяют жесткость и большие матрицы, и их индексы определяют коэффициенты, относительно которых они собраны. В каждой итерации Ньютона Гаусса нелинейный решатель собирает матрицы, соответствующие уравнениям

и затем производит аппроксимированный якобиан. Дифференцирования коэффициентов сделаны численно.

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

где коэффициенты $$\tilde{K}$$ и $$\tilde{F}$$ может зависеть от решения $$\tilde{U}$$. “Смешанный” подход аппроксимирует производное отображение невязки

Нелинейность граничных условий и зависимости коэффициентов на производных $$\tilde{U}$$ правильно не линеаризуются этой схемой. Когда такая нелинейность сильна, схема уменьшает до итерации фиксировать-точки и может сходиться медленно или нисколько. Когда граничные условия линейны, они не влияют на свойства сходимости схем итерации. В Неймановом случае они невидимы (H является пустой матрицей), и в случае Дирихле они просто утверждают, что невязка является нулем на соответствующих граничных точках.