Anderson-et-al-1 (1185923), страница 17
Текст из файла (страница 17)
Если разностная схема устойчива, то рост любого возмущения, вводимого на и-м шаге по времени, ограничен; для неустойчивых конечноразностных схем возмущение возрастает. Рассмотрим распределение погрешности на сетке в любой момент времени. Для удобства выберем момент времени 1=0. Схематически это распределение погрешности показано на е(х,О) Рнс.
3.)г. Начальное рас~ределенне погренгностгг. рис. 3.11. Предположим, что погрешность е(х,г) можно представить в виде суммы ряда Фурье е(х, 1)= ~ Ь (1)ег "', (3.104) причем период основной частоты (пг = 1) равен 2Е. Нас интересует решение в интервале длины /., поэтому волновые числа /г =тп//, т= О, 1, 2, ..., М, где М вЂ” число отрезков 'длины гхх, помещающихся в отрезке длины /.. Например, если интервал длины 2Л разбит на отрезки пятью узлами, то М= 2, а в сумму ряда входят лишь гармо- ники =й /2~=~/2/., /~=!/2/., =1, /о=О, т=О, /а=!//-, т=2. Напомним, что частота указывает на число волн, помещающихся в отрезке длины 2/..
Если пг =О, то /е — — О, а соответствующее слагаемое описывает стационарную составляющую решения. Так как погрешность удовлетворяет линейному уравнению, то поведение каждой гарглоники, входящей в (3.104), можно рассмотреть независимо. ! ассмотрим член е (х, г) = Ь (1) ег м". Гл. 3. Основы метода конечнык разоостей 90 в виде хе' ". При !=О «л| г=е, тогда 1) е«те~ага« Будем искать решение уравнения (и = 0) оно имеет вид е . Пусть «' а«лт «! х=е =е, а (х (3.105) причем /с«, вещественно, но а может быть и комплексным. Подставляя (3.105) в (3.!03), получим «и+ли и к «т и «» «с ы <«+л«у - «г ы «, «с ~а <«-л«!) е' ' л'е'"" — е е '" =г(е е '" — 2е е +е е '" где г = ссЛ»/(Лх)а.
Газделив на е'е "и использовав соотношение соз 6 = (е1а+ е-'а)/2, получим е«л' = 1+ 2г (соз р — !), где р = ам»ах. При помощи тригонометрического тождества яп'(р/2) = (1 — соз 6)/2 перепишем последнее соотношение в окончательном виде е"' = 1 — 4г я'п'(6/2)» (3.106) Так как для каждой гармоники з"+'=е'л'е", то погрешность !' округления не будет возрастать на каждом шаге по маршевой координате (времени), если !е«л'~ не превосходит единицы.
Следовательно, разностная схема устойчива при ! 1 — 4г э! па (()/2) ~ ~ (1. (3.107) Коэффициент 1 — 4г з1п«((!/2) (равный отношению з" +'/е") называют коэффициентом (или множителем) перехода и обозначают через 6. Отметим, что при анализе Фурье устойчивости конечноразностных схем реальные граничные условия не учитываются; вместо них для гармоник выставляют обычно периодические граничные условия. При решении неравенства (3.107) надо рассмотреть два возможных случая: 1.
Гсли (1 — 4гз!пт((!/2)) ) О, то 4гз!па((!/2)) О. 2. Гели (1 — 4г яп'((!/2) ) ( О, то 4» япа(6/2) — 1 < !. Первое неравенство выполняется для всех г ) О, а второе — лишь при г ( 1/2. Последнее неравенство и является условием устойчивости рассматриваемой конечно-разностной схемы; оно накладывает ограничение на соотношение шагов по времени и пространственной координате. Теперь можно легко объяснить, почему в примере, приведенном в конце п. 3.3.3, получались физи- $ З.б.
устойчивость нонечно-разностных схем 91 чески нереальные значения температуры. Шаг Л1 в этом примере был вдвое больше максимально допустимого условием устойчивости конечно-разностной схемы, поэтому решение резко росло. При сс(М/Ьхт) =!/2 устойчивость рассматриваемой разностной схемы легко проверить. Следует заметить, что выражение (3.106) для коэффициента перехода можно получить при помощи подстановки разложения погрешности в ряд Фурье (3.104) в конечно-разностное уравнение. Мы предлагаем читателю проделать это в качестве упражнения. Метод Неймана (или Фурье) применим и к анализу конечноразностных схем решения уравнений гиперболического типа. В качестве примера рассмотрим одномерное волновое уравненнс первого порядка — +с — =О, ди ди д1 дх (3.108) описывающего распространение волны, бегущей со скоростью с.
Это уравнение имеет только одно семейство характеристик, являющихся решением характеристического уравнения х~ = с. Общее решение уравнения (3.108) имеет вид и(х — с1) = сова(. Это решение требует, чтобы начальные данные, заданные при 1= 0, без изменения переносились вдоль характеристик. Лаке (1.ах, 1954] предложил конечно-разностную схему первого порядка для решения уравнений такого типа: л / л л К ил+'= ~+' ~ — с — ~ ~+' ' ').
(3.109) ! 2 ох 2 Первое слагаемое в правой части есть среднее значение неизвестной на предыдущем шаге по времени, а второе — конечноразностный аналог первой производной по пространству. Пусть ит — — е е м . Подставляя это выражение в разностное уравнелафИх ние, получим, что коэффициент перехода принимает вид еавг соз р — И з(п (1. Следовательно, схема Лакса устойчива прн 1 соз 11 — Ь з(п 51(1, где параметр т = сб1/Лх называется числом Куринта, Так как квадрат модуля комплексного числа равен сумме квадратов его вещественной и мнимой частей, то рассматриваемая конечноразностная схема устойчива при ~т(~(1. (3.110) Гл.
3. Основы метода конечных разностей Следовательно, и в этом случае устойчивость конечно-разностной схемы ойределяется соотношением шагов.по времени и пространственной координате. Условие (3.110) называется условием устойчивости Куранта — Фридрихса — Леви (КФЛ), которое подробно обсуждалось в связи со сходимостью и устойчивостью в исторически важной работе Куранта и др.
(Сопгап! е1 а1., 1928], которую обычно рассматривают как основополагающую для развития современных численных методов решения уравнений в частных производных. Коэффициент перехода, или, как его еще иногда называют, множитель перехода, для некоторой конечно-разностной схемы зависит от шагов сетки и волнового числа. Для конечно-разностной схемы Лакса коэффициент перехода имеет вид 6 — совр !ч знт 8 — ! 6 (сна —.1/сваей+ чз з(изб ееагс!н(-чзнв! (3.111) где З! — фазовый угол.
Из последнего соотношения ясно, как коэффициент 6 зависит от числа Куранта ч и параметра частоты р. Зависимость коэффициента 6 от этих параметров построена на рис. 3.12 при разных числах Куранта в фазовой плоскости. Тщательный анализ этих кривых позволяет сделать ряд интересных выводов.
Фазовый угол в методе Лакса меняется от 0 для малых частот до — и для больших, что легко проверить, вычислив этот угол для обоих предельных случаев. При числе Куранта, равном единице, все гармоники распространяются без затуханий. При числах Куранта, меньших единицы, низкочастотные и высокочастотные гармоники меняются слабо, тогда как гармоники со средней частотой затухают довольно сильно. По представленным на рис. 3.12 кривым можно оценить и затухание по величине фазового угла.
Очень важно понять физический смысл условия КФЛ (3.110) для гиперболических уравнений. Рассмотрим волновое уравнение второго порядка ии — с и„„=О. з (3. 1! 2) Его характеристики имеют вид х+ сг = сопи(= сы х — с! = сопз1 = с.. Решение в точке (х,1), как показано на рис. 3.13, зависит лишь от условий, заданных между пересекающимися в этой точке характеристиками. Следовательно, аналитическое решение в точке (х, г) зависит лишь от информации, содержащейся в области между характеристиками с~ и сь !.О О.в 0.7 О,В Ъ.
О.В 0.4 О.в 0,2 0.1 0.0 0.1 0.2 0,3 0.4 0.5 0.6 ОЛ О.В 0.8 !.0 -ф/тг Рнс. 3.!2 Зависимость модуля коэффипиента перехода от фазы при различных числах Куранта для схемы Лаков. Рис. 3.!3. Характеристики одномерного волнового ураннения второго порядка. Гл. 3, Основы метода конечных разностей Условием устойчивости большинства явных схем решения гиперболических уравнений в частных производных является условие КФЛ, имеющее вид (стэг/тах((1, т.
е. то же условие, что и (3.110). Вго можно представить в другой форме: (Ж/стх)а ( 1/са. Тангенс угла наклона характеристик с(г/с(х = .+1/с. Следовательно, условие КФЛ эквивалентно требованию, чтобы область зависимости аналитического решения лежала внутри области зависимости численного решения. Область зависимости численного решения может быть шире области зависимости аналитического решения, но не наоборот. Возможна и простая геометрическая интерпретация условия КФЛ: угол наклона прямых, соединяющих узлы (1 -+ 1, и) и (1, и + 1), по абсолютной величине должен быть меньше угла.
наклона характеристик. Условие КФЛ имеет физический смысл. На основе проведенных рассуждений можно ожидать, что численное решение будет ухудшаться, если слишком много излишней информации о течении учитывается в данной точке, т. е. если величина с(та1/Лх) сильно отличается от единицы. При численных расчетах так и получается. При применении явных схем для решения гиперболических уравнений наиболее точные результаты получены для близких к единице чисел Куранта. Напомним, что к этому же выводу мы пришли при анализе коэффициента перехода в методе Лакса (см.
рис. 3.12). Прежде чем перейти к описанию методов анализа устойчивости конечно-разностных схем для решения систем уравнений в частных производных, приведем пример применения метода Неймана к анализу устойчивости конечно-разностных схем для решения уравнений с большим числом независимых переменных. Пример 3.4. Определить условие устойчивости простой явной схемы для решения двумерного уравнения теплопроводности ди даи даи — = а — + а —. дт дка дуа ' Конечно-разностный аналог рассматриваемого уравнения имеет вид ци+1 ци 1 г .(ци, — 2ц" + ц" ~ «) + -1- г„(ц,с„, — 2ц," а + ц,",,).
95 й 9.6, устойчзвость конечно-разностных схем Здесь г, = сс(А//Ьх'), а г„= сс(А!/Луз). Будем искать решение в виде суммы гармоник нида з м га з га и иве=в е 'е х . Если р, =/г,йх, рз = Й„Лу, то е"' = 1 + 2г, (соз й, — 1) + 2гв (соз рз — 1). Используя тригонометрическое тождество з(пз(р/2) =(1— — совр)/2, для коэффициента перехода получим выражение 6 = 1 — 4г„яп — ' — 4гн яп — '. зй! 2~а Следовательно, разностная схема устойчива при [ 1 — 4г„япз ([)1/2) 4г з! пг(йз/2) [( 1 т.