Рассмотрим особенности реализации неявных многошаговых тодов. Перепишем схему (1.55) в виде

где величина gj зависит от уже найденных значений ш, ui-Hi известна.

Для решения уравнения (1.56) применяют два метода: метод пр стой итерации и метод Ньютона. Рассмотрим первый метод, g этом случае находим с помощью итерационного процесса

/(тт. {у-) + еи 5=1. 2, (1.57,

(«/+l)

где {ui-Y ~ начальное приближение к искомому значению функ. ции.

Для быстрой сходимости итерационного процесса нужно удачно выбрать начальное приближение [u-f. Обычно его получают с помощью явной многошаговой формулы. Тогда в целом алгоритм расчета на каждом шаге выглядит так; сначала предсказывается значение (ul-f, а затем проводится одно или несколько уточнений начального значения по формуле (1.57), т. е. получается предсказываю-ще-исправляющий метод, который часто называют методом предиктор-корректор.

Используют два способа реализации этого метода. В первом случае итерации проводят до тех пор, пока два последовательно полученных значения (и/+)*\ {ui-f не станут достаточно близки. Тогда число итераций может меняться от шага к шагу. Во втором случае на каждом шаге исправляющая формула применяется фиксированное число раз. Обычно на практике метод предиктор-корректор применяется с числом итераций не более трех.

Наиболее часто из линейных /г-шаговых схем используют схемы (1.49), (I.5I), называемые схемами Адамса. Можно доказать, что явная схема Адамса имеет порядок аппроксимации равный k, а неявная - (k + I). При использовании метода предиктор-корректор обычно применяют предсказывающую и исправляющую схемы одного порядка точности. В частности, широко применяется метод предиктор-корректор со схемами Адамса четвертого порядка, в котором предсказание делается по формуле (1-52), а уточнение - по (1.53)-

Выбор шага интегрирования и оценка погрешности численного решения. Обычно при реализации на ЭВМ большинства численных схем решения обыкновенных дифференциальных уравнений предусматривается автоматический выбор величины шага Дт для обеспечения определенной погрешности расчета. Этот выбор основан на оценке локальной погрешности численного решения на шаге, т- б-погрешности численного решения в точке Tj-i, оцениваемой в предположении, что в начале шага в момент времени ту значение искомой функции было известно точно.

Япя одношаговых методов существуют два основных способа и локальной погрешности на шаге. При первом - интегриро-0 от точки Tj до точки Tyj осуществляется по схеме одного и то-вани JJдJa точности р дважды: с шагом Дт и с шагом Дт/2. Можно 01$зать, что в этом случае для локальной погрешности e/+J реше-найденного с шагом Дт/2, справедлива следующая приближенная оценка:

е/+1 = [«/ +1 (Дт/2) -(Дт)]/(2Р- 1), (1.58)

/+1 (Дт/2), и/+(Дт) - решения в точке найденные на

нове одного и того же решения Ш в точке ту с шагами Дт/2 и Дт соответственно.

При втором способе интегрирование на отрезке 1ту, ту! также выполняется дважды: с помощью схемы порядка точности р и схемы порядка точности (/7 + 1) при одном и том же шаге Дт, и находятся соответственно два решения в точке ту: (р) и (р + 1). В этом случае локальная погрешность оценивается по формуле

8/+I лок

ui+ {р) - и + {р + \).

(1.59)

Обычно при расчетах указывается допустимое значение абсолютной нли относительной локальной погрешности, в соответствии с которым автоматически подбирается величина шага. Например, при первом способе оценки погрешности это делается так. В ходе расчетов сначала из некоторой точки ij делается шаг Дт, затем два шага длиной Дт/2, и на основе двух решений (Дт) и и (Дт/2) оценивается погрешность e+J в точке ту по формуле (1.58). Если погрешность не превосходит допустимого значения, то ui (Дт/2) принимается в качестве найденного значения сеточной функции, и из точки Tj-i продолжается расчет с шагами Дт и Дт/2. Если оценка локальной погрешности превосходит допустимое значение, то из точки Ту делают два шага длиной Дт/4 и оценивают погрешность в точке Ту-- Дт/2. При достижении требуемого уровня погрешности расчет из точки ту -- Дт/2 продолжается с шагами Дт/2 и Дт/4. В противном случае проводится дальнейшее измельчение шага. В том случае, когда при первой проверке погрешность оказывается существенно ниже допустимого значения, обычно шаг увеличивают путем его удвоения.

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

При использовании многошаговых схем по методу предиктор- Рректор оценку локальной погрешности иа шаге обычно удается

выразить через вычисляемые прн реализации схемы предсказываемое и исправленное значения искомой функции. Таким образом, в этом случае для оцеикн погрешности не требуется дополнительных вычнс-леинй функции / (т, и), что является преимуществом по сравнению с одношаговыми методами. Однако у многошаговых методов сложнее реализуется процедура изменения шага. При уменьшении шага необходимо либо путем интерполяции определить значения искомой функции в новых предшествующих точках, либо получать по одно-шаговой схеме новые «разгонные» точки.

При любой стратегии, основанной на оценке локальной погрешности, не учитывается накопление погрешности в ходе всего расчета н, следовательно, фактическая погрешность численного реи1еиия остается иензвестной.

Более надежный путь оцеикн погрешности - сравнение двух решений (и/У и («)" полученных во всей расчетной области 10, tnaJ прн двух сетках разной густоты - с шагами Дт и Дт". В этом случае полная погрешность е = Т - оценивается по следующей формуле;

8/-И1/[(4т/АЛ-и. (1-60>

где р - порядок точности схемы. Наиболее часто сравнивают решения иа сетках, шаги которых отличаются вдвое.

Подчеркнем, что несмотря па внешнее сходство с оценкой (1.58) для локальной погрешности, соотношение (1.60) справедливо для полной погреншостн численного решения, так как сравниваются решения, полученные путем расчета всего процесса от начального момента Tft " О до момента т на двух разных сетках.

Рассмотренный способ оценки полной погрешности не удается использовать прн автоматическом выборе шага в процессе счета иа основе Оценок локальной погрептостн.

Решение систем обыкновенных дифференциальных уравнений. Большинство рассмотренных выше схем решения одного дифференциального уравнения первого порядка может быть легко обобщено для решения системы уравнений:

АТг

(1.61)

Разностные схемы записывают для системы уравнений так же, как и для одного уравнения. Отличие заключается лишь в том,, что-вместо одного значения сеточных функций ш н f (Xj, ui) в кждр узле вычисляют векторные сеточные функции с Л компонентами:

W = \u{, ut]. Ji=\f{.....fl

Де fi - fi (Tj, ..... u).

Например, явная схема Эйлера записывается для системы двух уравнений в виде

и{, иО,

K+-i -f Дт/з(ту, и{, и{у

Таким образом, запись какой-либо схемы для системы обыкновенных дифференциальных уравнений ие вызывает дополнительных затруднений. Сложности могут возникнуть при ее численной реализации. Они обусловлены двумя обстоятельствами. Первое связано с необходимостью при применении неявных схем решения на каждом шаге систем алгебраических уравнений. Этот вопрос рассматривается в § 1.6 на примере неявной схемы Эйлера для системы уравнений теплового баланса.

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

Жесткие системы уравнений. Поясним понятие жестких систем и специфические трудности нх решения на примере все той же задачи расчета нестационарного теплового режима системы тел. Пусть имеется следующая задача расчета нагрева двух тел внутренними источниками:

ill.

Точное аналитическое решение этой задачи имеет вид « (т) = ехр iih т) - D.n ехр (цзт) + ТГ, п = \,2, Где Тгаи стационарный перегрев над средой «-го тела; гп - постоянные, определяемые нз начальных условий; (д.], (,ц - собственные числа задачи, равные

- + OieVCi, m, = (Oi, -b ОзеУС, - параметры, имеющие

смысл темпов охлаждения каждого нз тел.

Из структуры зависимости для ti,-: видно, что может возникнуть ситуация, когда отрицательные показатели экспонент 1Л1 н \iz сильно отличаются по абсолютной величине, т. е. рИ <С Это неравенство будет выполнено, например, если т. mi, т. е. тело 2 обладает су1дественно меньшей тепловой ииерцнен; или в случае т -- тз, 1 = Cj. сг, н (Tic <С cTia (тепловая связь между телами на-

ого сильнее связи со средой), тогда [ii = -OjJCi, \i.2 = -(cfla +

При условии 1 < llAgl в поведении решения Т„ (т) можно выделить два временных интервала: начальный и основной, причем длительность первого значительно меньше, чем длительность второго. На начальном интервале происходит быстрое изменение второго слагаемого в решении, имеющего экспоненциальный множитель с большим по абсолютной величине показателем. К началу основного этапа «быстрая» экспонента практически полностью затухает и процесс определяется первым слагаемым с большой постоянной времени.

Проблема численного решения системы уравнений с сильно различающимися собственными числами (х„ возникает в том случае, когда необходимо получить искомые функции на временном интервале, длительность которого значительно превышает начальный этап. Действительно, При использовании какой-либо из рассмотренных выше явных схем илн схем прогноза-коррекции нельзя выбирать шаг по времени большим, чем это допустимо по условию устойчивости. Типичное ограничение на шаг по времени для решения системы уравнений имеет вид

где постоянная b зависит от используемой схемы, а ртах - максимальное по модулю собственное число.

При реализации чпсленой схемы на начальном этане процесса ограничение на шаг Дт не является обременительным, поскольку для расчета быстрого процесса, определяемого членом ехр ([Хп,н\Т). из условия получения требуемой погрешности расчета значение uia-га Ат, как правило, все равно должно быть меньше, чем Дтща;?:- Неудобства возникают после выхода на основную стадию. «Быстрый экспоненциальный множитель в точном решении затухает и возникает потребность увеличения шага по времени для «отслеживания» медленно меняющегося процесса. Однако из-за свойств разностной схемы шаг увеличивать нельзя, поскольку сразу же начинает развиваться неустойчивость. В результате вся длительная основная стадия процесса рассчитывается с малым шагом по времени. Это приводит к недопустимому увеличению затрат машинного времени и накоплению погрешностей округления, которое может существенно исказить окончательные результаты.

Рассмотренный пример позволяет отметить две особенности жестких систем. Во-первых, в решении таких систем входят два вида функций: бьЕстроубывающне функции, которые полностью «затухают» в начальной стадии, и медленно изменяющиеся функции, которые определяют процесс на основной стадии с длительностью, намного превышающей начальную. Во-вторых, интервал, на котором требуется получить решение, значительно превышает начальную стадик»;

Таким образом, жесткость систем определяется как структуро решения, так и интервалом, на котором необходимо его получить-Строгое математическое определение понятия жесткости дано в кий1 гах [22. 29].

Нами был рассмотрен важный и довольно типичный, но отнюдь единственный пример жестких систем. Из этого примера вытекает "нбоее распространенный критерий, позволяющий выявить жест-ю систему в случае, когда ее собственные числа вещественные и отрицательные.

Кпитерий заключается в сопоставлении по модулю максимального Umax И минимального собствснных чисел линейной или линеа-пизованной системы. Если этн собственные числа сильно различаются по абсолютной величине (pniaxPmiu > 0. э промежуток времени, на котором нужно получить решение, значительно превосходит интервал затухания ехр (ртах), то система является жесткой и при ее численном решении по явным схемам возникнут проблемы с выбором шага.

Для решения жестких систем применяют специальные неявные разностные схемы, при реализации которых возникает ряд трудностей. Во-первых, при решении систем нелинейных разностных уравнений применение метода простой итерации неэффективно. Необходимо применять метод Ньютона. Во-вторых, линейные системы разностных уравнений, получающиеся либо непосредственно в процессе peuieHMH линейных систем обыкновенных ди()х})еренциальных уравнений, либо при реализации метода Ньютона в случае нелинейных систем, являются часто плохо обусловленными. Для их решения применяют специальные приемы.

В заключение отметим, что здесь мы лишь обратили внимание на существование класса жестких систем, прн численном решении которых возникают трудности, и дали весьма не строгий качественный анализ причин этих трудностей. Более подробно с методами решения жестких систем можно ознакомиться в книгах [22, 29].

Программное обеспечение решения систем уравнений. Для численного решения задачи Коши для обыкновенных дифференциальных уравнений и систем таких уравнений имеется достаточно большое число стандартных подпрограмм, реализующих различные одноша-говые и многошаговые методы П5]. При применении этих подпрограмм 1Юльзователь должен составить подпрограмму, в которой производится вычисление правых частей конкретной системы уравнений, а также Организовать вывод результатов - значений искомых функ-ии при интересующих значениях аргумента Ту. Особенности использования стандартных подпрограмм разберем на примере под-[Рограммы R KGS нз математического обеспечения ЕС ЭВМ, ко-рая реализует схему Рунге-Кутта четвертого порядка для системы "новенных дифференциальных уравнений с автоматическим lOopoM шага интегрирования. Пример применения этой подпрограммы приведен в следующем параграфе для решения задачи рас-

знестационарцого теплового режима системы тел.

-ращение к подпрограмме RKGS имеет внд

*LL RKGS(PRMT, Т, DERY, THLF, FCT, OUTP, AUX).