Численный метод решения жестких задач на неоднородных схемах второго порядка точности * Текст научной статьи по специальности «Математика»

УДК 519.622 Е.А. Новиков

ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ НА НЕОДНОРОДНЫХ СХЕМАХ ВТОРОГО ПОРЯДКА ТОЧНОСТИ*

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

Во многих приложениях возникает проблема численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Основные тенденции при построении численных методов связаны с расширением их возможностей для решения задач все более высокой размерности. Математические постановки практических задач постоянно уточняются, что приводит как к росту размерности, так и к усложнению правой части системы дифференциальных уравнений. В ряде случаев расчеты требуется проводить с так называемой инженерной точностью - порядка 1% и ниже. Это связано с точностью измерения констант, входящих в правую часть системы дифференциальных уравнений. Иногда такая точность расчетов является удовлетворительной с точки зрения поставленной цели. Известно (см., например, [1]), что порядок аппроксимации численной схемы следует сочетать с требуемой точностью расчетов. Здесь на основе явных методов типа Рунге-Кутты и L-устойчивого (2,1)-метода второго порядка точности построен комбинированный алгоритм интегрирования и приведены результаты расчетов.

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

У = I(У). У(0 = Уо. г0 ^^ Ч. (1)

где у и / - вещественные N -мерные вектор-функции;

? - независимая переменная, предложен класс (т,к )-методов.

В отличие от традиционных методов для описания данных численных схем требуется две постоянные: т - число стадий и к - количество вычислений правой части системы (1) на шаге интегрирования. Для решения задачи (1) рассмотрим (2,1)-схему вида

Уп+1 = Уп + РхК + Р2к2. ВпК = ¥(Уп I Впк2 = К (2)

где к1 и к2 - стадии метода, Оп = Е — акАп;

Е - единичная матрица;

к - шаг интегрирования;

Ап - некоторая матрица, представимая в виде Ап = /'п + кВп + 0(к2);

/' = д/ (Уп )/Эу - матрица Якоби системы (1);

* Работа поддержана грантами Президента НШ-3428.2006.9 и РФФИ №06-08-00920.

Bn - не зависящая от шага интегрирования произвольная матрица; а, p1 и p2 - числовые коэффициенты.

Использование матрицы Ап, представимой в виде Ап = /'п + кВп + 0(к2), позволяет применять (2) с замораживанием как аналитической, так и численной матрицы Якоби [3]. В частности, при использовании матрицы Якоби /п—к, вычисленной к шагов назад, имеем Вп =—к/п/п, где

/1 = д 7 (Уп )/Эу 2.

Набор коэффициентов а = 1 — 0,5л/2, р1 = а и р2 = 1 — а обеспечивает второй порядок точности и L-устойчивость схемы (2), а ее локальная ошибка имеет вид 8п = (к3 /6)[(6а — 2)/'2/ + //2 — 3Вп/] + 0(к4). Контроль точности (2) построим по аналогии [4].

Для этого введем обозначение у(]п) = ■’п (к2 — к1), где к1 и к2 вычисляются по формулам (2).

Тогда для контроля точности вычислений на каждом шаге нужно проверять неравенство

Ь(]„)!!<£, 1 < }п < 2, (3)

где £ - требуемая точность расчетов, !! •!! - некоторая норма в И1*, а целочисленная переменная ]п выбирается наименьшей, при которой выполняется неравенство (3). Оценку максимального собственного числа wn0 = кК,таХ матрицы Якоби системы (1), необходимую для перехода на явную формулу, оценим

через ее норму wn 0 = к!! д//Эу !!. Ниже данная оценка будет применяться для автоматического выбора численной схемы.

Теперь для решения задачи (1) рассмотрим явный двухстадийный метод типа Рунге-Кутта [5]

Уп+1 = Уп + 0,5(к1 + к2), к = к/ (Уп),, к2 = к/ (Уп + к1), (4)

который имеет второй порядок точности, а его локальная ошибка дп имеет

вид дп = (к3/12)[2//2/ — //2] + 0(к4). Построим неравенство для контроля точности вычислений. Для этого рассмотрим вспомогательную схему уп+11 = уп + к1 первого порядка точности. С помощью идеи вложенных методов оценку ошибки £п 2 метода второго порядка можно вычислить по формуле £п 2 = Уп+1 — Уп+1 1 = 0,5(к2 — к1). Построим неравенство для контроля устойчивости (4) предложенным в [6] способом. Для этого рассмотрим вспомогательную стадию к3 = к/(уп+1). Заметим, что к3 совпадает со стадией к1 , которая применяется на следующем шаге интегрирования, и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии к1, к2 и к3 применительно к задаче у = Ау, где А есть матрица с постоянными коэффициентами. В результате получим к1 = Хуп, к2 = (X + X2)уп и к3 = (X + X2 + 0,5X3)уп, где X = кА. Легко видеть, что к2 — к1 = X 2 уп и 2(к3 — к2) = X3 уп. Тогда согласно [6] оценку максимального собственного числа wп 2 = кК, тах матрицы Якоби системы (1) можно вычислить по формуле

w 2 = 2тах | к3 — к21 1 к2 — к! |. Интервал устойчивости схемы (4) второго порядка точности прибли-

п'2 1<;^ 3 2 2 1

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

методе применяется мало итераций и дополнительные искажения вносят нелинейность задачи (1). Поэтому контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг кп+1 будем вычислять следующим образом. Новый шаг кас по точности определим по формуле кас = дкп, где кп есть последний успешный шаг интегрирования, а д, учитывая соотношение к2 — к1 = 0(к2), задается уравнением д2!! к2 — к1 !!= £. Шаг к по устойчивости зададим формулой к = йкп, где й определяется из равенства dwn2 = 2. Тогда прогнозируемый шаг кп+1 вычисляется по формуле

кп+1 = тах[кп, тт(к°с, к")]. (5)

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

Теперь для численного решения задачи (1) рассмотрим схему вида

Уп+1 = Уп + Г1к1 + Г2 к2 , к1 = к/ (Уп) , к2 = к/ (Уп + к1) . (6)

Заметим, что при г1 = г2 = 0,5 численная формула (6) имеет второй порядок точности и совпадает с (4). Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (6) для решения скалярного тестового уравнения у = Лу , у(0) = у0, Яе(Л) < 0. Получим уп+1 = Q(х) уп,

где функция устойчивости Q(х) имеет вид Q(х) = 1 + (г + г2)х + г2х2, х = кЛ. Требование первого порядка точности приводит к соотношению г + г2 = 1, которое ниже будем считать выполненным. Теперь выберем г2 таким образом, чтобы метод (6) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Т2(г) = (2г2 — 1) на промежутке [-1,1]. Проведем замену переменных, полагая г = 1 — 2х/у. Получим Т2(х) = 1 — 8х/у+ 8х2/у2, при этом отрезок [у,0] отображается на отрезок [-1,1]. Нетрудно показать, что среди всех многочленов вида Р2(х) = 1 + х + с2х2 для Т2(х) неравенство !Т2(х) !< 1 выполняется на максимальном интервале [у,0], у =—8. Потребуем совпадения коэффициентов Q(х) и Т2(х) при у =—8. Это приводит к соотношениям г1 + г2 = 1 и г2 = 1/8. В результате имеем коэффициенты г1 = 7/8 и г2 = 1/8 метода первого порядка точности с максимальным

интервалом устойчивости, локальная ошибка 8п которого имеет вид 8п = (3/8)к2// + 0(к3). Для контроля точности численной формулы первого порядка будем использовать оценку локальной ошибки. Учитывая, что к2 — к1 = к2/П/п + 0(к3) и вид локальной ошибки, неравенство для контроля точности

записывается в виде !! к2 — к1 !!< 8£3, где !! •!! - некоторая норма в , £ - требуемая точность расчетов.

Построим неравенство для контроля устойчивости метода первого порядка. Для этого рассмотрим вспомогательную стадию к3 = к/(уп+1). Запишем стадии к1, к2 и к3 применительно к задаче

y = Ay, где A есть матрица с постоянными коэффициентами. В результате получим к1 = Xyn,

к2 = (X + X2)yn и k3 = (X + X2 + 0,125X3)yn, где X = hA. Легко видеть, что k2 - кх = X2yn

и 8(к3 — к2) = X3yn. Тогда согласно [6] оценку максимального собственного числа wn1 = hAn max матрицы Якоби системы (1) можно вычислить по формуле w , = 8max | k3 — к21 /1 к2 — к,' |. Интервал

11,1 1<i<N 3 2 2 1

устойчивости численной схемы (6) равен восьми. Поэтому для ее контроля устойчивости можно применять неравенство wn 1 < 8.

На основе построенных явных методов первого и второго порядков точности легко сформулировать алгоритм переменного порядка и шага. Расчеты всегда начинаются методом второго порядка как более точным. Переход на схему первого порядка осуществляется при нарушении неравенства wn 2 < 2. Обратный

переход на метод второго порядка происходит в случае выполнения неравенства wn 1 < 2. При расчетах

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

В случае использования схемы (2) формулировка алгоритма интегрирования также не вызывает

трудностей. Нарушение неравенства wn 1 < 8 вызывает переход на схему (2). Передача управления явным методам происходит в случае выполнения неравенства wn 0 < 8. Численную формулу (2) без потери порядка точности можно применять с замораживанием матрицы Dn. Отметим, что при замораживании матрицы Якоби величина шага интегрирования остается постоянной. Попытка замораживания матрицы Dn

осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях: 1) нарушение точности расчетов; 2) если число шагов с замороженной матрицей достигло заданного максимального числа iqh; 3) если прогнозируемый шаг больше последнего успешного в qh раз. Числами

iqh и qh можно влиять на перераспределение вычислительных затрат. При iqh = 0 и qh = 0 замораживания не происходит, при увеличении iqh и qh число вычислений правой части возрастает, а количество обращений матрицы Якоби убывает. Нижепостроенный алгоритм переменного порядка и шага, а также с автоматическим выбором явной или L-устойчивой численной схемы будем называть RKMK2.

Расчеты проводились на IBM PC Athlon(tm)XP 2000 с двойной точностью. В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Расчеты проводились с точностью £ = 10—2. Это связано с тем, что в построенном алгоритме применяются схемы низкого порядка точности, и поэтому данным методом осуществлять расчеты с более высокой точностью нецелесообразно. Сравнение эффективности проводилось с известным методом Гира в реализации А. Хиндмарша DLSODE из коллекции ODEPACK [7]. Ниже через ifu и ija обозначены, соответственно, суммарное число вычислений правой части и количество обращений (декомпозиций) матрицы Якоби задачи (1), которые позволяют объективно оценить эффективность алгоритма интегрирования.

Тестовый пример описывается системой двух дифференциальных уравнений в частных производных с начальными и граничными условиями. Лаборатория Akzo Nobel Central Research сформулировала задачу исследования проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма [8]. Исследование проводилось в диагностических и терапевтических целях. Рассматривается система одномерных уравнений реакции-диффузии

du d 2 u , dv ,

— = —т- — Kuv , — = — Kuv , (7)

dt Эх2 dt

к

которые возникают из химической реакции A + B ^ C, где A - антитело с радиоактивной меткой, реагирующее с субстратом B - тканью, пораженной опухолью; к - константа скорости реакции. Концентрации A и B обозначены через и и v соответственно. При выводе уравнений (7) предполагалось, что кинетика

реакции описывается законом действующих масс, причем реагент А подвижен, тогда как реагент В неподвижен.

Изучается полубесконечная пластина, внутри которой субстрат В равномерно распределен. Реагент А, попадая на поверхность пластины, начинает проникать в нее. Для моделирования проникновения уравнения (7) рассматриваются в полосе 8Т = {(х,£): 0 < х <^,0 < £ < Т} с начальными

и(х,0) = 0,у(х, 0) = У0, х > 0 и граничными и(0, £) = ф(£), 0 < £ < Т условиями, где У0 константа. Для численного решения переменная х преобразуется так, чтобы полубесконечная пластина преобразовалась в конечную. Такое преобразование обеспечивает специальное семейство преобразований Мебиуса £ = х/(х + с), С > 0. Каждое преобразование из этого семейства преобразует БТ в прямоугольник {(<£, £): 0 <^< 1, 0 < £ < Т}. С использованием £ задача (7) переписывается в виде

ди (С -1)4 д2и 2(С -1)3 ди ду

— = --------Т + о-------------------киу, — = - киу,

д£ с2 д^2 С2 д£ д£

с начальными и (£, 0) = 0, у(£, 0) = У0 ,С> 0 и граничными и (0, £) = ф(£), ди (1, £ )/д^ = 0, 0 < £ < Т условиями. Последнее граничное условие получено из соотношения ди(^,£)/дх = 0. Дискретизация производных по пространственным переменным с использованием метода прямых приводит к задаче Коши для системы обыкновенных дифференциальных уравнений. Для дискретизации применяется

равномерная сетка {^}, С'$ = 1 'АС, АС = 1^ , 1 — 1 — N . Через и ■ и у. обозначены аппроксимации и(£1, £) и у(£1, £) соответственно. Очевидно, что и$ и у$ являются функциями от £. Дискретизации производных первого и второго порядков по пространственной переменной соответственно имеют вид

ди■ и;+, -, д2и■ и, , -2и■ + и;+,

____з_ = 1+1 1-1 _____з_ = 1-1______1 1+1 1 — , — N

дС 2а£ ’ дС~ (а£)2 ’ •

Значения и0 и м№+1 получены из граничных условий, они имеют вид и0 = ф(£) и иИ+1 = и№. Полагая у = (м1, У1, и2, У2 ,..., им, У№ )Т и Т = 20, эта полудискретизированная задача записывается в виде

^ = /(£,у), у(0) = g, уе Я2\ 0 — £ — 20, (8)

Ж

где N - задаваемый пользователем параметр. Функция / определяется формулами

* =а У21+1 - У2 _-3 + в у2 _-3 - 2 у21-1 + у21+1 - ку у / 21-1 а1 2А£ 1 (АС)2 1 -1 1 ’

/21 =-ку2 .Л 1 -1,

где

а = 2(1А^-1)3 С2, р. = 0-Д£- 1)4с2, 1 — 1 — N,

АС = 1/ N, у-1(£) = ф(£), у2^+1 = У2,N -1 ,

g е Я2, g = (0,У0,0,У0, ...,0,У0)Т.

Функция p(t) = 2 при t е (0,5] и p(t) = 0 при t е (5,20], то есть р имеет разрыв первого рода в точке t = 5. Согласно [8] подходящими значениями для параметров k, v0 и с являются к = 100, v0 = 1 и

с = 4. В [9] приведены подпрограммы на фортране для вычисления правой части и матрицы Якоби задачи (8). Там же приведены результаты расчетов с высокой точностью и расширенной разрядной сеткой.

Ниже расчеты проводились с численной матрицей Якоби при N = 200, то есть система (8) состоит из 400 уравнений. Задача о нахождении разрыва функции p(t) при t = 5 возлагалась на алгоритм управления шагом. Решение данной задачи алгоритмом RKMK2 вычислено с затратами ifu = 14 106 и ija = 38. При расчетах только по L-устойчивой схеме (2) затраты ifu = 15 954 и ija = 43. При расчетах по явным методам переменного порядка и шага затраты ifu = 51 014. При расчетах программой DLSODE требуемая точность

10-2 достигается при задаваемой точности 10-3 с затратами ifu = 25 358 и ija = 62. При более высокой точности расчетов DLSODE эффективнее построенного алгоритма. Это естественно, потому что построенный алгоритм основан на схемах низкого порядка точности и проводить с его помощью высокоточные расчеты нецелесообразно. Заметим, что алгоритм на основе явных схем с переменным порядком и шагом по времени счета более чем в 2 раза эффективнее других методов, что является следствием достаточно большой размерности задачи (8). С ростом N преимущество явных методов возрастает.

Построенный алгоритм RKMK2 предназначен для расчетов с небольшой точностью - порядка 1% и ниже. В этом случае достигается его максимальная эффективность. В RKMK2 с помощью признака можно задавать различные режимы расчета: 1) явными методами первого или второго порядков точности с контролем или без контроля устойчивости; 2) явными методами с переменным порядком и шагом; 3) L-устойчивым методом с замораживанием или без замораживания как аналитической, так и численной матрицы Якоби; 4) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования.

Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f.

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

Литература

1. Об оптимизации параметров методов численного решения задачи Коши для обыкновенных дифференциальных уравнений / С.С. Артемьев [и др.] // Численные методы механики сплошной среды. - 1984. -№2. - С. 5-14.

2. Новиков, Е.А. Одношаговые методы решения жестких систем / Е.А. Новиков, Ю.А. Шитов, Ю.И. Шокин // ДАН СССР. - 1988. - Т. 301. - №6. - С. 1310-1314.

3. Новиков, Е.А. Алгоритм интегрирования жестких систем на основе (іті,к)-метода второго порядка точности с численным вычислением матрицы Якоби: Препринт №20 I Е.А. Новиков, Ю.А Шитов. - Красноярск: ВЦ СО РАН, 19ВВ. - 23 с.

4. Новиков, Е.А. О повышении эффективности алгоритма интегрирования на основе формулы типа Ро-зенброка второго порядка точности за счет замораживания матрицы Якоби: Препринт №592 I Е.А. Новиков, В.А. Новиков, Л.А. Юматова. - Новосибирск: ВЦ СО АН СССР, 19В5. - 26с.

5. Кнауб, Л.В. Алгоритм интегрирования переменного порядка и шага на основе явного двухстадийного метода Рунге-Кутты I Л.В. Кнауб, Ю.М. Лаевский, Е.А. Новиков II СибЖВМ. - 2007. - Т. 10. - №2. -С. 177-1В5.

6. Новиков, Е.А. Явные методы для жестких систем I Е.А. Новиков. - Новосибирск: Наука, 1997. - 197с.

7. http:IIwww.netlib.orgIodepackIindex.html.

В. Mazzia, F. Test Set for Initial Value Problem Solvers I F. Mazzia, F. lavernaro II Department of Mathematics,

University of Bari, August. - 2003.

9. http:IIpitagora.dm.uniba.itI testsetIsrcIproblemsImedakzo.f.

--------- ^-----------

УДК 514.7 А.Г. Рогачевский

ДВУМЕРНАЯ МЕХАНИКА ЭЛЕКТРОМАГНИТНОГО ПОЛЯ НА ПЛОСКИХ СЛОЕНИЯХ ПРОСТРАНСТВА r4 . ЧАСТЬ I

Рассматривается голономное векторное поле специального типа, являющееся полем базисного вектора ортогональной системы координат (векторное ОК-поле). В дальнейшем (часть II) это поле интерпретируется как векторный потенциал электромагнитного поля. В данной части развивается необходимый для приложения к электродинамике математический аппарат. В частности, изучается векторное ОК-поле на плоских двумерных слоениях псевдоевклидова пространства R14 .

Векторные ортогонально-коодинатные (ОК) поля

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

Векторное ОК-поле - это голономное поле А(г), которое может служить базисным вектором некоторой

тэ4

ортогональной криволинейной системы координат в R1 . Такие координаты up назовем А-координатами

(р = 0, 1, 2, 3). Соответствующие базисные векторы обозначим L(p) = L(p), причем пусть А = L(0) , u0 = а . Индексы координат тензоров в А-пространстве будут писаться в скобках, если нужно отличие от тензорных

индексов в координатах R4. Таким образом, L((p) = dxk /dup, xk- координаты в R4.

Далее будут использоваться те же обозначения, что и в [1-5]. Приведем некоторые из этих обозначений и определений. Интегральная линия Іа ОК-поля А(г) дается уравнением r = r(a). При заданном А(г) координата а - инвариантный параметр, так как А(1(а)) = r'a . Единичная касательная к La n = r'ш, где ы -натуральный параметр на La, одновременно является единичной нормалью к координатным гиперповерхностям f(r) = а, соответствующим А(г) как голономному полю [6]. Так как А(г) будет интерпретироваться как векторный потенциал, считаем, что (А, А) > 0; а также используем обозначение (А, А) = а2 . Как следствие, имеем (n'o , n'o) < 0.