Качественный анализ нелинейных систем

Постановка задачи

Описание модели

модель описывается системой дифференциальных уравнений 2-ого порядка;

задается область определения значений переменных в виде системы неравенств:

Исследование стационарных точек

в заданной области значений переменных определяются стационарные точки системы путем решения системы алгебраических уравнений:

выполняется линеаризация системы уравнений в стационарных точках.

Точки классифицируются путем нахождения собственных значений матрицы (якобиана) линеаризованной системы. Возможны шесть типов точек.

Обозначения:

— собственные значения якобиана ().

Устойчивый узел:

Собственные значения действительные (мнимая часть равна нулю), меньшие нуля:

Неустойчивый узел:

Собственные значения действительные, большие нуля:

Устойчивый фокус:

Собственные значения комплексно-сопряженные, действительная часть меньше нуля:

Неустойчивый фокус:

Собственные значения комплексно-сопряженные, действительная часть больше нуля:

" Центр«:

Собственные значения комплексно-сопряженные, действительная часть равна нулю:

" Седло«:

Собственные значения действительные, с разными знаками:

Построение фазовых траекторий

Фазовые траектории (портреты) строятся в окрестностях стационарных точек с помощью метода численного интегрирования.

Практическая реализация задачи

Задача реализована в среде системы Тополог.

Описание модели

Описание системы уравнений выполняется с использованием стандартных типов системы Тополог: функция и вектор. В виде функции описывается только правая часть уравнения.

Пример: пусть мы имеем систему из двух нелинейных уравнений вида:

тогда описание данной системы на языке системы Тополог будет выглядеть следующим образом:

данные

x: вещ

y: вещ

F1: функция (x, y) = x²/100 + y²/25 — 1

F2: функция (x, y) = y — 10 * Sin (X)

VF: вектор = (F1, F2) |

J: якобиан (VF)

Исследование стационарных точек

Решение системы уравнений в заданной прямоугольной области выполняется стандартным операционным блоком системы Найти_Все_Решения. Он имеет следующий синтаксис:

Найти_Все_Решения (<�вектор функций>, <�погрешность>,

<�граница 1>, <�граница 2>,<�граница 3>, <�граница 4>,

<�шаг сетки>, <�массив решений>)

Решения ищутся путем запуска метода Ньютона из узлов сетки, наложенной на прямоугольную область поиска, и накопления найденных различных решений. Прямоугольная область поиска задается параметрами: <�граница 1>, <�граница 2>,<�граница 3>, <�граница 4>. Шаг сетки равномерный по обеим координатам и задается параметром <�шаг сетки>. После завершения работы блока все решения скапливаются в массиве точек, описанным как:

Решения: массив из Точка_вещ

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

Для получения собственных значений якобиана системы в найденных стационарных точках используется стандартный операционный блок системы Тополог Собственные_Значения. Он имеет следующий синтаксис:

Собственные_Значения (<�матрица>,

<�максимальное число итераций>, <�погрешность>) — >

(<�массив реальных частей>,

<�массив мнимых частей>)

На выходе данный блок формирует два массива, содержащих реальные и мнимые части собственных значений квадратной матрицы (в нашем случае якобиана).

Численные методы нахождения собственных значений матриц делятся на два типа: прямые и итерационные. Прямые методы находят собственные значения путем решения характеристического уравнения матрицы с помощью одного из численных методов (метод Ньютона (касательных), метод секущих, метод хорд и др.). Итерационные методы при помощи преобразований, не меняющих собственных значений матрицы, приводят ее к диагональному (в случае действительных собственных значений) или блочно-диагональному (при наличии комплексно-сопряженных собственных значений) виду. Значения из диагонали полученной матрицы и являются собственными значениями исходной матрицы.

В основу блока Собственные_Значения заложен итерационный метод нахождения собственных значений по методу Якоби с понижением нормы для действительных матриц.

Основная суть метода Якоби:

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

Для реализации такого преобразования матрицы можно выбирать в виде произведения двух матриц, где матрица определяет так называемое вращение, а — сдвиг или комплексное вращение. Элементы этих двух матриц удовлетворяют соотношениям:

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

.

Параметр сдвига удовлетворяет соотношению

где

; ;

;

.

По окончанию итерационного процесса матрица будет нормальной, причем, и либо, либо для всех значений (,), а матрица будет блочно-диагональной: блоки размера () будут содержать действительные собственные значения, а () — комплексно-сопряженные.

Построение фазовых траекторий

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

Алгоритм метода может быть представлен в следующем виде:

,

,

;

гдечисло уравнений, -номер шага,

.

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

Примеры.

Пример задачи из электротехники (линейная система)

Пусть необходимо исследовать переходные процессы и выявить стационарные состояния электротехнической цепи:

Якобиан имеет вид:

Варианты возможных переходных процессов представлены ниже.

Устойчивый узел

Начальные условия:

Собственные значения:

Параметры интегрирования:

нелинейная система анализ уравнение

Неустойчивый узел

Начальные условия:

Собственные значения:

Параметры интегрирования:

Устойчивый фокус

Начальные условия:

Собственные значения:

Параметры интегрирования:

Неустойчивый фокус

Начальные условия:

Собственные значения:

Параметры интегрирования:

" Центр"

Начальные условия:

Собственные значения:

Параметры интегрирования:

" Седло"

Реальной схема, описываемая исследуемой системой уравнений, не может иметь стационарную точку типа «седло».

Приведенный пример был реализован в системе Тополог с помощью следующего описания модели:

ИМЯ «Решение системы дифференциальных уравнений»

ДАННЫЕ

{параметры модели}

E0: вещ R1: вещ R2: вещ L: вещ C: вещ k: вещ

{функциональные переменные и их рабочие значения}

Il: вещ Uc: вещ Iln: вещ Ucn: вещ

{параметры вычислительного эксперимента}

t: вещ tt: вещ tнач: вещ tкон: вещ h: вещ

{рабочие переменнные метода интегрирования}

di1: вещ di2: вещ di3: вещ di4: вещ

du1: вещ du2: вещ du3: вещ du4: вещ

{оптимизированный вариант}

dIldt: функция (Il, Uc) = (Il* (k — R1) — Uc + E0) / L

dUcdt: функция (Il, Uc) = (Il* (k + R2) — Uc + E0) / (C*R2)

{вектор производных}

V: вектор = (dIldt, dUcdt) |

{матрица Якоби}

J: якобиан (V)

{рабочие массивы}

Re: массив из вещ

Im: массив из вещ

{ Описание переменных и объектов для графического вывода }

i: цел n: цел

x: вещ y: вещ

mi: вещ mu: вещ Точка: точка_вещ

Il_t: ломаная_вещ

Uc_t: ломаная_вещ

IlUc: ломаная_вещ МОДЕЛЬ

{ параметры цепи }

E0 = 10.0

R1 = 1.0

R2 = 1.0

L = 1.0

C = 1.0

k = 1.5

{ начальные условия }

tнач = 0.0

Il = 0.0

Uc = 0.0

t = tнач

{ масштабные коэффициенты }

mi = 0.2

mu = 0.2

Собственные_Значения (J, 300, 1E-15) — > (Re, Im)

{ параметры вычислительного эксперимента }

h = 0.1

tкон = 10

{ установка начальной точки графиков Il (t) и Uc (t) }

Изменить_Абсциссу (t) — > (Точка) Изменить_Ординату (Il*mi) — > (Точка) Добавить_Элемент (Точка) — > (Il_t)

Изменить_Ординату (Uc*mu) — > (Точка) Добавить_Элемент (Точка) — > (Uc_t)

Обновить (графика)

t = tнач

{ расчетный цикл интегрирования }

пока t < tкон начало

{ расчет методом Рунге-Кутта 4-го порядка }

Iln = Il Ucn = Uc

{ 1-й подшаг }

tt = t

di1 = h * dIldt du1 = h * dUcdt

{ 2-й подшаг }

tt = t + h/2.

Il = Iln + di½. Uc = Ucn + du½.

di2 = h * dIldt du2 = h * dUcdt

{ 3-й подшаг }

tt = t + h/2.

Il = Iln + di2/2. Uc = Ucn + du2/2.

di3 = h * dIldt du3 = h * dUcdt

{ 4-й подшаг }

tt = t + h

Il = Iln + di3 Uc = Ucn + du3

di4 = h * dIldt du4 = h * dUcdt

Il = Iln + (di1 + 2. *di2 + 2. *di3 + di4) /6.0

Uc = Ucn + (du1 + 2. *du2 + 2. *du3 + du4) /6.0

t = t + h

Изменить_Абсциссу (t) — > (Точка) Изменить_Ординату (Il*mi) — > (Точка) Добавить_Элемент (Точка) — > (Il_t)

Изменить_Ординату (Uc*mu) — > (Точка) Добавить_Элемент (Точка) — > (Uc_t)

конец Число_Элементов (Il_t) — > (n)

цикл i = 1 до n

начало Ордината (Il_t [i]) — > (y)

Ордината (Uc_t [i]) — > (x)

Добавить_точку (y, x) — > (IlUc)

конец Обновить (ГРАФИКА) КОНЕЦ Пиримеры построения фазовых портретов для нелинейных систем Дана система дифференциальных уравнеий:

Модель для системы Тополог:

ИМЯ «Решение системы дифференциальных уравнений»

ДАННЫЕ

x1: вещ

x2: вещ

x1s: вещ

x2s: вещ

x1_мин: вещ = - 1.0

x1_макс: вещ = 3.5

x2_мин: вещ = - 1.0

x2_макс: вещ = 1.5

Re₁: вещ

Re₂: вещ

R: вещ Шаг_Сетки: вещ = 0.2

{ параметры вычислительного эксперимента }

t: вещ

tt: вещ

tнач: вещ

tкон: вещ

h: вещ

{ рабочие переменнные метода интегрирования }

X1n: вещ X2n: вещ

di1: вещ di2: вещ di3: вещ di4: вещ

du1: вещ du2: вещ du3: вещ du4: вещ

dX1: функция (x1, x2) = (3 — x1 — 5 * x2) * x1

dX2: функция (x1, x2) = (2 — 4 * x1 — 2 * x2) * x2

{ вектор производных }

V: вектор = (dX1, dX2) |

{ матрица Якоби }

J: якобиан (V)

{ рабочие массивы }

Re: массив из вещ

Im: массив из вещ Решения: массив из точка_вещ

{ Описание переменных и объектов для графического вывода }

Число_Решений: цел Число_Собственных_Значений: цел Текущее_Решение: цел Фазовый_Портрет: массив из ломаная_вещ Текущая_Ломаная: ломаная_вещ Устойчивые_Узлы: массив из точка_вещ Неустойчивые_Узлы: массив из точка_вещ Устойчивые_Фокусы: массив из точка_вещ Неустойчивые_Фокусы: массив из точка_вещ Центры: массив из точка_вещ Седла: массив из точка_вещ МОДЕЛЬ

{ начальные условия }

tнач = 0.0

{ параметры вычислительного эксперимента }

h = 0.05

tкон = 2

Найти_Все_Решения (V, 1E-8, x1_мин, x1_макс,

x2_мин, x2_макс, Шаг_Сетки, Решения) Число_Элементов (Решения) — > (Число_Решений) цикл Текущее_Решение = 1 до Число_Решений начало Абсцисса (Решения [Текущее_Решение]) — > (x1)

Ордината (Решения [Текущее_Решение]) — > (x2)

Собственные_Значения (J, 300, 1E-15) — > (Re, Im)

Число_Элементов (Re) — > (Число_Собственных_Значений) если Число_Собственных_Значений = 2 то начало

Re₁ = Re [1]

Re₂ = Re [2]

R = Re₁ * Re₂

если Im = 0.0 то { вещественные СЗ }

если R < 0.0 то Добавить_Элемент (Решения [Текущее_Решение]) — >

(Седла) иначе если Re₁ < 0.0 то Добавить_Элемент (Решения [Текущее_Решение]) — >

(Устойчивые_Узлы) иначе Добавить_Элемент (Решения [Текущее_Решение]) — >

(Неустойчивые_Узлы) иначе { комплексные СЗ }

если Re₁ < 0.0 то Добавить_Элемент (Решения [Текущее_Решение]) — >

(Устойчивые_Фокусы) иначе если Re₁ = 0.0 то Добавить_Элемент (Решения [Текущее_Решение]) — >

(Центры) иначе Добавить_Элемент (Решения [Текущее_Решение]) — >

(Неустойчивые_Фокусы) конец конец

x1s = 0

пока x1s < x1_макс начало

x2s = 0

пока x2s < x2_макс начало гудок

x1 = x1s

x2 = x2s

Очистить_Массив (Текущая_Ломаная) Добавить_Точку (x1, x2) — > (Текущая_Ломаная)

t = tнач

{ расчетный цикл интегрирования }

пока t < tкон начало

{ расчет методом Рунге-Кутта 4-го порядка }

X1n = x1 X2n = x2

{ 1-й подшаг }

tt = t

di1 = h * dX1 du1 = h * dX2

{ 2-й подшаг }

tt = t + h / 2.0

x1 = X1n + di½.0×2 = X2n + du½.0

di2 = h * dX1 du2 = h * dX2

{ 3-й подшаг }

tt = t + h / 2.0

x1 = X1n + di2/2.0×2 = X2n + du2/2.0

di3 = h * dX1 du3 = h * dX2

{ 4-й подшаг }

tt = t + h

x1 = X1n + di3 x2 = X2n + du3

di4 = h * dX1 du4 = h * dX2

x1 = X1n + (di1 + 2.0*di2 + 2.0*di3 + di4) /6.0

x2 = X2n + (du1 + 2.0*du2 + 2.0*du3 + du4) /6.0

t = t + h

Добавить_Точку (x1, x2) — > (Текущая_Ломаная) конец Добавить_Элемент (Текущая_Ломаная) — >

(Фазовый_Портрет)

x2s = x2s + Шаг_Сетки конец

x1s = x1s + Шаг_Сетки конец Обновить (ГРАФИКА) КОНЕЦ

Фазовый портрет:

Дана система:

Модель для системы Тополог взята из предыдущего примера с параметрами:

x1_мин: вещ = - 1.0

x1_макс: вещ = 5.5

x2_мин: вещ = - 1.0

x2_макс: вещ = 2.5

dX1: функция (x1, x2) = (4 — 2.5 * x2) * x1

dX2: функция (x1, x2) = (-2 + x1) * x2

h = 0.01

tкон = 2

Получен фазовый портрет: