Fletcher-1-rus (1185917), страница 21
Текст из файла (страница 21)
Например, если используется схема ВВЦП, определяемая алгоритмом (3.41), то сказанное выше означает, что фактическое вычисление 'Т~+~ осуществляется на основании данных о 'Т1" и 'Те и *Т~е.ь так что 'Т(+' = и ('Т",,) + (1 — 2а) ("Т,") + з ('Тг+ ~).
(4.16) Вспоминая, что точное решение алгебраических уравнений Т1 удовлетворяет уравнению (3.41), и пользуясь формулой (4.15), 8 К. Флетчер, т. 1 114 Гл. 4. Теоретические основы с учетом (4.16) получим однородное алгебраическое уравнение $";+ = з$1 1+ (1 — 2з) $1 + аозт"ч. ь (4.177 Если предположить, что заданы начальные и граничные значения, то все начальные ошибки $1О /=2, 3, ..., У вЂ” 1, а также граничные ошибки 5ч, и $", а=О, 1, 2, ..., в соответствии с уравнением (4.17) будут равны нулю. Если при расчете значения Т~ в каком-то внутреннем узле не будут внесены те или иные ошибки округления, то и результирующие ошибки решения останутся равными нулю. Двумя наиболее общеупотребительными методами анализа устойчивости являются матричный метод и метод Неймана В основе обоих методов лежит возможность предсказания, будет ли иметь место рост ошибки как разности между истинным решением численного алгоритма и фактически вычисляемым решением, т.
е. разности, учитывающей погрешности округления. Альтернативный вариант анализа устойчивости связан с предположением о том, что начальные условия представляются с помощью рядов Фурье. Каждая гармоника, или мода, ряда Фурье будет возрастать или затухать в зависимости от вида дискретизированного уравнения, которое, как правило, само порождает конкретную форму выражения для показателя роста (или затухания) соответствующей моды. Если данная мода может неограниченно возрастать, то дискретизированное уравнение имеет неустойчивое решение. Эта интерпретация понятия устойчивости [К!сЫшуег, Мог1оп, 1967) непосредственно используется в методе анализа устойчивости, предложенном Нейманом (см. п.
4.3.4 и 4.3.5). Неограниченный рост некоторой конкретной моды остается возможным, даже если дискретизированные уравнения решаются точно, т. е. если ошибки округления отсутствуют. Если же ошибки округления будут внесены, то сама неустойчивая природа дискретизированных уравнений будет вызывать неприемлемый рост ошибок. Следовательно, методика анализа устойчивости дискретизированных уравнений будет оставаться одной и той же независимо от проявления присущей им устойчивости или неустойчивости. 4.3Л. Матричный метод: схема ВВЦП Данный метод состоит в том, что система уравнений, определяющих распространение ошибки представляется в матричной форме, после чего мы анализируем собственные значения соответствующей матрицы. Для начала этот метод будет про- $4.3.
Устойчивость иллюстрирован на примере схемы ВВЦП, реализуемой согласно (3.41). Подставляя в (4.17) значения 1=2, 3, ..., У вЂ” 1 и замечая, что на границах ошибки равны нулю, т. е. что при всех п будет ь", = 5 = О, получим $з"" = (1 — 2з) $з" + зйз", ~з зьз + (1 2з) вз + з~~' в!'+ = зз1 ~ + (! — 2з) $,". + з$";.ьь $~+~ = з$г з + (! — 2з) $г (4.18) В матричной форме зто можно записать как $"+' = А$", и = О, 1, ... (4.19) где А — квадратная матрица порядка У вЂ” 2, а й" — вектор дли- ны У вЂ” 2; то и другое определяются выражениями (1 — 2з) з з (1 — 2з) з з (1 — 2з) Можно показать, что при увеличении и ошибки в" остаются ограниченными, если все собственные значения Х матрицы А различны и имеют абсолютные значения, меньшие или равные единице; иначе говоря, если ! в !(1 для любых Х.
(4.20) Условие устойчивости (4.20) позволяет использовать только такие значения з, которые удовлетворяют двойному неравенству — 1~(1 — 4ЗЗ1П (о(у 11) ь1. (4.22) Далее формула (9.48) позволяет найти собственные значения трехдиагональной матрицы А в виде Х =1 — 4зз(пз( ), т=1, 2, ..., У вЂ” 2. (4.21) Гл. 4. Теоретические основы 116 Правая часть этого неравенства удовлетворяется при всех т и з, тогда как для выполнения его левой части необходимо, чтобы было «!и '! 1 зз!п (,з!т В)~~я а это справедливо для всех т, если з ( 1/2. Поэтому алгоритм ВВЦП является устойчивым, если з ( 1/2. Форма матричного уравнения (4.19) остается той же самой, если в формулах (4.18) будет фигурировать не вычислительная ошибка, а само решение. Следовательно, матричному методу можно дать и другую интерпретацию как методу, исследующему рост индивидуальных мод в представлении начальных условий с помощью ряда Фурье.
Как будет указано в п. 4.3.3, изменение характера граничных условий при х = О и/или х = 1 приводит лишь к незначительной модификации метода. 4.3.2. Матричный метод; двухслойная схема общего вида В данном пункте мы покажем, как работает матричный метод в случае двухслойной схемы общего вида, примененной к уравнению диффузии. Эта схема, данная в форме (7.24)„ имеет вид а+1 н т — т — й1л „„Т1 — а(1 — (!) Т,„Т1 = О, (4.23) и+! н где Т.„„Т1 — — (Т,, — 2Т; + Т„!))бх'. В уравнении (4.23) величина а — коэффициент тепловой диффузии, тогда как параметр р контролирует степень неявности (см.
$ 7.2). Численная ошибка (4.15) определяется уравнением, эквивалентным (4.23). Это уравнение записывается в форме — з В!'-+ ' + (1 + 2зй) 0+' — з К+ !' = з (1 — Р) М- ! + + (1 — 2з (1 — ())] $1 + з (1 — Р) Ву+!. (4.24) Уравнение (4.24) пригодно для внутренних узлов. При наличии граничных условий Дирихле уравнения в концевых точках не требуются. Если уравнение (4.24) повторить для всех узловых точек, где 1' = 2, 3, ..., ) — 1, то в матричной форме можно написать ~ца+! ~цн (4.25) ф 4.3.
Устойчивость где 1+ 2зр — зр — зр 1+ 2зр — зр — з[3 1+2з5 — зб — з[3 1+ 2зр 1 — 2з (1 — р) з (1 — р) з(1 — р) 1 — 2з(1 — 6) з(1 — р) В= з(1 — 6) 1 — 2з(1 — [3) з(1 — 6) з (1 — р) 1 — 2з (1 — 0) Алгоритм (4.23) остается устойчивым, если величины собственных значений матрицы А-'В не превышают единицы.
С учетом структуры матриц А и В это эквивалентно ограничению [(Л ) У(Л„) [< !.0. (4.26). Принимая во внимание симметричный трехдиагональный характер матриц А и В, можно получить аналитические выражения для собственных значений [М!!с[!е!1, бг![[![[тз, 1980) Лл=[+2зр — 2зрсоз(~~ ) = =1+4зй~ш'(2(~' !)). Лв — — 1 — 2з(! — 5)+ 2з(1 — й) сов ( ! ) = =1 — 4 (1 — р)э!пт( ~ ).
Следовательно, условие устойчивости имеет вид $Т ~ ! — 44 (! — Р) в)п [)п/2 (1 — !)1 [ < ! ! + 44Р в)пт [)п/2(У вЂ” !)1 причем )п т. е. необходимо, чтобы 1 — 4з(1 — р) < 1+ 4зр, что удовлетворяется, и чтобы 1 — 4з(1 — р) ) — 1 — 4зр, иначе говоря. чтобы 2 )'4з(1 — 2р), или з < 0.5/(1 — 2р), если р < 0.5. 418 Гл. 4.
Теоретические основы Следовательно, если р ( 0.5, то для устойчивости требуется, чтобы з ( 0.51'(! — 2р). Если р ) 0.5, то алгоритм (4.23) является безусловно устойчивым. Весьма желательное свойство безусловной устойчивости будет подробнее обсуждаться в $ 7.2.
4.З.З. Матричный метод: граничные условия для производных Матричный метод может применяться и тогда, когда граничные условия имеют форму Неймана, т. е. когда они форму.лируются для производных. Предположим, что — =б при х=О. дТ дх (4.27) Это условие можно реализовать путем введения фиктивной точки Те, такой, что — = б = дТ 1 Те — То Т, = Те — 25 Лх. дх (1 2 ах Следовательно, если рассмотреть все подобные уравнения, то А$"+' = В$" + С, (4.29) где 1 + 2в(! — 2зб — зр 1+ 2вб — зр — з(1 1+ 2зр — зр — 2з (1 — 8) 2в (1 — 8) з (1 — 8) 1 — 2з (1 — 8) 0 1 — 2з (1 — р) з (! — 8) з (1 — 11) 1 — 2з (1 — В) Если это соотношение использовать для функции распределения ошибок, то в сочетании с (4.24) такой прием даст следующее соотношение, центрируемое в узле с номером 1: 41+ 2зй) ь",+ — 2зЯ,"+ =(! — 2в(1 — р)($",+2з(1 — 8)$," — 2зббх.
(4.28) 119' э 4.3. Устойчивость Как и прежде, для обеспечения устойчивости требуется, чтобы !(2,„ ,,) ) < 1.О. В общем случае, когда нельзя получить явных формул для собственных значений матрицы А-'В, (4.29) удобно переписать в виде 6"+' = 06" + А 'С, (4. 30). повторяемых до тех пор, пока Хн+' не совпадет с Х'. Если сходимость достигнута, то Х"+' является собственным вектором, соответствующим максимальному собственному значению,.
которое определяется максимальным элементом вектора Х"+' Обычно на каждой итерации масштаб подбирается так, чтобы максимальный элемент вектора Х" равнялся единице. Стартовое значение Х' выбирается произвольно. Карнаган, Лютер и Уилкс предлагают программу такого расчета и обсуждают возможные трудности, связанные с применением метода. Применительно к двухслойной схеме общего вида (4.23) максимальные собственные значения матрицы 0 = А — 'В для граничных условий Дирихле и Неймана, поставленных в точке х = О, приводятся в табл. 4.2 для некоторого набора значений р и 3. Эти результаты были получены с помощью степеннбггр Таблица 4.2. Данные о численной устойчивости дли двухслойной схемы общего вида Граничные услааия Дирикле при к О Граничные услааня Неймана нри к Е Максимальные ссбстаеиные зиачени» Максимальные ссбстаеииме значения 0.00 0.50 0.5! 0.52 0.30 !.20 1.30 0.60 1.20 1.80 0.9899 — !.0289 — !.0686 0.9753 — !.0273 0.9742 0.9721 0.9660 — !.0053 — 1.0446 — 0.9534 — 1.0161 0.9221 0.9!59 где 0 =А 'В.
Максимальные собственные значения матрицы 0 можно найти, пользуясь степенным методом [Сагпа)!ап е! а!., 1969]. Это потребует расчетов по формуле Х"+' = 0Х", (4. 31) Гл. 4. Теоретические осковы Я 20 метода. В случае граничного условия Дирихле в точке х = 0 результаты согласуются с требованием з( — ', если б (0.5, 0.5 1 — 2р ' обеспечивающим устойчивость (Х,„( 1). Ясно, что необходимость выполнения граничного условия Неймана (4.27) несколько уменьшает устойчивость схемы при одинаковых значениях р и з. 4.3.4. Метод Неймана: схема ВВЦП При определении критериев устойчивости наиболее общепринятым является метод Неймана, так как в общем случае его легче всего применять и так как он самый непосредственный и надежный.
К сожалению, однако, этот метод можно применять только для установления необходимого и достаточного условия для линейных задач с начальными условиями и постоянными коэффициентами. Как правило, практические задачи охватывают и переменные коэффициенты, и различного рода нелинейности, и усложненные варианты граничных условий. А если что-то из этого имеет место, то метод может применяться только локально и при замораживании нелинейных членов. Для такой более общей ситуации метод Неймана обеспечивает необходимые, но не всегда достаточные условия устойчивости. Строго говоря, метод применим только для внутренних точек.