ЖЕСТКАЯ ДИФФЕРЕНЦИАЛЬНАЯ СИСТЕМА

- система обыкновенных дифференциальных уравнений, при численном решении к-рой явными методами типа Рунге - Кутта или Адамса, несмотря на медленное изменение искомых переменных, шаг интегрирования обязан оставаться малым. Попытки уменьшить время вычисления решения Ж. д. с. за счет увеличения шага интегрирования приводят к резкому возрастанию погрешности (взрыв погрешности).

Автономная система обыкновенных дифференциальных уравнений

наз. Ж. д. с, если для любых начальных значений на заданном отрезке [0, Т], принадлежащем интервалу существования решения (1), выполнены условия:

а) максимальный модуль собственных значений матрицы Якоби (спектральный радиус) ограничен вдоль решения z(t):

б) существуют такие числа t _п, N,v, что при

справедливо неравенство

здесь

X(t) - фундаментальная матрица уравнения в вариациях для системы (1),

t _п - длина пограничного слоя. К Ж. д. с. относятся все системы вида (1), для к-рых условия а) и б) выполняются совместно после масштабирования компонент вектора z(t)на каждом решении.

Неавтономная нормальная система обыкновенных дифференциальных уравнений порядка тназ. Ж. д. с, если жесткой является равносильная ей автономная система порядка m+1. Примером жесткого неавтономного уравнения может служить скалярное уравнение:

j(х)О С ^p, -заданная функция. Другой пример Ж. д. с.-

линейная однородная система с постоянной (т m)-матрицей А, имеющей различные собственные значения, разделенные на две группы:

Условие а) для системы (3) очевидно выполняется. Норма ограничена сверху величиной при где С ₁, С ₂ - некоторые постоянные. При N=L/(C₁be^aT). выполняется условие б). Вектор d²z(t)/dt² при ограничен сверху:

Поэтому для вектора dz(t)/dt приможно построить интерполяционный алгебраич. многочлен нулевой степени с равномерно распределенными узлами и выбрать шаг интерполяции в виде Н=С ₃/С ₁bе ^{a Т}, где С ₃- некоторая постоянная (см. [1]), зависящая от заданной погрешности. С другой стороны, использование метода ломаных Эйлера для системы (3) требует ограничения сверху на шаг интегрирования на всем отрезке решения системы (см. [1]):

При этом составляющие приближенного решения системы (3) методом ломаных Эйлера, соответствующие первой группе собственных значений (4), будут представлены с достаточной точностью (см. [1]):

Ограничение вида (5) на шаг интегрирования характерно для экстраполяционных методов типа Рунге - Кутта и Адамса. Отношение H/h, к-рое может рассматриваться как качественная мера жесткости системы (3), достигает во многих случаях величин порядка 10⁴-10¹⁰. Математич. описание динамич. процессов и явлений в физике, химии, биологии, технике и экономике, связанное с учетом все большего числа факторов, повышающих порядок дифференциальных систем, приводит к увеличению жесткости. Ж. д. с. требуют специальных методов решения.

В некоторых случаях исходную систему (1) можно преобразовать, используя теорию и асимптотич. методы для дифференциальных уравнений с малым параметром при старшей производной, поведение решения к-рых в пограничном слое описывается экспоненциально затухающими функциями (см. [2]). Однако Ж. д. с. обычно трудно привести к подобному виду, и, кроме того, после асимптотич. преобразования жесткость не всегда существенно уменьшается. Поэтому для решения систем общего вида используются численные методы, основным свойством к-рых является обеспечение подавления быстро затухающих составляющих решения системы (1) вне пограничного слоя при величине применяемого шага интегрирования, близкой к величине шага алгебраич. интерполяции.

Для решения Ж. д. с. целесообразно использование неявных методов (см. [9], [17]). Для построения неявных схем можно использовать способ неопределенных коэффициентов на основе формулы Тейлора, записанной в виде:

где п- целое положительное число, H>0, nH=t_n. Напр., такие методы описаны в [9]:

r=1, а ₀=1, а ₁ =-1 - неявный метод ломаных; r=2, а ₀= a₁=-2, а ₂= -неявный метод 2-й степени; r=3, а ₀= a₁=-3, a₂= a₃=- неявный метод 3-й степени; r=4, а ₀= a₁=-4, а ₂ = 3, а ₃= а ₄= -неявный метод 4-й степени.

Под степенью метода понимается наибольшая степень Н ^к в разложении (7) по степеням Н, множитель при к-рой совпадает с соответствующим множителем в (6).

Применение неявного метода ломаных к системе (3) приводит к разностным уравнениям

Пусть система (3) асимптотически устойчива по Ляпунову. Тогда матрица ||Е-НА|| - неособенная при любом Н. При использовании формулы Лагранжа - Сильвестра решение (8) представляется в виде

Для неявного метода ломаных (8) условие асимптотич. устойчивости выполняется при любом Н, и в (9), с ростом псоставляющие решения, соответствующие второй груши собственных значений (4), будут быстро убывать по модулю. Величина Нограничивается только условиями требуемой точности приближенных решений. Стремление повысить степень линейных многошаговых методов вступает в определенное противоречие с их устойчивостью (см. [11]). r-Шаговый метод

наз. А-устойчивым, если все решения (10) стремятся к нулю при с фиксированным положительным Ив случае применения (10) к скалярному уравнению

где q- комплексная постоянная с отрицательной действительной частью. Явный r-шаговый метод не может быть А -устойчивым. Степень линейного многошагового А-устойчивого метода не может превышать двух. Наименьшая константа ошибки c_d=получается для неявного метода трапеций (см. [11]):

Требование А(a)-устойчивости накладывает менее жесткие ограничения на методы в отличие от А-устойчивости. r-Шаговый метод (10) при фиксированном H>0 наз. А(a) -устойчивым, если все решения (10) стремятся к нулю при в случае применения (10) к уравнению (11), q- комплексная постоянная, к-рая принадлежит множеству

С - комплексная плоскость. Явный r-шаговый метод не может быть А(О)-устойчивым. Существует только один А(О)-устойчивый r-шаговый метод со степенью равной r+l - метод (12). Для всех существуют А(a)-устойчивые методы со степенью, равной rпри r=3, 4 (см. [13]).

Следующее ослабление требований A-устойчивости содержится в так наз. определении Гира (см. [14]): метод (10) наз. жестко устойчивым, если все решения (10) стремятся к нулю при в случае применения (10) к уравнению (11). Hq- комплексная постоянная, принадлежащая множеству где

причем для обеспечивается заданная точность.

Методы (7) жестко устойчивы и, следовательно, А(a)-устойчивы. Для них

Можно построить согласно (6) неявные аналоги явных методов Рунге - Кутта любой степени, обладающие A-устойчивостью и жесткой устойчивостью. Напр., метод 2-й степени

Применение метода (13) к системе (3) приводит к разностным уравнениям

которые доказывают его А-устойчивость. Для метода 3-й степени:

аналогично получаются разностные уравнения:

Так же строятся А-устойчивые и жестко устойчивые методы более высоких степеней. Методы (13) и (14) и методы, опубликованные в [5], принципиально отличаются от так наз. неявных методов Рунге - Кутта (см. [16]), не получивших распространения из-за большой вычислительной трудоемкости. Применение неявных методов требует большей вычислительной работы на одном шаге, чем применение явных методов, но это оправдывается в Ж. д. с. за счет резкого увеличения шага интегрирования. При решении задачи (3) необходимо обращать матрицу или решать систему линейных уравнений. Здесь может возникнуть проблема плохой обусловленности, так как число обусловленности матрицы ||Е- НА|| увеличивается с увеличением H. В общем случае (1) на каждом шаге интегрирования необходимо решать нелинейную систему уравнений относительно вектора y_n+1. Обычно при этом пользуются модифицированным методом Ньютона с начальным условием, вычисленным по любой экстраполяционной формуле, использующей вектор у _п. Метод экстраполирования при этом наз. предиктором, а неявный метод - корректором. Метод простой итерации в Ж. д. с. неприменим из-за больших значений LH. Поскольку в неявных схемах применяется метод Ньютона для решения уравнений F(y_n+1)=0 и, следовательно, необходимо вычислить матрицу Якоби системы (1), то иногда включают эту матрицу непосредственно в формулы методов, к-рые при решении линейных систем также обладают свойством А- устойчивости (см. [12], [15]). Для решения Ж. д. с. широко применяется процедура Гира (метод Гира) с автоматич. контролем погрешности на шаге и изменением по этому критерию как степени методов, так и шага интегрирования [14]. В процедуре Гира в качестве корректора используются и методы (7) (см. [9]). Другой подход к созданию методов интегрирования жестких систем уравнений связан с учетом в формулах методов решений соответствующих линейных систем (см. [4] - [8], [10]). В первых работах этого направления рассмотрены жесткие системы уравнений с известными собственными значениями матрицы дf(z)/дz, по к-рым строились матричные коэффициенты методов. В связи с трудностью решения проблемы собственных значений это направление длительное время не развивалось. В [6], [7], [8] предложены способы построения: матричных коэффициентов без решения проблемы собственных значений для матрицы дf(z)/дz системы (1), Методы этого направления можно строить на основеравенства (см. [7]):

где - неособенная (т m)-матрица, С- матрица, не зависящая от т.

Различный выбор матриц j(t_n+t) и Сприводит к разностным уравнениям, соответствующим тем или иным методам численного интегрирования, если пренебречь правой частью в (15). При С=0 получаются явные методы, а при - неявные. Пусть j(t_n+t)=Е. Тогда неявный метод ломаных (7) соответствует С=-НЕ;. метод трапеций (12) явный метод ломаных С=0. При j(t_n+t)=e^At, где А- постоянная (т m)-матрица, получается обобщенный метод ломаных (см. [7]):

Формула (16) переходит в явный метод ломаных при А=0. Метод (16) дает точное решение системы уравнений

для дискретных значений t_n=nH. Если известна матрица

то использование краз рекуррентной формулы

приводит к матрице

применяемой в (16). В качестве начального приближения для (18) при достаточно малом целесообразно употреблять приближенную формулу:

или, если собственные значения матрицы действительные, то

Формула (19) дает соответствие по устойчивости в смысле Ляпунова решений системы дифференциальных (17) и разностных (16) уравнений в случае комплексных собственных значений Ас малой действительной частью. Если область существования решений является замкнутой, выпуклой по z, и приближенные решения принадлежат той же области, то погрешность метода (16) удовлетворяет разностному уравнению (см.[7]):

где z_n=z_n-y_n, a z_n=z(t_n), y_n=y(t_n)- решения (1) и (16)соответственно. Пусть -норма вектора (норма матриц подчинена принятой норме вектора) и в области G: .

По матрице Ас действительными элементами вычисляется число

Тогда справедлива оценка

Если погрешность e_n можно оценивать, считая В=0. В оценках возможны другие нормы векторов, соответствующие им нормы матриц и логарифмич. нормы (см. [3]). Приведенные оценки показывают, что при решении системы (1) шаг интегрирования Нможно использовать значительно большим, чем в классич. методах. Матрицу Аследует выбирать близкой по всем элементам к матрице Якоби системы (1). В пограничном слое, когда переменные изменяются быстро, грубо оценивая m₁, m₂, l, R по приближенному решению, можно менять матрицу А, добиваясь необходимой точности. Так как за пограничным слоем переменные системы (1) меняются медленно, на практике часто оказывается достаточно одной матрицы А, чтобы вычислить все решение при Проверку точности можно производить по правилу Рунге (см. [1]).

С целью увеличения точности на основе метода (16) в [7] предложен класс системных методов численного интегрирования. Для методов этого класса выполняются требования: 1) системный метод s-й степени должен быть точным при интегрировании алгебраич. многочленов степени s-1 при A=0; 2) системный метод любой степени должен быть точным при решении уравнений (17).

Одношаговый метод 2-й степени имеет вид:

Там же построены явный одношаговый метод 3-й степени и многошаговые системные методы. Даны асимптотич. оценки их погрешностей. При A=0 эти методы переходят в формулы типа Рунге - Кутта или Адамса. При аппроксимации интегралов от е ^At дробно рациональными матричными полиномами от A с учетом неособенности соответствующих матриц системные методы переходят в явную форму соответствующих неявных методов, получающихся после итераций по методу Ньютона. Неявным системным методом 1-й степени является метод

Он получается из (15), если выбрать

Уравнения (21) получаются способом приведения дифференциальных уравнений к интегральным (см. [8]). Интеграл от е ^Аt в методе (21) вычисляется по формулам (18) с начальным условием (20). Исследован способ коррекции этого интеграла по ходу решения (см. [8]).

Лит.:[1] Бахвалов Н. С, Численные методы, 2 изд.. М., 1975; [2] Васильева А. В., "Матем. сб.", 1960, т. 50, в. 1, с. 43-58; [3] Былов Б. Ф. и др., Теория показателей Ляпунова..., М., 1966; [4] Гавурин М. К., в сб.: Методы вычислений, 1963, в. I, с. 45-51; [5] Ракитский Ю. В., "Докл. АН СССР", 1970, т. 193, № 1, с. 40 - 42; [6] его же, там же, 1972, т. 207, № 4, с. 793-95; [7] его же, "Тр. Ленингр. политехи, ин-та", 1973, № 332, с. 88-97; [8] Павлов Б. В., Повзнер А, Я., "Ж. вычисл. матем. и матем. физ.", 1973, т. 13, №4, с. 1056-59; [9]Curtiss С. F.,Hirschfelder J. О., "Ргос. Nat. Acad. Sci. USA", 1952, v. 38, p. 235-43; [10] Mah R. S. H., Michaelson S., Sargent R. W., "Chem. Engng Sci.", 1962, v. 17, p. 619-39; [11] Dahlquist G., "Nord. tidskr. Informationsbehandling", 1963, b. 3, s. 27-43; [12] Rosenbrock H. H., "Comput. J.", 1963, v. 5, p. 329- 30; [13] Widlund О. В., "Nord. tidskr. Informationsbehandling", 1967, b. 7, s. 65-70; [14] Gear C. W., в кн.: Information Processing 68, v. 1, Amst., 1969, p. 187-93; [15] Lambert J. D., Sigurdsson S. Т., "SIAM J.Numer Anal.", 1972, y. 9, p. 715-33; [16] Lambert J. D., Computational methods in ordinary differential equations, N.Y., 1973; [17] Stiff differential systems, [The IBM research simposia series], N.Y.- K., 1974.

Ю. В. Ракитский.