Для устойчивости самого метода проведем выбор шага интегрирования h исходя из условия:
(60).
где — собственное значение матрицы Якоби.
Для комплексного значения условие имеет вид:
(61).
Собственными значениями матрицы Якоби порядка n называют корни, где, ее характеристического уравнения, определяемого по формуле:
(62).
где, А — матрица Якоби динамической модели;
Е — единичная матрица.
Произведем расчет матрицы Якоби по формуле (55), подставляя начальные значения фазовых координат:
Тогда характеристическое уравнение имеет вид:
(64).
Вычислим корни характеристического уравнения с помощью программы MathCad, тогда собственные значения матрицы Якоби имеют вид:
Корни характеристического уравнения имеют отрицательные и нулевые значения, что говорит об устойчивости системы.
Для гидравлической системы рекомендуемый шаг интегрирования h=0.5с. Выполним проверку устойчивости численного метода Эйлера при данном шаге.
При л=0: =1;
При л=-0.760: =1.38.
Проверка условий выполняется, следовательно, шаг h=0.5 обеспечит устойчивость метода и приемлемую точность вычислений.
Решение систем дифференциальных уравнений методом Эйлера
Формула численного интегрирования неявного метода Эйлера имеет вид:
Совместное преобразование двух последних выражений приводит к записи:
(67).
где — модифицированная матрица Якоби на k+1 шаге, которая формируется по следующему правилу:
Диагональные элементы матрицы Якоби на k-ом шаге пересчитываются по формуле:
(68).
Остальные элементы не изменяются. Для матрицы размерности 3×3 получаем:
(69).
— модифицированный вектор входных воздействий на k+1 шаге, определяемый по формуле:
(70).
Решение системы уравнений (66) дает значение фазовых координат на k+1 шаге, то есть в момент времени tk+1.
Алгоритм неявного метода Эйлера с постоянным шагом интегрирования h:
- 1) задание шага интегрирования h;
- 2) задание начальных значений фазовых переменных при t0=0;
- 3) вычисление времени tk+1=tk+h, где k=0,1,2… ;
- 4) вычисление модифицированных матриц и на k+1 шаге;
- 5) решение системы уравнений (66) с целью определения в момент времени tk+1;
- 6) переход к этапу (3) до тех пор, пока в случае устойчивой системы фазовые координаты не достигнут состояния конечного значения .
Начальные значения вектора определяются на основании входных воздействий системы. В качестве начальных значений фазовых переменных берем вектор начальных значений .
Рисунок 7 — Графики фазовых координат M (n)0, M (n)1, M (n)2, M (n)3.
Рисунок 8 — Переходный процесс гидросистемы.