Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
^п+1 = О,. + cO + т/0„+1)-
В-аппроксимация. В стандартной теории численного интегрирования установление аппроксимации, использующее предположение
о гладкости искомого решения, является тривиальным упражнением, решаемым простыми разложениями в ряд Тейлора. Ho установление В-аппроксимации не тривиально, и некоторые схемы, аппроксимирующие уравнение в обычном смысле, свойством Б-аппро-
230
ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
[Ч. II
ксимации просто не обладают. Поясним это парадоксальное на первый взгляд обстоятельство.
Вычисление погрешности аппроксимации состоит в подстановке в разностные уравнения ограничения на сетку точного решения (предполагаемого гладким). Возьмем явную схему Эйлера и вычислим невязку а = (Хп+1 — Хп)/ъ ~ /Oxn)- Так как
мент времени, а эта величина по предположению есть 0(1).
То же самое относится и к неявной схеме Эйлера. Однако известная схема «трапеция» допускает две в обычном случае равноценные модификации. В одном варианте
В этом случае можно использовать интегральное тождество
где Xn+n2 = X(tn + т/2). Обозначение 0(т2) му употребляем в
дальнейшем только для величин, допускающих оценку типа Ст2, где С = 0( 1) и не зависит от жесткости системы. В данном случае эта величина зависит только от гладкости кривой X(t).
Оценим теперь
*в+і = *(*» + *) = + т Hq + у x(tn + е*)> коп) = /(*„),
то здесь все в порядке: в оценке участвуют только X в какой-то мо-
Gt
X
T
2
/п + Х
=S /[*(01 dt.
Так как f[X0) \ = X(t) — гладкая функция, то
и здесь все в порядке: а = 0(т2). Однако в другом варианте:
аналогичные оценки не проходят. В самом деле,
= *(*» + I) + = Хп+т + 0( т2),
ЖЕСТКИЕ СИСТЕМЫ ОДУ
231
где
(в силу гладкости Х(()). Дальнейшие оценки дают
/(*„+,/2 + г) = f(Xn+U2) + fx(Xn+m + Br) г.
Последнее слагаемое не есть 0(г2), хотя г — 0(г2); величина /х вычисляется не на траектории, да и вообще она уже зависит от жесткости системы. Оценки типа /хО(г2) в 5-теории не принимаются. Считается, что такие величины могут быть сколь угодно большими. Конструирование разностных схем, обладающих свойством 5-аппроксимации, не так просто, как могло бы показаться.
Схемы интерполяционного типа. Рассмотрим пример схемы, для которой удается сравнительно просто доказать требуемые В-теорией свойства. Опишем стандартный шаг численного интегрирования — переход от хп в момент tn к хп+1 в момент In 4- т. Основу схемы составляет аппроксимация очевидного тождества (для точного решения)
<„+г
x(tn +г) = x(tn) + \ f[x(t)]dt. (34)
tn
Обозначим f[t\ = f\x(t)] и вычислим интеграл по какой-нибудь хорошей квадратурной формуле, считая пока f[t] известной функцией.
Сделаем замену переменных t = tn + \х (Ё, Є [0, 11) и введем на [0, 1 ] сетку узлов 0 < < ... < < 1. Обозначим
tl = tn + ?'т, /' = /[*'], Ii(X) — интерполяционный базис на сетке {!'}, состоящий из полиномов степени S- 1, определяемых свойствами Ii(X) = 6{. Построим интерполяционный полином Лагранжа:
Z(X) = JfiIi(X), (35)
/=1
и вычислим приближенно
t +I I
п *
\ т dt^ т ( j?(X) dx
t о
я
Если f[t\ — гладкая функция, погрешность квадратуры есть
0(ts+1). (Сделав замену t = tn + |х, мы строим полином для функ-
X(t)+X(tn + x) і т \ *
Г =...."...------------x[ta + ± ) = 0(х2)
232 ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ [Ч. II
ции /[!=] = f[tn + ?т]. В оценку точности интерполяции полиномом степени s — 1 входит s-я производная f[\] по которая, очевидно, есть s-я производная по /, умноженная на т*. Еще одна степень т добавляется при интегрировании.) Теперь формула (35) переписывается в виде стандартной механической квадратуры:
*n+i = хп + т 2 (36)
/=і
где У зависят только от сетки
Пока этой формулой мы воспользоваться не можем, так как значения fJ не известны. Ho /у = /[х(/у)], а промежуточные значения
x(tJ) (обозначим их Yi) могут быть приближенно вычислены точно так же:
УІ = Хп+ S /[*1 dL (37)
Используя для вычислений интегралов тот же интерполяционный полином, получаем серию соотношений:
S
Yi = Xn +т 2 atJf(YJ), /=1,2, ...,s. (38)
J = і
В результате мы имеем следующий алгоритм численного интегрирования:
а) зная хп, решаем систему s dim х нелинейных уравнений (38)
относительно неизвестных Y1; здесь же вычисляются и f(Y‘) (/= 1, 2, ..., s);
б) вычисляем хп+1 по формуле (36).
Подчеркнем, что коэффициенты схемы {агу}, {АД не зависят от вида системы и ее размерности. Они зависят только от выбора узлов ?/.
Почти очевидна следующая теорема.
Теорема 2. Схема интерполяционного типа (36), (38) обладает свойством 5-аппроксимации порядка s.
Доказательство. Вычислим невязку в (36) при подстановке в это соотношение сеточной функции Xn:
S
а=*»+,-*„-¦* 2 ^/[*(*')]•
J=1
ЖЕСТКИЕ СИСТЕМЫ ОДУ
233
Ho
(п + Х
хп+1-хп= \ ладі dt,
*п
а г ^ Ы /[X(Ii)] есть интеграл от интерполяционного полинома степени 5—1, построенного для гладкой функции /[X(t)] = X(t). Таким образом, эта сумма отличается от Хп+1 — Хп на величину 0(т5+1) (погрешность интерполяции 0(х4) интегрируется по интервалу длиной т). Итак , невязка a = 0(xs+l). Точно таким же образом оцениваем невязки в соотношениях (38). Обозначая Xі = X(tj), имеем
