Численное интегрирование дифференциальных уравнений. Задачи электротехники, приводящие к решению дифференциальных уравнений.

ЛЕКЦИЯ №4 Численное интегрирование дифференциальных уравнений Задачи электротехники, приводящие к решению дифференциальных уравнений Обыкновенные дифференциальные уравнения (ОДУ) используются во всех областях электротехники чрезвычайно широко. Это объясняется тем, что большинство основных соотношений между электрическими и магнитными величинами устанавливают связь не между самими величинами, а скоростью их изменения во времени и пространстве. Как известно, физический смысл первой производной функции f(x) - - это скорость изменения функции в каждой ее точке. При написании уравнения математической модели электромагнита (глава 1, формула (1.3)) мы воспользовались соотношением . Это значит, что напряжение на катушке индуктивности зависит от скорости изменения тока, и чем она выше, тем напряжение больше. Здесь неизвестный ток входит в уравнение в виде своей производной. В электромеханике, кроме уравнений, связывающих только электрические величины, есть уравнения, связывающие электрические и механические величины, которые тоже являются дифференциальными, так как в основе их лежит второй закон Ньютона: , где ускорение, как производная от скорости по времени; F - механическая сила электрического происхождения, например, от взаимодействия двух токов. Как составляется дифференциальное уравнение простой электрической цепи, описано в главе 1. Приведем пример составления системы дифференциальных уравнений более сложного устройства - трансформатора, имеющего две обмотки, намотанные на сердечник из электротехнической стали, по которому замыкается магнитный поток, пронизывающий обе обмотки. Этот магнитный поток называется потоком взаимоиндукции, а параметр трансформатора, связанный с ним - коэффициентом взаимоиндукции М (рис.4.1). По второму закону Кирхгофа для каждой из обмоток трансформатора можно записать уравнения равновесия напряжений: , (4.1) где R1 и R2 - активные сопротивления обмоток, L1 и L2 - индуктивности обмоток, М - коэффициент взаимоиндукции обмоток, u1 и u2 – напряжения на обмотках. Обозначив производные токов обмоток и , перенеся слагаемые без производных в правую часть, упорядочивая, получим: (4.2) или в матричной форме , где ; ; ; ; . Умножив обе части на обратную матрицу (зная, что - единичная матрица, имеющая тот же смысл в матричных операциях, что и единица), получим . (4.3) Мы получили в матричной форме систему дифференциальных уравнений трансформатора, записанную в форме Коши, то есть разрешенную относительно производных. Сравнив матричное уравнение (4.3) с уравнением (1.3), мы видим полное совпадение их структуры. Интегрирование системы исходных дифференциальных уравнений электротехнического или электромеханического устройства, как правило, происходит с применением численных методов. Задача Коши для обыкновенных дифференциальных уравнений В наиболее общем виде любое ОДУ n-го порядка можно записать так: (4.4) Из теории ОДУ известно, что уравнение (4.4) можно заменить системой из n уравнений первого порядка путем подстановки y(m)=ym [1]: (4.5) где к = 1,2,3...,n. Уравнения (4.4) и (4.5) имеют бесчисленное множество решений. Единственные решения выделяют с помощью дополнительных условий, которым должны удовлетворять искомые решения. В пособии рассматриваются задачи с начальными условиями. Для таких задач кроме исходного уравнения (4.4) в некоторой точке x0 , должны быть заданы начальные условия, т.е. значения функции y(x) и ее производных в начальной точке: Для системы ОДУ начальные условия задаются в виде: . (4.6) Задача с начальными условиями носит название задачи Коши. К численному решению ОДУ приходится обращаться, когда не удается найти аналитическое решение через известные функции. Для некоторых задач численные методы оказываются более эффективными даже при наличии аналитических методов. В пособии рассматриваются некоторые численные методы решения задачи Коши. Выбор методов проведен так, чтобы дать понятие о явных и неявных, одношаговых и многошаговых методах, их достоинствах и недостатках. Методы Эйлера Формой Коши для системы ОДУ является форма, в которой уравнения системы разрешены относительно производных: (4.7) где k = 1,2,3,....,n. Система (4.7) для решения задачи Коши дополняется начальными условиями. В пособии для упрощения выкладок будем решать задачу Коши для одного уравнения: . (4.8) Из теории рядов известно, что любую функцию y(x) в малой окрестности точки х0 можно разложить в ряд по степеням х, который носит название ряда Тейлора: (4.9) Ограничив ряд при малом шаге двумя членами, получим (4.10) где O(h2) - бесконечно малая величина порядка h2. Заменим производную входящую в формулу (4.10), на правую часть уравнения (4.8). Отбрасывая бесконечно малую О(h2), получим соотношение, позволяющее найти приближенное решение дифференциального уравнения в точке , зная значение подынтегральной функции в точке , величину производной в данной точке и значение шага интегрирования h: (4.11) Рассчитав значение функции в точке , принимаем его за начальное для следующего шага и повторяем расчет. Для любого i-го шага интегрирования будем иметь: или (4.12) Такая формула носит название рекуррентной. Метод численного интегрирования, использующий два члена разложения ряда Тейлора, носит название явного метода Эйлера. Явным этот и другие подобные методы называются потому, что в уравнении используется уже известное значение производной - в точке xi , то есть производная задается явно. При этом решение находится за конечное, заранее определенное число шагов. Как известно, уравнение прямой имеет вид y=ax+b, где а - тангенс угла между осью абцисс и прямой, b - значение y при x, равном нулю. Также известно, что геометрический смысл производной - величина производной f’(x) равна тангенсу угла наклона касательной, проведенной к графику кривой функции в данной точке. Принимая за аргумент шаг h , можно построить на графике прямую линию из точки x0,y0 в точку x1,y1 под углом, тангенс которого равен производной в точке x0 - f(x0,y0) (рис.4.2). Из точки х1 под углом, тангенс которого равен f(x1,y1), проводим прямую в точку х2 и т. д. В качестве решения получим ломаную линию. В связи с этим метод Эйлера еще носит название метода ломаных. Подставив в формулу (4.12) вместо производной в исходной точке f(xi,yi) значение еще неизвестной производной в искомой точке xi+h - f(xi+1,yi+1), получим формулу неявного метода Эйлера: (4.13) Следовательно, неявным называется такой метод, в исходные уравнения которого входят производные еще не рассчитанных точек интегрируемой функции. В общем случае решение уравнений вида (4.13) ищется итерационными методами со свойственными им недостатками (зацикливание, расхождение процесса). Однако, если производную удается представить в линейном виде (линеаризовать) , (4.14) что в задачах электротехники, как правило, имеет место (см. уравнение (1.3) в главе 1 - ), от итерационного метода можно перейти к прямому методу интегрирования неявным методом Эйлера. Подставив выражение для производной из (4.14) в (4.13), получим: . Разрешив данное уравнение относительно yi+1 ,получим прямую формулу для интегрирования ОДУ неявным методом Эйлера: . (4.15) Геометрическая интерпретация неявного метода Эйлера изображена на рис.4.3. Неявные методы обладают большой устойчивостью при интегрировании жестких ОДУ. Они позволяют проводить интегрирование при больших значениях шага интегрирования h , чем явные методы. Более подробно об этом будет сказано ниже. На каждом шаге метода Эйлера решение y(x) определяется с погрешностью за счет отбрасывания членов ряда Тейлора, пропорциональных h в степени выше первой. Это означает, что локальная погрешность (погрешность на шаге интегрирования) метода Эйлера имеет второй порядок. Глобальная погрешность (погрешность на всем интервале интегрирования) имеет первый порядок. При постоянном шаге h для оценки погрешности применима формула Рунге , (4.16) где yh(x) - приближенное решение ОДУ в точке х, полученное с шагом h; ykh(x) - приближенное решение того же уравнения, полученное с шагом kh; p - порядок метода. Процедуры, реализующие явный и неявный методы Эйлера, приведены в конце главы (ПРОГРАММЫ 4.1 и 4.2). Сами процедуры дают решение на шаге интегрирования. Для поиска решения на всем интервале изменения аргумента в вызывающей программе должен быть организовано циклическое изменение аргумента х. Как видно, метод Эйлера как в явной, так и в неявной форме достаточно просто реализуется на ЭВМ. Коэффициенты А и В неявного метода рассчитываются в вызывающей программе. Для повышения точности решения разработаны одношаговые методы, позволяющие учесть большее количество членов ряда Тейлора, чем в методе Эйлера. Методов достаточно много [9]. Рассмотрим некоторые из них. 4.4. Методы Рунге-Кутты В зависимости от старшей степени h, с которой учитываются члены ряда Тейлора, построены вычислительные схемы Рунге-Кутты разных порядков точности. Рассмотрим методы второго порядка. Для них получено однопараметрическое семейство схем вида [5]. где - свободный параметр; . Локальная погрешность данной схемы имеет третий порядок, глобальная - второй, т.е. решение ОДУ полученное по этой схеме, равномерно сходится к решению с погрешностью О(h2) . Для параметра L наиболее часто используются значения L=0.5 и L=1. В первом случае формула (4.17) приобретает вид: , (4.17) во втором: . (4.18) Геометрическая интерпретация метода в обоих вариантах представлена на рис 4.4.В случае а) сначала делается прогноз по методу Эйлера (точка yэ), затем рассчитывается значение производной в этой точке, берется полусумма этих значений и из начальной точки проводится прямая под полученным углом. Ее пересечение с прямой x0+h дает искомый результат y1. В случае б) сначала находится значение функции на половине шага по методу Эйлера yh/2. Далее из начальной точки проводится прямая под полученным углом. Ее пересечение с прямой x0+h дает искомый результат y1. Процедура, реализующая приведенные методы, приведена ниже (ПРОГРАММА 4.3). Как и процедура, реализующая явный метод Эйлера, процедура R_Kutt2 имеет дальний способ вызова Far , что необходимо для реализации вызова функции расчета производной как параметра. С точки зрения точности решения обе приведенные схемы эквивалентны. На рис 4.5 изображены абсолютные погрешности (по отношению к точному аналитическому решению) явного метода Эйлера и метода Рунге - Кутты второго порядка. Как видно из графика на рис.4.5, точность метода второго порядка при одном и том же шаге аргумента значительно выше, чем первого. Но при этом нужно провести больше расчетов: три вызова функции F при L=0.5 вместо одного в методе первого порядка. На практических занятиях предлагается оценить соотношение ‘’величина шага - число операций” для данных методов и сделать соответствующие выводы. Наряду с одношаговыми методами в электротехнике применяют и многошаговые методы численного интегрирования ОДУ [8,9]. Многошаговые методы интегрирования В многошаговых методах для получения решения дифференциального уравнения используются результаты нескольких предыдущих шагов интегрирования путем использования различных алгоритмов экстраполяции (extra - вне , pole - узел). Многошаговые методы позволяют сократить вычисления за счет использования результатов расчета предыдущих точек, но их недостаток в том, что они требуют «разгонки», т.е. первые шаги необходимо делать одним из одношаговых методов, что усложняет алгоритм. Другим недостатком многошаговых методов может служить их относительно меньшая устойчивость при интегрировании ОДУ, описывающих быстропеременные процессы, о чем будет сказано ниже. В качестве примера рассмотрим современный многошаговый неявный метод ФДН, не требующий привлечения одношаговых методов для «разгонки», получивший распространение при расчете математических моделей в электротехнике, особенно в электромеханике при расчете процессов в электрических машинах. Метод ФДН получил название от формулы записи аппроксимирующего интегрируемую функцию многочлена [7]: , (4.19) где r - порядок многочлена. Соотношение (4.19) получило название формулы дифференцирования назад - ФДН. Как видно, в соотношении (4.19) присутствует производная искомой точки , т.е. метод ФДН является неявным. Метод ФДН, как и другие неявные методы, достаточно сложен при программной реализации в сравнении с явными, требует больше оперативной памяти ЭВМ [8]. Как и все неявные методы, ФДН в общем случае для нахождения решения требует применения итерационных процедур с присущими им свойствами. Но, если производную удается представить в линейном виде , можно перейти к прямым способам поиска коэффициентов, входящих в формулы метода. Для многих задач электротехники это, как правило, имеет место (см. параграф 1.1, формулу (1.3) и параграф 4.1, формулу (4.3)). Ниже приводятся основные соотношения, с помощью которых реализуется метод ФДН и его программная реализация для одного уравнения (ПРОГРАММА 4.3). Порядок метода r зависит от количества предшествующих рассчитываемой точке точек (узлов), информация о которых используется при расчете (рис 4.6). Если r=1 , то мы имеем одношаговый метод, известный как неявный метод Эйлера (см. параграф 4.3). Это может служить проверкой правильности программной реализации алгоритма метода ФДН. Расчетная формула метода ФДН имеет вид , (4.20) где , а r+1 коэффициентов вычисляются из системы уравнений следующего вида: . (4.21) Система (4.21) относительно неизвестных решается методом Гаусса. Ниже (ПРОГРАММА 4.4) приведены процедуры расчета коэффициентов и самого метода ФДН. При "разгонке" ФДН не нужно прибегать к другим методам интегрирования, нужно только, увеличивая порядок метода от 1 до принятого, пересчитывать коэффициенты . Пересчет их нужно вести и при изменении шага интегрирования h . Если и то и другое неизменно, расчет ведется по формуле (4.20) при постоянных параметрах. Как и другие методы численного интегрирования, ФДН легко обобщается на системы дифференциальных уравнений. Возможность оценки погрешности на шаге интегрирования (см. формулу(4.16)) позволяет автоматизировать выбор шага интегрирования как для явных, так и для неявных методов, и тем самым сократить время решения ОДУ [8]. Как говорилось ранее, неявные методы обладают большей устойчивостью при решении жестких уравнений. Дадим наглядное понятие о жестких ОДУ и об устойчивости решения ОДУ. Понятие о жесткости и устойчивости систем дифференциальных уравнений Переходные процессы (см. главу 2) в электрических и электромеханических устройствах описываются с помощью систем дифференциальных уравнений. Их порядок определяется сложностью системы. В частном случае линейной электрической цепи, состоящей из активных сопротивлений и индуктивностей (R-L цепь), после исключения токов других ветвей, дифференциальное уравнение для тока k-й ветви будет иметь вид [4]: , (4.22) где n - число ветвей. Уравнение (1.3) - его частный случай при n = 1. Система уравнений трансформатора (4.3) также сводится к уравнению (4.22) , к уравнению второй степени, в чем можно убедиться самостоятельно, используя прием повторного дифференцирования [9]. Из курса математики известно , что полный интеграл дифференциального уравнения (4.22) состоит из частного решения, определяемого видом функции F(t), и полного решения однородного уравнения: . (4.23) Для тока ветви ik это обозначает, что он состоит из двух составляющих - , где называется принужденной составляющей, определяемой ЭДС после окончания переходного процесса, а - свободная составляющая, затухающая (по экспоненте для R-L цепи) во время переходного процесса. В уравнении модели катушки электромагнита принужденная составляющая , а свободная составляющая , где - постоянная времени цепи. При числе ветвей n свободная составляющая имеет вид: . (4.24) Свободные токи других ветвей также будут состоять из суммы экспонент с различной скоростью затухания, определяемой постоянными времени при другой величине и знаках коэффициентов А. При численном интегрировании уравнения (4.22) шаг h выбирают таким, чтобы он был соизмерим с наименьшей постоянной времени . Для явного метода Эйлера это определяется следующим ограничением на шаг [6]: . (4.25) Если , шаг интегрирования будет определяться наименьшей постоянной времени, а время расчета - наибольшей. Такие уравнения называются жесткими. При численном интегрировании уравнение (4.22) заменяется системой из n уравнений первого порядка [1], и ее решение сводится к решению системы линейных алгебраических уравнений. Если исходная система дифференциальных уравнений жесткая, то матрица системы становится плохо обусловленной (см. пример в главе 1). Матрицу системы в связи с этим еще называют матрицей жесткости. Интегрирование такой системы при шаге больше граничного приводит к неустойчивости, быстрому накоплению погрешностей и аварийному завершению интегрирования. Примером жесткой системы дифференциальных уравнений могут служить уравнения трансформатора (4.2) c матрицей системы L, у которого , где L1 и L2 - индуктивности обмоток, М - коэффициент взаимоиндукции. Причем L1 L2 в реальном трансформаторе всегда больше М2. Соотношение имеет место в промышленных трансформаторах большой мощности. Легко убедиться, что определитель матрицы L в уравнении (4.2) равен L1L2-M2, то есть при он будет равен нулю и система не будет иметь решения. Все методы интегрирования по влиянию шага h на устойчивость процесса интегрирования делятся на условно устойчивые и абсолютно устойчивые [6]. Условно устойчивыми считаются методы, имеющие ограничение на максимальный шаг интегрирования h, абсолютно устойчивыми - не имеющие такого ограничения. Все явные методы относятся к условно устойчивым, среди неявных методов есть абсолютно устойчивые. Подводя итог, можно сказать, что жесткость уравнений математической модели является свойством самого объекта исследования (электрической цепи, электрической машины, трансформатора), его параметров. При этом принято говорить о коэффициентной устойчивости. Скорость протекания процессов в системе зависит также и от скорости изменения воздействующих на нее внешних сил. В случае модели электромагнита замка это будет частота ЭДС источника. При выборе шага интегрирования h следует иметь в виду, что шаг должен быть согласован как с параметрами математической модели, так и со скоростью изменения внешних сил, воздействующих на модель. Для электрических цепей ими являются ЭДС. Как отмечалось ранее, многошаговые методы при этом уступают одношаговым. Для иллюстрации того, как влияют параметры модели и скорость протекания процесса, обусловленная приложенной ЭДС, на устойчивость и погрешность интегрирования различными методами, были проведены расчеты на модели электромагнита замка. Результаты расчетов приведены ниже. Шаг интегрирования h при исследовании влияния параметров модели на устойчивость процесса интегрирования принимался одинаковым и малым по сравнению с периодом ЭДС (T=40h), чтобы погрешности интегрирования, связанные с изменением тока в результате действия ЭДС были минимальны. Для иллюстрации интегрирования быстро изменяющихся процессов многошаговыми и одношаговыми методами шаг выбирался соизмеримым с частотой питающей ЭДС, а постоянная времени - много большей, чем шаг (=2500h). На рис.4.7 и 4.8 приведены результаты решения уравнения модели катушки электромагнита явным и неявным методами Эйлера при различных постоянных времени цепи катушки электромагнита. При шаге h = 0.01(рис.4.7) оба метода дают одинаковую погрешность, но с разными знаками. При (рис.4.8) явный метод становится неустойчивым, а неявный дает результат, близкий к аналитическому. . При интегрировании более точными методами Рунге- Кутты и ФДН с шагом h = 0.01 получился результат практически совпадающий с аналитическим (рис.4.9,а и 4.10,а). При шаге расчет явным методом Рунге-Кутты ведется с возрастающей погрешностью, а неявным методом ФДН с той же точностью, что и с малым шагом по отношению к постоянной времени LH/RH (рис.4.9,б и 4.10,б). . Анализируя сказанное выше, можно сделать следующий вывод. Шаг интегрирования системы ОДУ должен быть выбран так, чтобы быть соизмеримым с постоянной времени процесса, скорость которого в моделируемом устройстве будет максимальной. Жесткой же называется система ОДУ, которая описывает одновременно протекающие динамические процессы, например в электрической цепи, скорости которых отличаются на несколько порядков. Скорость изменения интегрируемых величин - токов, напряжений и т.д. зависит не только от параметров объекта, но и от частоты действующих электродвижущих сил, и шаг интегрирования нужно выбирать, учитывая, что одношаговые методы при этих условиях точнее многошаговых. Примером могут служить графики на рис.4.11, где сравниваются методы Рунге - Кутты и ФДН. Как видно из рисунка, при одинаковых условиях погрешность расчета методом ФДН получилась значительно выше. Из сказанного выше можно сделать вывод о том, что не существует методов численного интегрирования, которые удовлетворяли бы всем условиям. Метод выбирается исходя из конкретной задачи, стоящей перед исследователем и возможностей ЭВМ (быстродействие, объем оперативной памяти). Для интегрирования уравнений систем, в которые входят электрические машины, следует отдавать предпочтение неявным методам.