Качественный анализ нелинейных систем
Постановка
задачи
Описание модели
1
модель
описывается системой дифференциальных уравнений 2-ого порядка;
1
задается
область определения значений переменных в виде системы неравенств:
Исследование стационарных точек
1
в
заданной области значений переменных определяются стационарные точки системы
путем решения системы алгебраических уравнений:
1
выполняется
линеаризация системы уравнений в стационарных точках.
Точки классифицируются путем нахождения собственных значений
матрицы (якобиана) линеаризованной системы. Возможны шесть типов точек.
Обозначения:
- собственные
значения якобиана ().
1.
Устойчивый
узел:
Собственные значения действительные (мнимая часть равна нулю),
меньшие нуля:
1.
Неустойчивый
узел:
Собственные значения действительные, большие нуля:
1.
Устойчивый
фокус:
Собственные значения комплексно-сопряженные, действительная
часть меньше нуля:
1.
Неустойчивый
фокус:
Собственные значения комплексно-сопряженные, действительная
часть больше нуля:
1.
"Центр”:
Собственные значения комплексно-сопряженные, действительная
часть равна нулю:
1.
"Седло”:
Собственные значения действительные, с разными знаками:
Построение фазовых траекторий
Фазовые траектории (портреты) строятся в окрестностях
стационарных точек с помощью метода численного интегрирования.
Практическая реализация задачи
Задача реализована в среде системы Тополог.
Описание модели
Описание системы уравнений выполняется с использованием
стандартных типов системы Тополог: функция и вектор. В виде функции описывается
только правая часть уравнения.
Пример: пусть мы имеем систему из двух нелинейных уравнений
вида:
тогда описание данной системы на языке системы Тополог будет
выглядеть следующим образом:
данные
x: вещ
y: вещ
F1: функция (x, y) = x^2/100 + y^2/25 - 1
F2: функция (x, y) = y - 10 * Sin (X)
VF: вектор = (F1, F2) |
J: якобиан (VF)
Исследование стационарных точек
Решение системы уравнений в заданной прямоугольной области выполняется стандартным
операционным блоком системы Найти_Все_Решения. Он имеет следующий
синтаксис:
Найти_Все_Решения (<вектор функций>, <погрешность>,
<граница 1>, <граница 2>,<граница 3>,
<граница 4>,
<шаг сетки>, <массив решений>)
Решения ищутся путем запуска метода Ньютона из узлов сетки,
наложенной на прямоугольную область поиска, и накопления найденных различных
решений. Прямоугольная область поиска задается параметрами: <граница 1>,
<граница 2>,<граница 3>, <граница 4>. Шаг сетки равномерный
по обеим координатам и задается параметром <шаг сетки>. После завершения
работы блока все решения скапливаются в массиве точек, описанным как:
Решения: массив из Точка_вещ
Если при поиске решения из некоторого узла сетки метод выходит за
границу прямоугольной области или якобиан системы уравнений на текущей итерации
становится неопределенным, то итерации прекращаются и метод запускается из
следующего по очереди узла.
Для получения собственных значений якобиана системы в найденных
стационарных точках используется стандартный операционный блок системы Тополог Собственные_Значения.
Он имеет следующий синтаксис:
Собственные_Значения (<матрица>,
<максимальное число итераций>, <погрешность>) - >
(<массив реальных частей>,
<массив мнимых частей>)
На выходе данный блок формирует два массива, содержащих реальные и
мнимые части собственных значений квадратной матрицы (в нашем случае якобиана).
Численные методы нахождения собственных значений матриц делятся на
два типа: прямые и итерационные. Прямые методы находят собственные значения
путем решения характеристического уравнения матрицы с помощью одного из
численных методов (метод Ньютона (касательных), метод секущих, метод хорд и
др.). Итерационные методы при помощи преобразований, не меняющих собственных
значений матрицы, приводят ее к диагональному (в случае действительных
собственных значений) или блочно-диагональному (при наличии
комплексно-сопряженных собственных значений) виду. Значения из диагонали полученной
матрицы и являются собственными значениями исходной матрицы.
В основу блока Собственные_Значения заложен итерационный
метод нахождения собственных значений по методу Якоби с понижением нормы для
действительных матриц.
Основная суть метода Якоби:
Пусть задана действительная матрица размера и матрица (или ) формируется как последовательность
преобразований . При определенном выборе матрица будет диагональной с элементами, равными собственным значениям
исходной матрицы ( и представляют соответственно правую и
левую систему собственных векторов ).
Для реализации такого преобразования матрицы можно выбирать в виде произведения двух
матриц , где матрица определяет так называемое вращение, а - сдвиг или комплексное вращение.
Элементы этих двух матриц удовлетворяют соотношениям:
Пары индексов (,) выбираем на каждом шаге итерации путем циклического перебора.
Параметр вращения можно определить из соотношения , причем нужно выбрать такое значение , чтобы после очередного преобразования
вращения евклидова норма -ого столбца преобразованной матрицы была
не меньше нормы -ого столбца. Параметр сдвига выбирается так, чтобы на каждом шаге
преобразований уменьшалась евклидова норма матрицы , равной
.
Параметр сдвига удовлетворяет соотношению
; ;
;
.
По окончанию итерационного процесса матрица будет нормальной, причем , и либо , либо для всех значений (, ), а матрица будет блочно-диагональной: блоки размера
() будут содержать действительные
собственные значения, а () - комплексно-сопряженные.
Построение
фазовых траекторий
Фазовые траектории образуются в процессе решения задачи
интегрирования системы дифференциальных уравнений, описывающих исследуемую
модель. Для этого использован метод численного интегрирования Рунге-Кутта
четвертого порядка.
Алгоритм метода может быть представлен в следующем виде:
,,
,,
;
где -число уравнений, -номер шага,
,
,
,
.
При незначительном, по сравнению с методом Эйлера, усложнении
расчетных формул метод Рунге-Кутта обеспечивает значительно более высокую
точность, имеющую порядок .
Примеры.
Пример задачи из электротехники (линейная система)
Пусть необходимо исследовать переходные процессы и выявить
стационарные состояния электротехнической цепи:
Якобиан имеет вид:
Варианты возможных переходных процессов представлены ниже.
1.
Устойчивый
узел
Начальные условия:
Собственные значения:
Параметры интегрирования:
нелинейная система анализ уравнение
1.
Неустойчивый
узел
Начальные условия:
Собственные значения:
Параметры интегрирования:
1.
Устойчивый
фокус
Начальные условия:
Собственные значения:
Параметры интегрирования:
1.
Неустойчивый
фокус
Начальные условия:
Собственные значения:
Параметры интегрирования:
1.
"Центр”
Начальные условия:
Собственные значения:
Параметры интегрирования:
1.
"Седло”
Реальной схема, описываемая исследуемой системой уравнений,
не может иметь стационарную точку типа "седло”.
Приведенный пример был реализован в системе Тополог с помощью
следующего описания модели:
ИМЯ "Решение системы дифференциальных уравнений"
ДАННЫЕ
{параметры модели}
E0: вещ R1: вещ R2: вещ L: вещ C: вещ k: вещ
{функциональные переменные и их рабочие значения}
Il: вещ Uc: вещ Iln: вещ Ucn: вещ
{параметры вычислительного эксперимента}
t: вещ tt: вещ tнач: вещ tкон: вещ h: вещ
{рабочие переменнные метода интегрирования}
di1: вещ di2: вещ di3: вещ di4: вещ
du1: вещ du2: вещ du3: вещ du4: вещ
{оптимизированный вариант}
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= h * dIldt du1 = h * dUcdt
{ 2-й подшаг }= t + h/2.= Iln + di1/2. Uc = Ucn +
du1/2.= h * dIldt du2 = h * dUcdt
{ 3-й подшаг }= t + h/2.= Iln + di2/2. Uc = Ucn +
du2/2.= h * dIldt du3 = h * dUcdt
{ 4-й подшаг }= t + h= Iln + di3 Uc = Ucn + du3=
h * dIldt du4 = h * dUcdt= Iln + (di1 + 2. *di2 + 2. *di3 + di4) /6.0= Ucn +
(du1 + 2. *du2 + 2. *du3 + du4) /6.0 = 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)
конец
Обновить (ГРАФИКА)
КОНЕЦ
Пиримеры построения фазовых портретов для нелинейных систем
Модель для системы Тополог:
ИМЯ "Решение системы дифференциальных уравнений"
ДАННЫЕ: вещ: вещs: вещs: вещ_мин: вещ = - 1.0_макс: вещ = 3.5_мин:
вещ = - 1.0_макс: вещ = 1.5_1: вещ_2: вещ: вещ
Шаг_Сетки: вещ = 0.2
{ параметры вычислительного эксперимента }: вещ: вещнач: вещкон:
вещ: вещ
{ рабочие переменнные метода интегрирования }n: вещ X2n: вещ: вещ
di2: вещ di3: вещ di4: вещ: вещ du2: вещ du3: вещ du4: вещ: функция (x1, x2) =
(3 - x1 - 5 * x2) * x1: функция (x1, x2) = (2 - 4 * x1 - 2 * x2) * x2
{ вектор производных }: вектор = (dX1, dX2) |
{ матрица Якоби }: якобиан (V)
{ рабочие массивы }: массив из вещ: массив из вещ
Решения: массив из точка_вещ
{ Описание переменных и объектов для графического вывода }
Число_Решений: цел
Число_Собственных_Значений: цел
Текущее_Решение: цел
Фазовый_Портрет: массив из ломаная_вещ
Текущая_Ломаная: ломаная_вещ
Устойчивые_Узлы: массив из точка_вещ
Неустойчивые_Узлы: массив из точка_вещ
Устойчивые_Фокусы: массив из точка_вещ
Неустойчивые_Фокусы: массив из точка_вещ
Центры: массив из точка_вещ
Седла: массив из точка_вещ
МОДЕЛЬ
{ начальные условия }нач = 0.0
{ параметры вычислительного эксперимента }= 0.05кон = 2
Найти_Все_Решения (V, 1E-8, x1_мин, x1_макс,_мин, x2_макс,
Шаг_Сетки, Решения)
Число_Элементов (Решения) - > (Число_Решений)
цикл Текущее_Решение = 1 до Число_Решений
начало
Абсцисса (Решения [Текущее_Решение]) - > (x1)
Ордината (Решения [Текущее_Решение]) - > (x2)
Собственные_Значения (J, 300, 1E-15) - > (Re, Im)
Число_Элементов (Re) - > (Число_Собственных_Значений)
если Число_Собственных_Значений = 2 то
начало
Re_1 = Re [1]_2 = Re [2]= Re_1 * Re_2
если Im [1] = 0.0 то { вещественные СЗ }
если R < 0.0 то
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Седла)
иначе
если Re_1 < 0.0 то
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Устойчивые_Узлы)
иначе
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Неустойчивые_Узлы)
иначе { комплексные СЗ }
если Re_1 < 0.0 то
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Устойчивые_Фокусы)
иначе
если Re_1 = 0.0 то
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Центры)
иначе
Добавить_Элемент (Решения [Текущее_Решение]) - >
(Неустойчивые_Фокусы)
конец
конецs = 0
пока x1s < x1_макс
началоs = 0
пока x2s < x2_макс
начало
гудок= x1s= x2s
Очистить_Массив (Текущая_Ломаная)
Добавить_Точку (x1, x2) - > (Текущая_Ломаная)= tнач
{ расчетный цикл интегрирования }
пока t < tкон
{ расчет методом Рунге-Кутта 4-го порядка }n = x1 X2n = x2
{ 1-й подшаг }
tt = t= h * dX1 du1 = h * dX2
{ 2-й подшаг }= t + h / 2.0= X1n + di1/2.0 x2 = X2n +
du1/2.0= h * dX1 du2 = h * dX2
{ 3-й подшаг }= t + h / 2.0= X1n + di2/2.0 x2 = X2n +
du2/2.0= h * dX1 du3 = h * dX2
{ 4-й подшаг }= t + h= X1n + di3 x2 = X2n + du3= h *
dX1 du4 = h * dX2= X1n + (di1 + 2.0*di2 + 2.0*di3 + di4) /6.0= X2n + (du1 +
2.0*du2 + 2.0*du3 + du4) /6.0
t = t + h
Добавить_Точку (x1, x2) - > (Текущая_Ломаная)
конец
Добавить_Элемент (Текущая_Ломаная) - >
(Фазовый_Портрет)s = x2s + Шаг_Сетки
конецs = x1s + Шаг_Сетки
конец
Обновить (ГРАФИКА)
КОНЕЦ
Фазовый портрет:
Дана система:
Модель для системы Тополог взята из предыдущего примера с
параметрами:_мин: вещ = - 1.0_макс: вещ = 5.5_мин: вещ = - 1.0_макс: вещ = 2.5
: функция (x1, x2) = (4 - 2.5 * x2) * x1: функция (x1, x2) = (-2 +
x1) * x2
= 0.01кон = 2
Получен фазовый портрет: