Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Установление сходимости, а лучше, и порядка сходимости (еще лучше, с хорошей оценкой константы С) есть основная цель теоретического обоснования метода приближенного решения. Этому вопросу посвящен § 7.
Схема Адамса. Опишем общую конструкцию схем численного интегрирования, достоинством которой является ее экономичность. Каждый шаг интегрирования требует только одного вычисления правой части /, в то же время порядок точности метода может быть (формально) любым желаемым. В методах Рунге—Кутты (они описаны в § 7) число вычислений / на шаг равно порядку точности метода.
Итак, пусть задача решается на равномерной сетке, значения Xn (и все предшествующие Xn..і, хп_2, ..., х0) уже найдены. По значениям ft = f(xn_ +l)' для г = 0, 1, ..., р (р определяет порядок точности метода) в узлах Iі — tn_p+i построим интерполяционный полином степени рг.
Lp{t, {/'},{/,» = 2/'(0/і-
г=0
Его можно применить и для экстраполяции функции f[x(t)\ на интервале [tn, tn + x = tn+l].
Теперь используем очевидное тождество
‘n + t
x(tn + x) = x(tn) + \f[x(t)]dt. (6)
t
п
Заменяя /[*(/)] интерполяционным полиномом и вычисляя интеграл, получаем формулу
*„+1 = *» + xOoZp + aIfp-I + ••• + «р/о)> (?)
где а0, ait ..., ар — некоторые универсальные (не зависящие от шага т) числа. Они очевидным образом вычисляются через интегралы от базисных интерполяционных полиномов l‘p(t). При вычислении а і делается замена переменных t = tn + ijx и рассматривается стандартный шаг по ?, равный единице.
Оценим погрешность аппроксимации, предполагая /(х), а следовательно, и решение х(г) достаточно гладкими. Погрешность экстраполяции \\f[x(t)} — Lp(t)\\ — 0(хР+1) (см. § 3). При интегриро-
62
ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ
[Ч.І
вании по [tn, tn+l] в оценке погрешности вычисления интеграла появляется еще множитель т. Переписывая (6), (7) в форме, дающей в пределе дифференциальное уравнение, получаем
X<t~+X\ —--=ClfJ(Xn) + A1ZOn..,) + ••• + apf(xn-p) + 0(xP + l).
Итак, порядок погрешности аппроксимации равен числу используемых в (7) точек р+ I.
Неудобством метода является необходимость помнить некоторое число прошлых значений Это, конечно, мелочь, если не считать самого начала процесса интегрирования, когда прошлого просто нет. Приходится несколько первых шагов выполнять нестандартно, например методом Рунге—Кутты (см. § 7).
Метод Адамса является характерным примером схемы, формальный порядок которой превышает порядок дифференциального уравнения, Стандартный алгоритм начинает работать лишь при задании, кроме начальных данных х0, еще и значений X1, х2, ..., хр. Таким образом, общее решение разностного уравнения содержит больше, чем нужно, произвольных постоянных и, следовательно, какие-то лишние «решения».
Полезно иметь представление о том, во что переходят лишние решения в пределе при г -*• 0. Рассмотрим простейшее уравнение
X = ах и две схемы второго порядка — примитивную схему (4) и квалифицированную схему Адамса второго порядка:
-+1~- = I /(хп) - \ /(*„_,), /(*) S ах.
В этом простом случае можно вычислить и проанализировать общие решения разностных уравнений. Они ищутся в стандартной для однородных разностных уравнений с постоянными коэффициентами форме хп = ClQl + C2q%, где Q1, Q2 — корни характеристического уравнения, C1, C2 — постоянные, определяющиеся в данном случае начальными данными X0 и X1.
Характеристические уравнения получаем, подставляя дп в уравнение. Для простейшей схемы (4) имеем
д2 — 2axq —1 = 0, т.е. Q1 г = ах ± Vl + (at)2.
Первый корень (при |ах| <«1):
Qi = I + ах + і (ax)2 + 0(х3) = еах + 0(х3),
дії _ еапх(1 _|_ о(х2)), П 5» 0(1/х).
«5]
ИНТЕГРИРОВАНИЕ ЗАДАЧИ КОШИ ДЛЯ СИСТЕМ ОДУ
63
Это решение в пределе дает решение дифференциального уравнения. Второй корень:
<?2 » — I + ах, q% « (—1)пе~апх,
т.е. сеточная функция никакого разумного предела не имеет (из-за множителя (—1)" = (—1)'/х). Это есть паразитическое решение, появившееся из-за превышения порядка разностного уравнения над порядком дифференциального.
В общем решении хп = C1^i + C2q\, и для того чтобы схема имела второй порядок точности, нужно обеспечить соотношения C1 = X0 + 0(г2), C2 = 0(г2). '
Аналогичные выкладки для схемы Адамса приводят к характеристическому уравнению
q2 = q +1 aiq — ~ ах.
Его корни:
<7i,2 = i + f ax±V| + \ ах + ^ а2х2.
Предоставим читателю убедиться, что
q2 ах/2, = 1 + ах + а2х2/2 + 0(х3) = еах (1 + 0(х3)).
Таким образом, q{lx = eat(\ + 0(х2)), а паразитическое решение у2 ^ (ах)" очень быстро стремится к нулю (мы, конечно, считаем [ ах I -«sc 1, например ахгкО.1).
Итак, выбор X1 в схеме Адамса должен обеспечить соотношение
C1 = X0 + 0(х2). Полезно проверить, что выбор X1 = х0 + xf(x0) обеспечивает требуемые соотношения. Более высокий интеллектуальный уровень схемы Адамса (по сравнению с примитивной схемой с центральной разностью) сказался в том, что паразитическое решение этого метода очень быстро убывает — как (ax)tlx.
Численное интегрирование на ЭВМ. Представленный выше анализ погрешностей приводит к выводу, что точность численного интегрирования тем выше, чем меньше шаг х. Это верно только до известного предела — до тех пор, пока погрешности округления, связанные с конечной разрядностью машинной арифметики, остаются пренебрежимо малыми.
