1612725465-542b3179b36a4849700e0b2ecf7da111 (828844), страница 3
Текст из файла (страница 3)
. . , M,u0j = u0 (xj ), j = 0, . . . , N,j = 0, . . . , N − 1,(1.37)которая аппроксимирует дифференциальную задачу (1.36) с порядкомO(τ + h). Используя принцип максимума и спектральный признак Неймана, можно показать, что схема (1.37) при a < 0 будет устойчива привыполнении условияτ−a ≤ 1.hС другой стороны, при a > 0 схема (1.37) будет абсолютно неустойчивой(см.
задачу 1.2).13Таким образом, мы построили две условно устойчивые явные схемыс направленными против потока разностями для уравнения переносас постоянным коэффициентом aun+1− unjunj − unj−1j+a= fjn ,τhun+1− unjunj+1 − unjj+a= fjn ,τhеслиa > 0,(1.38)еслиa < 0.Они устойчивы при выполнении неравенства|a|τ≤ 1.h(1.39)Во внутренних узлах сетки противопоточную схему (1.38) можнозаписать в виде одного уравненияun+1− unja + |a| unj − unj−1a − |a| unj+1 − unjj++= fjn .τ2h2h(1.40)Аналогично выглядит явная противопоточная схема и в случае знакопеременного коэффициента a(x, t).
Например, если на границах отрезка [0, l] выполнены условияa(0, t) > 0,a(l, t) < 0,0 ≤ t ≤ T,то получим такую противопоточную схемуun+1− unjunj − unj−1unj+1 − unjj+ a++ a−= fjn ,τhhun0 = µn0 ,unN = µnl ,n = 0, . . . , M,0uj = u0 (xj ),j = 0, . . . , N,где anj + anj a =,2j = 1, . . . , N − 1, anj − anj a =,2−+(1.41)(1.42)которая аппроксимирует с порядком O(τ + h) начально-краевую задачуut + a(x, t)ux = f (x, t), 0 < x < l, 0 < t ≤ T,u(0, t) = µ0 (t), u(l, t) = µl (t), 0 ≤ t ≤ T,u(x, 0) = u0 (x), 0 ≤ x ≤ l.14(1.43)С помощью принципа максимума можно доказать (см.
задачу 1.10),что для устойчивости противопоточной схемы (1.41) с переменным коэффициентом a(x, t) достаточно выполнения условияmax |a(x, t)|x,tτ≤ 1.h(1.44)1.5. Схема Лакса. Далее для простоты изложения будем рассматривать начально-краевую задачу (1.27) с однородным уравнением переносаut + aux = 0.(1.45)В схеме Лакса разностное уравнение, аппроксимирующее уравнение переноса (1.45), записывается так()un+1− 0, 5 unj+1 + unj−1unj+1 − unj−1j+a= 0,τ2hj = 1, . . . , N − 1. (1.46)Для локальной погрешности аппроксимации имеем выражениеnψh,j=τh2utt −uxx + . . .
,22τпоэтому при τ = O(h2 ) схема Лакса не будет аппроксимировать уравнение переноса, а при законе предельного переходаr=aτ= consth(1.47)будет аппроксимировать с порядком O(τ +h). Таким образом, аппроксимация имеет место лишь при определенной связи между шагами τ и h,т. е. схема Лакса принадлежит к классу условно аппроксимирующихсхем.Для множителя перехода получаем формулуλ(φ) = cos φ − ir sin φ.Следовательно, при законе предельного перехода (1.47) необходимоеусловие устойчивости схемы Лакса заключается в выполнении неравенства r ≤ 1, т. е.aτ≤ 1.(1.48)h151.6. Схема Лакса – Вендроффа. Разностные уравнения этой схемы выглядят так()u∗j+1/2 − 0, 5 unj+1 + unjunj+1 − unj+a= 0,τ /2h(1.49)u∗j+1/2 − u∗j−1/2un+1− unjj+a= 0.τhСхема Лакса – Вендроффа относится к семейству двухшаговых схем.В этой схеме сначала в полуцелых узлах xj+1/2 = xj +h/2 по схеме Лакса вычисляются вспомогательные величины u∗j+1/2 , относящиеся к моменту времени tn + τ /2.
Затем, на втором шаге, вычисляются значенияискомой сеточной функции un+1на (n + 1)-м слое по времени.jДля исследования аппроксимации и устойчивости двухшаговых схемпредварительно производится исключение из схемы вспомогательныхвеличин u∗ . В результате исключения получим одношаговую схему Лакса – Вендроффаun+1− unjunj+1 − unj−1a2 τ unj+1 − 2unj + unj−1j+a=·,τ2h2h2(1.50)которая, как легко проверить, аппроксимирует уравнение переноса (1.45)со вторым порядком по τ и h.Для множителя перехода имеем такое выражениеλ = 1 − ir sin φ − 2r2 sin2φ.2Поэтому необходимое условие устойчивости |λ| ≤ 1 будет равносильновыполнению неравенства(φ )21 − 2r2 sin2+ r2 sin2 φ ≤ 1,2илиφφφ(φ)1 − 4r2 sin2 + 4r4 sin4 + 4r2 sin21 − sin2≤ 1.2222Последнее неравенство равносильно условию r2 ≤ 1. Таким образом,необходимое условие устойчивости схемы Лакса – Вендроффа совпадает с необходимым условием (1.48) устойчивости схемы Лакса.1.7. Диссипация и дисперсия.
Наряду с уравнением переносаut + aux = 0,16a = const(1.51)рассмотрим еще два уравненияut + aux = µuxx ,µ = const > 0,ut + aux + νuxxx = 0,ν = const.(1.52)(1.53)Пусть начальная функция в задаче Коши для этих уравнений представляется в виде ряда Фурье∑u(x, 0) =bm eimx .(1.54)mБудем искать решение каждого из этих уравнений методом разделенияпеременных∑∑bm λt eimx =bm um (x, t),(1.55)u(x, t) =mmгде um (x, t) – гармоника с волновым числом mum (x, t) = λt eimx ,(1.56)λ подлежит определению. Действительная и мнимая часть гармоникипредставляют собой m-волны, длина l которых связана с волновым числом формулой2πl=.(1.57)mТак как уравнения (1.51)–(1.53) линейны, то поведение каждой изгармоник можно рассматривать независимо.
Подставляя гармоникус волновым числом m в уравнение переноса (1.51), получаемln(λ) + aim = 0илиλ = e−aim .Следовательно, если гармоника (1.56) является решением уравненияпереноса, то она имеет видum (x, t) = eim(x−at) .(1.58)Обозначая ξ = x − at, получаемum (x, t) = eimξ = um (ξ, 0).17(1.59)Таким образом, в любой момент времени t > 0 гармоника um получается сдвигом начальной гармоники на величину at. Следовательно,уравнение переноса описывает движение m-волн, которые независимоот их длины распространяются с постоянной скоростью vm = a без искажения своей формы.Легко проверить, что гармоника (1.56) будет решением второго уравнения (1.52), еслиln(λ) + aim = −µm2илиλ = e−aim e−µm ,2т.
е. гармоника в этом случае имеет видum (x, t) = e−µm t eim(x−at) .2Следовательно, для всех гармоник происходит затухание амплитудыволн (диссипация волн). Поскольку m = 2π/l, то короткие волны затухают быстрее длинных. Скорость vm распространения волн не зависитот длины волн и равна по-прежнему a. За диссипацию волн отвечаетчлен µuxx со второй производной от решения.Наконец, подстановка гармоники в уравнение (1.53) даетln(λ) + aim + ν(im)3 = 0илиλ = e−im(a−νm ) ,2откуда получаем, что2um (x, t) = eim(x−(a−νm)t).Таким образом, третье уравнение описывает движение волны без изменения ее амплитуды (без диссипации).
Но скорость ее распространениязависит от длины волныvm = a − νm2 .(1.60)Из этой формулы видно, что волны разной длины распространяютсяс различными скоростями (волны диспергируют). Более значительнымизменениям подвергается скорость распространения коротковолновыхвозмущений (большие m). За дисперсию волн отвечает член νuxxx с третьей производной от решения.18Рассмотрев поведение отдельных гармоник, мы теперь сможем предсказать качественное поведение решения (1.55) задачи Коши для этихуравнений.
Пусть, например, начальная функция u(x, 0) имеет вид ступеньки{1, x ≤ 0,u(x, 0) =(1.61)0, x > 0и a > 0. Разложение такой функции в ряд Фурье (1.54) будет содержатьвесь набор гармоник.Решение задачи Коши для уравнения переноса (1.51) представляетсяв таком виде∑∑bm eim(x−at) =bm eimξ = u(ξ, 0),u(x, t) =(1.62)mmт. е. решением задачи будет движущийся со скоростью a начальныйпрофиль.Решение∑∑22bm e−µm t eim(x−at) =bm e−µm t eimξu(x, t) =(1.63)mmзадачи Коши для уравнения (1.52) с диссипативным членом, в которомкороткие волны сильно затухают, будет иметь вид размазанной ступеньки.Наконец, решение∑2u(x, t) =bm eim(x−(a−νm )t)(1.64)mзадачи Коши для уравнения (1.53), в котором волны разной длины движутся с разными скоростями, имеет немонотонный, осциллирующий характер. Согласно формуле (1.60) при ν > 0 волны малой длины будутиметь скорость ме́ньшую, чем волны большой длины, а при ν < 0 –наоборот.
Поэтому осцилляции будут отставать от основного решения(описываемого первыми гармониками) при ν > 0 и, соответственно, перемещаться впереди при ν < 0.1.8. Дифференциальное приближение разностной схемы.Вернемся к численному решению задачи Коши для уравнения переноса (1.51). В качестве начального профиля возьмем ступеньку{1, x ≤ x0 ,u(x, 0) =(1.65)0, x > x019и проведем расчет по явной противопоточной схемеun+1− unjunj − unj−1j+a= 0,τh(1.66)a = const > 0.В результате получим решение в виде размазанной ступеньки (рис.
3),т. е. решение будет качественно таким же как и решение уравнения (1.52)с диссипативным членом. В чем дело? Ведь мы хотели решить уравнение переноса, в котором диссипативного члена нет. Дело в том, что мыискали численно решение не уравнения переноса, а решение разностнойсхемы. Таким образом, свойства решений аппроксимируемого дифференциального уравнения и аппроксимирующего разностного уравнениямогут не совпадать. Как же в таком случае предсказать свойства решения разностного уравнения?y1.00.003210.51020x30Рис.
3. Графики точного решения (штриховые линии) и численногорешения (сплошные линии), полученного с помощью противопоточнойсхемы (1.66) в моменты времени t = 1 (1); t = 8 (2); t = 15 (3). a = 1;x0 = 10; aτ /h = 0, 5Это можно сделать с помощью метода дифференциального приближения, с которым мы сейчас кратко познакомимся.
Суть этого метода заключается в замене исходного разностного уравнения специальным дифференциальным уравнением, которое обладает всеми свойствами исследуемого разностного уравнения. Поэтому вместо исследованияразностного уравнения исследуют это дифференциальное уравнение,что во многих случаях сделать гораздо проще. Получение дифференциального уравнения, соответствующего разностному уравнению, начинается с записи этого разностного уравнения в виде так называемой теоретической разностной схемы, в которой разностные операторы действуют в том же функциональном пространстве, что и аппроксимируемые ими дифференциальные операторы. Например, разностное уравнение (1.66) записывается в виде следующей теоретической разностной20схемыu(x, t + τ ) − u(x, t)u(x, t) − u(x − h, t)+a= 0.(1.67)τhРешением такой схемы является функция u(x, t) непрерывных аргументов x и t в то время, как решением уравнения (1.66) является сеточнаяфункция uh , определенная только в узлах сетки.Пусть достаточно гладкая функция u(x, t) является решением теоретической разностной схемы (1.67).
Подставим ее в эту схему и выразимu(x, t + τ ) и u(x − h, t) через значения функции и ее производных в точке (x, t) по формуле Тейлора. В результате получим дифференциальноеуравнение, эквивалентное разностной схеме (1.67)ut + aux +ττ2hh2utt + uttt − a uxx + a uxxx + . . . = 0.2626(1.68)Определение. Дифференциальное уравнение бесконечного порядка (1.68), полученное после разложения по формуле Тейлора решенияu(x, t) теоретической разностной схемы (1.67), называется дифференциальным представлением разностной схемы (1.66).Некоторые свойства разностной схемы можно исследовать уже с помощью этого дифференциального представления, но для наших целейбудет удобнее использовать другую форму дифференциального представления, получающуюся в результате исключения из (1.68) всех производных по времени кроме той, которая входит в аппроксимируемоеуравнение (1.51), т.