ББК 31.21 Р24 УДК 621.3.011/014:681.3 Авторы: М. Г. Александрова, А. Н. Белянин, В. Брюкнер, X. Вернер, А. И. Гой, Л. В. Данилов, В. М. Золотницкий, Р. Зюссе, Г. Ульман, Е. С. Филиппов, Г. Шайнерт, Ю. Шмидт Расчет электрических цепей и электро- Р24 магнитных полей на ЭВМ/ М. Г. Александ- рова, А. Н. Белянин, В. Брюкнер и др.: Под ред. Л. В. Данилова и Е. С. Филиппо- ва.—М.: Радио и связь, 1983.—344 с., ил. В пер.: 1 р. 10 к. Книга содержит теоретические сведения, алгоритмы, программы и задачи по расчету электрических цепей и стационарных электромагнитных полей на ЭВМ. В ней от- ражены результаты, полученные в рамках совместной науч- ной и учебно-методической работы специалистов СССР и ГДР. Программы написаны на языках ФОРТРАН и АЛГОЛ. Для инженерно-технических работников и студентов электротехнических и радиотехнических специальностей. 2402020000-036 Р——————————26-82 046(01)-83 ББК 31.21 РЕЦЕНЗЕНТЫ: д-р техн. наук, профессор А. И. Князь, канд. техн. наук Л. П. ГАВРИЛОВ Редакция литературы по радиоэлектронике © Издательство «Радио и связь», 1983 ПРЕДИСЛОВИЕ Одна из важнейших задач в подготовке современно- го инженера—обучение работе на ЭВМ. Это в полной мере относится к подготовке специалистов в области радиоэлектроники, автоматики, энергетики. В то же вре- мя книг, из которых читатель мог бы одновременно по- черпнуть сведения как о теории вопроса, так и об алго- ритмах и программах расчета соответствующих устройств, мог бы на примерах и задачах понять осо- бенности применения ЭВМ в конкретных ситуациях, сравнительно немного. Цель настоящей книги хотя бы в малой степени восполнить этот пробел. По своему научному направлению и плану построения она близка к книге Л. О. Чуа и Пен-Мин Лина «Машинный анализ электронных схем» (Москва, «Энергия», 1980 г.), одна- ко не повторяет содержания последней. Книга написан": коллективом авторов—специалистов по теоретической электротехнике Советского Союза и Германской Демократической Республики. Многолетнее научное сотрудничество между Ленинградским электро- техническим институтом им. В. И. Ульянова (Ленина) и Высшей Технической Школой г. Ильменау (ГДР) по- зволило получить результаты, используемые как в на- учной, так и в учебно-методической работе вузов обеих стран. Было признано целесообразным изложить неко- торые из этих результатов в книге, полезной как для студентов и аспирантов электротехнических и радиотех- нических специальностей, так и. для специалистов. Это и определило содержание настоящей книги. Авторы ясно сознавали невозможность сколь-нибудь полного охвата в рамках одной книги всего накопленного к настоящему времени огромного материала по расчету электрических цепей и стационарных электромагнитных полей на ЭВМ. Поэтому при отборе материала в первую очередь 3 отдавалось предпочтение либо полностью оригинальным результатам, либо оригинальным изменениям и уточне- ниям известных алгоритмов и программ, либо результа- там, слабо освещенным в современной литературе. В то же время для удобства использования книги в учебных целях авторы сочли целесообразным включить в нее материалы чисто учебного характера. Книга состоит из двух частей. Часть I посвящена рас- чету электрических цепей. Главы 1—3, в которых рас- сматриваются вопросы формирования уравнений узло- вых напряжений, расчет линейных цепей на основе построения дискретных резистивных схем и расчет ли- нейных цепей лестничной структуры, носят учебный ха- рактер (за исключением § 2.3, 2.4). Во всех остальных главах ч. I содержится материал, представляющий не только учебный, но и методический и научный интерес. Описанный в § 2.3, 2.4 метод расчета плохо обусловлен- ных линейных систем является оригинальным, а его эф- фективность подтверждается приведенными примерами. В гл. 4 излагается один из подходов к анализу нелиней- ных резистивных цепей, основанный на введении так на- зываемых функций состояния и использующий обобще- ние принципа минимума мощности. Большое внимание уделено численным методам расчета переходных про- цессов в нелинейных цепях на основе явных алгоритмов численного интегрирования и сравнению эффективности этих алгоритмов (гл. 5). В гл. 6 рассматривается при- менение кубической сплайн-аппроксимации для расчета периодических режимов в нелинейных цепях. Там же изложен другой эффективный метод анализа периодиче- ских режимов в нелинейных цепях, являющийся разно- видностью метода точек. Две последние главы ч. I посвящены моделированию электрических цепей. В гл. 7 изложен подход к построе- нию моделей нелинейных цепей, имеющих меньший по- рядок по сравнению с прототипом. Глава 8 посвящена построению с помощью ЭВМ моделей нелинейных цепей, осуществляющих отображение заданного вектора вход- ных сигналов в заданный вектор выходных сигналов. Часть II книги посвящена вопросам расчета элек- тромагнитных полей и применению методов теории поля для расчета емкостей и индуктивностей. Рассматрива- ются только стационарные поля. При отборе материала для этого раздела авторы стремились к т^му, чтобы рассматриваемые задачи при всем их многообразии сво- ДилНсь к использованию небольшого числа простых вы- числительных алгоритмов — вычислению рядов и инте- гралов (гл. 9), решению линейных интегральных урав- нений Фредгольма первого и второго рода (гл. 10 и 12), решению систем линейных алгебраических уравнений (гл. 11 и 13). Несмотря на простоту применяемых алго- ритмов, обоснование некоторых из них в отношение точ- ности и численной устойчивости представляет трудную задачу, недостаточно освещенную в литературе (особен- но это относится к алгоритму решения интегрального уравнения Фредгольма первого рода). Однако накоп- ленный в настоящее время опыт применения указанных алгоритмов позволяет считать их достаточно апробиро- ванными и надежными при решении электротехнических задач. Форма изложения материала в ч. II несколько отли- чается от ч. I. Если в теории цепей обычно удается раз- работать алгоритмы и программы, обслуживающие до- статочно широкие классы задач, то в теории электро- магнитного поля это сделать значительно труднее, а иногда вообще невозможно. Поэтому в каждой главе ч. I вначале излагаются теоретические сведения и опи- сание алгоритма, затем приводится программа, обслу- живающая достаточно широкий класс цепей, и даются примеры решения конкретных задач. В ч. II в каждой главе наряду с изложением теоретических сведений рас- сматривается конкретная задача и приводятся алгорит- мы и программы расчета, специально предназначенные для решения этой задачи. Расчет индуктивностей и емкостей рассматривается в гл. 9, 10. Остальные главы ч. II (и частично гл. 9) посвящены расчету полей постоянных магнитов различ- ной конфигурации и высоковольтных разрядников ме- тодами вторичных источников, сеток и эквивалентных зарядов. Языки программирования ФОРТРАН и АЛГОЛ в настоящее время наиболее популярны. Специалисты ГДР чаще применяют язык АЛГОЛ, в то время как в СССР наибольшее распространение получил ФОРТРАН. Это и определило использование в книге двух алгоритмических языков: ФОРТРАНА в советских программах и АЛГОЛА в немецких. Главы ч. 1 написаны Л. В. Даниловым, В. М. Золот- ницким, Е. Филипповым, Р. Зюссе, Г. Ульманом, В. Брюкнером, Г. Шайнертом; главы ч. II—М. Г. Алек- 5 сандровой, А. Н. Беляниным, А. И. Гоем, X. Вернером, Ю. Шмидтом. Перевод рукописи, написанной немецки- ми авторами, выполнен М. Г. Александровой и В. М. Зо- лотницким. Авторы выражают глубокую благодарность Д. Д. Сотникову, М. С. Жихареву, Б. Н. Макееву и С. И. Коннику за помощь в составлении и отладке про- грамм. При подготовке рукописи к печати большую по- мощь оказали В. А. Прохорова, Л. П. Коновалова, Т. А. Латынова, которым авторы также выражают свою искреннюю признательность. Значительному улучшению содержания книги способствовали замечания, высказан- ные в отзывах рецензентов: д-ра техн. наук, профессора А. И. Князя и канд. техн. наук доцента Л. П. Гаврилова. Часть I. РАСЧЕТ ЭЛЕКТРИЧЕСКИХ ЦЕПЕЙ НА ЭВМ Глава 1. ФОРМИРОВАНИЕ УРАВНЕНИЙ УЗЛОВЫХ НАПРЯЖЕНИЙ Анализ электрических цепей начинается с формиро- вания уравнений, среди различных форм которых широ- кое распространение получили уравнения узловых на- пряжений. Формирование этих уравнений отличается простотой и универсальностью алгоритмов. При анализе периодических режимов в линейных электрических це- пях методом комплексных амплитуд, а также при ана- лизе переходных процессов в линейных и нелинейных цепях на основе дискретных схем замещения расчет сво- дится к исследованию некоторой вспомогательной рези- стивной цепи, для которой уравнения узловых напря'--^- ний представляют собой систему линейных алгебраиче- ских уравнений. В данной главе приводятся алгоритм и программа формирования уравнений узловых напряже- ний цепей, содержащих элементы R, L, С и работающих в установившемся синусоидальном режиме. Решение полученных уравнений в данной главе не рассматривается. Существует ряд стандартных про- грамм, основанных на алгоритмах Гаусса, Краута и дру- гих, позволяющих решить уравнения узловых напряже- ний. Однако если система уравнений плохо обусловлена, что характерно, например, для систем уравнений боль- шой размерности, то стандартные алгоритмы могут дать неверные результаты. В этом случае применяются дру- гие методы. Некоторые из них рассмотрены в гл. 2. 7 1.1. АЛГОРИТМ ФОРМИРОВАНИЯ УРАВНЕНИЙ УЗЛОВЫХ НАПРЯЖЕНИЙ Рассмотрим линейную цепь, состоящую из двухпо- люсных элементов, зависимых источников тока, управ- ляемых напряжениями, и независимых источников тока. Для общности примем, что каждая ветвь цепи является составной и включает источник тока. Введем составные ветви двух видов: 1) параллельное соединение элемента Gk и независимого источника тока Ik (рис. 1.1,й); 2) па- 4 tth о- Рис. 1.1 (1.1) (1.2) раллельное соединение зависимого источника тока, управляемого напряжением ветви I (ii=gmUi), и неза- висимого источника тока Ik (рис. 1.1,6). Уравнение составной ветви первого вида l'&==G&t4—Ik, где Gk — в общем случае комплексный параметр. Уравнение составной ветви второго вида lk=gmUl—Ik. После установления всех Л^в составных ветвей нуме- руются ветви от 1 до NB и все узлы от 0 (базисный узел) до N=Ny—\. Наиболее рациональный, легко программируемый путь формирования уравнений узловых напряжений с помощью ЭВМ состоит в последовательном заполне- нии элементов матрицы Gy узловых проводимостей и вектора 1у узловых токов. Вначале принимаются Gy==0 и 1у=0. При поступлении перфокарты с данными вет- вей, начиная с первой, их параметры поочередно вводят- ся в ЭВМ в качестве соответствующих элементов Gy 1у. Формирование матрицы узловых проводимостей и век- тора узловых токов завершается при поступлении пер- фокарты с данными последней ветви. Параметры составных ветвей обоих видов войдут в элементы Gy и 1у следующим образом: 1. Составная ветвь с элементом G, направленная от узла / к узлу k. Ток ветви, выходящий из узла / и вхо- дящий в узел k, согласно (1.1) равен l,-= G (Uy,—Uyft) —/==—t'ft, (1.3) где Uy, — потенциал /-го узла, /=1, 2, ..., N. Поэтому параметр G ветви войдет в матрицу узловых проводимостей так: (1.4) г, - '[•• / . G. k .—G. . .' "У- т .—О . . G . . . Ток / источника тока войдет элементом строк / и k вектора 1у: / * 1у= [.../...-/...]. (1.5) Если ветвь направлена от узла / к базисному узлу (k=0), то проводимость ветвей G войдет только в эле- мент с индексом jj, а ток источника — только в состав- ляющую /. 2. Составная ветвь с источником тока, управляемым напряжением, направленная от узла /' к узлу k и управ- ляемая ветвью =2. Решение. В результате расчета получены следующие урав- нения: (2,7 + j0,75) 0, — (2,5 + j) Us + j0,25[/„ = 3 -j4; 12 — 0,5(/, + (0,833 + /•) Us — 0,333[/, == 3 + /8, ( — 2+/0.75) t/, + (1,667 + ДОг + (0,583 — /0,25)t/, = 0. 1.2. Составить уравнения узловых напряжений для цепк, содер- жащей шесть независимых узлов (рис. 1.3), при <о=2. Параметры цепи: /?i=l, R2=2, Rs==3, ^4=2, Rs=5, Р^=5, ^=1; Li=l, Ls=2, Ls=4; Ci=l; C2=0,5, Сз=0,2; C4=0,4; I(ii=3, /02=—/2, /оз=3+/4; Ую=2+/3, Уо2=6-Н6. ® Рис. 1.2 Рис. 1.3 Решение. Расчет привел к следующим уравнениям: (4,5 +;4,5) Ut — (2 + ;5) 0,, - \U, — 0,5^ + /0,5^ =3,2; — /2(/i + (0,3333 + ;3) Us — JU, — 0,3333^6 = 0,2; (3 — ;3) У, - (4 + /4)С/2 + (7,5 + /6,75)^3 + /0,25С/5 — 0,5У, = 0,2 ; -0.5[/,+(0,7+J04)^-0,2t/6=- -;2; + /0,56'i — 0,3333^ + /0,25t/3 — 0,2^ + (0,7333 — /0,875) 0^ — —0,2^=0,2; — 05U,— 0,20,, + (0.7 + /0,8)t/. = 3 + /4. 13 1.3. Составить уравнения узловых напряжений для электриче- ской цепи, изображенной на рис. 1.4. Параметры цепи приведены на схеме, <в=2. Рис. 1.4 Решение. Система уравнений имеет следующий вид: ;0,75^-^=4-H2; (0,667 — /0,417) t/a — 0,677^ = 1 — j2; (0,333 + !)0г +.0,333^ —0,667^ = — 2 + /; - ^1 + ( - 1 - Л U, + (0,333 + /•) t/з + (0,667 + /2) ^ = - 3 - /. Глава 2. АНАЛИЗ ЛИНЕЙНЫХ ЦЕПЕЙ С ПОМОЩЬЮ ДИСКРЕТНЫХ РЕЗИСТИВНЫХ СХЕМ ЗАМЕЩЕНИЯ В данной главе рассматривается составление дис- кретных резистивных схем замещения, которые являют- ся моделью численного алгоритма решения дифференци- альных уравнений, описывающих цепь. При этом отпа- дает необходимость в составлении дифференциальных уравнений исходной цепи. Анализ составленной дискрет- чой резистивной цепи сводится к решению системы ли- нейных алгебраических уравнений, поэтому эффектив- ность анализа переходных процессов во многом опреде- ляется эффективностью алгоритма решения линейных алгебраических уравнений. Проблема эффективности указанных алгоритмов существенно усложняется при 14 увеличении размерности цепи, а также для жестких си- стем, характеризуемых большим разбросом постоянных времени. В данной главе рассмотрен алгоритм расчета систем линейных алгебраических уравнений, отличаю- щийся большей точностью расчета по сравнению с наи- более употребительными алгоритмами современных программ. 2.1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ Численный анализ переходных процессов в динами- ческой цепи, содержащей наряду с резистивными индук' тивные и емкостные элементы, можно проводить с по- мощью приближенной резистивной дискретной схемы замещения. Для всех шагов вычислений получаются резистивные схемы одинаковой структуры, но с различными значе- ниями параметров. Анализ резистивной дискретной схе- мы проводится многократно последовательно на каждом шаге по уравнениям узловых напряжений. Рассмотрим простейшую дискретную схему замеще- ния линейной динамической цепи на основе неявного алгоритма Эйлера. Неявная (обратная) формула алго- ритма Эйлера для численного решения дифференциаль- ного уравнения первого порядка x=f(x, t) (2.1) Xn+l=Xn+hf(Xn+l, tn+\)-- порядковый номер шага. значение производной на где h — величина шага; п — Отсюда приближенное (п+1)-м шаге x,i+\=(xn+\—Xn)/h. (2.2) Используя это выражение, получаем дискретные схе- мы замещения емкостного и индуктивного элементов. Для линейного емкостного элемента с характеристикой i=Cdu/dt имеем с учетом (2,2) in+i w (C/h) un+i- (C/h) Un. (2.3) Этому выражению соответствует дискретная схема для (п+1)-го шага (рис. 2.1,а), состоящая из парал- лельно соединенных резистивного элемента с постоянной (не зависящей от п) проводимостью G==C/h и источни- ка тока с током Ciin/h, зависящим от напряжения на емкости в начале шага. На рис. 2.1,6 показана последо- вательная схема с источником напряжения, полученная в результате преобразования источника тока. 15 Для линейного индуктивного элемента с характери- стикой u==-Ldl/dt с помощью (2.2) получаем прибли- женное выражение, дуальное (2.3): Un+i w (L/h) in+i- (L/h) in. (2.4) Соответствующая дискретная схема для (/г+1)-го шага из последовательно соединенного сопротивления I^=L/h и источника напряжения Lin/h показана на рис. 2.1,в. Преобразовав источник напряжения, получим схему ((РИС. 2.1,г) из параллельно соединенной проводи- о-» и /г+/ а). /?=/,//? и П*1 \G=fl/L в) о~ г) Рис. 2.1 мости G==h/L и источника тока с током in, зависящим от тока через индуктивность в начале шага. При замене всех емкостей и индуктивностей приве- денными приближенными схемами получим дискретную резистивную цепь для (?г+1)-го шага, соответствующую неявному алгоритму Эйлера. При анализе цепи по урав- нениям узловых напряжений используются параллель- ные резистивные схемы с источниками тока. Для иллюстрации составления дискретной схемы рассмотрим простейшую цепь второго порядка (рис.2.2,а) с численно заданными значениями элементов, подклю- ченную к источнику ступенчатого тока. Приняв шаг /г==0,1, получим дискретную схему замещения (рис. 2.2,6). 16 В соответствий Си схемой записываем уравнения для дискретных значений узловых напряжений г 13 -ПГа^(„+1)1П+10^1(„)1 [-1 22 J L ^ (,.,.,) J-[ 20^ („) J- Характерной особенностью дискретных схем замеще- ния на основе неявного алгоритма Эйлера является не- зависимость параметров резистивных элементов от но- мера шага, что приводит к постоянству матрицы узло- (Т) ЕЗ ® ОФ isW Ми,, ~гп а) 6) Рис. 2.2 вых проводимостей, т. е. левой части уравнений узловых напряжений. Изменяется лишь правая часть уравнений, определяемая вектором узловых токов, который обра- зуется токами источников. Последние зависят от значе- ния узловых напряжений в конце предыдущего интер- вала '(шага), а также от времени (при наличии источ- ников внешнего сигнала). Поэтому вектор узловых токов должен вычисляться к началу анализа цепи на каждом очередном шаге. Итак, дискретная схема замещения линейной дина- мической цепи состоит из постоянных резистивных эле- ментов и источников, токи (или напряжения) которых изменяются от шага к шагу. Анализ дискретных рези- стивных схем удобно проводить по уравнениям узловых напряжений, которые можно записать в матричной форме: Gyiiy=iy, (2.5) где Gy—матрица узловых проводимостей; iy— вектор узловых токов. Для цепей с большим числом узлов и ветвей уравне- ния формируются с помощью ЭВМ. Для этого можно применить алгоритм и программу, описанные в гл. 1. Выбор шага расчета h при построении дискретных схем замещения является в общем случае достаточно слож- ной проблемой, особенно при исследовании жестких си- 2—252 17 стём с большим разбросом собственных чисел матрицы узловых напряжений. Рассмотрение этого вопроса выхо- дит за рамки книги. Описание возникающих здесь труд- ностей и путей их преодоления можно найти, например, в [24]. 2.2. ПОЯСНЕНИЯ К ПРОГРАММЕ АНАЛИЗА ЭЛЕКТРИЧЕСКИХ ЦЕПЕЙ ВО ВРЕМЕННОЙ ОБЛАСТИ С ИСПОЛЬЗОВАНИЕМ ДИСКРЕТНЫХ МОДЕЛЕЙ НА ОСНОВЕ НЕЯВНОГО АЛГОРИТМА ЭЙЛЕРА Приведем пояснения к подпрограмме DIS, которая позволяет рассчитать переходный процесс во временной области при воздействиях i(t)=Asmi, ступенчатом и в виде прямоугольного импульса тока. Перед подготов- кой входных данных должна быть произведена нумера- ция узлов, последний номер должен иметь опорный узел. При гармоническом воздействии необходимо про- нормировать параметры цепи по частоте (угловая часто- та воздействия должна быть равна 1). Входными данными подпрограммы являются масси- вы R, L, С, дающие информацию о топологии и пара- метрах цепи, массивы U и IJ4, дающие информацию о воздействующих функциях и точках их подключения, число узлов N и число независимых узлов Nl==N—1, шаг расчета Н и время расчета Т. Массивы R, L, С за- полняются одинаково, для характеристики соответству- ющего элемента и точек его подключения используются три ячейки: в первых двух указываются номера узлов, в третьей—параметр в омах, генри, фарадах соответст- венно. Для характеристики воздействия в массивах U и TJ4 также используются три ячейки: в первых двух указаны номера узлов (номер узла, от которого направ- лен ток, указывается первым), в каждой третьей ячейке массивов U и U4 записываются соответственно амплиту- да тока или его постоянное значение (массив U) и 0 или длительность импульса тока (массив U4). Ввод данных в программе 2.1 осуществляется с по- мощью оператора NAMLIST. В качестве примера при- водится карта входных данных для задачи 2.2: ^=1,3,1,2,3.100; L=1,3,10; 0=1,2,20,2,3,20; U=3,\,2; [/4=3,1,5; Г=5; Я=0,5. Подпрограмма DIS выводит на печать массив узло- вых напряжений схемы (массив U2). 18 ПРОГРАММА 2 Л REAL a(i50)/i50»0.0/,L(150)A50*0.&/,C(l50)/150»0.0/, *UC30)/30»0.0/,U4(30)/30»Q.O/,Z1(30,30),M100),M2(30) NAMELIST/D/R,I.,C,U,H,-U4 ,M,T BEAD(5,D) WRITE(b,D) Ш^-З. CALL MS(Zl.,Ml,M2,Ml,R,b,C)U,H,WiNiT) STOP END SUBROUTINE DIS(Z1 ,M1 ,M2 ,N1 ,R,Ii,C,U,H,U4 ,M.T) DIMENSION U(30),M(30),Z1(N1,NI),M1(N1),M2(N1.) BEAL R(i50) ,L(150) ,C(150) ,Ul(30)/30»t>.0/,U2(30)/30*0.0/, »U3(30)/30*0.0/,Z(30,30)/900*0.0/ K4=0 DO 1 1=1,150,3 IF((R(I)+C(I)+L(I)).EQ.O,0) GO TO 6 DO 1 J=1,3 IF(J-2)2,3,4 2 R1=R(I+2) IF(Rl.EQ.O) GO TO 1 B1=1/R1 . . ' K=K(I) K1=R(I+1) 00 TO 5 3 B1=L(I+2) .- ' IF(Rl.EQ.O) GO TO 1 : B1=H/R1 K=L(I) K1=L(I+1) GO. TO 5 4 S1=C(I+2)/H IF(Rl.EQ.O) GO TO 1 K=C(I) Kl=0(1+1) 5 Z(K,K)=Z(K,K)+R1 Z(K1,K1)=Z(K1,K1)+R1 Z(K,K1)=Z(K,K1)-R1 Z(K1,K)=Z(K1,K)-R1 1 CONTINUE 6 N=N-1 DO 32 J=1,N DO 32 1=1 ,N 32 Zl(I,J)=Z(I,J) CALL MIMV(Z1,N,D,M1,M2) 10 DO 8 1=1,10 J=I*3 IF(U(J).EQ.O) GO TO 29 K=U(J-2) K1.=U(J-1) K3=U4(J) IP(K3)12,11,12 11 R1=U(J)*SIN(K4»H) GO TO 9 12 IF(K4*H-K3)l4,14,13 13 Rl=0 ' GO TO 9 2* 14 R1=U(J) 9 U1(K)=U1(K)-R1 UI(K1)=U1(K1)+R1 8 CONTINUE 29 DO 7 1=1,150,3 IF((L(I).+C(I)).BQ.O) GO TO 15 K=L(I) K1=L(I+1) J=C(.I) K3=C(I+1) IF(K.EQ.O) GO TO 16 R1=(U2(K)-U2(K1))»H/L(I+2)+U3(I) U3(I)=H1 •U1(K)=U1(K)-R1 U1(K1)=U1(K1)+R1 16 IF(J.EQ.O) GO TO 7 R1=(U2(J)-U2(K3))*C(I+2)/H •Ul(J)=Ul(J)+Rl •U1(K3)=U1(K3)-RI 7 CONTIOTJE 15 K4=K4+1 K=T/H IF(K+,OT.K) GO TO 17 TO 20 X=l,M •S=0.0 DO 21 Jsl.N 21. S=Z1(I,J)»UI(J)+S 20 TJZ(I)=S TO 23 Xel,M 23 Ut(I)=Q,0 WBI'CB(6,305(V2(I),I=l,K) 30 FORMAT(8F9.4) GO 'TO 10 17 BETURK END Ограничение для программы—воздействующие функ- ции должны быть функциями тока. Таким образом, источники напряжения должны быть преобразованы в источники тока. Если ветвь с источником напряжения не содержит сопротивления, то в нее следует включить очень малое сопротивление и выполнить преобразование, получив дополнительный узел с неизвестным напряже- нием. Например, в задаче 2.4 в ветвь с источником на- пряжения достаточно включить сопротивление 0,001 Ом. Правильность выбора этой величины подтверждает рас- чет: рассчитанное напряжение узла должно быть равно заданному напряжению источника. 2.3. АНАЛИЗ ЛИНЕЙНЫХ РЕЗИСТИВНЫХ ЦЕПЕЙ ПРИ ПЛОХО ОБУСЛОВЛЕННЫХ МАТРИЦАХ Анализ линейных резистивных цепей представляет интерес как один из этапов исследования описанных вы- ше дискретных схем замещения, а также нелинейных 20 резистивных и инерционных цепей. Линейные резистив- ные цепи описываются системой линейных алгебраиче- ских уравнений, для численного решения которых в на- стоящее время разработано много методов [4, 26]. Тем не менее эта проблема остается актуальной, так как усложнение схем приводит к увеличению числа уравне- ний системы и, как следствие, к ухудшению обусловлен- ности, т. е. к увеличению разброса абсолютных значе- ний собственных чисел матрицы линейной системы. При анализе же плохо обусловленных систем возникают вы- числительные трудности, связанные с тем, что при обра- щении плохо обусловленных матриц ошибки расчета могут стать недопустимо большими. Методы, не требу- ющие обращения матриц (итерационные методы), схо- дятся лишь при довольно жестких условиях, налагаемых на матрицу системы. В данном параграфе представлен метод, относящийся к итерационным и не требующий обращения матриц и в то же время налагающий значительно более слабые условия на матрицу системы. Рассмотрим систему линейных алгебраических урав- нений, описывающих линейную резистивную цепь Ах+В=0, где А—квадратная неособая матрица; В (2.6) постоянная матрица-столбец. Самым простым методом решения такой системы яв- ляется метод простой итерации ^(А+^-Ч+В, где х^)—k-я итерация решения; I—единичная матрица. Однако условия сходимости метода очень жестки. Для сходимости итераций к решению необходимо и до- статочно, чтобы все собственные числа матрицы (A+I) были по модулю меньше единицы [26]. Суть предлагаемого подхода заключается в такой модификации метода простой итерации, которая позво- ляет существенно расширить область сходимости мето- да, а именно необходимо лишь, чтобы все собственные числа матрицы А лежали строго в левой полуплоскости. Для пассивных линейных резистивных цепей это усло- вие всегда выполняется, если, например, рассматривать систему (2.6) как систему уравнений узловых напря- жений. 21 (2.7) Рассмотрим вначале вспомогательную систему ли- нейных дифференциальных уравнений х=Ах+В. (2.8) Как известно из теории линейных дифференциальных уравнений, система (2.8) дает решение, сходящееся к решению уравнения (2.6) при t—>-oo, в том и только в том случае, когда все собственные числа матрицы А лежат строго в левой полуплоскости. Нетрудно видеть, что уравнение (2.?) представляет собой не что иное, как k-ю итерацию решения уравне- ния (2.8) с помощью явного метода Эйлера с шагом h, равным единице [2]. К сожалению, при достаточно ма- лых собственных числах матрицы метод Эйлера, как известно [2], дает расходящиеся результаты. Однако, как показано далее, и расходящееся решение несет в се- бе полную информацию о точном решении системы (2.8). Рассмотрим систему линейных дифференциальных уравнений (2.9) где А—квадратная постоянная матрица гаХ"; i(t}— заданная вектор-функция. Формула Эйлера численного интегрирования систе- мы (2.9) x(й)==(/гA+I)x(й-l)+/гf(/г-l), (2.10) где h—шаг расчета; х^—k-я итерация, дающая при- ближенное значение неизвестной x(t) в точке t=kh; f при t=kh; /г=0, 1,...; [ 0 при остальных t, где xW определяется формулой (2.10). Очевидно, что составляющие х^) являются импульс- ными функциями. Применяя к (2.12) и (2.13) преобразование Лапласа, получаем соответственно: /ЗД=АХ(/0+Р, (2.14) X^^e-^A+^X^+fte-^F, (2.15) где Х(р), Хэ(р) и F—изображения по Лапласу соответ- ственно вектор-функций \(t), Хэ(<) и fo(0. Таким образом, F=[^oi, xo2, ..., ^о"]7. (2.16) Решая (2.14) и (2.15), получаем соответственно: (2.17) (2.18) (2.19) Х(р) ==[/?! i-ePft- I - А F. ад=[—г Сравнивая (2.17) и (2.18), видим, что ад=х(^1). Пользуясь свойствами преобразования Лапласа, можно на основании (2.19) установить связь между х(<) и Хэ(0. Пусть в результате решения уравнения (2.13) мы определим Хэ(0. Таким образом, в силу (2.19) X^-^-)-.^^). (2.20) Пусть х(р)=[х,(р); Х2(р), ..., Xn(p)V; Хэ(р)=[Хэ1(р); Хэ2(р), ..., Х,п(р)Г; (2.21) Xh(p)-^fc(0; Xsh{p)^x^(t), k=\, 2, ..., п. 23 Следовательно, из (2.26) получаем xk[ep^гl•}—^(t),k==l.2....,n. (2.22) Будем последовательно применять следующие пре- образования: 1. Заменим в левой части (2.22) р на p/h. По теоре- ме подобия [10] получим (2.23) 2. Из теории преобразования Лапласа [10] известно следующее соотношение: если Ф (/»)-> у (О, то со 0(ln/^J-^-y(,)^, (2.24) где Г (т) —гамма-функция. Заменим в левой части (2.23) р на \пр и применим (2.24): fP-l \ _. Г ^-1 ^ h ) J~^ x,flL-L\^?_ 3. Заменим в левой части (2.25) р на /5+1. По тео- реме смещения [10] получим (2.25) оо x^-i•)-^^'j-^^lx„(Ь)hd.. (2.26) 4. Заменим в левой части р на р+1. По теореме по- добия получим оо ^ (Р) -> ^1(0 = е-^ J -^д ^ (A,) d,. (2.27) Эта формула и дает требуемое соотношение между Xk(t) И X3k(t). В (2.27) можно перейти от интеграла к степенному ряду. Для этого воспользуемся тем, что для принятого 24 вида вектор-функции t(t} решение уравнения (2.13) представляет собой «частокол» из импульсных функций: 1"*^ x,k{t}-= Va^(t-hi}, k- со У а^.8 (t - hi}, k=l, 2,..., п, (2.28) t=i где a.ki — постоянные коэффициенты. Суммирование в (2.28) берется от г'=1, а не от г'===0, так как по пред- положению Xyk (0) =0. Если подставить (2.28) в (2.27), то, учитывая свой- ства импульсной функции, получаем aki (l-i)\ о ' ' г=1 Известно, что если (2.28) есть решение уравнения (2.13), то коэффициенты a&, растут не быстрее, чем по экспоненциальному закону. Это свойство позволяет на основе теории обобщенных функций заключить, что ряд, стоящий в последнем выражении под знаком интеграла, можно интегрировать почленно. В результате, пользуясь свойствами интегралов, содержащих импульсные функ- ции, получаем ^O—r^ST^^)'". <2-29) г=1 В силу указанного выше свойства коэффициентов а-ki ряд (2.29) сходится даже в том случае, когда результа- ты эйлеровских итераций являются расходящимися. Пусть теперь в (2.9) f(f)—постоянный вектор и на- чальные условия—нулевые. Таким образом, (2.30) i(t)=[b„ 6„...ДГ; F(p)=[b,fp, b,Jp,...,b,]pY; /еР'«. (pp;i ^ k-я составляющая вектор-функции F —т, ляется формулой W(eP'l—l)==6ftfte-PV(l—e--P'l), k=\, 2, ..., п опреде- и ее оригинал ^hb^(t -hi), k==l, 2, ... , п. Ы 25 Следовательно, в правую часть (2.13) вместо i,{t—h) подставляется вектор-функция [s hb,8 (t - hi), | /^8 [t - АО, ..., S ^,8 (/ - hi)]. (2.31) L =0 (=0 i=0 ] Затем, решая (2.13), применяем формулу (2.27) или (2.29). Если собственные числа матрицы А лежат стро- го в левой полуплоскости, то lmx{t) дает реше- f-»00 ние уравнения (2.6). Как указывалось выше, линейные пассивные рези- стивные цепи обладают тем свойством, что собственные числа матрицы А в (2.6) лежат строго в левой полупло- скости. Ограничение на расположение собственных чисел матрицы А можно снять с помощью предварительных преобразований. Пусть в уравнении Ах+В==0 матрица А (вообще говоря, комплексная) имеет собст- венные числа, расположенные в произвольных точках комплексной плоскости, за исключением начала коорди- нат. Умножим это уравнение слева на A'r*—транспони- рованную, сопряженную матрицу: А^Ах+С^, С^А^В. Матрица А^А является эрмитовой положительно опре- деленной, и все ее собственные числа вещественны и по- ложительны. Следовательно, необходимо рассматривать вспомогательное дифференциальное уравнение в виде х == - [А^А] х - С. Алгоритм, использованный при решении систем линей- ных алгебраических уравнений предложенным методом, имеет следующие особенности: 1. От системы '(2.6) делаем переход к системе (2,8), которая решается методом Эйлера. 2. Не следует выбирать шаг расчета h сразу же слишком большим, так как расходимость метода Эйлера в этом случае может быть столь сильной, что либо про- изойдет переполнение разрядной сетки, либо сходимость ряда (2.29) будет очень медленной. В то же время прак- тика расчетов показала, что шаг эйлеровских итераций 26 h удобно выбирать таким, чтобы он был значительно больше (на порядок) шага во времени при расчете Xk(t) по формуле (2.29). Более конкретные указания по вы- бору шага приведены в комментариях к программе, при- веденной в конце параграфа. 3. Найденные итерации подставляем в (2.29), огра- ничиваясь таким числом итераций, чтобы слагаемое ря- да (2.29), соответствующее первой отброшенной итера- ции, было по абсолютной величине в е раз меньше модуля предыдущей частичной суммы ряда при значе- нии переменной /, определяемом выбранным шагом по аргументу t (шаг расчета по аргументу обозначим че- рез А^). Практика показала, что обычно достаточно взять e=10-5. 4. Значение частичной суммы ряда (2.29) при t=^t принимаем в качестве новых начальных условий и по- вторяем весь расчет, определяя значение вектор-функ- ции \(t) при t=2At и т. д. 5. Расчет заканчивается, когда значения составляю- щих вектора х, найденные в результате применения двух последовательных циклов расчета по пп. 1—4, от- личаются друг от друга не более чем на 0,001%. 2.4. ПОЯСНЕНИЯ К ПРОГРАММЕ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Программа 2.2 SUBROUTINE, PROG31 (N, EPS, PSI, TO, H, HI, A, U, G, D) написана на языке ФОРТРАН и предназначена для решения системы ли- нейных алгебраических уравнений Ах+В==0, где А—квадратная матрица п-го порядка; В—вектор- столбец свободных членов системы; х—вектор-столбец неизвестных. На самом деле программа решает систему линейных дифференциальных уравнений вида х = Ax + В, решение которой совпадает на бесконечности с реше- нием системы Ах+В=0 "в том случае, если собственные числа матрицы А лежат в левой комплексной полуплоскости. Частичная сумма 27 •ПРОГРАММА 2.2 SUBROUTINE PBOG3'l(N,EPS,PSI,'rO ,H,H1 ,A,TT,&,D) KMBNSION .A(N,N)-,B(N) ,VW1 »У(К) , »G(N);Z(N,2), *C(H,2).R(H,2);, *Xj, <=/'+ 1, ... , П; k=\ I / п; (2.32; d, --a, ^Cikdk,, i