Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 14

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 14)

Явная схемаБудем рассматривать следующую краевую задачу для уравнения теплопроводности с постоянными коэффициентами. В области {0 < x < 1, 0 < t 6 T } требуется найти решение уравнения∂u∂2u=∂t∂x2(87)u(x, 0) = u0 (x)(88)удовлетворяющее начальному условиюи граничным условиямu(0, t) = 0,u(1, t) = 0.(89)Примером явной разностной схемы для данной задачи может служить следующая схема:un − 2unm + unm−1un+1− unmm= m+1.τh2(90)При этом краевые уловия дают дополнительные уравнения:un0 = unM = 0,u0m = u0 (mh).Спектральный признак устойчивости для указанной схемы даёт следующее условие:1 − 4τ sin2 πkh 6 1 ⇒ τ 6 1 .2h2 h22(91)Итак, построенная схема даёт первый порядок аппроксимации по времени и второй по пространству (это будетдоказано позже), однако, для использования этой схемы шаг по времени следует брать достаточно малым(τ 6 12 h2 ).5.5.2.

Неявная схемаИспользование неявной схемы позволяет снять условие на зависимость шагов по времени и по пространству.Для этого нужно взять схему+ un+1un+1 − 2un+1un+1− unmmm−1m= m+1.(92)τh2Эта схема, имеет тот же порядок аппроксимации, что и предыдущая (явная) схема (это будет доказано позже),однако условие устойчивости для неё (согласно спектральному признаку) имеет вид−11 + 4τ sin2 πkh 6 1,(93)h22 откуда немедленно следует, что эта схема устойчива при любых τ и h.5.5.3.

Схема с весамиСхема с весами является обобщением уже рассмотренных схем. А именно, вводя обозначениеΛunm =unm+1 − 2unm + unm−1,h2введём схему вида(94)un+1− unmm= σΛun+1+ (1 − σ)Λunm ,(95)mτгде σ — фиксированная константа из отрезка [0, 1]. Таким образом, при каждом значении константы σ получаемнекоторую разностную схему; в частности, при σ = 1 имеем чисто неявную схему, а при σ = 0 — явную.

Отдельного упоминания заслуживает случай σ = 12 . Такая схема (она носит название «схема Кранка – Николсон»)является абсолютно устойчивой и имеет второй порядок аппроксимации по времени и по пространству (это будет доказано позже). При вычислениях по этой схеме требуется решить систему с трёхдиагональной матрицей,в частности, это можно можно сделать методом прогонки.Исследуем погрешность аппроксимации схемы (95) для произвольной константы σ из отрезка [0, 1]. Представим решение задачи (95) в видеyin = u(xi , tn ) + zin ,51где u(xi , tn ) — точное решение исследуемой дифференциальной задачи. Тогда для погрешности получим системууравнений:zin+1 − zin= σΛzin+1 + (2 − σ)Λzin + ψin ,(96)τгде сеточная функция ψ, входящая в правую часть уравнения, равнаψin = σΛun+1+ (1 − σ)Λuni −iun+1− uniiτ(97)и называется погрешностью аппроксимационной схемы (95). Разлагая выражение для ψin по формуле Тейлорав точке (xi , tn + 12 τ ), получим1h2ψin = (u′′ − u̇) + σ −τ u′′ + uIV + O(τ 2 + h4 ).212Учитывая исходное уравнение, окончательно можем записать1h2 ′′nψi =σ−τ+u̇ + O(τ 2 + h4 ).212(98)(99)Отсюда, в частности, следуют указанные порядки аппроксимации для схем при σ = 0, 1, 21 .

Кроме того, следуетh2отметить, что при σ = 12 − 12τсхема имеет повышенный порядок аппроксимации — второй по времени ичетвёртый по пространству.5.6. Сеточные теоремы вложенияСейчас мы докажем простой аналог одной теоремы из курса УрЧП, относящейся к теории пространстваСоболева.Рассмотрим две нормы на пространстве сеточных функций на отрезке [0, X].2kun k1,h := hX un+1 − un 222kun k0,h := hX.|un |2 .(100)(101)Теорема 5.3 (вложения).

Пусть u0 = uN = 0. Имеют место следующие неравенства:22kun k0,h 6 X kun kC ,(102)Xkun k21,h .2(103)kun k2C 6 Первое неравенство очевидно (оценка интеграла максимумом модуля функции, помноженной на длинуотрезка), а второе тривиально. Пусть k — точка, в которой достигается максимум модуля функции un (можносчитать, что k 6 N2 ). Тогдаk−1k−1XX √ un+1 − un√uk =(un+1 − un ) =h·.(104)hn=0n=0Применим к произведению в правой части неравенство Коши, Буняковского, Шварца и ещё многих товарищей:k−1X√|uk |2 6 k · ( h)2 ·h| {z }X/2n=0un+1 − un√hНу вроде мы именно этого и добивались.

522=X2kun k1,h .2(105)5.7. Методы стрельбы и прогонки5.7.1. Метод прогонкиДля решения трёхдиагональных систем, часто возникающих при построении разностных схем, можно использовать метод прогонки.Пусть нам дана трёхдиагональная система Ax = b, где α1 , . . . , αn — числа на диагонали, β1 , . . . , βn−1 —числа над диагональю, а γ2 , . . . , γn — числа под диагональю. Обратите внимание на индексы — они именнотакие, поскольку мы хотим, чтобы в каждой строке они были одинаковые.Попробуем рекуррентно выразить xi друг через друга. А именно, будем искать решение в виде(106)xi−1 = Ai xi + Bi .Числа Ai и Bi называются прогоночными коэффициентами. Попробуем найти формулы для Ai и Bi .Запишем систему, используя эти соотношения: α1 x1 + β1 x2 = b1 ,...γi xi−1 + αi xi + βi xi+1 = bi ,...γn xn−1 + αn xn = bn .(107)Перепишем общее уравнение через прогоночные коэффициенты:γi (Ai xi + Bi ) + αi xi + βi xi+1 = bi .(108)Ура! Нам повезло, и мы избавились от третьей переменной.

Выразим xi :βibi − γi Bixi = −xi+1 +.γi Ai + αiγi Ai + αi{z}|| {z }Ai+1(109)Bi+1А теперь мы видим, что нам ещё больше повезло, поскольку дроби содержат только переменные с индексами i.Обозначая их соответственно Ai+1 и Bi+1 , получаем как раз формулу нужного вида:xi = Ai+1 xi+1 + Bi+1 .(110)Итак, теперь уже ясно, как надо решать систему. Вычислим сначала первые два коэффициента: A2 = − αβ11и B2 = αb11 .

Далее вычисляем все остальные коэффициенты по формуламAi+1 = −βi,γi Ai + αiBi+1 =bi − γi Bi,γi Ai + αi(111)не забывая складывать нажитое непосильным трудом в массивы. Кстати говоря, можно смело использоватьпамять, в которой у нас лежат коэффициенты, поскольку она нам больше не пригодится. Когда все Ai и Biвычислены, можно приступать к нахождению xi .Из последних двух уравнений выразим xn . Имеем(xn−1 = An xn + Bn(112)γn xn−1 + αn xn = bn .Поэтомуxn = −Bn −An +bnγnαnγn.(113)Все остальные переменные вычисляем по формуле xi−1 = Ai xi + Bi .Вот такой простенький метод.

Работает он за линейное время, точнее говоря, за 8n операций. Остаётся толькодоказать, что он работает. А именно, есть проблема с делением на нуль при вычислении Ai и Bi .Теорема 5.4. Пусть α1 = αn = 1 (этого всегда можно добиться). Пусть имеется диагональное преобладание |αi > |γi | + |βi |, причем γi , βi 6= 0. Пусть |β1 | + |γn | < 2. Тогда в алгоритме не возникнет деления нануль, и |Ai | 6 1, то есть погрешность не будет расти.53 Индукция. База: |A2 | = αβ11 6 1.

Пусть уже доказано, что |Ai | 6 1. Тогда|Ai γi + αi | − |βi | > |αi | − |Ai | · |γi | − |βi | > |γi |(1 − |Ai |) > 0.| {z }(114)>0Значит, знаменатель дробей не меньше числителя (по модулю), и потому дробь меньше 1.Далее, рассмотрим два случая.1◦ Если |β1 | < 1, то |A2 | < 1, а потому предыдущее неравенство строгое, и все |Ai | < 1. А если β1 > 1, то поусловию |γn | < 1.

Поэтому αγnn + An 6= 0. 2◦ Если |β1 | > 1, то по условию |γn | < 1. Опять неравенства строгие, и αγnn < 1, а потому αγnn + An 6= 0. 5.7.2. Метод стрельбы2Напомним, что δ (z) := zn+1 − 2zn + zn−1 .Для решения краевых задач вида ∆y − py = f,y(0) = a,y(X) = b(115)можно использовать метод стрельбы. Он заключается в том, что мы строим обычную разностную схему дляисходного уравнения, и для однородного (то есть когда f ≡ 0). А именно,δ2y− pn y n = f n ,h2(116)δ2 z−pz=0.n nh2Тогда зададим краевые условия для y и z: y0 := a, z0 := 0. А теперь зададим (от фонаря) значения функцииво втором узле: y1 := q, z1 := r 6= 0.

Из уравнений, написанных выше, мы найдём все остальные yi и zi .Естественно ожидать, что краевое условие для y будет нарушено. Тогда, исходя из общей теории дифуров,можно найти константу C из уравнения yN + CzN = b. После этого остаётся лишь вычислить «правильное»решение по формуле yen := yn + Czn .Замечание. Насчёт уточнения крайне абстрактного математического термина «от фонаря»: чтобы схемуне разнесло, лучше брать y1 = a + O(h).

Да и z1 лучше взять порядка h. Естественно, условие z1 6= 0 нужно,чтобы не получить тривиальное решение.5.8. Повышение порядков аппроксимации. Метод балансаОсновная идея повышения порядка заключается в получении некоторых соотношений на погрешности в силусамой системы.5.8.1. Пример номер разРассмотрим пример: ′′ y (x) − p(x)y(x) = f (x),y(0) = a,(117)y(X) = b.Легко видеть, что самая тупая аппроксимация этой схемы даст 2-й порядок. А мы сейчас сделаем 4-й.Имеемδ 2 y y (4) (xn ) 2y ′′ (xn ) = 2 +h + O(h4 ).(118)h12А вот теперь применяем трюк: в силу системы′′δ 2 p(xn )y(xn ) + f (xn )(4)y = p(x)y(x) + f (x) =+ O(h2 ).(119)h2Итогоδ 2 yn1− pn yn − δ 2 (pn yn + fn ) = fn .(120)h212δ 2 yn11− pn yn − δ 2 (pn yn ) = fn + δ 2 (fn ).(121)2h1212Заметим, что мы не вылетели из класса трёхдиагональных систем, но получили схему 4 порядка.

Весь фокус —1во «вкусовых добавках» с коэффициентом 12.545.8.2. Пример номер два ′′ y (x) − p(x)y(x) = f (x),y ′ (0) − αy(0) = a,y(X) = b.(122)Здесь нам кроме самой системы ещё и краевое условие нужно аппроксимировать. Самая простая схема даётпервый порядок:y1 − y0− αy0 = a.(123)hФигово. А можно сделать вот что: разложить y(0) по Тейлору до 2 порядка:y1 − y0hh− αy0 − a = y ′ (0) + y ′′ (0) − αy(0) − a + O(h2 ) = y ′′ (0) + O(h2 ).h22(124)В силу системы имеем y ′′ (0) = p(0)y(0) + f (0). Поэтому схемаy1 − y0h− αy0 − (p0 y0 + f0 ) = a.h2(125)имеет второй порядок.5.8.3.

Характеристики

Список файлов лекций