2.5. Явные и неявные схемы (Математическое моделирование задач газодинамики и тепло-массообмена. Молчанов А.М. 2013)

5. Явные и неявные схемы.5.1. Явная, Кранка-Николсона и полностью неявная схемы.Ограничения на шаг по времени, связанныес устойчивостью(см.формулы (4.31), (4.34) и т.п.), могут быть очень неудобны и приводить кбольшим затратам времени на численное решение задач.Введем понятие числа Куранта (иногда используется название: числоКуранта-Фридрихса-Леви):CFL =u ∆t∆x(5.1)По условию (4.34) для обеспечения устойчивости необходимо:CFL ≤ 1(5.2)Эти ограничения можно снять (частично или полностью), если перейти отявных разностных схем, которые рассматривались ранее, к неявным схемам.Неявные схемы используют уравнения, которые выражают данные на( n + 1) -ом шаге по времени на через несколько соседних точек сетки. Длянахождения результата решается система линейных уравнений.Рассмотримследующуюразностнуюсхему длярешениязадачитеплопроводности (4.25)Φ nm+1 − Φ nmΦ n +1 − 2Φ mn +1 + Φ mn +−11− ν m +1=0,∆t∆x 2(5.3)ν =Γ/ρ(5.4)где коэффициентдля простоты полагается постоянным (на анализ устойчивости, исходя изпринципа «замороженных» коэффициентов, это не влияет).Подставим решение в форме (4.18) - Φ nm = λ n eiα m в задачу (5.3) и найдемспектр λ (α ) :1λ=1−ν∆t∆x 2(eiα=−2+e− iα)14ν∆t 2 α1+sin∆x 22(5.5)Очевидно, что λ (α ) ≤ 1 при любом положительном значении ν , т.е.необходимое условие устойчивости выполняется при любом шаге повремени.Можно сформулировать следующее правило: неявные схемы чащебывают устойчивыми по сравнению с явными.

Но при этом усложняетсяпроцесс решения системы линейных алгебраических уравнений, которыеполучаются после дискретизации исходных уравнений.Кроме чисто явных и чисто неявных схем, бывают смешанные схемы.Например, для уравнения теплопроводности этоΦ nm+1 − Φ nmν− 2  f ( Φ nm++11 − 2Φ nm+1 + Φ nm+−11 ) + (1 − f ) ( Φ mn +1 − 2Φ mn + Φ mn −1 )  = 0 ,∆t∆x(5.6)где параметр f принимает значения от 0 до 1. При f = 0 получается чистоявная схема, при f = 1 - полностью неявная, при f = 0.5 получается такназываемая схема Кранка-Николсона.Легко показать, что при f ≥ 0.5 схема удовлетворяет необходимомуусловию устойчивость при любых соотношениях шагов по времени и попространственной координате.Такие схемы называются безусловно устойчивыми.Различные значения f можно интерпретировать как характеристикуизменения Φ m при переходе во времени от момента t n к t n +1 , показанного нарис.

4. Явная схема по существу предполагает, что старое значение Φ nmсуществует в пределах всего временного шага, за исключением точкиt n +1 = t n + ∆t . Неявная схема предполагает, что в моментtnΦmрезкоуменьшается от Φ nm до Φ nm+1 , а затем остается равной Φ nm+1 на всем временномшаге и функция в пределах временного шага характеризуется новымзначением Φ nm+1 .

Схема Кранка-Николсона предполагает линейное изменениеΦ m . На первый взгляд линейное изменение должно быть более разумным,чем две другие альтернативы.Рис.4. Изменение искомой функции Φ m температуры во времени дляявной схемы (1), схемы Кранка-Николсона (2) и полностью неявной схемы(3).Тем не менее, чисто неявная схема является предпочтительной. Причиныдве.1) В случае схемы Кранка-Николсона могут иметь место колеблющиесярешения.

Устойчивость в математическом смысле просто гарантирует, чтоэти колебания будут, в конечном счете, затухать, но это не обеспечиваетфизически правдоподобного решения. Несколько примеров подобныхрешений, полученных с помощью схемы Кранка-Николсона, можно найти в[1].2) Точные решения рассматриваемой задачи имеют экспоненциальныйхарактер, поэтому линейная интерполяция хуже отражает истинный характеризменения функции.УпражнениеДоказать, что схема Кранка-Николсона действительно может даватьнефизичные результаты, поскольку для нее не всегда выполняютсятребования, описанные в разделе 2.3.5.2. Метод прогонкиКак уже говорилось, в случае использовании неявной схемы получаетсясистема уравнений относительно неизвестных на ( n + 1) -ом шаге по временизначений функции Φ .

Например, уравнение (5.3) преобразуется к видуν∆t  n +1 ν∆t n +1 ν∆t n +1n1 + 2 2  Φ m = 2 Φ m +1 + 2 Φ m −1 + Φ m ,∆x ∆x∆xm = 1, 2,..., N(5.7)Для решения этой системы используется так называемый методпрогонки.Методпрогонкииспользуетсядлярешениясистемылинейныхалгебраических уравнений, в которой каждая m -ая строка представлена ввидеamΦ m = bm Φ m +1 + cmΦ m −1 + d m(5.8)где m = 1, 2,..., N − 1, N .Для случая (5.7) коэффициенты равныν∆t am =  1 + 2 2  ,∆x ν∆tν∆tbm = 2 , cm = 2 , d m = Φ nm∆x∆x(5.9)Для m = 1, m = N коэффициентыc1 = 0, bN = 0(5.10)Матрица такой системы называется трехдиагональной.

В английскойлитературе метод называется алгоритмом Томаса или TDMA (Tri-diagonalMatrix Algorithm- трехдиагональный матрицы алгоритмом).Если на границах области заданы значения искомой функции, то:a1 = 1, b1 = c1 = 0, d1 = Φ1 ,(5.11)aN = 1, bN = cN = 0, d N = Φ N(5.12)Записанные условия (5.8), (5.10) означают, что Φ1 известна в зависимостиот Φ 2 . Уравнение для m = 2 представляет собой соотношение междуΦ1 , Φ 2 , Φ 3 . Но поскольку Φ1 может быть выражена через Φ 2 , это соотношениеприводится к соотношению между Φ 2 и Φ 3 .

Другими словами, Φ 2 можновыразить через Φ 3 . Процесс подстановки можно продолжать до тех пор, показначение Φ N не будет выражено через Φ N +1 . Но поскольку Φ N +1 несуществует, мы в действительности на данном этапе получим численноезначение Φ N . Это позволит нам начать процесс обратной подстановки, вкотором Φ N −1 получится из Φ N , Φ N − 2 — из Φ N −1 , ..., Φ 2 — из Φ 3 и Φ1 — из Φ 2 .Это и составляет существо алгоритма трехдиагональной матрицы.Предположим, что при прямой подстановке имеем зависимостьΦ m = Pm Φ m +1 + Qm(5.13)Φ m −1 = Pm −1Φ m + Qm −1(5.14)после того, как полученоПодставляя (5.14) в (5.8), получаем следующее соотношение:Φm =bmc Q + dmΦ m +1 + m m −1( am − cm Pm −1 )( am − cm Pm −1 )Сравниваем это выражение с (5.13) и получаем(5.15)формулы дляпрогоночных коэффициентов Pm , Qm :Pm =bmc Q + dm, Qm = m m −1( am − cm Pm −1 )( am − cm Pm −1 )Эти рекуррентные соотношения определяют Pm , Qm(5.16)через Pm −1, Qm −1 .Заметим, что в начале рекуррентного процесса уравнение (5.8) для m = 1 поформе почти совпадает с (5.13).

Таким образом, P1 , Q1 определяются вследующем виде:P1 =b1d, Q1 = 1a1a1(5.17)На другом конце последовательности Pm , Qm имеем bN = 0 . Это дает PN = 0 ,и из (5.13)получаем(5.18)QN = Φ NС этого момента осуществляется обратная подстановка с помощьюуравнения (5.13).Краткое описание алгоритма1. Рассчитываем P1 , Q1 из уравнения (5.17).2.

Используя рекуррентные соотношения (5.16), получаем Pm , Qm дляm = 2,3,..., N3. Полагаем QN = Φ N4. Используя уравнение (5.13) дляm = N − 1, N − 2,...,3, 2,1 ,получаемтрехдиагональностиматрицы,Φ N −1 , Φ N − 2 ,..., Φ 3 , Φ 2 , Φ1 .Алгоритм,использующийсвойствоявляется мощным и удобным методом решения алгебраических уравнений,которые можно представить в виде (5.8). В отличие от общих матричныхметодов метод прогонки требует машинной памяти и машинного времени,пропорциональных N , а не N 2 или N 3 ..