РАБОТА 2. МОДЕЛИРОВАНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Цель работы: освоение методики моделирования линейных дифференциальных уравнений.

1. ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ

Многие физические процессы, такие как колебания маятника, движение стрелки гальванометра, изменение высоты при посадке самолета, процессы в электрическом колебательном контуре могут быть описаны линейным однородным дифференциальным уравнением второго порядка

y''+a1y'+a0y=0,

здесь a0, a1 - постоянные коэффиценты, определяющие характер процесса.

Компьютерное моделирование следует из представления коэффициентов правой и левой частей уравнения (передаточной функции) параметрами: N=1, D=[1 a1 a0]. Амплитуда переменной y(t) зависит от начальных условий, например, от начального отклонения маятника от положения равновесия и его начальной скорости. Иными словами, от вектора начального состояния x=[y0 ý0].


%toolbox control, N=1, D=[1 2 2], 
t=time(5 1000), u=zero(t), x=[1 1], y=lsims(N D x u t), [t y]=??

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

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

p2+a1p+a0=0
.

Вещественные корни λ1=(√D-a1)/2, λ2=-(√D+a1)/2 зависят от дискриминанта D=a12-4a0. Если он отрицателен D=-1|D|, то в формуле появляются сомножители j=√-1 и √|D|.

Если корни вещественные и различные λ1, λ2, то решение имеет вид

y(t)=с1eλ1t2eλ2t.

Этот случай отвечает у маятников их верхнему положению. Постоянные c1 и c2 находят, подставляя начальные условия в выражениях для y(t) и y'(t) при t = 0. Получаем систему уравнений c1+c2 = y(0) и λ1c1– λ2c2 = y'(0).

Для комплексных показателей λ1=α+jω, λ2=α-jω (нижнее положение у маятников) решение сводится к виду:

y(t)=eαt1sin(ωt)+с2cos(ωt)).

Зависимость постоянных от начальных условий c1=(y'(0)-αy(0))/ω, c2=y(0) следует из тех же соображений, что и выше.

Пример. Дано уравнение

ÿ + 2ý + 2y = 0, y0 = ý0 = 1.

Его характеристическое уравнение p2 + 2p + 2 = 0

имеет корни p1,2 = –1 ± j. Следовательно, общее решение будет следующим:

y(t) = c1 e–t sin t + c2 e–t cos t.

Дифференцируя, находим выражение для производной:

ý(t) = – (c1 + c2) e–t sin t + (c1 – c2) e–t cos t.

При t = 0, с учетом начальных условий, получаем из системы уравнений значения коэффициентов с1 = 2, с2 = 1. Следовательно, y(t) = e–t (2sin t + cos t). График кривой строится ниже.


t=time(5), y=fun(exp(-t)*(2*sin(t)+cos(t))), [t y]=??

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

В нашу прагматичную эпоху аналитические задачи помогают решать сетевые роботы, эти возможности заложены в аватор сайта, см. меню. Аналитические построения выполняет и скриптовый исполнитель:


%toolbox control, y=solve("y''+2y'+2y=0"), 
puts(y), t=time(5), c1=2, c2=1, y=fun(eval(y)), plot(t y)

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

Моделирование системы уравнений. Многие технические объекты можно описать системой n линейных дифференциальных уравнений первого порядка

dx/dt=Ax+bu, y=cx+du,

пусть

A =
3
-2
-6
-1

b =
1
0

c =
1
0

где x=x(t) - вектор состояния (пространства Rn, n=2), u=u(t) - управление (входной сигнал), y=y(t) - реакция на него, A – квадратная матрица коэффициентов, b – вектор-столбец, c - вектор-строка определяет, какая из переменных вектора состояния является измеряемой величиной.

Матрица S=[A b ; c d] часто используется для описания таких систем, неинерционная связь обычно отсутствует d=0. При моделировании летательного аппарата составляющими вектора состояния являются текущие координаты самолета и скорости их изменения, матрица A характеризует динамику канала тангажа, курса или крена, а управление описывает воздействия, формируемые летчиком или автопилотом.

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

Для рассматриваемых выше параметров запишем

y'1=3y1-2y2+5, y'2=-6y1-y2.

После дифференцирования первого уравнения, получаем y''1=3y'1-2y'2, подставляя выражения для производных, имеем y''1=21y1-4y2+15. Исключим вторую переменную, используя первое уравнение, получим y''1-2y'1-15y1=5.

Еще один метод построения такого дифференциального уравнения состоит в нахождении передаточной функции системы Q(p)=c(pI-A)-1b, где

pI-A =
p-3
2
6
p+1

Найдем обратную матрицу делением матрицы алгебраических дополнений (алгебраические дополнения матрицы второго порядка получают перестановкой диагональных элементов и инверсией знака элементов вне главной диагонали) на определитель d(p)=(p-3)(p+1)-12

(pI-A)-1 =
(p+1)/d(p)
-2/d(p)
-6/d(p)
(p-3)/d(p)

Умножение на вектор-строку с=[1 0] оставляет от инверсной матрицы только первую строку, домножение ее на вектор входа дает Q(p)=(p+1)/d(p)=(p+1)/(p2-2p-15), отсюда

y''1-2y'1-15y1=u'+u=5,

так как на входе константа u(t)=5, то u'(t)=0. Как видно, результат обоих методов составления дифференциального уравнения одинаков.

Дифференциальное уравнение решается известным нам способом, однако следует учесть, что в качестве начальных условий заданы y1(0) и y2(0), а не y1(0) и y'1(0). Начальное значение производной первого сигнала y'1(0) отыскивается по заданным условиям из первого дифференциального уравнения исходной системы, в данном случае y'1(0)=3y1(0)-2y2(0)+5.

Для построения фазового портрета динамической системы расчет ведется по двум выходам системы, отвечающим двум ее состояниям. При нулевых начальных состояниях численное моделирование учитывает только вход u(t)=5


A=[3 -2; -6 -1], b=[1 0], 
% СОБСТВЕННЫЕ ЗНАЧЕНИЯ МАТРИЦЫ A (КОРНИ),
L=eig(A), diags(L), tr(L), L[0]=?, L[1]=?,     
t=time(0.5 1000), u=one(t), u={5*u}, 
c=[1 0], S=ss(A b c), y1=lsim(S u(t)), 
c=[0 1], S=ss(A b c), y2=lsim(S u(t)), plot(y1 y2)

Полезно проверять также собственные значения матрицы A, отвечающие корням ее характеристического уравнения α+jω, α-jω. Положительные вещественные части α>0 свидетельствуют о неустойчивости системы. Для устойчивых систем время моделирования берется значительно большим, чтобы дать развиться процессу, не 0.5, а больше, чтобы получить график, как на рисунке.

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

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

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

y'' + 2y' + 3y = 0, y(0) = 2, y'(0) = 4.

Для построения схемы моделирования воспользуемся методом понижения производной (методом Кельвина). В нем можно выделить четыре шага.

Шаг 1. Разрешаем уравнение относительно старшей производной. В частности для этого однородного дифференциального уравнения получаем y'' = – 2y' – 3y.

Шаг 2. Полагаем старшую производную y'' известной и выполняем ее последовательное интегрирование цепочкой интеграторов, получая все низшие производные и саму переменную y. В случае данного уравнения для этого потребуется два последовательно включенных интегратора, на выходе которых получим сигналы ý и y.

Шаг 3. Замыкаем цепочку интеграторов, формируя старшую производную на основе уравнения, полученного на первом шаге. В нашем примере для этого потребуется сумматор, складывающий сигналы y и y', домноженные на коэффиценты -2 и -3.

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

Итоговая схема содержит два интегратора, два масштабных усилителя и сумматор (обозначается кружочком).



Rambler's Top100