Консалтология Общая теория консалтинга (стр. 9 )

Следует отметить, что формулы, аналогичные выражениям (8.52) и (8.50), можно вывести и для общего случая, когда проис­ходит изменение шага h и порядка r формулы, если интерполя­ционный полином в уравнении (8.46) выразить через разделенные разности.

С учетом выражения (8.48), выбирая значения для р, q и r, получаем много других формул численного интегрирования. В ча­стности, для линейных многошаговых методов в общем случав справедливо выражение

(8.58)

которое соответствует как явным (b-1=0), так и неявным (b-1≠0) методам.

На базе уравнения (8.58) получена известная неявная формула Гира для жестких дифференциальных уравнений:

(8.59)

из которой легко определяется аппроксимация для производной:

(8.60)

коэффициенты которой в предположении, что вычисления ведутся с переменным шагом,

Ей сопутствует формула прогноза, аналогичная выражению (8.50) при использовании раздельных разностей,

(8.61)

где

Неявный метод интегрирования с обратными разностями пер­вого порядка Брайтона, или формула дифференцирования назад (ФДН) отличается от выражения (8.60) и (8.61) введением в соот­ветствующие выражения обратной разности

∆xn-i= хn+1-i — xn-i, і = l, 2, ..., k.

При этом получаем:

(8.62)

(8.63)

где

Неявный метод численного интегрирования жестких систем диф­ференциальных уравнений, основанный на применении разностей высших порядков, предусматривает прогнозируемое значение вектора неизвестных математиче­ской модели консультируемой проблемы хп+1 определять с помощью интерполяцион­ного полинома Ньютона для неравных промежутков:

(8.64)

где

При этом для алгебраизации производной по времени исполь­зуется формула

(8.65)

где

(8.66)

(8.67)

— разделенные разности, вычисленные в предыдущий момент вре­мени tn, откуда

Формулы (8.60), (8.62) и (8.65) могут использоваться для све­дения дифференциально-алгебраической системы уравнений (8.1) к системе нелинейных алгебраических уравнений (8.2), соответ­ствующей временной точке tn+1. При оценке начального значения искомого вектора неизвестных х(0)n+1 система нелинейных алгебраи­ческих уравнений решается методом Ньютона с использованием формул (8.61), (8.63) и (8.64).

Формулы (8.60), (8.62) и (8.65) эквивалентны друг другу, если их вывод проводить при одинаковых исходных предположениях (постоянном или переменном шаге h), и между их коэффициентами существует однозначная связь. Например, формула второго по­рядка

(k =2) для каждого из перечисленных методов приводится к виду, аналогичному формуле Шихмана:

(8.68)

Об эквивалентности методов (8.60), (8.62) и (8.65) свидетель­ствует и их одинаковая численная устойчивость при k ≤ 6 (рис. 8.2, а и б).

Рис. 8.2. Годографы устойчивости неявных многошаговых методов чис­ленного интегрирования Гира—Брайтона (а, б) и Энрайта (в); (А) — зона устойчивости

Известно, что при k≤2 эти методы обладают А-устойчивостью, т. е. область их абсолютной устойчивости включает всю левую полуплоскость Re(hλ)<0, где λ — собственные значения матрицы Якоби для правой части уравнения (8.45). Для k = 3, ..., 6 они, с одной стороны, А(α)-устойчивы, т. е. область их аб­солютной устойчивости включает бесконечный клин |arg(—hλ)|<α, а с другой,— жестко устойчивы, так как область их абсо­лютной устойчивости ограничена частью левой полуплоскости, при этом в окрестности начала координат обеспечивается не только устойчивость, но и точность. В табл. 8.1 приведены значения α для формул различных порядков k≤6 при относительном изменении шага вычислений q = hj/hj-1.

Таблица 8.1

А(α)-устойчивость многошаговых неявных формул численного интегрирования

Среди часто используемых на практике методов вычисления при­ближения хп+1, кроме рассмотренных выше, следует отметить сле­дующие.

1. Одношаговый А-устойчивый метод трапеций, для которого

(8.69)

Для повышения точности вычислений рекомендуется этот метод использовать совместно со схемой экстраполяции

(8.70)

где xn+1(h) и х2n+2(h/2) — решения, полученные по формуле (8.69) соответственно с шагом h и h/2 для каждой точки tn+1.

2. Многошаговый неявный метод Энрайта, в котором исполь­зуется вторая производная,

(8.71)

годографы устойчивости которого показаны на рис. 8.2, в для k≤7. При решении нелинейных задач для вычисления второй про­изводной

(8.72)

необходимо вычислять матрицы Якоби на каждой итерации и на каждом шаге, поэтому данный метод слишком трудоемкий. Он обладает большей точностью и устойчивостью, чем методы, выра­женные уравнением (8.60), (8.62) и (8.65).

3. Композиционный многошаговый неявный метод Скила, соответствующий по смыслу методу Энрайта и представляющий собой комбинацию метода Гира — Брайтона, или ФДН, и метода Адамса — Мултона (8.51):

(8.73)

где матрица Якоби для уравнения (8.45).

Выражение (8.73) имеет порядок k+1 для каждого шага вы­числений k. Параметр γ(k) в уравнении (8.73) выбирают исходя из улучшения численной устойчивости метода, которая в данном случае охватывает методы до 12-го порядка включительно.

Данные о значениях γ(k) и параметре α для характеристики A(α) устойчивости сведены в табл. 8.2.

Таблица 8.2

A(α)-устойчивость композиционных методов неявного численного интегрирования

4. Блочные или циклические методы, основанные на идее одно­временного вычисления «блока» приближения xn+1, xn+2, ..., хп+l. Блок просчитывается с постоянным шагом на выбран­ном интервале

tn = t0 + nh, где п = ml + j (j = 1,2,...,l и m= 0, 1, 2, ...), при этом сразу используется несколько независимых многошаговых методов:

(8.74)

Например, при выборе исходных многошаговых методов тре­тьего порядка (k = 3) получаем:

Следовательно, если блок приближений [xml-2, xml-1, хml]t известен, то по приведенным формулам определяется следующий блок

[xml+1, xml+2, хml+3]t, a затем процесс приближений повторяется.

Неявные циклические мето­ды характеризуются большей точностью и устойчивостью, чем методы, выраженные уравне­ниями (4.48), (4.50) и (4.53), что следует из табл. 8.3.

Таблица 8.3

А(α)-устойчивость циклических неявных методов численного интегрирования

Их целесообразно применять для моделирования консультируемых проблем, ком­плексные собственные значения моделей которых расположены вблизи мнимой оси плоскости hλ. Установлено, что эти методы жестко устойчивы при порядке k≤25.

Оценка локальной погрешности. Каждый из методов численного интегрирования характеризуется локальной погрешностью, или погрешностью на данном шаге вычислений. Для k-шагового метода, в котором приближение хп+1 к точным значениям x(tn+1) вычис­ляется через ранее вычисленные приближенные значения xn-k, хп— k +1, ..., хп—1, можно записать обобщенную формулу вычислений в виде

(8.75)

Раскладывая функцию в уравнении (8.75) в ряд Тейлора в окре­стности точки хп+1, получаем оценку для локальной погрешно­сти дискретизации как оценку отсечения. Если метод имеет k-й порядок, то приближенно