1626435587-55f52a4de97976f3c6215fa7c103f544 (Смирнов 2017 - Основы вычислительной физики ч2), страница 14

Схематически данную связь удобно изображать с помощьючисленной схемы, показанного на рис. 11 (б). Точки на рисунке соответствуютчетырём узлам сетки, связь значений численного решения в которыхопределяется соотношением (91).Обратим внимание на чрезвычайную простоту явной схемы (89).Действительно, чтобы получить схему (89), нам потребовалось лишьзаменить частные производные в уравнении (86) соответствующимиконечными разностями. Сделав это вполне очевидным способом, мыполучили формулу (91), позволяющую нам совместно с граничнымиусловиями вычислять искомые значения численного решения слой заслоем по времени.Однако использование явной схемы на практике для построениячисленного решения оказывается не столь простым, как вывод и программирование формул (90), (91).

Для понимания основной проблемы,связанной с использованием явной схемы, полезно посмотреть на систему разностных уравнений (89) как на схему Эйлера для системы ОДУпервого порядка на ( − 1) функцию времени 1 (), . . . , −1 (). Ранеемы уже видели [1, с. 101], что численное решение, полученное с использованием явных схем и недостаточно малого шага сетки, оказываетсянеустойчивым.

В этой связи исследуем вопрос устойчивости схемы (89).Для простоты рассмотрим задачу Дирихле (0, ) = (1, ) = 0 для однородного уравнения ( (, ) ≡ 0). Общее решение такой задачи можетбыть получено методом разделения переменных:шаблона(, ) =∞∑︁ sin() exp(− 2 2 ),=1где — коэффициенты Фурье начального распределения температурыпри = 0. Каждый член ряда представляет собой гармонику Фурье,равную нулю на границах отрезка [0, 1], с амплитудой, экспоненциальноубывающей во времени.

Посмотрим, сохранится ли это свойство базисных решений при переходе к разностной схеме (89). Для этого рассмотрим вначале действие разностного оператора+1− 2 + −1(92)ℎ2на сеточную функцию = sin( ) = sin(ℎ), где —амплитуда гармоники с номером на -м шаге по времени, ℎ — шаг70ˆ =сетки по : (︀)︀ˆ = sin(ℎ + ℎ) − 2 sin(ℎ) + sin(ℎ − ℎ) .ℎ2Раскрывая синус суммы и разности, приводя подобные члены и используя тождество cos(ℎ) − 1 ≡ −2 sin2 (ℎ/2), имеем:)︂)︂(︂(︂ˆ = − 4 sin2 ℎ sin(ℎ) = − 4 sin2 ℎ · .(93)ℎ22ℎ22Таким образом, мы получили, что = sin(ℎ) является собственной функцией разностного оператора (92). Подставляя полученˆ (89), получаемный результат в численную схему ( +1 − )/ = выражение для амплитуды синуса на ( + 1)-м слое:(︂(︂)︂)︂4ℎ+1 = · 1 − 2 sin2.(94)ℎ2Поскольку квадрат синуса не превосходит единицы, из полученного соотношения следует, что амплитуда численного решения будет затухать в геометрической прогрессии (экспоненциально по ) независимоот номера гармоники при условии<ℎ2.2(95)условно устойчивойСледовательно, явная схема (89) является.

В случае, если условие (95) не выполнено, происходит развитие неустойчивости в численном решении. Как видно на рис. 12 (а), вначале численноерешение может хорошо согласовываться с точным, однако начиная снекоторого момента в численном решении становятся заметны высокочастотные осцилляции (см. рис. 12 (б)), амплитуда которых возрастает экспоненциально по времени (в геометрической прогрессии сувеличением числа шагов). В соответствии с (94), быстрее всего возрастает амплитуда гармоники Фурье с максимальной частотой, близкой к частоте Найквиста: ∼ sin( −1 ).

Это кардинально противоречит физическим ожиданиям, в соответствии с которыми высокочастотные компоненты решения должны затухать быстрее всего. Дальнейшийрост неустойчивости приводит к появлению очень больших значенийчисленного решения и скорому выходу за пределы разрядной сетки —в выводе программы появляются значения inf и nan.Почему при нарушении соотношения (95) численное решение, построенное по явной схеме (89), оказывается неустойчивым и не аппроксимирует точное решение уравнения (86) даже в пределе , ℎ → 0?71(б)xtРешениеv(а)0.20.1000.51xРис.

12.(а) Динамика развития и (б) проявление неустойчивостиЧтобы ответить на этот вопрос, сравним распространение возмущенийв уравнении теплопроводности и численной схеме (89). В непрерывномслучае распространение возмущений описывается функцией Грина(︂)︂( − ′ )2( − ′ )′′( − , − ) = √︀,(96)exp −4( − ′ )4( − ′ )где — функция Хевисайда.

Формально () > 0 ∀ при > ′ , чтоозначает наличие дальнодействия в уравнении теплопроводности (86).Однако, хотя функция Грина нигде не обращается в ноль при > ′ ,она быстро убывает с ростом | − ′ |, что с физической точки зрения ограничиваетобласть влияния возмущения масштабом расстоя√ний ℎ* ≈ 2 . Тем не менее, даже несмотря на конечный размер области возмущений ℎ* , скорость распространения возмущений остаётсянеограниченной: ℎ* / → ∞ при → 0.

Принципиально иначе происходит распространение возмущений в численном решении, построенном по явной схеме (89). Если в некоторой точке ( , ) произошловозмущение (например, выделилось тепло), то в соответствии с (89)через время это приведёт к изменению температуры всего в трёх узлах сетки: −1 , и +1 . На каждом шаге область возмущения(показана крестиками на рис. 13 (а)) расширяется ровно на один шагсетки ℎ влево и вправо, что соответствует постоянной и конечной скорости распространения возмущений 2ℎ/ . Следовательно, чтобы схема(89) корректно воспроизводила точное решение с неограниченно высокой скоростью распространения возмущений, необходимо положить = (ℎ) = (ℎ2 ), что согласуется с условием устойчивости (95) явнойсхемы (89).Наконец, исследуем вопрос порядка аппроксимации численной схемы (89).

Запишем невязку точного решения на численной схеме (89):=+1− +1 − 2 + −1−− .ℎ272(а)tt0Рис. 13.(б)0xxРаспространение возмущения в численном и точном решенииРаскладывая (, ) в ряд Тейлора в точке ( , ), получаем:= + 2 /2 + . . . 2 ℎ2 /2 + 2 ℎ4 /24 + . . .−− ,ℎ2где нижними индексами обозначены производные по и . В силу (86) − − = 0, откуда = ℎ2− + . . . = ( + ℎ2 ).212(97)Следовательно, явная схема (89) имеет первый порядок точности повремени и второй — по пространственной координате .Таким образом, платой за простоту явной схемы является её низкаяэффективность, связанная с условной устойчивостью и первым порядком аппроксимации по . Как следствие, использование явной схемытребует большого количества шагов для построения численного решения.

В следующих двух пунктах мы покажем, как можно исправитьуказанные недостатки.3.3. Неявная схемаЕсли заменить частные производные и в уравнении (86) конечными разностями не в точке ( , ), как это было сделано в п. 3.2,но в ( , +1 ), получим схему+1+1+1− 2+1 + −1+1 − =+ +1 .ℎ2(98)Разностные соотношения (98) при = 1, . . . , − 1 образуют систему+1уравнений на ( − 1) неизвестную величину 1+1 , . . . , −1 .

Обратимвнимание, что количество уравнений и неизвестных в системе на два+1меньше количества узлов сетки: входящие в (98) величины 0+1 и 73определяются условиями на границе. Значения температуры на следующем, ( + 1)-м, слое по времени получаются в результате решениясистемы уравнений и уже не могут быть выписаны явно, как это было сделано в предыдущем параграфе (см. формулу (91)). В этой связисхему (98) относят к классучисленных схем.В каждое уравнение системы (98) входят три неизвестных величины+1(+1 , ±1) на слое = +1 и известное значение , что соответствует шаблону на рис. 11 (в), с. 68.Система уравнений на величины +1 ( = 1, .

. . , − 1) имееттрёхдиагональную матрицу и может быть решена методом прогонки,см. п. 1.5, с. 15. Коэффициенты трёхдиагональной матрицы (18) прирешении задачи Дирихле с нулевыми условиями на границе имеют вид:{︂0, = 1; =− /ℎ2 , = 2, . . . , − 1;2 = 1 + 2 ;{︂ ℎ 2− /ℎ , = 1, . . . , − 2; =0, = − 1;неявных= + +1 .В случае решения задачи Неймана с условием = 0 на грани+1цах значения 0+1 и следует подставлять в (98) из соотношений(87,88), что приведёт к изменению матричных коэффициентов системыв первом и последнем уравнениях:⎧ = 1,⎨ 0,− /ℎ2 , = 2, .

. . , − 2, =⎩ 2− 3 /ℎ2 , = − 1;{︂21 + 2 /ℎ , = 2, . . . , − 2, =1 + 32 /ℎ2 , = 1 и = − 1;⎧ 22 = 1,⎨ − 3 /ℎ ,− /ℎ2 , = 2, . . . , − 2, =⎩0, = − 1;= + +1 .Заметим, что в обоих рассмотренных случаях матричные коэффициенты удовлетворяют строгому условию диагонального преобладания(| | ≥ | | + | |), что гарантирует устойчивость решения системы линейных уравнений: малое изменение распределения температуры 74при = будет приводить к малому изменению решения +1 наследующем слое при = +1 .Решение задачи с периодическими граничными условиями (0) == (1), (0) = (1) потребует применения соответствующей модификации метода прогонки, рассмотренной в п. 1.6 на с. 16.Исследуем устойчивость неявной схемы (98).