Численное моделирование RC-цепей L-устойчивым методом первого порядка точности Текст научной статьи по специальности «Математика»

УДК 519.622 Е.А. Новиков, Б.У. Уатай

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ RC-ЦЕПЕЙ L-УСТОЙЧИВЫМ МЕТОДОМ ПЕРВОГО ПОРЯДКА ТОЧНОСТИ

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

Ключевые слова: электрическая цепь, дифференциально-алгебраическая система, L-устойчивый метод, контроль точности.

Е.А. Novikov, B.U. Uatay

NUMERICAL MODELING OF RC-CHAINS BY L-STABLE TECHNIQUE OF THE FIRST ORDER ACCURACY

The technique of the first order accuracy for solving differential and algebraic tasks is constructed. An algorithm of differential step is formulated. The results of RC-schemes modeling are given.

Key words: electric circuit, differential and algebraic system, L-stable technique, accuracy control.

При моделировании кинетики некоторых химических реакций, электрических цепей и в ряде других важных приложениях возникает необходимость численного решения жестких задач, неразрешенных относительно производной [1-6]

Р{х,х') = 0, х(70) = х0, ?0 < ? < ^, (1)

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

х' = /(х), х(70) = х0, (2)

Однако сведение (1) к виду (2) требует больших дополнительных затрат на шаг интегрирования, потому что это связано с обращением матрицы Ру=дР(х,у)/ду, которая может быть вырожденной. Разрешенная задача (2), как правило, жесткая, что приводит к необходимости применения методов, которые, в свою очередь, нуждаются в обращении матрицы Якоби. При решении жестких задач широкое распространение получили методы типа Розенброка [7] благодаря простоте реализации и достаточно хорошим свойствам точности и устойчивости. Наибольшее распространение получили схемы, в которых при вычислении каждой стадии применяется одна и та же матрица Якоби [8].

Здесь построен одностадийный ^-устойчивый метод решения задачи (1), который отличается от классического метода типа Розенброка приближенным нахождением производной решения. Сформулирован алгоритм интегрирования с автоматического выбора величины шага интегрирования. Приведены результаты расчетов РС-схем.

Для решения задачи (2) методы типа Розенброка записываются в виде

т 1-1

*Л+1=*Л+ЕМ. /), = Н - а!Ц , 1)пк1 = /1/(хп + £ ДД ), (3)

2=1 2=1

где h - шаг интегрирования; Е - единичная матрица, /п-д/(хпуду - матрица Якоби системы (2); к, 1</<М - стадии метода, а, р, в, 1</<т, 1</</-1 - коэффициенты, определяющие свойства точности и устойчивости (3). В настоящее время методы типа Розенброка трактуются более широко. Под ними понимаются все численные схемы, в которых матрица Якоби или ее аппроксимация вводятся непосредственно в формулу интегрирования [9].

Рассмотрим одностадийный метод (3) вида

хп+1 = хп+рА . /},А = ¥ 00, (4)

где матрица йп определена в (3). При коэффициентах р1=а=1 схема (4) имеет первый порядок точности и является ^-устойчивой. Для контроля точности численной формулы (4) можно применять неравенство [10]

IIК \\<е,

где £ - требуемая точность интегрирования; ||^| - некоторая норма в №.

В случае большой размерности задачи (2) основные вычислительные затраты связаны с обращением матрицы йп. Обычно вместо обращения решается линейная система алгебраических уравнений йпЬ=Ь^П) с применением Ш - разложение матрицы йп с выбором главного элемента по строке или столбцу, а иногда и по всей матрице. Это означает порядка М3 арифметических операций. Обратный ход метода Гаусса стоит порядка N арифметических операций. Таким образом, при большой размерности исходной задачи общие вычислительные затраты фактически полностью определяются временем декомпозиции матрицы йп.

Используя обозначение х-у, задачу (1) можно переписать в виде

х' = у, Г{х,у) = 0,

к'

(5)

с начальными условиями х(Ь)=хо и у(?о)=уо. Дополнительное условие у(/о)=уо можно вычислить, например, решая задачу Р(х0,у)=0 методом установления.

Ниже будем предполагать существование и единственность решения задачи (1) [11]. Будем так же предполагать, что функция Р нужное число раз дифференцируема на шаге интегрирования, а матрица йп=Рпх+аНіРпу не вырождена, что, как правило, выполняется. При записи йп использованы обозначения Рпх=дР(хп,уп)/дх и Рпу=дР(хп,уп)/ду. Тогда т-стадийный метод типа Розенброка применительно к решению задачи (5) можно записать в виде

х , = х

и+1 п

+Емх. Уп+і=Уп+Т.р>к?'

2=1 2=1 2-1 2-1 2-1

ДА'=ЧР„,-ІУ*+Е +2>А')>'

7=1

п I]

7=1

п I]

7=1

(6)

1 г 1 к=-^\к;]}.

т

т

где рI и вц, 1</<т, 1</</-1, - числовые коэффициенты. Отметим, что при применении (6) для решения задачи (2) получим методы типа Розенброка вида (3). При решении задачи (5) они отличаются от (3) приближенным нахождением производной решения. Так как уо вычисляется приближенно, будем предполагать выполненным соотношение ||Яо||=СЛ™ где Яо=Я(хо,уо), |Н| - некоторая норма в №, С - не зависящая от шага интегрирования постоянная. Заметим, что методы типа Розенброка применительно к неавтономной неявной задаче (1) приведены в [4].

Рассмотрим одностадийную формулу типа Розенброка, то есть

*л+і = *„ + МГ. уп+і =Уп + РхЧ.

ак

Б =F -\-ahF

пу

Для определения коэффициентов а и рі разложим стадии кіх и кіу в ряды Тейлора в окрестности точки (хп,уП) до членов с Н2 включительно в предположении, что матрица Рпу невырожденная. В результате получим

к* = Ну - К Лр 1 - аІЇК-\Г у - ^ ] + Шг3),

1 іу п пу п* тгу I- пху п пх пу п -I V

^ = --F ^ + -НР\Р у Р^Л + ОІИ2).

1 пу п пу пху п пх пу п J V /

а а

Подставляя полученные соотношения в (7), запишем

Хп+1=Хп +РМУп -Ку1']^-СР^Ку\1'п,Уп -РщРщРЛ+ОЦ?}.

Ниже будем предполагать невырожденность матрицы Ру=дР(х,у)/ду в некоторой окрестности решения х(<) задачи (5). Разложим точное решение х^) в ряд Тейлора в окрестности точки П Учитывая очевидное равенство Рху+Ру/=0, имеем у'=Ру~]Рху, где у=х'. Теперь окончательно можно записать

х(К+1) = х(гп)+ЫК) - о, 5/г/<; 7<;_у(/(;)+<9(/г3),

где матрицы Рх и Ру вычислены на точном решении (х(Ц у(Ц). Пусть хп=х(ц и уп-1=у(?п-1). Сравнивая ряды Тейлора для приближенного х„+1 и точного х(^+1) решений, получим, что численная схема (7) будет иметь первый порядок точности, если р1=1 и Рп=0(Л), причем а есть свободный коэффициент. Значение параметра а определим из условия Р„+1=0(Л2). Разлагая Р„+1 в ряд Тейлора в окрестности точки (хп, уп), имеем

/7 + Р к?+Р Ц+ООг2).

п+1 п пх I пу 1 \ у

Подставим сюда ряды Тейлора для к1х и к1у. Учитывая очевидные соотношения к1х =0(Л) и к1у =0(Л), запишем Р„+1=а-1(а-1)Р„+0(Л2). Отсюда видно, что из условия Р„=0(Ь) следует равенство Р„+1=0(Л). Выберем а=1, тогда Рп+1 есть величина 0(Л2). Следовательно, она не будет вносить вклад в главный член локальной ошибки - в первый член при разложении ошибки в ряд Тейлора. Кроме того, при а=1 схема (7) является ^-устойчивой, если ее применять для решения скалярного тестового уравнения х'=Лх, где А - произвольная комплексная постоянная, Ре(А)<0.

Контроль точности схемы (7) можно организовать по аналогии [10], то есть на каждом шаге следует проверять неравенство

кхх \\<є,

где к1х задается формулой (7), £ - требуемая точность расчетов, ||^|| - некоторая норма Ям. В расчетах использовалась норма вида

\\¥„\\= шах- ¥п

\ X \+г

Если |х'„|<г, то по /-й компоненте решения контролируется абсолютная ошибка ге, в противном случае контролируется относительная ошибка е. Аналогичная норма применялось для контроля точности схемы (4).

В отличие от методов типа Розенброка (3), применительно к решению задач вида (2), производная решения в (7) вычисляется приближенно. Поэтому при выборе величины шага интегрирования дополнительно проверяется неравенство

В расчетах использовалась норма вида

||<£.

11=тах|

1<г<№

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

системы алгебраических уравнений Р(х0,у0)=0, например, методом установления. Для этих целей можно применить ^-устойчивый (2,1)-метод, описанный в [13]. При решении линейных систем алгебраических уравнений, возникающих при вычислении стадий метода (4), применяется Ш - разложение матрицы йп с выбором главного элемента по столбцу. Так как речь идет о жестких задачах, то расчеты предполагается проводить с двойной точностью. Это означает, что в мантиссе содержится 14 значащих цифр.

Пусть приближенное решение уп в точке и вычислено с шагом Л„. При описании алгоритма учитываются соотношения к1х=0(Л„) и Р„=0(Л„2).

Шаг 1. Вычисляется вектор Рп=Р(х„, у„).

Шаг 2. Вычисляются матрицы РПх и РПу по формулам

= д17(ХтУп) /дх>Рпу= д17(Хп>Уп) / ЗУ .

Шаг 3. Формируется матрица по формуле йп=Рпх+ЬпРпу.

Шаг 4. Выполняется декомпозиция матрицы Оп. Если йп вырожденная матрица, то Ьп полагается равным 0,9Л„ и происходит переход на шаг 3.

Шаг 5. Вычисляется вектор к1х из линейной системы алгебраических уравнений

I) к* = ЫР у -Р).

П 1 V пуу П ПУ

Шаг 6. Вычисляется число д1 из соотношения

^ || &* ||= £ ,

где е есть требуемая точность интегрирования. При записи данного соотношения учтено равенство

к1х=0(Ьп).

Шаг 7. Если ф<1, то есть требуемая точность не выполнена, Лп, то полагается равным дА и происходит переход на шаг 3 (возврат).

Шаг 8. Вычисляется вектор к1у по формуле

к(=\(К-Ну„)■

ап

Шаг 9. Вычисляются П+1, приближенное решение хп+1 и производная приближенного решения уп+1 по формулам

К+1 = К + К Хп+х = Хп + РхК ■ Уп+1 =Уп+ рАу ■

Шаг 10. Вычисляется число д2 из соотношения

Я,2||£Г\р \\=£.

12 11 п п II

При записи данного уравнения учтено равенство Рп=0(Лп2). Отметим, что декомпозиция матрицы йп производилась ранее, и поэтому для задач большой размерности вычисление Ц2 приводит к незначительным вычислительным затратам на фоне декомпозиции, которая (декомпозиция) фактически полностью определяет общие вычислительные затраты.

Шаг 11. Вычисляется новый шаг интегрирования Лп+1 по формуле

/2и+1=1ШП(^2)/2и.

Шаг 12. Если ^+1+Ь1п+1>^, то Лп+1 перевычисляется по соотношению Лп+1=М„+1 (точное попадание в конечную точку интервала интегрирования).

Шаг 12. Если | Лп+1|>10-14, то выполняется следующий шаг интегрирования (переход на шаг 1). В противном случае окончание расчетов.

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

где А и В - квадратные матрицы размерности (№Ы) соответственно с элементами вц и Ьц, д=(д1,...д)т

- М-мерный вектор с постоянными коэффициентами, v(t) - скалярная функция. Во всех приведенных ниже примерах д=(-1, 0, ... , 0)т; матрица А диагональная, а матрица В трехдиагональная с элементами

Решение системы (8) с нулевыми начальными условиями осуществлялось на промежутке [0,к с начальным шагом Ь0=10-5.

Для решения задачи (8) применялись методы (4) и (7). В разрешенном виде задача (8) приведена в [11]. Расчеты проводились с задаваемой точностью £=10-1 и £=10-2. Проводить вычисления методом первого порядка с более высокой точностью нецелесообразно. Разрывы функции v(t) явно не выделялись, ответственность за точность вычислений возлагалась на алгоритм управления величиной шага интегрирования. Для проверки точности вычисления решения по построенному алгоритму оно было вычислено несколькими методами высокого порядка точности.

Задача 1. N=2, ^=5, а11=1, а22=1. Решение в точке ^=5 принимает значения х1=х2=0,99326.

Задача 2. N=4, tk=5, вц=1, в22=10-3, взз=10-6 и в44=10-9. Решение в точке tk=5 принимает значения

х1=0,99323, х2=х3=х4=0,99322.

Задача 3. N=4, tk=0,05, ац=в22=азз=а44=10-2. Решение в точке tk=0,05 принимает значения х1=0,76192, х2=0,55453, хз=0,40285, х4=0,32318.

Задача 4. N=8, tk=0,05, а,-,=10-/, 1</<М. Решение в точке tk=0,05 принимает значения х1=0,98847, х2=0,98720, хэ=0,98708, Х4=Х5=Хб=Х7=Х8=0,98706.

Задача 5. N=8, tk=0,05, а„=10-2, 1</<N. Решение в точке tk=0,05 принимает значения х1=0,75091, х2=0,52607, ха=0,34268, х4=0,20723, х5=0,11645, х6=0,61298Ч0-\ х7=0,31475Ч0-\ ха=0,18768^10-1.

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

1. Брайтон Р.К., Густавсон Ф.Г., Хачтел Г.Д. Новый эффективный алгоритм решения алгебраических систем дифференциальных уравнений, основанный на использовании формул численного дифференцирования в неявном виде с разностями назад // Тр. Института инженеров по электронике и радиоэлектронике. - 1972. - 60. - №1. - С. 136-148.

2. Численные методы решения сингулярных систем / Ю.А. Бояринцев, В.А. Данилов, А.А. Логинов [и др.]. - Новосибирск: Наука, 1989. - 223 с.

3. Gear C.W. Differential-algebraic equations index transformations // SIAM J. Sci. Stat. Comput. - 1988, 9. -№1. - P. 39-47.

4. Новиков Е.А., Юматова Л.А. Некоторые методы решения обыкновенных дифференциальных уравнений, неразрешенных относительно производной //ДАН СССР. - 1987, 295. - №4. - С. 809-812.

5. Левыкин А.Н., Новиков Е.А. Класс (m,k)-методов решения неявных систем // Докл. РАН. - 1996. - №4.

6. Петренко А.И., Тимченко А.П., Слюсарь П.Б. Сравнение программ схемотехнического проектирования на базе набора тестовых задач // Изв. вузов. Радиоэлектроника. - 1980, 23. - №2. - С. 5-12.

Ax' + Вх + v(t)g = 0, 0<t<tk,

(8)

Ъг= 2, l<i<N-l, Ьш=-1, bJ+hJ=bu+1—!. l<i<N-l.

Функция v(t) определяется по формуле

Литература

- С. 442-445.

7. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. - 1963. - №5. - Р. 329-330.

8. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685c.

9. Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений. - М.: Мир, 1988. - 340с.

10. Демидов Г.В., Новиков Е.А. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Числ. мет. мех. сплошной среды. - 1985, 16. - №1. - С. 27-39.

11. Эльцгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. - М.: Наука, 1965. -424с.

12. Deuflhard P., Hairer E., Zugck J. One-step and extrapolation methods for differential-algebraic systems // Numer. Math. - 1987, 51. - P. 501-516.

13. Новиков E.A. (2,1)-метод решения жестких неавтономных задач // Системы управления и информационные технологии. - 2008. - №2(32). - С. 12-15.

УДК 582.28 (571.51) М.Р. Ратова, А.С. Шишикин, А.Н. Борисов

СТРУКТУРА БАЗЫ ДАННЫХ МАКРОМИЦЕТОВ

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

Ключевые слова: грибы съедобные, база данных, условия местообитания.

M.R. Ratova, A.S. Shishikin, A.N. Borisov MACROMYCETE DATABASE STRUCTURE

Issues of database structure development oriented on fungi data analysis and ecological factor research which influence on their growth directly are considered in the article. Possibilities for database application for fungal community information research are given.

Key words: edible fungi, database, habitat conditions.

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

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

Видовое разнообразие грибов обусловлено его симбиотическими связями с древесными породами. Известно, что в сообщество сосны обыкновенной входят белый гриб (Boletus edulis f. pinicola (Vassilk)), масленок зернистый (Suillus granulatos Fr. O.Kuntze), масленок поздний (Suillus luteus Fr.), рыжик сосновый (Lactarius deliciosus Fr.) и др. Микологический комплекс березы представлен белым грибом (березовая форма) (Boletus edulis f. beticola (Vassilk)), подберезовиком (Leccinum scabrum (Bull. ex Fr.) S. F. Gray), груздем настоящим (Lactarius resimus (Fr.)); лиственницы - масленком лиственничным (Suillus grevillei (Klotsch) Sing.), ели - рыжиком еловым (Lactarius deliciosus Fr.), мокрухой еловой (Gomphidius glutinosus Fr.).

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