Анализ устойчивости фон Неймана

В численном анализе анализ устойчивости фон Неймана (также известный как анализ устойчивости Фурье) — это процедура, используемая для проверки устойчивости конечно -разностных схем применительно к линейным уравнениям в частных производных . Анализ основан на Фурье разложении числовой ошибки и был разработан в Лос-Аламосской национальной лаборатории после того, как он был кратко описан в статье 1947 года британскими исследователями Джоном Крэнком и Филлис Николсон . Этот метод является примером явного интегрирования по времени , когда функция, определяющая управляющее уравнение, вычисляется в текущий момент времени.Позже метод получил более строгое рассмотрение в статье. в соавторстве с Джоном фон Нейманом .

Численная стабильность

Устойчивость численных схем тесно связана с численной ошибкой . Конечно-разностная схема устойчива, если ошибки, допущенные на одном временном шаге расчета, не приводят к увеличению ошибок по мере продолжения вычислений. — Нейтрально устойчивая схема это схема, в которой ошибки остаются постоянными по мере продолжения вычислений. Если ошибки затухают и в конечном итоге затухают, численная схема называется устойчивой. Если, наоборот, ошибки растут со временем, то численную схему называют неустойчивой. Устойчивость численных схем можно исследовать путем выполнения анализа устойчивости фон Неймана. Для задач, зависящих от времени, стабильность гарантирует, что численный метод дает ограниченное решение всякий раз, когда решение точного дифференциального уравнения ограничено. В целом устойчивость может быть трудно исследовать, особенно если рассматриваемое уравнение является нелинейным .

В некоторых случаях устойчивость фон Неймана необходима и достаточна для устойчивости в смысле Лакса – Рихтмайера (как это используется в теореме эквивалентности Лакса ): модели PDE и конечно-разностной схемы линейны; УЧП имеет постоянный коэффициент с периодическими граничными условиями и имеет только две независимые переменные; и в схеме используется не более двух временных уровней. Устойчивость фон Неймана необходима в гораздо более широком круге случаев. Его часто используют вместо более детального анализа устойчивости, чтобы получить хорошее представление об ограничениях (если таковые имеются) на размеры шагов, используемых в схеме, из-за ее относительной простоты.

Иллюстрация метода

Метод фон Неймана основан на разложении ошибок в ряды Фурье . Чтобы проиллюстрировать процедуру, рассмотрим одномерное уравнение теплопроводности $${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}$$определяется на пространственном интервале $L$, с обозначением $$u_{j}^{n}=u(x_{j},t^{n})$$

Мы можем дискретизировать уравнение теплопроводности как

$$u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)$$ ( 1 )

где $$r={\frac {\alpha \,\Delta t}{\left(\Delta x\right)^{2}}}$$

Тогда решение $u_{j}^{n}$ дискретного уравнения аппроксимирует аналитическое решение $u(x,t)$ PDE на сетке.

Определить ошибку округления $\epsilon _{j}^{n}$ как $$\epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}$$где $u_{j}^{n}$ является решением дискретизированного уравнения ( 1 ), которое будет вычислено при отсутствии ошибки округления, и $N_{j}^{n}$ — численное решение, полученное в арифметике конечной точности . Поскольку точное решение $u_{j}^{n}$ должно точно удовлетворять дискретизированному уравнению, ошибка $\epsilon _{j}^{n}$ также должно удовлетворять дискретизированному уравнению. Здесь мы предполагали, что $N_{j}^{n}$ также удовлетворяет уравнению (это справедливо только для машинной точности). Таким образом

$$\epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)$$ ( 2 )

является рекуррентным соотношением для ошибки. Уравнения ( 1 ) и ( 2 ) показывают, что и ошибка, и численное решение имеют одинаковый характер роста или затухания во времени. Для линейных дифференциальных уравнений с периодическими граничными условиями пространственное изменение ошибки можно разложить в конечный ряд Фурье по $x$, в интервале $L$, как

$$\epsilon (x,t)=\sum _{m=-M}^{M}E_{m}(t)e^{{i}k_{m}x}$$ ( 3 )

где волновое число $k_{m}={\frac {\pi m}{L}}$ с $m=-M,\dots ,-2,-1,0,1,2,\dots ,M$ и $M=L/\Delta x$. Зависимость ошибки от времени учитывается в предположении, что амплитуда ошибки $E_{m}$ является функцией времени. Часто делается предположение, что ошибка растет или затухает со временем экспоненциально, но для анализа устойчивости это не является необходимым.

Если граничное условие не является периодическим, то мы можем использовать конечный интеграл Фурье по $x$:

$$\epsilon (x,t)=\int _{-{\frac {\pi }{\Delta x}}}^{\frac {\pi }{\Delta x}}E_{k}(t)e^{ikx}dk$$ ( 4 )

Поскольку разностное уравнение ошибки линейно (поведение каждого члена ряда такое же, как и самого ряда), достаточно рассмотреть рост ошибки типичного члена:

$$\epsilon _{m}(x,t)=E_{m}(t)e^{ik_{m}x}$$ ( 5а )

если используется ряд Фурье или

$$\epsilon _{k}(x,t)=E_{k}(t)e^{ikx}$$ ( 5б )

если используется интеграл Фурье.

Поскольку ряд Фурье можно рассматривать как частный случай интеграла Фурье, мы продолжим разработку, используя выражения для интеграла Фурье.

Характеристики устойчивости можно исследовать, используя именно эту форму погрешности без потери общности. Чтобы узнать, как ошибка меняется с течением времени, подставьте уравнение ( 5b ) в уравнение ( 2 ), заметив, что $${\begin{aligned}\epsilon _{j}^{n}&=E_{m}(t)e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=E_{m}(t+\Delta t)e^{ik_{m}x}\\\epsilon _{j+1}^{n}&=E_{m}(t)e^{ik_{m}(x+\Delta x)}\\\epsilon _{j-1}^{n}&=E_{m}(t)e^{ik_{m}(x-\Delta x)},\end{aligned}}$$давать (после упрощения)

$${\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}=1+r\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right).$$ ( 6 )

Представляем $\theta =k_{m}\Delta x\in [-\pi ,\pi ]$ и используя тождества $$\sin \left({\frac {\theta }{2}}\right)={\frac {e^{i\theta /2}-e^{-i\theta /2}}{2i}}\qquad \rightarrow \qquad \sin ^{2}\left({\frac {\theta }{2}}\right)=-{\frac {e^{i\theta }+e^{-i\theta }-2}{4}}$$уравнение ( 6 ) можно записать как

$${\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}=1-4r\sin ^{2}(\theta /2)$$ ( 7 )

Определить коэффициент усиления

$$G\equiv {\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}$$ ( 8 )

Необходимым и достаточным условием ограниченности ошибки является то, что $|G|\leq 1.$ Таким образом, из уравнений ( 7 ) и ( 8 ) условие устойчивости определяется выражением

$$\left|1-4r\sin ^{2}(\theta /2)\right|\leq 1$$ ( 9 )

Обратите внимание, что термин $4r\sin ^{2}(\theta /2)$ всегда положительный. Таким образом, чтобы удовлетворить Уравнению ( 9 ):

$$4r\sin ^{2}(\theta /2)\leq 2$$ ( 10 )

Чтобы вышеуказанное условие выполнялось для всех $m$ (и поэтому все $\sin ^{2}(\theta /2)$). Наивысшее значение, которое может принять синусоидальный член, равно 1, и для этого конкретного выбора, если условие верхнего порога удовлетворено, то то же самое будет и для всех точек сетки, таким образом, мы имеем

$$r={\frac {\alpha \Delta t}{\left(\Delta x\right)^{2}}}\leq {\frac {1}{2}}$$ ( 11 )

Уравнение ( 11 ) дает требование стабильности для схемы FTCS применительно к одномерному уравнению теплопроводности. Он говорит, что для данного $\Delta x$, допустимое значение $\Delta t$ должно быть достаточно малым, чтобы удовлетворять уравнению ( 10 ).

Подобный анализ показывает, что схема FTCS для линейной адвекции безусловно неустойчива.