Портал | Содержание | О нас | Авторам | Новости | Первая десятка | Дискуссионный клуб | Чат Научный форум --> Первая десятка "Русского переплета" Темы дня:

Президенту Путину о создании Института Истории Русского Народа. |Нас посетило 40 млн. человек | Чем занимались русские 4000 лет назад?

| Кому давать гранты или сколько в России молодых ученых?

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

ПРИ ЧИСЛЕННОМ РЕШЕНИИ ЭВОЛЮЦИОННЫХ ЗАДАЧ

В. В. БОБКОВ

Белорусский государственный университет, Минск

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

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

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

u'(t) = lu(t),

связывающее во времени t искомую функцию u(t) с ее производной u'(t) посредством равенства (1) при заданном числовом множителе l ? 0. Непосредственной подстановкой легко убедиться, что решением такого дифференциального уравнения является функция вида

u(t) = e ltc

с произвольной постоянной c, конкретное значение которой, очевидно, может быть задано постановкой дополнительно к (1), скажем, начального условия u(t0) = u0 , фиксирующего значение неизвестной функции u(t) в момент времени t = t0 . Экспоненциальная функция e lt при этом может быть определена, например, посредством равенства

с бесконечно большим числом слагаемых в правой части, где k ! = 1 " 2 " _ " k. Не останавливаясь здесь (в силу вспомогательной роли подобного материала для основной цели нашего разговора) на доказательстве известного факта, что сумма в правой части (3) имеет конечное значение при любом фиксированном lt, и на других замечательных свойствах этой функции, отметим лишь некоторые из них, используемые в наших рассуждениях и почти очевидные:

e 0 = 1, e l(t + t) = e lt " e lt, y'(t) = ly(t),

где y(t) = e lt. Неплохое приближение для значения этой функции при фиксированном lt можно получить сохранив в правой части (3) достаточно большое число первых слагаемых. К примеру,

Представление о графике такой показательной функции e x (x = lt, e © 2,7) школьники уже имеют (рис. 1).

Полезным для наших целей примером может быть также обобщающее (1) линейное неоднородное уравнение

u'(t) = lu(t) + a(t),

где a(t) определяется по заданной функции j(t) посредством равенства

a(t) = j'(t) - lj(t).

С учетом (5) уравнение (4) приводится к однородному виду

n'(t) = ln(t),

где n(t) = u(t) - j(t). Принимая во внимание (1), (2), в качестве решения уравнения (6) можно взять функцию n(t) = e ltc, что позволяет считать решением линейного дифференциального уравнения (4) с неоднородностью (5) функцию вида

u(t) = e ltc + j(t)

с произвольной постоянной c. Так как на скорости изменения функции y(t) = e lt существенно сказывается значение числового параметра l, а выбор функционального параметра j(t) не привязан к l, то соответствующим их заданием можно сделать оба слагаемых суммы (7) изменяющимися на некоторых отрезках времени в разных временных масштабах. Понятны также роль и значимость при этом константы c, связанной с постановкой начального условия для уравнения (4), (5). Существенная разномасштабность составляющих решения (7) может быть, очевидно, причиной многих трудностей при его приближенном представлении.

Линейные уравнения (1), (4) являются частными случаями эволюционного уравнения вида

u'(t) = f (t, u(t)),

правая часть которого задается известной функцией двух переменных: независимой t и зависимой u(t). В общем случае нелинейной функции f решение уравнения (8) в отличие от (1), (4) не удается уже выписать явно. Тем более, естественно, не удается точно решить эту задачу в общем случае системы уравнений

которую часто записывают в виде (8), понимая при этом под u(t) векторную функцию с компонентами ui(t), i = 1, 2, _, n. Аналогичный смысл придается и функциям u'(t), f (t, u(t)). Еще более трудной подобная задача будет, конечно, в случае таких эволюционных уравнений, решения которых являются функциями не только временнЧй, но и пространственных переменных. Правые части уравнений типа (9) (или в более компактной форме типа (8)) в этом случае могут зависеть также от пространственных переменных и производных решения по этим переменным. Использование такого рода сложных эволюционных уравнений в приложениях связано, как правило, с разработкой способов их приближенного решения, ориентированных на современную вычислительную технику. Многие конструктивные недостатки таких численных методов могут быть выявлены уже при их применении даже в случае простых уравнений вида (1) или (4), (5), решения которых нам известны.

Аналогичную (1) роль для случая системы эволюционных уравнений может играть, например, система дифференциальных уравнений вида

которую, подобно (8), можно записать в компактной форме

u'(t) = Au(t),

понимая под матрицей A таблицу чисел

Искомая векторная функция u(t) при этом понимается в виде таблицы-столбца соответствующей размерности:

Правило умножения матрицы A справа на вектор u(t) в (11) очевидным образом может быть восстановлено сопоставлением (10) и (11).

Подобное же обобщение на случай системы легко может быть проведено, конечно, и для уравнения вида (4), (5).

По аналогии с (1), (2) решение системы (11) может быть записано в форме

u(t) = e Atc

с произвольным вектором c размерности n, при этом матрица e At определяется подобно (3):

Здесь буквой I обозначена единичная матрица вида (12) соответствующей размерности: ai j = 0, если i ? j, aii = 1. Под умножением матрицы на число понимается умножение на это число каждого элемента таблицы (12); столбцы матрицы произведения двух матриц определяются согласно упомянутому выше правилу умножения матрицы (12) на вектор-столбец типа (13) для каждого из соответствующих (по номеру) столбцов правой матрицы в этом произведении; при сложении матриц складываются их соответствующие элементы.

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

Сначала рассмотрим вопрос о том, как по известному в момент времени t значению u = u(t) решения уравнения (или системы) (8) можно найти (точно или приближенно) значение этого решения в момент времени t + t, где t > 0. В общем случае (8) приходится удовлетвориться, как уже отмечалось, лишь некоторым приближением . Наиболее простой способ нахождения дает, например, метод Эйлера

где y © u, построение которого может быть, очевидно, сведено к приближенной замене в (8) в момент времени t значения производной u'(t) простейшим разностным отношением

(u(t + t) - u(t))/ t.

Применительно к (1) метод Эйлера приводит к расчетным формулам вида

где

E (lt) = 1 + lt.

В соответствии с (2) взамен (18) мы должны были бы иметь

И аналогично в случае системы (11) вместо (18), (19), (20) можно записать (с прежним смыслом обозначений):

E (At) = I + At,

Заметим прежде всего, что использование расчетных формул (21), (22) (взамен точного соотношения (23)) в случае системы вида (11) по сравнению со случаем одного уравнения (1) сопряжено не только с порождаемым увеличением количества уравнений ростом объема вычислительной работы, но может приводить и к принципиальным отличиям в качестве приближений, обусловленным возможной разномасштабностью составляющих решения этой системы. При этом если в случае одного линейного уравнения эффект разномасштабности мог возникать (см. (7)) лишь за счет неоднородности (скажем, вида (5)), то для системы эволюционных уравнений подобный эффект оказывается возможным даже в однородном случае (11). И если разномасштабность в (7) легко прогнозировалась по виду исходного уравнения (4), (5), то для случая системы (11) подобная разномасштабность в (14) носит, вообще говоря, скрытый характер. Чтобы составить представление об источниках такой скрытой (и поэтому труднее учитываемой при аппроксимации) разномасштабности, будем рассматривать интересующий нас n-мерный вектор решения не в привычной координатной форме, а запишем его в специальном виде, привязанном к данной матрице A.

Любой вектор u в координатной форме записи (13), очевидно, может быть разложен в сумму вида

u = u1l 1 + u2l 2 + _ + unl n,

где вектор-столбец l i, i = 1, 2, _, n, имеет нулевые компоненты, за исключением компоненты номера i, равной единице. В случае диагональной матрицы A (ai j = 0 в таблице (12), если i ? j) операция Au с учетом (24) приводит, как легко видеть, к равенству

Au = u1a11l 1 + u2a22l 2 + _ + unannl n.

И аналогично (см. (3), (15), (22), (23)):

E (At)u = u1(1 + a11t)l 1 + u2(1 + a22t)l 2 + _

_ + un(1 + annt)l n,

В случае же произвольной матрицы A вида (12), подобной (25)-(27), простой картины мы не имеем. Чтобы получить близкую к ней по наглядности ситуацию, используем вместо (24) разложение вектора u не по универсальной координатной системе базовых векторов l i, i = 1, 2, _, n, а по специальной системе векторов xi, i = 1, 2, _, n, тесно привязанной к данной матрице A и обладающей аналогичными возможностями для разложения любого вектора в рассматриваемом n-мерном векторном пространстве. Решение этой задачи связано с построением подобного базиса из так называемых собственных векторов матрицы A, удовлетворяющих требованиям

Axi = li xi (xi ? 0).

Число li , при котором система однородных линейных алгебраических уравнений (28) имеет нетривиальное решение xi, называется обычно собственным значением матрицы A, соответствующим собственному вектору xi. Очевидно, что в случае диагональной матрицы в качестве xi можно взять l i, при этом li = aii , i = 1, 2, _, n. Существует достаточно широкий класс матриц вида (12), в случае которых в n-мерном векторном пространстве можно выбрать базис из собственных векторов матрицы A. В частности, существует он и в случае симметрической матрицы, для которой в таблице (12) ai j = aji . Для случая такой матрицы все собственные значения (ее спектр), оказывается, действительны. Для большей аналогии с (1) можно считать в дальнейшем матрицу A в (11) симметрической. Читатель, испытывающий затруднения с использованием базиса из собственных векторов матрицы A, может, конечно, в качестве A в системе (11) представлять себе и диагональную матрицу, при этом, однако, не следует пользоваться конкретной информацией о ее собственных векторах и собственных значениях, ограничиваясь лишь фактом их существования, так как в общем случае подобная информация скрыта через (28) в матрице A и ее явное использование обычно сопряжено с большими вычислительными затратами.

Итак, пусть наряду с (24) для вектора u существует разложение

по базису из собственных векторов xi, i = 1, 2, _, n, матрицы A, где ai , i = 1, 2, _, n, - некоторые константы, вообще говоря неизвестные нам, как и собственные векторы. Тогда с учетом (28) аналогично (25)-(27) можно записать

Собственные значения li , i = 1, 2, _, n, матрицы A, используемые в последних равенствах, нам также, вообще говоря, неизвестны. В общем случае разномасштабности эти значения могут быть сильно разбросаны (см. случай диагональной матрицы A, когда li = aii , i = 1, 2, _, n, и система (11) фактически распадается на n независимых уравнений вида (1), каждое из которых может описывать существенно разный по скорости эволюции процесс). Поэтому изменение на одном шаге интересующего нас вектора от значения u (см. (29)) до (см. (31)) может быть существенно разным по слагаемым соответствующей суммы. Это должно предполагать разный уровень тщательности приближения таких слагаемых. Однако метод (16) в силу (30), к сожалению, приближает такие составляющие единообразно, что характерно также и для других разностных методов и часто является основной причиной многих вычислительных трудностей.

Локально на одном шаге применения к системе (11) метода (16), когда можно считать y = u, для локальной ошибки метода в силу (30), (31) можно записать (см. также (18), (19))

где

e(lt) = e lt - E (lt).

Заметим, что с учетом (28) (и связанной с этим терминологии) функцию (19) можно рассматривать как спектральную функцию матрицы (22). Аналогичную роль играет и функция (3) для матричной экспоненты (15). В этом смысле (33) можно рассматривать в качестве функции спектральной ошибки метода на шаге t. Значения этой функции при каждом l = li , i = 1, 2, _, n, представляют собой ошибки в приближении посредством E (lit) спектральных множителей перехода в разложениях типа (30), (31).

Так как для существенно разных значений l = li при данном t величина ошибки вида (33) приближения e lt двумя первыми слагаемыми (см. (19)) суммы типа (3) может быть принципиально разной (см. рис. 1), то по-разному могут быть приближены посредством и соответствующие составляющие (см. (32)) интересующего нас вектора . При достаточно хорошем качестве приближения составляющих, связанных с малыми по абсолютной величине значениями li , уровень приближения других может оказаться и вовсе не приемлемым (см. рис. 1). В случае отрицательных значений li , i = 1, 2, _, n, например, это прежде всего обостряет проблему сохранения в процессе дискретизации (16) свойств устойчивости исходной системы (11). При таких li , i = 1, 2, _, n, в случае любого t > 0 справедливы, например, неравенства

Через условия типа (34), в частности, и может быть обеспечено свойство устойчивости к возмущениям решений системы (11): любые отклонения в решении (искажения коэффициентов ai , i = 1, 2, _, n, в разложении (29)) должны затухать во времени (см. (31)). Очевидно, что сохранить подобное (34) свойство для функции E (lt) в применении к системе вида (11) (при большом разбросе собственных значений такую систему обычно называют жесткой) не может не только метод Эйлера (16), но и любой другой метод, для которого характерно приближение суммой конечного числа первых слагаемых в представлении типа (3). Поиски выхода из этой ситуации привели к получившим сейчас широкое распространение неявным разностным схемам. Простейший пример такой схемы представляет собой так называемый неявный метод Эйлера

построение которого можно связать, например, с приближенной заменой в момент времени t + t значения производной u'(t + t) в уравнении (8) разностным отношением (17). В применении к (1) уравнение (35) относительно скалярной неизвестной легко приводится к виду (18), где

E (lt) = (1 - lt)- 1.

В случае же системы (11) нахождение компонент векторной неизвестной связано, очевидно, с решением системы линейных алгебраических уравнений

Решение этой системы формально можно записать в виде

где через (I - tA )- 1 обозначена матрица, обратная к I - tA. Не затрагивая вопросов существования обратной матрицы, заметим лишь, что элементы матрицы X, обратной к B, связаны с элементами исходной матрицы посредством соотношений, задаваемых матричным равенством BX = I, которое для нахождения каждого из столбцов x i, i = 1, 2, _, n, матрицы X дает систему линейных алгебраических уравнений вида

Bx i = l i.

По аналогии с изложенным ранее функцию (36) можно трактовать как спектральную функцию матрицы из правой части (38). Очевидно, что для случая такой аппроксимации e lt сформулированные выше требования типа (34) выполняются уже при любом t > 0, что является главным аргументом в пользу неявного метода Эйлера. Это же справедливо, конечно, и для других чисто неявных методов, в случае которых вместо спектральной функции e lt матричной экспоненты e At возникает приближение вида

Можно говорить и об иных неявных методах, приводящих к аппроксимациям e lt, как и в (40), основанным на разложениях типа (3) для числителя и знаменателя дроби

где 0 # a < 1. Правда, такие методы (со "взвешенной" через значения параметра a "явностью" и "неявностью") не вносят принципиальных отклонений в логику приводимых здесь рассуждений, поэтому мы и ограничиваемся сейчас, ориентируясь на читателя-школьника, только случаем простейшего неявного метода (35), который, как уже отмечалось, при любых t > 0 сохраняет в аппроксимации (36) свойства (34) спектральной функции e lt. Подобными свойствами неявных методов и объясняется повышенный интерес сегодня во всем мире к таким численным методам. Однако (если даже не принимать во внимание вычислительные затраты на решение систем вида (37) или (39) и забыть о других вычислительных проблемах, связанных с такими системами, например не поднимать вопроса о возможной высокой чувствительности решений этих систем к ошибкам входных данных) сам факт использования неявных методов порождает новые проблемы, особенно в случае существенной разномасштабности составляющих решения исходной эволюционной задачи.

Первая из этих проблем связана с тем, что, к примеру, для неявного метода Эйлера спектральная функция (36) при положительных l может (см. также рис. 1) оказаться не только отрицательной, но и быть меньше единицы по абсолютной величине вопреки естественному свойству

e lt > 1

точной спектральной функции для таких l. Связанная с этим "лишняя устойчивость" неявного метода Эйлера может быть, очевидно, просто опасной при численном моделировании быстрых переходных процессов с локальной неустойчивостью, когда, скажем, в линейном приближении системы (8) в спектре связанной с ним матрицы на короткое время могут появиться и положительные собственные значения. Естественно возникает задача о построении таких численных методов, которые в применении к (11) одновременно обеспечивали бы в случае любого t > 0 сохранение в аппроксимации спектральной функции e lt как условий типа (34), так и свойства вида (42) на соответствующих по знаку собственных значениях матрицы A, при этом не искажалось бы и принципиально важное свойство (спектральной монотонности): , если lj > lk . На традиционном пути конструирования численных методов, сопряженном с упомянутыми выше классическими аппроксимациями e lt (см. (3), (40) и (41)), решение данной задачи найти не удается.

Но даже и в том случае, когда все li < 0, i = = 1, 2, _, n, на который, казалось бы, и ориентированы неявные разностные методы типа (35), характерная для них излишняя абсолютизация свойства устойчивости также порождает серьезные трудности.

Заметим прежде всего, что поведением спектральной функции (33), (36) полностью предопределяются свойства локальной погрешности (32) неявного метода Эйлера применительно к системе (11). В рассматриваемом случае для этой функции можно записать: e(x) = e x - 1/(1 - x), x = lt < 0. Анализ поведения e(x) не вызывает затруднений. Так как e'(x) = e x - 1/(1 - x)2, то функция e(x) достигает своего экстремального значения в точке, являющейся корнем уравнения (1 + z)2 = e z, где z = - x. Школьники легко могут проверить, что положительный корень этого уравнения существует и единствен (z © 2,5). Тем самым можно утверждать, что для каждого собственного значения li < 0, i = 1, 2, _, n, матрицы A существует такое значение t, при котором функция - e(lit) как функция аргумента t (рис. 2) принимает максимальное значение (приближенно равное 0,2). В силу (32) соответствующим выбором начальных условий (предопределяющих значения коэффициентов ai , i = 1, 2, _, n) при t © - 2,5/ lk , 1 # k # n, локальная погрешность (32) рассматриваемого метода может быть сделана неприемлемой. Очевидно также (см. рис. 2), что с уменьшением t (в области изменения шага, расположенной правее левой границы точек указанных экстремумов) абсолютная величина погрешности вида (33) может (вопреки естественным ожиданиям) и возрасти (как и, скажем, евклидова мера вектора-ошибки (32)). При существенном разбросе собственных значений матрицы A подобная немонотонность локальной погрешности по шагу может порождать большие трудности, особенно при расчете быстрых переходных процессов. И лишь при малых значениях шага дискретизации (из области шагов, расположенной левее упомянутой границы точек экстремумов) процесс становится управляемым посредством изменения t, как и в случае явных методов. Разброс собственных значений, однако, может быть настолько значительным, что указанная граница будет очень близкой к нулю. Поэтому решать эту проблему за счет малости шага не просто недостойно науки, но и практически невозможно, вообще говоря, даже на самых мощных ЭВМ как из-за недопустимо большого времени счета, так и из-за роста вычислительных ошибок. В силу этого сформулированные выше требования к конструированию численных методов естественно дополнить условием монотонности по шагу соответствующей локальной погрешности.

Заметим к тому же, что и при малых значениях шага абсолютная величина значений функции (33) как функции спектральной переменной l, ln # l # l1 , может оказаться существенно разной для разных значений l = li , i = 1, 2, _, n, что является главной причиной многих трудностей, возникающих при численном моделировании быстрых переходных процессов (к примеру, погранслоя) в динамических системах, описываемых соответствующими эволюционными уравнениями. Поэтому естественно потребовать дополнительно, чтобы конструируемые численные методы допускали избирательное регулирование уровня приближения разномасштабных составляющих решения.

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

ЛИТЕРАТУРА

1. Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред. Дж. Холла, Дж. Уатта. М.: Мир, 1979. 312 с.

2. Беллман Р. Введение в теорию матриц. М.: Наука, 1969. 368 с.

3. Бобков В.В. О построении методов численного решения систем обыкновенных дифференциальных уравнений // Дифференц. уравнения. 1995. Т. 31, ╧ 7. С. 1174-1178.

4. Бобков В.В. Об одном способе построения методов численного решения жестких систем // Вестн. Белорус. ун-та. Сер. 1, Физика, математика, механика. 1995. ╧ 2. С. 41-45.

* * *

Владимир Васильевич Бобков, доктор физико-математических наук, профессор, зав. кафедрой вычислительной математики Белорусского государственного университета. Лауреат Государственной премии БССР. Область научных интересов - вычислительная математика. Автор более 150 научных публикаций, в том числе десяти книг.