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

(8.76)

где х (k+1)(τ) — (k + 1)-я производная от решения x(t) в проме­жуточной точке tп < τ <tn+1; 0(hk+2) — члены отсечения, содер­жащие h в степени

k + 2 и выше и которыми обычно пренебре­гают; Ck+1 — константа, характеризующая численный метод.

Например, для формулы Адамса—Мултона (8.51)

для формулы Гира (8.60) при постоянном шаге

для метода трапеций с учетом k = 1 С2 = 1/12 и др.

Следует отметить, что выражение (8.76) мало удобно для прак­тических применений, поэтому для оценки локальной погрешно­сти εт применяют различные косвенные методы.

Прежде всего, можно производные высоких порядков выразить приближенно через разности высоких порядков, используемые, например, в формуле (8.65):

(8.77)

Далее можно применить экстраполяцию значений переменных xn+1(h) и xn+1(h/2), полученных интегрированием от точки tn до точки tn+1 один раз с шагом h и другой раз с шагом h/2, при этом

(8.78)

где k — порядок используемого метола интегрирования.

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

определим оценку локальной погрешности

(8.79)

Здесь, как и при экстраполяции, имеется выбор: продолжать вычисления по значению, полученному методом порядка k или ме­тодом порядка k + 1.

Если учесть, что формулы прогноза (8.61), (8.63) и (8.64) имеют порядок на единицу больше, чем основные расчетные уравнения (8.60), (8.62) и (8.65), то к способу вложенных методов можно отне­сти следующую оценку локальной погрешности:

(8.80)

Следует отметить, что в выражение (8.80) входит то же зна­чение х(0)п+1, которое используется в итерациях Ньютона при реше­нии нелинейного уравнения (8.2). Справедливость уравнения (8.80) непосредственно следует из выражения (8.77), если учесть фор­мулу (8.64).

Рассмотрим взаимосвязь локальной погрешности εт с общей (глобальной) погрешностью εmах, определяемой консультантом при моделировании различных консультируемых проблем и равной раз­ности вычисляемых хп+1 и точных х*п+1 значений вектора пере­менных уравнений модели.

Чаще всего исходными при моделировании являются интервал интегрирования Т (t0, t0+T) и максимальная допустимая глобаль­ная погрешность Еmах на этом интервале. Тогда допустимая на еди­ничном шаге вычислений локальная погрешность определяется по формуле

(8.81)

Экспериментальные исследования показывают, что при исполь­зовании выражения (8.81) решение получается часто более точным, чем требуется, особенно в начале интервала интегрирования, в ре­зультате чего выполняется завышенное число шагов вычислений.. Определенную избыточность решения можно уменьшить, если вместо соотношения (8.81) ввести другое, учитывающее действитель­ный характер изменения переменных при решении уравнения мо­дели объекта.

Например, линеаризация уравнения (8.1)

с учетом возможных аппроксимаций для производной вида (8.60)

приводит к выражению:

(8.82)

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

где

— локальная погрешность вычислений в единицу времени, определяемая с учетом уравнений (8.80) или (8.76). Тогда по аналогии с предыдущим записываем

(8.83)

и, вычитая линеаризированные уравнения (8.82) и (8.83), находим искомое уравнение для общей (глобальной) погрешности

(8.84)

где

Следует отметить, что матрица А в полученном уравнении (8.84) совпадает с матрицей линеаризованного уравнения модели консультируемой проблемы, поэтому вычисления глобальной погрешности при извест­ных матрицах L и U выполняются только процедурой прямой и обратной подстановок согласно выражениям (8.16) и (8.17). В результате оценка величины En+1 по выбранной εn+1 (и наобо­рот) не требует ни коренной переделки программы решения, ни чрезвычайно большего времени работы процессора.

Найденное по уравнению (8.84) значение общей погрешности En+1 используется наряду с максимально допустимой ее величиной Emax для управляемого выбора допустимого значения локальной погрешности

(8.85)

где δ = 0,5 (0 < δ < 1) — коэффициент надежности, при этом вы­бирают εт, 0 = 0,1Еmах и ограничивают εт, min≠0.

Выражение (8.85) можно получить, если учесть из уравнения (8.84) ропорциональность Еn+1 и εn+1, а из выражения (8.76) зависимость h и εт, т. е. записать

(8.86)

Стратегии выбора шага (порядка). Подставляя формулу (8.76) в выражение (8.81), находим взаимосвязь максимально допусти­мой погрешности на единичном шаге вычислений ε с максимальным размером шага h:

(8.87)

которая изображена на рис. 8.3 для k = 1, 2, 3, и 4.

Рис. 8.3. Кривые зависимости допустимой погрешности на единичном шаге вычислений от из­менений порядка k алгоритма

Из рисунка следует, что при ε>εа алгоритм первого порядка (неявный ме­тод Эйлера) обеспечивает наибольший допустимый шаг вычисле­ний, при εb< ε < εа — алгоритм второго порядка (метод Ших-мана), при εс< ε < εb — алгоритм третьего порядка и при ε < εс — алгоритм четвертого и более высоких порядков. Так как положение кривых на рис. 8.3 зависит от текущих значений Ck+1 и x(k+1)(tn+1), то даже при выбранном постоянном ε оптимальные порядок k алгоритма и шаг вычислений h могут изменяться от од­ного временного шага к другому. Как следует из рис. 8.3, любое увеличение допустимой погрешности ε не приведет к увеличению порядка формулы для максимизации шага вычислений и любое уменьшение ε не может привести к уменьшению оптимального порядка k используемого метода. Для достижения максимальной вычислительной эффективности целесообразно в процессе интегри­рования одновременно изменять порядок алгоритмов и шаг с помощью многошаговых алгорит­мов.

С точки зрения программирова­ния изменение порядка алгоритма эквивалентно смене набора коэффи­циентов, определяющих ту или иную многошаговую формулу. При увели­чении порядка k алгоритма увели­чивается число коэффициентов в наборе. Это сопровождается ростом загрузки памяти ЭВМ для хранения «прошлых» значений вектора пере­менных, используемых для оценки приближения формулами высоких порядков. Как правило, указанные дополнения не требуют существен­ных затрат. Наоборот, изменение шага связано со значительными вычислительными затратами, так как часто хранимые «прошлые» значения вектора переменных, соответствующие шагу hj-1 не­обходимо интерполировать для получения преобразованных «прошлых» значений, соответствующих новому шагу вычислений hj. Процедура интерполяции упрощается, если «прошлые» значе­ния хранить в виде вектора обратных разностей

или в виде вектора Нордсика

старые и новые значения которых связаны матрицей, НЭ которой зависят только от коэффициента q = hj/hj-1.

Успешно используются разные стратегии изменения по­рядка и шага методов интегрирования. В соответствии с одной из них, предложенной Гиром, выбор шага осуществляется по формуле

(8.88)

получаемой на основании выражений (8.87) и (8.81).

На основании одного из выражений (8.76) или (8.80) оценивает­ся погрешность вычисления каждой составляющей вектора х сразу для методов трех ближайших порядков

(8.89)

Далее выбираются погрешности, соответствующие наименее точным составляющим вектора переменных х:

Порядок метода после каждых k шагов вычислений выбирается равным порядку, при котором величина еi минимальна (i = 1, 2, 3). Новое значение шага h определяется в соответствии с выраже­нием (8.88). При этом шаг вычислений при неизменном порядке метода корректируется не сразу, а только при достижении опре­деленных соотношений между действительной и заданной погреш­ностями вычислений между используемым и разрешенным значе­ниями шага интегрирования. Напомним, что в соответствии с табл. 8.1 частые и большие изменения шага не разрешаются, при этом