Моделирование во временной области (расчет переходных процессов)

Модели и методы анализа проектных решений Лекция 10 – Моделирование во временной области (расчет переходных процессов) Введение Переходный процесс определяется как реакция РЭУ на входное воздействие (изменение тока или напряжения) во времени. В качестве входного воздействия используют как реальные сигналы, так и идеализированные воздействия (тестовые сигналы) единичный скачок или единичную дельта-функцию, позволяющие проводить сравнение реакций различных цепей. Так, реакция цепи на единичный скачок называется переходной характеристикой, а реакция на единичную дельта-функцию носит название импульсной характеристики. Причем переходная и импульсная характеристики определяются при отсутствии других воздействий на схему. Реакцию РЭУ (цепи) на реальный сигнал (периодический или нет) будем называть обобщенно переходным процессом. 2 Общие положения Методы анализа РЭУ во временной области Нелинейные цепи Линейные цепи спектральный подход (частотные методы) методы гармонического баланса X-параметры метод дискретных моделей прямые методы численного дифференцирования методы на основе рядов Вольтерра метод переменных состояний численное преобразование Лапласа прямые методы формирования ММ РЭУ 3 Метод переменных состояния для линейных РЭУ При анализе переходных процессов необходимо выбрать искомые переменные таким образом, чтобы их число было минимальным, но достаточным для полного анализа электромагнитных процессов, протекающих в цепи. В качестве таких переменных удобно выбрать переменные, характеризующие энергетический запас в электрической цепи. Будем называть энергетическое состояние цепи просто состоянием цепи, а переменные, характеризующие это состояние – переменными состояния. Переменные состояния – потокосцепление или ток через индуктивность, заряд либо напряжение на емкости. Если во время коммутации емкость и индуктивность не меняется, то в качестве переменных состояния берут ток индуктивности iL и напряжение емкости UC, в противном случае - потокосцепление индуктивности ψL и заряд емкости QC. Метод переменных состояний применим для расчета как линейных, так и нелинейных цепей, поэтому его часто используют для расчетов переходных процессов в современных САПР электронных цепей и устройств. 4 Метод переменных состояния для линейных РЭУ Уравнения, составленные в нормальной форме относительно этих переменных, называются уравнениями переменных состояния (слева стоят первые производные от каких-то функций, справа – функции реакций и возмущений с некоторыми постоянными коэффициентами). К примеру, для цепи второго порядка уравнения переменных состояния имеют вид ⎧ dx1 ⎪⎪ dt = a11 x1 + a12 x2 + b11e1 , ⎨ ⎪ dx2 = a21 x1 + a22 x2 + b21e2 , ⎪⎩ dt (1) где х1, х2 – реакции; e1, e2 – возмущения; а11, а12, а21, а22 – некоторые постоянные коэффициенты, зависящие от параметров и конфигурации цепи. Система уравнений (1) в матричном виде примет вид ⎡ dx1 ⎤ ⎢ dt ⎥ ⎡ e1 ⎤ ⎡ x1 ⎤ ⎢ dx ⎥ = A ⋅ ⎢ ⎥ + B ⋅ ⎢ ⎥, ⎣e2 ⎦ ⎣ x2 ⎦ ⎢ 2⎥ ⎣ dt ⎦ где ⎡ a11 a12 ⎤ A=⎢ ⎥ ⎣a21 a22 ⎦ ⎡ e1 ⎤ B=⎢ ⎥ ⎣e2 ⎦ – матрица коэффициентов переменных состояния; – вектор коэффициентов источников сигналов. 5 Метод переменных состояния для линейных РЭУ Для составления уравнений переменных состояния (1) необходимо: для цепи после коммутации составить уравнения по законам Кирхгофа, затем исключить токи и напряжения резистивных ветвей, т. е. выразить их через переменные состояния. После этого необходимо найти напряжение индуктивности и ток емкости через переменные состояния и записать их через соответствующие производные: u L (t ) = L dU C (t ) di (t ) ; iC (t ) = C dt dt Так поступают, если схема простая по топологии, если схема «сложная», то для составления уравнений необходимо воспользоваться графом с емкостным деревом. При составлении графа необходимо, чтобы ветви дерева содержали только емкости, активное сопротивление и источники, т.е. ветви дерева не должны содержать индуктивности. Затем составляют уравнения по законам Кирхгофа для главных сечений и главных контуров и опять исключают токи резистивных элементов (выбор емкостного дерева позволяет избежать математических трудностей при исключении токов и напряжений резистивных ветвей). 6 Метод переменных состояния для линейных РЭУ Составим уравнения переменных состояния для схемы показанной на рисунке Составим уравнения по законам Кирхгофа для схемы после коммутации ⎧− iL (t ) + i2 (t ) + iC (t ) = 0, ⎪ ⎨uL (t ) + i1 (t ) R1 + i2 (t ) R2 = E , ⎪u (t ) − i (t ) R = 0; 2 2 ⎩ C (2) здесь i1(t) = iL(t), uC(t) – переменные состояния, так как они характеризуют запасы энергии в этой цепи. Также i3(t) = iC(t). Определим из системы (2) напряжение на индуктивности и ток через емкость, выраженные через переменные состояния, чтобы найти уравнения переменных состояния. Из третьего уравнения системы (2) имеем i2 (t ) = uC (t ) R2 7 Метод переменных состояния для линейных РЭУ Из второго уравнения системы (2) u L (t ) = L или diL (t ) E – i (t) R – i (t) R = E – i (t) R – u (t) = L 1 2 2 L 1 C dt diL (t ) R 1 1 = − 1 iL (t ) − uC (t ) + E L L L dt Из первого уравнения системы (2) iC (t ) = C u (t ) duC (t ) = iL (t ) − i2 (t ) = iL (t ) − C , R2 dt 1 duC (t ) 1 = iL (t ) − uC (t ). dt C CR2 Запишем эти выражения в виде системы уравнений R1 1 1 ⎧ diL (t ) = − ( ) − ( ) + i t u t E, L C ⎪⎪ dt L L L ⎨ du (t ) 1 ⎪ C = iL (t ) − 1 uC (t ) + 0. ⎪⎩ dt C CR2 Обозначим вектор переменных состояния через x(t) = [iL(t), uC(t)]T. Вектор нереактивных токов и напряжений обозначим как y(t). 8 Метод переменных состояния для линейных РЭУ Тогда полную систему уравнений состояний для рассматриваемой цепи можно записать в виде ⎧ dx (t ) = A ⋅ x (t ) + B ⋅ e(t ); ⎪ ⎨ dt ⎪⎩ y (t ) = C ⋅ x(t ) + D ⋅ e(t ), (3) где С – матрица коэффициентов для резистивной части схемы; D – вектор коэффициентов источников тока и напряжений, входящих в схему. Для нашей схемы можно записать или в матричной форме где ⎡ R1 ⎢− L А=⎢ 1 ⎢ ⎣ C ⎡ iL (t ) ⎤ ⎡ e1 (t ) ⎤ ⎡ i1 (t ) ⎤ = ⋅ + ⋅ C D ⎢u (t )⎥ ⎢i (t )⎥ ⎢e (t )⎥, ⎣2 ⎦ ⎣ 2 ⎦ ⎣ C ⎦ 1 ⎤ ⎡1 ⎡1⎤ L ⎥; B = ⎢ ⎥ ; C = ⎢ L 1 ⎥ ⎢0 ⎢0⎥ − ⎥ ⎣ ⎦ ⎣ CR2 ⎦ − ⎧i1 (t ) = 1 ⋅ iL (t ) + 0 ⋅ uC (t ), ⎪ 1 ⎨ ( ) = ⋅ ( ) + ⋅ uC (t ), i t i t L 2 ⎪ R ⎩ 2 0⎤ 1 ⎥; D = ⎡ ⎤. ⎢0 ⎥ ⎣ ⎦ R2 ⎥⎦ 9 Метод переменных состояния для линейных РЭУ Система уравнений вида (3) может быть решена: 1. Аналитически; 2. Численно, c помощью компьютера. Решение систем диф. уравнений В зависимости от числа независимых переменных дифференциальные уравнения делятся на две существенно различные категории: обыкновенные дифференциальные уравнения (как правило, t - время), содержащие одну независимую переменную, и уравнения в частных производных, содержащие несколько независимых переменных (чаще всего, x, y, z – пространственные координаты и t-время). Далее мы рассмотрим численные методы решения систем обыкновенных дифференциальных уравнений (задачи Коши) 10 Решение систем ОДУ Обыкновенным дифференциальным уравнением называют такое уравнение, которое содержит одну или несколько производных от искомой функций у = у(х). Его можно записать в виде F (x, y, y', y'‘, ... , y(n)) = 0, где х — независимая переменная. Наивысший порядок входящей в это уравнение производной называется порядком дифференциального уравнения. В частности, запишем уравнения первого: F (x, y, y' ) = 0 и второго порядков: F(x, y, y', y'' ) = 0. Решением дифференциального уравнения называется всякая n раз дифференцируемая функция у = ϕ(x), которая после ее подстановки в исходное уравнение превращает его в тождество. 11 Решение систем ОДУ Общее решение обыкновенного дифференциального уравнения n-го порядка содержит n произвольных постоянных С1, С2, ... , Сn: у = ϕ (x, С1,С2, ... ,Сn), Эта функция является решением диф. уравнения при любых значениях С1, С2, ... , Сn. Таким образом, частное решение дифференциального уравнения получается из общего, если этим произвольным постоянным задать определенные (конкретные) значения. Для уравнения первого порядка общее решение зависит только от одной произвольной постоянной: у = ϕ (x, С), Если постоянная в этом выражении принимает определенное значение С = С0, то получаем частное решение диф. уравнения вида у = ϕ (x, С0). В зависимости от способа задания дополнительных условий для получения частного решения дифференциального уравнения существуют два различных типа задач: задача Коши и краевая задача. 12 Решение систем ОДУ Если эти условия задаются в одной точке, то такая задача называется задачей Коши. Дополнительные условия в задаче Коши называются начальными условиями, а точка х = x0, в которой они задаются – начальной точкой. В общем виде задача Коши для дифференциального уравнения 1-го порядка выглядит следующим образом: dy = f ( x, y ) dx y(x0) = y0 Примеры задачи Коши: dy/dx = y2 cos x при x > 0, y(0) = 1; у'' = у'/х + х2 при х > 1, y(1) = 2, y' (1) = 0. 13 Решение систем ОДУ Если же для уравнения порядка n > 1 дополнительные условия задаются в более чем одной точке, т. е. при разных значениях независимой переменной, то такая задача называется краевой. Сами дополнительные условия называются при этом граничными (или краевыми) условиями. На практике обычно граничные условия задаются в двух точках х = а и х = b, являющихся границами отрезка, на котором рассматривается решение дифференциального уравнения. Постановка краевой задачи для дифференциального уравнения 1-го порядка имеет следующий вид: dy = f ( x, y ) dx y(a) = ya, y(b) = yb. Примеры краевых задач: у''+ 2у' – у = sin х, при 0 ≤ х ≤ 1, y(0) = 1, y(1) = 0; у"' = х + у у', при 1 ≤ х ≤ 3, y(1) = 0, y' (1) = 1, y(0) = 1, y' (3) = 2. 14 Решение систем ОДУ Методы решения обыкновенных дифференциальных уравнений можно разбить на следующие группы: - графические, - аналитические, - приближенные, - численные. Далее мы будем рассматривать только численные методы решения дифференциальных уравнений, которые в настоящее время являются основным инструментом при исследовании научно-технических задач, описываемых в виде дифференциальных уравнений. При этом нужно подчеркнуть, что данные методы особенно эффективны в сочетании с использованием современных средств вычислительной техники и САПР. Наиболее распространенным и универсальным численным методом решения дифференциальных уравнений является метод конечных разностей. 15 Решение систем ОДУ: метод конечных разностей Сущность метода состоит в следующем. Область непрерывного изменения аргумента (например, отрезок) заменяется дискретным множеством точек xi (i = 0, 1, … n), называемых узлами. Эти узлы составляют разностную сетку. Расстояние между двумя соседними узлами hi = xi +1 – xi, (i = 0, 1, … n –1) называют шагом сетки. (на практике часто используют равномерную сетку, для которой hi = h = const, т.е. расстояние между узлами одинаковое). Искомая функция непрерывного аргумента приближенно заменяется функцией дискретного аргумента на заданной сетке. Эта функция называется сеточной. Исходное дифференциальное уравнение заменяется разностным уравнением относительно сеточной функции. При этом для входящих в уравнение производных используются соответствующие конечно-разностные соотношения. Такая замена дифференциального уравнения разностным называется его аппроксимацией на сетке (или разностной аппроксимацией). Таким образом, решение дифференциального уравнения сводится к отысканию значений сеточной функции в узлах сетки. Значения сеточной функции зависит от шага сетки h, т.е. на самом деле мы имеем не одну разностную задачу, а семейство задач, зависящих от параметра h. Это семейство задач называют разностной схемой решения. 16 Решение систем ОДУ: метод конечных разностей Решение разностной задачи, в результате которого находятся значения сеточной функции yi в узлах xi, приближенно заменяет решение исходной дифференциальной задачи. Однако не всякая разностная схема дает удовлетворительное решение, т. е. получаемые значения сеточной функции yi не всегда с достаточной точностью аппроксимируют значения искомой функции y(xi) в узлах сетки. Здесь важную роль играют такие понятия, как устойчивость и сходимость разностной схемы. Устойчивость метода численного интегрирования характеризует поведение ошибки интегрирования во времени. Нарастание ошибки свидетельствует о неустойчивости метода. Устойчивость метода зависит от шага интегрирования и значения корней характеристического уравнения диф. уравнения, таким образом выбор величины шага интегрирования является результатом компромисса между требуемой точностью решения и устойчивостью метода. 17 Задача Коши: общие положения Требуется найти функцию y = y(x), удовлетворяющую уравнению y' = f (x, y(x)) и принимающую при х = х0 заданное значение y0: y(x0) = y0 При этом будем для определенности считать, что решение нужно получить для значений х > x0, т.е. начальные условия определены в начале отрезка интегрирования h x1 x2 x3 x4 … xi xi+1 … xn x Введем разностную сетку x0, x1, ... и шаги hi = xi +1 – xi, (i = 0, 1, ... , n –1). В каждой точке xi, вместо значений функции y(xi) вводятся числа yi, аппроксимирующие точное решение y(x) на данном множестве точек. Функцию у, заданную в виде таблицы {xi, yi} (i = 0, 1, ... , n), называют сеточной функцией. Если в правой части разностного уравнения отсутствует величина yi + 1, т. е. значение yi + 1 явно вычисляется по k предыдущим значениям yi, yi – 1, … , yi – k + 1 то такая разностная схема называется явной. Если в правую часть уравнения входит искомое значение yi + 1, то решение этого уравнения усложняется. В таких методах, называемых неявными, его приходится разрешать относительно yi + 1 с помощью итерационных методов 18 Задача Коши: одношаговые методы Метод Эйлера (метод ломанных). Рассмотрим диф.уравнение y'= f (x, y(x)). В окрестностях узлов х = xi (i = 0, 1, ... , n) и заменим в левой части производную y' аппроксимацией в виде сеточной функций yi следующего вида y −y i +1 hi i = f ( xi , yi ) Будем считать для простоты узлы равноотстоящими, т.е. hi = xi + 1 – xi = h = const. Тогда получим выражение вида yi + 1 = yi + h⋅f (xi, yi), (i = 0, 1, ... , n). Заметим, что из исходного диф. уравнения y'= f (x, y(x)) следует, что y' (xi, yi) = f (xi, y(xi)) = f (xi, yi), т.е. это выражение представляет собой приближенное нахождение значения функции y(x) в точке xi + 1 при помощи разложения в ряд Тейлора с отбрасыванием членов второго и более высоких порядков. Другими словами, приращение функции полагается равным ее дифференциалу. 19 Задача Коши: одношаговые методы Метод Эйлера (метод ломанных). Полагая i = 0, находим значение сеточной функции у1 при х = х1: y1 = y0 + h⋅f (x0, y0). Требуемое здесь значение функции у0 уже задано начальным условием, т. е. y0 = y(x0). Аналогично могут быть найдены значения сеточной функции в других узлах: y2 = y1 + h⋅f (x1, y1), . . . . . . yn = yn – 1 + h⋅f (xn – 1, yn – 1). Построенный нами алгоритм называется прямым методом Эйлера (или методом ломанных). Она имеет вид рекуррентной формулы, с помощью которой значение сеточной функции yi +1 в узле xi +1 вычисляется по ее значению уi в предыдущем узле xi : yn+1 = yn + h⋅y'n. В связи с этим прямой метод Эйлера относится к одношаговым методам. 20 Задача Коши: одношаговые методы Метод Эйлера (метод ломанных). Геометрическая интерпретация метода Эйлера дана на рис. Ниже. Изображены первые два шага, т. е. проиллюстрировано вычисление сеточной функции в точках x1 и x2. точное решение диф. уравнения для начального условия y0 = y(x0) ε1 = Δy 1 h Интегральные кривые 0, 1, 2 описывают точные решения диф. уравнения. При этом кривая 0 соответствует точному решению задачи Коши для начального условия y(x0) = y0 так как она проходит через начальную точку А(х0, у0). Точки В, С получены в результате численного решения задачи Коши методом Эйлера. Их отклонения от кривой 0 характеризуют погрешность метода. Отрезок АВ – отрезок касательной к кривой 0 в точке А и ее наклон характеризуется значением производной y'(x0) = f (x0, у0). Погрешность появляется потому, что приращение значения функции при переходе от точки х0 к х1 заменяется приращением ординаты касательной к кривой 0 в точке А. 21 Задача Коши: одношаговые методы Метод Эйлера (метод ломанных). Блок-схема алгоритма решения задачи Коши методом Эйлера приведена на рисунке Пример 1. Рассмотрим решение задачи Коши следующего вида: y' = (x + y), при y(0) = 1 на отрезке [0, 0.4] c шагом h = 0.1 Аналитическое решение имеет вид: ϕ(x) = 2ex – x – 1. Решение: 1 шаг: i = 1, x0 = 0, y0 = 1, f = (x0 + y0) = (0 + 1) = 1 y0 = y0 + h⋅ f = 1 + 0.1⋅1 = 1.1; (точное решение ϕ(0.1) = 1.110342); x0 = x0 + h = 0 + 0.1 = 0.1 2 шаг: i = 2, x0 = 0.1, y0 = 1.1, f = (x0 + y0) = (0.1+ 1.1) = 1.2 y0 = y0 + h⋅ f = 1.1 + 0.1⋅1.2 = 1.22; (ϕ(0.2) = 1.242805) x0 = x0 + h = 0.1 + 0.1 = 0.2 3 шаг: i = 3, x0 = 0.2, y0 = 1.22, f = (x0 + y0) = (0.2+ 1.22) = 1.42 y0 = y0 + h⋅ f = 1.22 + 0.1⋅1.42 = 1.362; (ϕ(0.3) = 1.399718) x0 = x0 + h = 0.2 + 0.1 = 0.3 4 шаг: i = 4, x0 = 0.3, y0 = 1.362, f = (x0 + y0) = (0.3+ 1.362) = 1.662 y0 = y0 + h⋅ f = 1.362 + 0.1⋅1.662 = 1.5282; (ϕ(0.4) = 1.58365) x0 = x0 + h = 0.3 + 0.1 = 0.4 начало Ввод x0, y0, h, N Вывод x0, y0 i =1…N F =F(x0,y0) конец y0 = y0 + h*F x0 = x0 + h Вывод x0, y0 22 Задача Коши: одношаговые методы Обратный метод Эйлера. Можно выразить значение функции y(x) в следующей точке xk+1 через ее значение в предыдущей точке хk и значение производной в искомой точке x‘k+1. Если, полагать, что наклон (т.е. производная) функции на интервале h остается неизменным. Тогда можно записать yn+1 = хn + h⋅y'n+1. Это выражение описывает обратную формулу Эйлера. Обратим внимание на то, что здесь значение функции yn+1 в точке xn+1 входит как в правую, так и в левую части уравнения, так как y'n+1 = f (yn+1, xn+1). Трансцендентный характер формулы обратного метода Эйлера требует для вычислений итерационные методы расчета, причем значение функции в предыдущей точке должно быть известно, а приближенное значение функции в искомой точке предсказано, например, по прямой формуле Эйлера. 23 Задача Коши: одношаговые методы Методы Эйлера. y(x) Приближенное решение диф. уравнения yn+1 Δy* yn*+1 Δy** yn yn**+1 h xn Точное решение диф.уравнения xn+1 xn+1 x Иллюстрация прямого (красный цвет линии) и обратного (зеленый цвет) методов Эйлера Наклон ломанной линии (численного решения) соответствует производной (т.е. наклону касательной) в точке xn для прямого метода, и xn+1 – для обратного. Ошибка метода зависит от величины шага h и «кривизны» функции y(x) 24 Задача Коши: одношаговые методы Погрешность метода дифференцирования Погрешность аппроксимации производной конечными разностями состоит из двух частей: εi = εi` + εi``. Составляющая εi` – определяется погрешностью начального значения ε0 = y0 – Y(x0). Как правило, начальное значение задается точно, т.е. y0 = Y(x0). и тогда ε0 = 0 и следовательно, равна нулю та часть погрешности решения εi`, которая связана с ε0. Погрешность εi``, обусловлена отброшенными членами ряда в разложении функции по формуле Тейлора. На каждом шаге эта погрешность имеет порядок O(h2), так как именно члены такого порядка отброшены. При нахождении решения в конечной точке хn, отстоящей на конечном расстоянии T от первой точки x0, погрешность суммируется. Суммарная погрешность, очевидно, будет равна nO(h2). Если учесть, что шаг равен h = T/n, то для суммарной погрешности получаем окончательное выражение: nO(h2) = T/h O(h2) = O(h). Таким образом, рассмотренные методы Эйлера имеют первый порядок точности относительно шага h. 25 Задача Коши: одношаговые методы Помимо методов Эйлера существуют и другие явные одношаговые методы. Так, рассмотренные методы Эйлера являются частными случаями методов первого порядка, относящихся к классу методов Рунге-Кутта. Эти методы используют для вычисления значения yi+1 (i = 0, 1, ... , n) значение уi, а также значения функции f (x, у) при некоторых других, специальным образом выбираемых значениях x∈[xi, xi+1] и y(x). На их основе могут быть построены разностные схемы разного порядка точности. Широкое распространение для решения задач научного и инженерно-технического характера нашел метод Рунге-Кутта четвертого порядка (РК-4) точности. Метод РК-4 эффективен, надежен и легко реализуется программными средствами. 26 Задача Коши: одношаговые методы Метод трапеций Прямую и обратную формулы Эйлера также можно записать в виде (yn+1 – yn)/h = y'n (yn+1 – yn)/h = y'n+1 Такая интерпретация прямой и обратной формул Эйлера подводит к выводу о возможности вывода других конечно-разностных формул интегрирования. Так, определив приращение искомой функции на каждом шаге интегрирования в виде линейной комбинацией производных y'n и y'n+1 : (yn+1 – yn) / h = b0 y'n + b1 y'n+1 и принимая b0 = b1 =1/2, получаем известную формулу трапеций (также называемую, модифицированной формулой Эйлера) yn+1 = yn + h ⋅ ( yn′ + yn′ +1 ) y + h ⋅ [ f ( xn , yn ) + f ( xn+1 , yn+1 )] = n 2 2 Формула трапеций является простейшей из многошаговых формул, когда для определения значения функции в текущей точке используются значения функции и ее производных в других точках (в данном случае, в двух). Можно показать, что метод трапеций имеет второй порядок точности относительно 27 величины шага h (т.е. ошибка пропорциональна O(h2)). Задача Коши: методы прогноза - коррекции В обратной формуле Эйлера и формуле трапеций для вычисления значений функции в текущей точке нужно знать значения производных в следующей точке. Для выхода из этой ситуации используются итерационные методы, которые требуют задания начального приближения для значений функции. В качестве начального значения можно либо использовать значение функции в предыдущей точке, либо определить ее по прямой формуле Эйлера. Таким образом, мы приходим к понятию прогноза, когда предсказывается значение искомой функции в следующей точке. На основе заданного предсказания значения функции по соотношению y'n+1 = f(yn+1, xn+1) теперь можно определить значение производной в этой точке и скорректировать («уточнить») – это значение, используя обратную формулу Эйлера или формулу трапеций. Отсюда следует понятие коррекции. Описанная методика численного интегрирования, когда для предсказания используются прямые методы, а для коррекции – обратные, получила название метода прогноза-коррекции. Итерации уточнения позволяют фиксировать ошибку и адаптировать шаг численного интегрирования по мере необходимости для получения решения с требуемой точностью. 28 Метод прогноза-коррекции Метод трапеций Блок-схема алгоритма решения задачи Коши методом трапеций (прогноз-коррекция) приведена на рисунке Пример 1. Рассмотрим решение задачи Коши следующего вида: y' = (x + y), при y(0) = 1 на отрезке [0, 0.4] c шагом h = 0.1 Аналитическое решение имеет вид: ϕ(x) = 2ex – x – 1. Решение: 1 шаг: i = 1, x0 = 0, y0 = 1, f0 = (x0 + y0) = (0 + 1) = 1 y1 = y0 + h ⋅ f0 = 1 + 0.1 ⋅1 = 1.1; (прогноз) x1 = x0 + h = 0 + 0.1 = 0.1; f1 = (x1 + y1) = (0.1 + 1.1) = 1.2 y0 = y0 + 0.5 ⋅ h ⋅ (f0 + f1) = 1 + 0.5 ⋅ 0.1 ⋅(1 + 1.2) = 1.11 (коррекция) (точное решение ϕ(0.1) = 1.110342); x0 = x1 = 0.1 2 шаг: i = 2, x0 = 0.1, y0 = 1.11, f0 = (x0 + y0) = (0.1 + 1.11) = 1.21 y1 = y0 + h ⋅ f0 = 1.11 + 0.1 ⋅1.21 = 1.231; (прогноз) x1 = x0 + h = 0.1 + 0.1 = 0.2; f1 = (x1 + y1) = (0.2 + 1.231) = 1.431 y0 = y0 + 0.5 ⋅ h ⋅ (f0 + f1) = 1.11 + 0.5 ⋅ 0.1 ⋅(1.21 + 1.431) = 1.24205 ϕ(0.2) = 1.242805) x0 = x1 = 0.2 начало Ввод x0, y0, h, N Вывод x0, y0 i =1…N F0 =F(x0,y0) конец прогноз y1 = y0 + h*F0 x1 = x0 + h F1 =F(x1,y1) коррекция y0 = y0 + 0.5*h*(F0 +F1) x0 = x1 Вывод x0, y0 29 Задача Коши: методы прогноза - коррекции По сравнению с одношаговыми методами методы прогноза-коррекции имеют рад особенностей: 1. Для реализации методов прогноза-коррекции необходимо иметь информацию о нескольких предыдущих точках: другими словами, они не относятся к числу «самостартующих» методов. Для получения исходной информации приходится прибегать к какому-либо одношаговому методу. Если в процессе решения дифференциальных уравнений методом прогнозакоррекции изменяется шаг, то обычно приходится временно переходить на одношаговый метод. 2. Поскольку для методов прогноза-коррекции требуются данные о предыдущих точках, то соответственно предъявляются и повышенные требования к объему памяти компьютера. 3. Одношаговые методы и методы прогноза-коррекции обеспечивают примерно одинаковую точность результатов. Однако вторые в отличие от первых позволяют легко оценить погрешность на каждом шаге. По этой причине, пользуясь одношаговыми методами, величину шага h выбирают несколько меньше, чем это, строго говоря, необходимо, и поэтому методы прогноза-коррекции требуют почти вдвое меньше машинного времени, чем методы Рунге-Кутта сравнимой точности. 4. Применяя метод Рунге-Кутта 4-го порядка точности, на каждом шаге приходится вычислять четыре значения функции, в то время как для обеспечения сходимости метода прогнозакоррекции того же порядка точности достаточно двух значений функции. Поэтому методы прогноза-коррекции требуют вдвое меньше затрат машинного времени, чем методы РунгеКутта сравнимой точности. 30 Расчет переходных процессов в РЭУ Полученные формулы численного интегрирования пригодны как для линейных, так и для нелинейных дифференциальных уравнений. Получим выражения этих формул интегрирования для систем линейных дифференциальных уравнений, которые описывают ММ эл. цепи. Пусть имеем систему линейных дифференциальных уравнений в нормальной форме Коши и векторно-матричном представлении: X' (t) = A⋅X(t) + B⋅E(t), где X(t) – вектор искомых функций времени t; A – матрица коэффициентов системы уравнений; B – вектор коэффициентов источников; E(t) – вектор воздействий (источников) в виде известных функций от времени. Применив прямую формулу Эйлера к этой системе уравнений, в конечно-разностном представлении получим: Xn+1 = Xn + h⋅X'n = Xn + h ⋅ (A⋅Xn + B⋅Еn). Откуда, приводя подобные, можно записать матричную форму прямой формулы Эйлера: Xn+1 = (1 + h⋅A)⋅Xn + h⋅B⋅Еn. Формула имеет итерационный характер, т.е. по известному на предыдущем шаге вектору Xn находим его текущее значение Xn+1, затем, подставляя в правую часть, находим следующее значение и т.д. 31 Расчет переходных процессов Аналогичным образом можно построить процесс решения на основе обратного метода Эйлера: Xn+1 = Xn + h⋅X'n+1 = Xn + h⋅(A⋅Xn+1 + В⋅Еn+1). После преобразований, получим: (1 – h⋅A)⋅Xn+1 = Xn + h⋅В⋅Еn+1. (это СЛАУ вида Т⋅x = w) Как видно, в случае линейных систем дифференциальных уравнений (т.е. если матрица А не зависит от предыдущего состояния ММ) трансцендентность обратной формулы Эйлера исчезает. Решение системы также имеет итерационный характер. Но, в данном случае, на каждом шаге нам придется решать СЛАУ: Xn+1 = (1 – h⋅A)-1⋅(Xn + h⋅В⋅Еn+1). Видно, что обратную матрицу (1 – h⋅A)-1 здесь можно вычислить только один раз, а затем менять только вектор правой части СЛАУ, который отвечает за воздействия на цепь. Для нахождения решения СЛАУ лучше пользоваться методами, основанными на LU- или QR-факторизации (разложении) исходной матрицы коэффициентов (1 – h⋅A). 32 Расчет переходных процессов В случае формулы трапеций матричное уравнение примет вид: Xn+1 = Xn + 0.5⋅h⋅(A⋅Xn + B⋅Еn + A⋅Xn+1 + B⋅Еn). Перегруппировывая компоненты этого выражения, получим новое уравнение: (1 – 0.5⋅h⋅A)⋅Xn+1 = (1 + 0.5⋅h⋅A)⋅Xn + 0.5⋅h⋅B ⋅(Еn + Еn+1). Это уравнение решается аналогично, т.е. путем решения СЛАУ при заданных начальных условиях Xn, подстановки найденного решения Xn+1 в правую часть на место Xn, нахождения нового решения и т.д. Аналогичным образом можно построить более сложные методы численного расчета дифференциальных уравнений, которые описывают ММ цепи, основанных на объединении одно- и многошаговых явных и неявных формул численного дифференцирования. 33 Расчет переходных процессов. Пример. Рассмотрим решение задачи расчета переходных процессов в RLC-цепях методом переменных состояний с помощью рассмотренных методов численного решения ОДУ. Для примера возьмем цепь второго порядка Система уравнений для расчета переменных состояния iL(t) и uC(t) имеет следующий вид 1 1 R1 ⎧ diL (t ) ⎪⎪ dt = − L iL (t ) − L uC (t ) + L E , ⎨ ⎪ duC (t ) = 1 iL (t ) − 1 uC (t ) + 0. ⎪⎩ dt C CR2 Эту систему необходимо дополнить двумя уравнениями для расчета токов через резисторы R1 и R2 ⎧i1 (t ) = 1⋅ iL (t ) + 0 ⋅ uC (t ), ⎪ 1 ⎨ ⎪i2 (t ) = 0 ⋅ iL (t ) + R ⋅ uC (t ). 2 ⎩ 34 Расчет переходных процессов. Пример. Таким образом, результирующая система уравнений примет вид где ⎡ R1 ⎢− L A=⎢ 1 ⎢ ⎢⎣ C 1 ⎤ ⎡1⎤ L ⎥ ; B = ⎢ ⎥; L 1 ⎥ ⎢0⎥ ⎥ − ⎣ ⎦ CR2 ⎥⎦ − ⎡1 C = ⎢0 ⎢ ⎣ 0⎤ 1 ⎥; D = ⎡ ⎤; ⎢0 ⎥ ⎣ ⎦ R2 ⎥⎦ x(t) = ⎡ iL (t ) ⎤ ⎢u (t )⎥ ; ⎣ C ⎦ ⎧ dx(t ) = A ⋅ x(t ) + B ⋅ e(t ); ⎪ ⎨ dt ⎪⎩ y (t ) = C ⋅ x(t ) + D ⋅ e(t ), ⎡ E (t )⎤ e(t) = ⎢ 0 .⎥ ⎣ ⎦ Первое уравнение – это система из двух ОДУ, а второе – СЛАУ из двух уравнений. Запишем итерационный процесс для этой системы на основе обратного метода Эйлера. Для первого уравнения ⎡ diL ⎤ ⎡− R1 ⎢ dt ⎥ ⎢ L ⎢ du ⎥ = ⎢ 1 ⎢ C⎥ ⎢ ⎣ dt ⎦ ⎣ C 1 ⎤ ⎡E⎤ L ⎥ ⋅ ⎡ iL (t ) ⎤ + ⎢ ⎥ 1 ⎥ ⎢⎣uC (t )⎥⎦ ⎢ L ⎥ − ⎥ ⎣0⎦ CR2 ⎦ − в соответствии с обратной формулой Эйлера, получим где ⎡ R1 ⎢− A = ⎢ 1L ⎢ ⎢⎣ C 1 ⎤ ⎡E⎤ L ⎥; B = ⎢ ⎥; L 1 ⎥ ⎢0⎥ ⎥ − ⎣ ⎦ CR2 ⎥⎦ − X= (1 – h⋅A)⋅Xn+1 = Xn + h⋅B⋅Еn+1 ⎡ iL (t ) ⎤ ⎢u (t )⎥ ⎣ C ⎦ 35 Расчет переходных процессов. Пример. Соответственно, учитывая, что (1 – h⋅A) = R1 h ⎤ ⎡ ⎢1 + h L L ⎥ ⎢ h h ⎥ ⎢ − 1+ ⎥ CR2 ⎦ ⎣ C получим следующую систему линейных уравнений (здесь верхний индекс определяет номер итерации): R1 h ⎤ ⎡ 1 + h ( k +1) (k ) ⎡1⎤ ⎢ L L ⎥ ⋅ ⎡ iL ⎤ = ⎡ iL ⎤ + h ⋅ ⎢ ⎥ ⋅ E ( k +1) ⎢ L h h ⎥ ⎢u ( k +1) ⎥ ⎢u ( k ) ⎥ ⎢0⎥ 1+ ⎢ − ⎥ ⎣ C ⎦ ⎣ C ⎦ ⎣ ⎦ CR2 ⎦ ⎣ C Для вычисления токов через резисторы R1 и R2, на каждом шаге нужно также рассчитывать второе уравнение по найденным переменным состояния, т.е. ⎡i1( k ) ⎤ ⎡1 ⎢ ( k ) ⎥ = ⎢0 ⎣i2 ⎦ ⎢⎣ 0 ⎤ ⎡ (k ) ⎤ 1 ⎥ ⋅ ⎢ iL ⎥ (k ) R2 ⎥⎦ ⎣uC ⎦ 36 Расчет переходных процессов. Пример. Приведем текст программы на MATLAB для данного примера function my_trans_1 clear all close all % -- задаем значения элементов цепи R1 = 1; R2 = 5; C = 1; L = 1; E = 2; % Формируем матрицы и вектор системы ОДУ A = [-R1/L -1/L; 1/C -1/(C*R2)]; B = [1/L; 0]; C = [1 0; 0 1/R2]; D = [0; 0]; % Задаем интервал интегрирования от 0 до 30 сек tspan = [0 30]; % ----- Метод Эйлера t1 = tspan(1); t2 = tspan(2); h = abs(t2-t1)/39; % - начальный шаг % вычисление обр. матрицы коэф. СЛАУ AA = inv(eye(size(A)) - h*A); E_T(1) = my_e(t1,E); % отображение графиков subplot(311); plot(T, X(1,:),'-r'); title('i_L(t) in A'); subplot(312); plot(T,X(2,:),'-b',T,E_T,'-m'); title('E(t) and u_C(y) in V'); subplot(313); plot(T,Y(1,:),'-b',T,Y(2,:),'-m'); title('i_R_1(t) and i_R_2(t) in A'); xlabel('Time in sec'); % --- зависимость источника сигнала от времени E(t) function y = my_e(t,E) y = E % y = -1*E*sign(mod(t,30)-15); % начальные условия в точке t0 X = [iL=0; uC=0] X(:,1) = [0;0]; % начальные токи через R1 и R2 Y = [iR1=0;iR2=0] Y(:,1) = C*X(:,1) + D*E_T(1); T = t1:h:t2; % -- вектор точек времени % цикл по времени for k=2:length(T) E_T(k) = my_e(T(k),E); X(:,k) = AA*(X(:,k-1) + h*B*E_T(k)); Y(:,k) = C*X(:,k) + D*E_T(k); end 37 Пример Результат работы программы для N=40 iL(t) 1 0.8 0.6 0.4 0.2 5 10 15 20 25 30 20 25 30 20 25 30 E(t) и uC(t) 2 1.5 1 0.5 5 10 15 iR1(t) и iR2(t)) 1 0.8 0.6 0.4 0.2 5 10 15 Time in sec 38 Прямые методы формирования ММ Аналогично можно получить вычислительный процесс расчета переходных процессов в РЭУ на основе ММ полученных прямыми методами. Основное требование, чтобы ММ была представлена в виде X' (t) = A⋅X(t) + B⋅E(t) Для этого в операторной форме система уравнений ММ (матрица коэффициентов) должна выглядеть следующим образом: T = G + p⋅C, G – матрица действительной части; где C – матрица мнимой части; р – оператор Лапласа. Для соблюдения этого условия необходимо емкостные ветви представлять в виде комплексной проводимости (адмитанса), а индуктивные ветви – в виде комплексного сопротивления (импеданса). В этом случае оператор Лапласа окажется во всех компонентах реактивной части матрицы коэффициентов в числителе. 39 Прямые методы формирования ММ Соответственно алгебраическая система уравнений может быть записана как (G + p⋅C)⋅Х = W Применяя к ней обратное преобразование Лапласа (переход во временную область), получим систему ОДУ следующего вида: G⋅Х + B⋅Х' = W. или C⋅Х' = W – G⋅Х . (это СЛАУ) Так как матрица мнимой части С в общем случае может быть вырожденной (т.е. det(C)=0), нам необходимо найти такой переход к системе дифференциальных уравнений в нормальной форме Коши, который не требует явного обращения матрицы С. 40 Прямые методы формирования ММ Рассмотрим обратную формулу Эйлера, предварительно умноженную на матрицу С: C⋅Xn+1 = C⋅Xn + h⋅C⋅Х'n+1. и, подставляя вместо C⋅Х'n+1 значения из дискретизованного для этого случая уравнения, получим C⋅Xn+1 = C⋅Xn + h⋅(Wn+1 – G⋅Хn+1) или, сгруппировав компоненты, окончательно можно записать (C + h⋅G)⋅Xn+1 = C⋅Xn + h⋅Wn+1. Это уравнение представляет собой обратную формулу Эйлера для интегрирования систем дифференциальных уравнений с матрицей коэффициентов, заданной в комплексной форме. Это уравнение также показывает, что нет необходимости описывать цепь дифференциальными уравнениями в нормальной форме Коши, если известны действительная G и мнимая C части матрицы коэффициентов Т. При этом матрицы G и С могут быть вырожденными, в то время как матрица (C + h⋅G) конечной системы уравнений будет не вырожденная. 41 Прямые методы формирования ММ Таким образом, методы интегрирования или вычисления переходных процессов в линейных цепях сводятся к последовательному решению систем линейных алгебраических уравнений. В случае нелинейных цепей матрицы коэффициентов G, C и вектор свободных членов W являются функциями от вектора переменных Х, т.е. система алгебраических уравнений C⋅Х' = W – G⋅Х становится нелинейной. В этом случае применение преобразования Лапласа для перехода к системе дифференциальных уравнений становится проблематичной. Однако с определенной степенью условности его формальное применение приведет к системе нелинейных дифференциальных уравнений. Для решения нелинейных систем дифференциальных уравнений также можно применять те же самые численные методы интегрирования. Однако, если в случае линейных систем достаточно один раз сформировать исходную алгебраическую систему уравнений и на каждом шаге переформировывать только вектор свободных членов в формуле интегрирования с учетом предыдущего решения, то в случае нелинейных цепей необходимо на каждом шаге интегрирования переформировывать как исходную алгебраическую систему, так и систему дифференциальных уравнений с учетом решения на предыдущем шаге. 42 Прямые методы формирования ММ Рассмотрим расчет переходных процессов в данной цепи модифицированным методом узловых потенциалов Сформируем ММ в операторной форме ⎡0 ⎤ ⎢ ⎥ Yn1 = 0 1 / R1 ⎢ ⎥ ⎢⎣0 0 1 / R2 + pC ⎥⎦ ⎡− 1 1 ⎤ A2 = ⎢ 1 0⎥ ⎢ ⎥ ⎢⎣ 0 − 1⎥⎦ ⎡ Yn1 ⎢Y ⋅ A t ⎣ 2 2 ⎡1 0⎤ Y2 = ⎢ ⎥ ⎣0 1 ⎦ A 2 ⎤ ⎡U n ⎤ ⎡ J n ⎤ ⋅ = Z 2 ⎥⎦ ⎢⎣ I 2 ⎥⎦ ⎢⎣ W2 ⎥⎦ 0 ⎤ ⎡0 Z2 = ⎢ ⎥ ⎣0 − pL ⎦ ⎡0 ⎤ J n = ⎢0 ⎥ ⎢ ⎥ ⎢⎣0⎥⎦ ⎡V ⎤ W2 = ⎢ 1 ⎥ ⎣0⎦ ⎡1 0⎤ ⎡− 1 1 0 ⎤ ⎡− 1 1 0 ⎤ Y2 ⋅ A t2 = ⎢ ⎥ ⋅ ⎢ 1 0 − 1⎥ = ⎢ 1 0 − 1⎥ 1 ⎦ ⎣ ⎦ ⎦ ⎣ ⎣ Получим ММ в виде: −1 1 ⎤ ⎡ U1 ⎤ ⎡ 0 ⎤ ⎡0 ⎢ 0 1/ R 1 0 ⎥ ⎢U 2 ⎥ ⎢ 0 ⎥ 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢0 0 1 / R2 + pC 0 − 1 ⎥ ⋅ ⎢U 3 ⎥ = ⎢ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − 1 1 ⎢ ⎥ ⎢ IV 1 ⎥ ⎢V1 ⎥ ⎢⎣ 1 −1 0 − pL ⎥⎦ ⎢⎣ I L ⎥⎦ ⎢⎣ 0 ⎥⎦ 43 Прямые методы формирования ММ Представим систему в виде (G +pС)*X = W ⎛⎡ 0 ⎜⎢ ⎜ ⎢ 0 1 / R1 ⎜⎢ 0 0 1 / R2 ⎜⎢ ⎜ ⎢− 1 1 ⎜⎢ −1 ⎝⎣ 1 −1 1 1⎤ ⎡0 0 0 0 0 ⎤ ⎞ ⎡ U1 ⎤ ⎡ 0 ⎤ ⎢0 0 0 0 0 ⎥ ⎟ ⎢U ⎥ ⎢ 0 ⎥ 0⎥ ⎥ ⎢ ⎥⎟ ⎢ 2⎥ ⎢ ⎥ − 1⎥ + p ⋅ ⎢0 0 C 0 0 ⎥ ⎟ ⋅ ⎢U 3 ⎥ = ⎢ 0 ⎥ ⎥ ⎢ ⎥⎟ ⎢ ⎥ ⎢ ⎥ 0⎥ ⎢0 0 0 0 0 ⎥ ⎟ ⎢ IV 1 ⎥ ⎢V1 ⎥ ⎢⎣0 0 0 0 − L ⎥⎦ ⎟⎠ ⎢⎣ I L ⎥⎦ ⎢⎣ 0 ⎥⎦ 0 ⎥⎦ Соответственно, обратная формула Эйлера имеет вид (C + h⋅G)⋅Xn+1 = C⋅Xn + h⋅Wn+1. Искомые переменные можно найти как VC= U3; IR1= -U2/R1 ; IR2= U3/R2; 44 Расчет переходных процессов. Пример. Приведем текст программы на MATLAB для данного примера function my_transient_2 clear all close all % -- задаем значения элементов цепи R1 = 1; R2 = 5; C = 1; L = 1; E = 2; % Формируем матрицы и вектор системы ОДУ G = [0 0 0 -1 1; 0 1/R1 0 1 0; 0 0 1/R2 0 -1; -1 1 0 0 0; 1 0 -1 0 0]; C = [0 0 0 0 0; 0 0 0 0 0; 0 0 C 0 0; 0 0 0 0 0; 0 0 0 0 -L]; % Задаем интервал интегрирования от 0 до 12 сек tspan = [0 12]; X(:,1) = [0 0 0 0 0]'; % начальные условия в точке t0 Y(:,1) = [0 0 0]'; % ----- T = t1:h:t2; % -- вектор точек времени for k=2:length(T) E_T(k) = my_e(T(k),E); W = [0 0 0 E_T(k) 0]'; X(:,k) = AA*(C*X(:,k-1) + h*W); Y(:,k) = [X(3,k); -X(5,k)/R1; X(3,k)/R2]; end % отображение графиков subplot(311); plot(T, X(5,:),'-r'); title('i_L(t)'); axis([0 t2 -1.5 0]) subplot(312); plot(T,Y(1,:),'-b',T,E_T,'-m'); title('E(t) и u_C(t) '); axis([0 t2 0 2]) subplot(313); plot(T,Y(2,:),'-r',T,Y(3,:),'-g'); title('i_R_1(t) и i_R_2(t))'); xlabel('Time in sec'); axis([0 t2 -0.5 1.5]) Метод Эйлера t1 = tspan(1); t2 = tspan(2); h = abs(t2-t1)/199; % - начальный шаг AA = inv(C + h*G); % вычисление обр. матрицы коэф. СЛАУ E_T(1) = my_e(t1,E); Y(:,1) = [X(3,1); -X(5,1)/R1; X(3,1)/R2]; % --- зависимость источника сигнала от времени E(t) function y = my_e(t,E) y = E % y = -1*E*sign(mod(t,30)-15); 45 Расчет переходных процессов. Пример. iL(t) -0.5 -1 -1.5 2 4 6 E(t) и uC(t) 8 10 12 2 ток IL получился со знаком «минус», так как индуктивность L «включена» неправильно в схеме в Micro-CAP 1.5 1 0.5 2 4 6 iR1(t) и iR2(t)) 8 10 12 2 4 6 Time in sec 8 10 12 1.5 1 0.5 -0.5 46 Расчет переходных процессов. Пример в Micro-CAP. ток IL получился со знаком «минус», так как индуктивность L «включена» неправильно в схеме в Micro-CAP 47 Прямые методы формирования ММ Прямое применение формул численного интегрирования для нелинейных систем дифференциальных уравнений не нашло широкого применения в реальных САПР в силу малой устойчивости и низкой точности получаемых решений. Более корректный подход к вычислению переходных процессов нелинейных цепей основан на том, что система нелинейных дифференциальных уравнений, представленная в нормальной форме Коши, на каждом шаге интегрирования представляется как система нелинейных алгебраических уравнений, и уже к ней применяются известные итерационные методы решения типа метода НьютонаРафсона и др. Формулы интегрирования в этом случае представляют собой результат применения итерационных алгоритмов к традиционным численным методам. Такой подход обеспечивает необходимую устойчивость и точность решения. 48 Заключение. По итогам лекции 8 студенты должны знать: 1. Методы построения ММ для линейных цепей для расчета переходных процессов методом переменных состояния 2. Построение ММ для линейных цепей на основе прямых методов для расчета переходных процессов 3. Базовые алгоритмы для моделирования РЭУ во временной области 4. Уметь оценивать сложность вычислительных задач при анализе цепей во временной области Должны уметь решать задачи на: 1. Формирование ММ линейных цепей для расчета переходных процессов 49 Спасибо за внимание! Ваши вопросы … 50 Устойчивость методов решения диф. уравнений Понятие устойчивости метода Одним из важных устойчивость. свойств методов интегрирования является его Устойчивость метода численного интегрирования характеризует поведение ошибки интегрирования при изменении времени расчета. Нарастание ошибки свидетельствует о неустойчивости метода. Устойчивость метода, как будет показано ниже, зависит от шага интегрирования и значения корней характеристического уравнения. Поскольку точность зависит от величины шага интегрирования, выбор его размера часто является результатом компромисса между точностью и устойчивостью. Изложим проблему устойчивости метода интегрирования на примере дифференциального уравнения вида x' = λ⋅x, (*) которое имеет аналитическое выражение решение в виде x = х0⋅exp(λ⋅t), где λ – корень характеристического уравнения (некоторое действительное или комплексное число). 51 Устойчивость методов решения диф. уравнений Понятие устойчивости метода Исследование устойчивости методов интегрирования начнем с прямой формулы Эйлера: хn+1 = хn + h⋅x'n. Используя формулу (*), получим выражение для производной x'n = λ⋅x n. Подставив ее в формулу Эйлера, получим хn+1 = (1 + λ⋅h)⋅хn. На первом шаге интегрирования: х1 = (1 + λ⋅h)⋅х0. На следующем шаге, соответственно, получим х2 = х1 + h⋅x'1 = (1 + λ⋅h)⋅х1 = (1 + λ⋅h)2⋅х0. Продолжив процедуру интегрирования в пределах n шагов, можно записать хn = (1 + λ⋅h)n⋅х0. 52 Устойчивость методов решения диф. уравнений Понятие устойчивости метода Предположим, что время интегрирования t → ∞ и, соответственно, n → ∞. Тогда, для того чтобы хn было ограниченным, для устойчивого дифференциального уравнения (когда Reλ < 0) необходимо выполнение условия | 1 + λ⋅h | < 1, где (**) h – размер шага (действительное число); λ – корень характеристического уравнения (возможно комплексный). Найдем область устойчивости, удовлетворяющей данному неравенству. Введем обозначение λ⋅h = q = u + j⋅v и подставив его в неравенство (**), получим | 1+ u + j⋅v | ≤ 1 или (1+ u)2 + v2 ≤ 1. 53 Устойчивость методов решения диф. уравнений Понятие устойчивости метода Это выражение описывает область устойчивости решения внутри круга с центром в точке (– 1, 0), изображенную на рис. Область устойчивости (ОУ) для метода Эйлера Im q ОУ –2 –1 Re q Полученный результат используется следующим образом: при Reλ<0 шаг h выбирается исходя из условия q = λ⋅h < 1, что соответствует точке внутри области устойчивости. В этом случае прямая формула Эйлера дает устойчивые решения, т.е. значения функции x(t) при увеличении n будут оставаться конечной величиной. При этом величина шага h определяет как устойчивость, так и достижимую точность интегрирования (погрешность решения). 54 Устойчивость методов решения диф. уравнений Понятие устойчивости метода Пример: применим прямую формулу Эйлера к простому дифференциальному уравнению вида x' = – x, при x0 = 1, соответствующему уравнению (*), при λ = – 1. Результаты численного интегрирования при разных h, сведены в таблицу n Значения x(t) h = 0,1 h=2 h=3 1 1 1 1 0,9 –1 –2 2 0,81 1 4 3 0,729 –1 –8 4 0,6561 1 16 5 0,59049 –1 –32 6 0,53144 1 64 q=λ⋅h -1 -2 -3 Внутри ОУ На границе ОУ Вне ОУ 55 Устойчивость методов решения диф. уравнений Решение ОДУ вида dx/dt = -x, 1 h=2 0.8 0.6 0.4 1 0.2 0.9 h = 0,1 0.8 -0.2 -0.4 -0.6 0.7 -0.8 0.6 -1 2 4 6 8 10 12 14 16 18 300 0.5 200 0.4 100 0.3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -100 красный – точное решение ОДУ синий – численное решение ОДУ -200 -300 h=3 -400 -500 -600 5 10 15 20 25 30 56 Устойчивость методов решения диф. уравнений Выводы по примеру Первая колонка решений, при h = 0,1, соответствует q = λ⋅h = – 1, т.е. области внутри единичной окружности и устойчивому решению ОДУ. Вторая колонка соответствует решению при h = 2 (q = – 2), т.е. на границе области устойчивости. Решение x(t) в этом случае колеблется между значениями 1 и – 1, но не нарастает. В третьей колонке решение при h = 3 (q = – 3) находится вне области устойчивости. В этом случае решение осциллирует и нарастает, хотя в соответствии с аналитическим решением, равным x(t) = exp(– t), должно уменьшаться. 57 Устойчивость методов решения диф. уравнений Продолжим исследование устойчивости методов интегрирования для обратной формулы Эйлера хn+1 = хn + h⋅x'n+1. Используя формулу (*), можно записать x'n+1 = λ⋅xn+1. Подставив это выражение в обратную формулу Эйлера, получим хn+1 = хn + h⋅λ⋅хn+1. На первом шаге интегрирования: х1 = х0 + h⋅λ⋅х1 или х1 = х0 /(1 – h⋅λ) = х0 /(1 – q). На следующем шаге, соответственно, получим х2 = х1 + h⋅x'2 = х1 /(1 – h⋅λ) = х0 /(1 – λ⋅h)2. Продолжив процедуру интегрирования в пределах n шагов, можно записать хn = х0 /(1 – λ⋅h)n = х0 /(1 – q)n. Для устойчивого решения, когда Reλ < 0, при t → ∞ и, соответственно, n → ∞, необходимо выполнение условия | 1/( 1 – q) | < 1. 58 Устойчивость методов решения диф. уравнений Найдем область устойчивости (ОУ), удовлетворяющую данному неравенству. Раскрывая обозначение q, получим | 1/(1 – u + j⋅v ) | ≤ 1 или (1 – u)2 + v2 ≥ 1. Этому выражению соответствует окружность с центром в точке (1,0). Неравенство выполняется вне окружности (см. рис.). Таким образом, обратная формула Эйлера устойчива для всех ОДУ, когда Reλ < 0 при любом шаге интегрирования h. Область устойчивости (ОУ) для обратного метода Эйлера Im q ОУ 1 2 Re q Размер шага h при этом можно выбирать лишь из соображений необходимой точности интегрирования (погрешности решения). 59 Устойчивость методов решения диф. уравнений Решение ОДУ вида dx/dt = -x, 1 0.8 0.6 1 h=2 0.4 0.2 0.9 h = 0,1 0.8 -0.2 -0.4 0.7 -0.6 0.6 -0.8 0.5 -1 0.4 300 0.3 2 4 6 8 10 12 14 16 18 200 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100 красный – точное решение ОДУ синий – численное решение ОДУ (прямой метод Эйлера) розовый - численное решение ОДУ (обратный метод Эйлера) -100 -200 -300 h=3 -400 -500 -600 5 10 15 20 25 30 60 Устойчивость методов решения диф. уравнений . В заключение, рассмотрим формулу трапеций: Можно показать (см. лит-ру),что границей устойчивости в этом случае является мнимая ось Im q, а областью устойчивости – вся левая полуплоскость (см. рис.). Таким образом, формула трапеций устойчива для всех схем, для которых выполняется условие Reλ < 0, при любом шаге интегрирования h. xn +1 = xn + h ⋅ ( xn′ + xn′ +1 ) 2 Область устойчивости (ОУ) для метода трапеций Im q ОУ Re q Замечание: формула трапеций дает устойчивый отклик для устойчивых цепей (для которых выполняется условие Reλ < 0) и нестабильный отклик для неустойчивых цепей (когда Reλ > 0). Это ценное свойство метода интегрирования, когда поведение решения для диф. уравнения заранее нам неизвестно. 61 Устойчивость методов решения диф. уравнений 1) Следует понимать, что выполнение условий устойчивости численного метода не подразумевает правильности выполнения расчетов. Это лишь означает, что любая ошибка не будет увеличиваться при последующих шагах. Одним из путей обеспечения точности при выбранном шаге h является использование формул с большим порядком интегрирования. Наиболее популярны на практике методы интегрирования 4-го порядка (методы Рунге-Кутты) 2) Следует также отметить, что полученные для простых методов интегрирования границы областей устойчивости справедливы лишь для рассматриваемого дифференциального уравнения вида x' = λ⋅x . Для других дифференциальных уравнений и других методов интегрирования границы областей устойчивости будут иными. 62 Терминология Формулы интегрирования, основанные на значениях функции и ее производных в предыдущие моменты времени, подобные прямой формуле Эйлера, называют иногда явными. Явные формулы в силу их слабой устойчивости используются главным образом для предсказания начальных значений при использовании других, неявных формул интегрирования. Неявными формулами интегрирования (аналогичными обратной формуле Эйлера и формуле трапеций) называются формулы, связывающие значение функции в следующей точке с ее производной в этой же точке, а также значениями функции и ее производными в предыдущих точках. Неявные формулы, как правило, более устойчивы и используются, в основном, для коррекции решения. В силу трансцендентного характера, когда неизвестное значение находится в правой и в левой частях уравнения, неявные формулы решаются итерационными методами, используемыми при решении нелинейных алгебраических уравнений. Объединение явных и неявных формул интегрирования приводит к методам прогноза-коррекции. 63 Терминология Формулу интегрирования называют А-устойчивой, если она дает ограниченное решение тестового дифференциального уравнения х' = λ⋅x для произвольных размеров шага h и любого числа шагов n при Reλ < 0. Обратная формула Эйлера и формула трапеций обладают этим свойством. Большинство радиоэлектронных цепей, встречающихся на практике, являются А-устойчивыми. Для линейной цепи это означает, что действительные части корней характеристических уравнений должны быть отрицательными (т.е. удовлетворяют условию Reλ < 0). Однако встречаются ситуации (например автоколебательные системы), когда схема в рабочей точке неустойчива и в установившемся режиме в ней наблюдаются периодические колебания. При численном интегрировании дифференциальных уравнений таких цепей следует применять методы, учитывающие эти свойства. Обратная формула Эйлера, например, в этой ситуации не подходит, т.к. дает затухающие решения. В этом случае лучше использовать формулу трапеций или другие методы интегрирования жестких систем диф. уравнений (например, многошаговые методы). 64