Основной целью компьютерного моделирования ДС, очевидно, является оценка временной эволюции системы, т.е. выявление характера изменения ее переменных состояния во времени. Таким образом, происходит восстановление динамики системы по ее математической модели. В случае непрерывных ДС этот процесс фактически выливается в процедуру численного интегрирования дифференциальных уравнений при конкретных начальных условиях на заданном временном интервале. Для дискретных систем организуется итерационная процедура, позволяющая вычислять последующее состояние системы на основании текущего состояния по заданной рекуррентной форме. Результаты вычислений необходимо представить в наглядной и удобной для анализа форме. Наибольшей иллюстративностью в этом плане обладают соответствующие графики – фазовый портрет системы (для систем 2-го и 3-го порядка) и графики изменения переменных состояния во времени (графики переходных процессов).
Отмеченные особенности моделирования ДС вызывают необходимость введения в структуру файла-сценария трех основных блоков: блока формирования исходных данных, блока вычисления и блока визуализации. В блоке исходных данных задаются значения параметров системы, ее начальное состояние и время моделирования, в вычислительном блоке на основании численной процедуры формируется массив переменных состояния на заданном временном интервале, а блок визуализации организует построение необходимых графиков.
Модель ДС формируется в отдельном файле-функции, который вызывается основной программой по его имени. Модель непрерывной системы представляется в виде модели в переменных состояния, т.е. в виде правых частей системы обыкновенных дифференциальных уравнений и записывается в соответствующем m-файле, т.н. ODE функции.
Рассмотрим пример создания ODE функции для простого математического маятника, динамика которого описывается следующей моделью в переменных состояния:
|
где x 1 – угол отклонения маятника от вертикальной оси, x 2 – угловая скорость вращения маятника относительно точки подвеса, g – ускорение свободного падения, l – длина маятника, k – коэффициент трения.
Файл-функция, соответствующий этой модели, организуется следующим образом.
% pendulum.m function y = pendulum(t, x)
g=9.81; l=1; k=0.01;
y = [x(2); -g/l*sin(x(1))-k*x(2)];
Помимо переменных состояния модель ДС часто содержит параметры, т.е. условно постоянные величины. В представленном выше примере параметрами модели маятника являются ускорение свободного падения, длина маятника и коэффициент трения. Параметры системы могут задаваться в том же файле, что и модель или в файле-сценарии. В последнем случае они передаются в файл-функцию, если они объявлены глобальными переменными. Для этого в начале файла сценария, и каждого из используемых файлов-функций записывается команда:
global a1 a2 a3 a4 …….
Целью моделирования системы может являться так называемый параметрический анализ, заключающийся в оценке влияния изменения параметров на общую динамику системы. Эта возможность должна быть предусмотрена при организации вычислительной процедуры. Значения параметров могут быть заданы и изменены как в самой программе путем присвоения им численного значения (например, l=3), так и с помощью команды интерактивного ввода input. Например, для интерактивного ввода значения длины маятника можно воспользоваться командой:
l=input(‘Введите значение длины маятника l=’)
После выполнения этой команды в командном окне MATLAB появится соответствующий запрос. Рационально задавать все параметры ДС в начале файла-сценария и при необходимости изменять их значение в теле этого файла с помощью соответствующих команд присвоения или интерактивного ввода.
Численное интегрирование ДС в MATLAB осуществляется с помощью решателей, каждый из которых реализует определенный алгоритм численного интегрирования системы обыкновенных дифференциальных уравнений. Под численным интегрированием системы обыкновенных дифференциальных уравнений (ОДУ) понимают нахождение приближенного решения задачи Коши, т.е. приближенных значений функций xi (t) на заданном временном интервале. Применение того или иного решателя определяется свойствами ДС и требованиями к точности полученных результатов и быстродействию вычислительной процедуры. В MATLAB используются различные решатели. Для начальной пробы можно порекомендовать решатель ode45, реализующий одношаговый явный метод Рунге-Кутта 4-го и 5-го порядка и обычно демонстрирующий хорошие результаты.
Вызов решателя, осуществляющего выполнение процесса численного интегрирования системы ОДУ, производится командой
[t,x] = solver(@Fname,tspan,x0)
Здесь: solver – имя решателя (например, ode45), Fname – имя файла-функции, содержащего модель динамической системы; tspan – интервал времени интегрирования, задаваемый в форме вектора-строки [tmin tmax]; x0 – вектор начальных условий, число элементов которого должно совпадать с размерностью системы уравнений; t – выходной массив времени (вектор-столбец); x – выходной массив переменных состояния системы, в общем случае представляющий собой матрицу с количеством столбцов равным порядку системы.
Например, для численного интегрирования маятника на временном интервале от 0 до 20 секунд при единичных начальных условиях следует воспользоваться командой:
[t,x] = ode45(@pendulum,[0,20],[1,1])
Результатом выполнения этой команды будут вектор столбец t и матрица x с двумя столбцами, в каждом из которых расположены значения соответствующей переменной на заданном временном интервале. Необходимо заметить, что число элементов массива времени и число строк массива переменных состояния заранее неизвестно, а определяется самой процедурой численного интегрирования. Соответствие времени и значений переменных задается построчно, т.е. в i-ой строке массива x находятся значения переменных в момент времени t[i].
Интервал времени интегрирования tspan и вектор начальных условий x0 могут быть заданы непосредственно в команде вызова решателя или перед ней с помощью соответствующих команд присвоения или интерактивного ввода. Для двухмерных систем начальные условия удобно задавать в поле фазового портрета с помощью команды графического ввода данных ginput(1). Эта команда позволяет вводить данные с помощью мыши из активного графического окна. Для этого сначала создается графическое окно, определяется диапазон изменения переменных в виде пределов осей абсцисс и ординат графика. После этого непосредственно используется команда графического ввода. После ее вызова в текущих осях появляется перекрестие, управляемое движением мыши. Щелчок левой кнопки мыши в нужной точке плоскости, формирует вектор с координатами [x y], т.е. вектор начальных условий. Данный алгоритм действий реализуется следующей последовательностью команд:
figure(1)
axis([-x1max x1max -x2max x2max]); x0 = ginput(1)
Двумя основными формами графического представления динамики системы являются графики переходных процессов переменных состояния, т.е. графики функций xi (t), полученные при определенных начальных условиях, и ее фазовый портрет как совокупность фазовых траекторий, полученных для набора начальных условий. Построение графиков переходных процессов не требуется ничего большего, кроме как правильного использования команды plot. Например:
plot(t,x(:,1)) – график функции x 1(t); plot(t,x(:,2)) – график функции x 2(t); plot(t,x(:,1),’r’ t,x(:,2),’b’) – графики функций x 1(t) и x 2(t) в одних координатных осях красным и синим цветом соответственно;
Одна фазовая траектория также строится весьма просто:
plot(x(:,1),x(:,2)) – график траектории на фазовой плоскости; plot3(x(:,1),x(:,2),x(:,3))) – график траектории в трехмерном фазовом пространстве.
Для получения фазового портрета необходимо производить численное интегрирование системы на множестве начальных условий и строить соответствующие фазовые траектории. Таким образом, в файле-сценарии должна быть организована процедура многократного вызова решателя и применения команды plot. Разумеется многократный запуск программы на решение – процесс весьма утомительный. Его можно успешно избежать, если использовать в теле файла-сценария операторы цикла.
Любая фазовая траектория характеризуют движение состояния системы во времени. Очень наглядно видеть это движение при визуализации результатов. Для отображения движения точки на плоскости или в трехмерном пространстве в MATLAB есть команды comet и comet3. Движущая точка представляется ядром кометы с выделенным цветом хвостом. Синтаксис этой команд для построения фазовых траекторий:
comet(x(:,1),x(:,2)) comet3(x(:,1),x(:,2),,x(:,3))).
Вот фактически все минимальные сведения для создания программ моделирования ДС. Структура файла-сценария задается разработчиком в соответствии с целями вычислительного эксперимента.