Матлаб. К. Ю. Петрова введение в matlab учебное пособие

• Название: К. Ю. Петрова введение в matlab учебное пособие
• Анкор: Матлаб
• Дата: 15.09.2022
• Размер: 2.57 Mb.
• Имя файла: mironovsky_petrova_matlab.pdf
• Тип: Учебное пособие; #678847
• страница: 14 из 18

bilinear. Ее синтаксис имеет вид:
[bz, az]= bilinear(b, a, F), где F – частота дискретизации в герцах.
Цифровые фильтры можно строить непосредственно, не опираясь на аналоговый прототип. Для этого имеются следующие команды:
butter – синтез цифрового фильтра Баттерворта.
cheby1 – синтез цифрового фильтра Чебышева первого типа.
cheby2 – синтез цифрового фильтра Чебышева второго типа.
ellip – синтез цифрового эллиптического фильтра.
Число входных аргументов этих команд колеблется от трех (у фильтра Баттерворта) до пяти
(у эллиптического фильтра).
Перечисленные команды позволяют рассчитывать не только цифровые, но и аналоговые фильтры, для них в качестве последнего входного параметра надо указать опцию 's'
. Например, по команде butter(n, w0, type, ' s')
будет построен аналоговый фильтр Баттерворта порядка n с частотой среза

0
Широкий набор средств для моделирования линейных и нелинейных дискретных систем предоставляет SIMULINK. Он содержит специальную библиотеку дискретных элементов, которая поддерживает все виды описания дискретных систем, а также допускает возможность построения и моделирования гибридных систем, содержащих непрерывные и дискретные элементы.
Решение дифференциальных уравнений
Решение задачи Коши
Задачей Коши называется задача о решении обыкновенного дифференциального уравнения с известными начальными условиями. В MATLAB имеются три возможности для решения задачи
Коши, не считая моделирования в SIMULINK.
Первая из них касается численного решения линейных дифференциальных уравнений с известной правой частью или систем таких уравнений. Оно может быть выполнено с помощью команды lsim (для решения однородных уравнений достаточно команды initial).
Вторая возможность – аналитическое решение линейных и простых нелинейных дифференциальных уравнений с помощью решателя dsolve тулбокса SYMBOLIC.
Пример. Пусть требуется решить линейное уравнение второго порядка
0 2


y
a
y

В командной строке набираем
:
>>y=dsolve('D2y+a^2*y=0'),
MATLAB выдает ответ:
y=C1*sin(a*t)+C2*cos(a*t).
Это означает, что общее решение данного дифференциального уравнения имеет вид cos sin
2 1
at
C
at
C
y


Задав начальные условия
,
,
2 0
0
a
y
y



получаем задачу Коши. Для ее решения набираем:
>>y=dsolve('D2y+a^2*y=0, y(0)=2, Dy(0)=a')
Получаем ответ: y = sin(a*t)+2*cos(a*t)
Аналогично решатель dsolve применяют для систем дифференциальных уравнений.
Пример. Требуется решить систему линейных дифференциальных уравнений

92
u
w
w
v
v
u
2
,
2
,
2







с начальными условиями
3
,
0 0
0 0



w
v
u
В командной строке набираем:
>> S=dsolve('Du=2*v, Dv=2*w, Dw=-2*u','u(0)=0, v(0)=0, w(0)=3').
MATLAB выдает сообщение о структуре решения:
S = u: [1x1 sym] v: [1x1 sym] w: [1x1 sym]
Для доступа к полям структуры S набираем:
>>u = S.u, v = S.v, w = S.w
MATLAB выдает ответ: u=exp(-2*t)+3^(1/2)*exp(t)*sin(3^(1/2)*t)-exp(t)*cos(3^(1/2)*t) v =-exp(-2*t)+3^(1/2)*exp(t)*sin(3^(1/2)*t)+exp(t)*cos(3^(1/2)*t) w =exp(-2*t)+2*exp(t)*cos(3^(1/2)*t).
Для записи этих формул в более удобном виде можно воспользоваться командой pretty. Еще лучше перенести их через клипборд в пакет MAPLE, где они принимают вид:

u


e
(
)

2 t
3 e
t
(
)
sin
3 t
e
t
(
)
cos
3 t

v



e
(
)

2 t
3 e
t
(
)
sin
3 t
e
t
(
)
cos
3 t

w

e
(
)

2 t
2 e
t
(
)
cos
3 t
Решим ту же систему, заменив начальные условия краевыми
:
)
1
(
,
1
)
0
(
,
1
)
0
(
2





e
w
v
u
>> S=dsolve('Du=2*v, Dv=2*w, Dw=-2*u', 'u(0)=1, v(0)=-1, w(1)=exp(-2)');
>>u = S.u, v = S.v, w = S.w
Получаем ответ:
u =exp(-2*t), v =-exp(-2*t), w =exp(-2*t).
Третья возможность – численное решение нелинейных дифференциальных уравнений с помощью команд типа ode23 и ode45. Буквенная часть названия этих команд – сокращение от
Ordinary Differential Equation, цифры указывают порядок используемой версии метода Рунге-
Кутта. Команда ode45 дает более точное решение, но требует больше времени. Основная модификация команды ode23 имеет вид
[t, X] = ode23('fun', [T0 T1]
, Х0).
Она обеспечивает решение системы дифференциальных уравнений, записанных в форме Коши
0 0
)
(
),
,
(
X
T
X
t
X
F
X



на интервале времени
1 0
T
t
T


Результатом ее выполнения является массив отсчетов времени t и соответствующий им массив значений
X
. Для того чтобы команда была выполнена, надо предварительно составить программу для вычисления вектор-функции
)
,
(
t
X
F
, стоящей в правой части дифференциального уравнения. Эта программа должна быть оформлена в виде m-файла, которому присваивается любое имя, например,
'fun'
Пример. Воспользуемся командой ode23 для моделирования нелинейного уравнения
0.
(0)
,
1
(0)
,
0,8
sin
3




y
y
y
y
y



Перепишем его в виде системы двух уравнений, обозначив х
1
= y,
y
x


2
:

93
,
sin
0,8
,
1 3
1 2
2 1
x
x
x
x
x












0 1
X
0
Функцию для вычисления правых частей оформляем в виде m-файла с именем fun: function F = fun (t, X)
F = [X(2); .8*X(1)^3-sin(X(1))];
Далее задаем численные значения параметров
>>
Т0 = 0, Т1 = 20, X0 = [1; 0];
после чего уже может быть выполнена основная команда
>> [t, X] = ode23('fun', [T0 T1], X0).
Ее результатом будут одномерный массив времени t на интервале от 0 до 20 секунд и двумерный массив Х, содержащий значения x
1
(t),
x
2
(t). Как правило, шаг времени – переменный.
График результата может быть получен командой plot(t, X).
Для последующего перехода к равномерной сетке времени (если это необходимо), можно использовать команду interp1:
>>tt = 0:0.1:20; XX = interp1(t, X, tt, 'spline');
результатом которой будет массив переменных
ХХ
для равноотстоящих моментов времени tt
Другая возможность получения равномерного шага связана с использованием команды
deval . Для этого слегка изменим синтаксис команды ode23:
>> SOL = ode23(@fun,[T0 T1], X0)
Теперь результатом будет структура
SOL
, о чем свидетельствует сообщение:
SOL = solver: 'ode23' extdata: [1x1 struct] x: [1x90 double] y: [2x90 double]
Структура
SOL имеет поля x и
y.
Поле x содержит набор отсчетов времени, а поле y
–
массив значений вектора
)
(t
X
Доступ к полям производится командами
SOL.x и
SOL.y
Первый аргумент команды deval – структура
SOL
, а второй – вектор точек, в которых нужно вычислить аппроксимацию решения.
Задаем равномерную сетку по времени и строим график:
>>t=0:.1:20; y = deval(SOL,t); figure(2), plot(t,y,
’o')
Результат приведен на рис. 5.4.
0 5
10 15 20
-1
-0.5 0
0.5 1
t
X
Рис. 5.4
В качестве дополнительного аргумента в команде ode23 можно указать требуемую точность решения eps (по умолчанию eps = 0.001). Команда ode45 выполняется аналогично и имеет такие же модификации. Обе команды можно использовать для моделирования нелинейных систем

94 автоматического управления, если в правой части дифференциальных уравнений учесть управляющее воздействие, задавая его как явную функцию времени.
В MATLAB существуют и другие решатели дифференциальных уравнений, например,
ode15s, ode23s, ode23t. Команды с буквой S предназначены для моделирования жестких (Stiff) дифференциальных уравнений, буква Т указывает на использование метода трапеций.
Большинство из них способны решать и дифференциально-алгебраические системы уравнения вида
),
,
(
)
,
(
t
X
F
X
t
X
M


где M – матрица, F – вектор-функция.
Решение краевых задач
В предыдущем пункте речь шла о решении задачи Коши, когда заданы дифференциальные уравнения и начальные условия. Значительно большую трудность представляют краевые задачи, когда задаются начальные и конечные условия или некоторые их комбинации.
Решение краевых нелинейных дифференциальных уравнений в MATLAB выполняется с помощью команды bvp4c (читается boundary value problem for continous). Она решает двухточечную краевую задачу для обыкновенных дифференциальных уравнений методом коллокаций.
Вызов в формате sol=bvp4c (F, G, solinit)
интегрирует систему обыкновенных дифференциальных уравнений вида
)
,
( Y
x
F
Y


на интервале
 
b
a,
, с краевыми условиями вида
0
))
(
),
(
(

b
Y
a
Y
G
. В качестве входных параметров команда bvp4c
принимает три аргумента – имя функции F для формирования правых частей дифференциального уравнения
)
,
( Y
x
F
, имя функции
G для вычисления вектора краевых условий
))
(
),
(
(
b
Y
a
Y
G
и структуру solinit, содержащую начальное приближение решения.
Пример. Продемонстрируем смысл входных параметров решателя bvp4c и вид возвращаемого результата на стандартной демо-функции twobvp. В ней требуется найти решения дифференциального уравнения
0



y
y
, которые удовлетворяют краевым условиям
2
)
4
(
,
0
)
0
(



y
y
Представим исходное уравнение в виде
)
,
( Y
x
F
Y


Введем обозначения


T
y
y
Y
y
y
y
y
2 1
2 1
,
,




. Тогда уравнение можно записать как






1 2
2 1
y
y
y
y


, т.е.








1 2
y
y
Y
Правая часть этого уравнения формируется в функции twoode, она не зависит явно от x .
На каждом шаге итерационного алгоритма, решающего краевую задачу, проверяется выполнение краевых условий (boundary conditions). Для этого значения переменных в начальной и конечной точках
)
(a
Y
и
)
(b
Y
подаются на вход специальной функции, вычисляющей вектор
))
(
),
(
(
b
Y
a
Y
G
. В нашем случае это вектор


T
y
y
2
)
4
(
)
0
(

. Он зависит только от первых компонент векторов
)
(a
Y
и
)
(b
Y
, так как в краевые условия не входят производные, и вычисляется в функции twobc
Тексты функций twoode и twobc приведены ниже: function dydx = twoode(x,y) dydx = [ y(2); -abs(y(1)) ]; function res = twobc(ya,yb) res = [ ya(1); yb(1) + 2 ];
Структура solinit имеет поля x
и y.
Поле x
содержит набор отсчетов независимой переменной, а поле y
– массив значений вектора
)
(x
Y
. Их начальные значения можно задавать вручную, или при помощи специальной команды bvpinit. Первым аргументом этой команды является набор

95 отчетов x
Второй аргумент – либо строка констант, которыми заполняется весь массив y
, либо имя функции fun
, возвращающей строку значений в точке x
. Проиллюстрируем обе эти возможности: function y= fun (x) y=[x^2, x-1];
>> solinit = bvpinit([0 1 2 3 4],[1 0]) solinit = x: [0 1 2 3 4] y: [2x5 double]
>> solinit.y ans = 1 1 1 1 1 0 0 0 0 0
>> solinit = bvpinit([0 1 2 3 4],@ fun) solinit = x: [0 1 2 3 4] y: [2x5 double]
>> solinit.y ans = 0 1 4 9 16
-1 0 1 2 3
Теперь сформированы все вспомогательные функции и решатель bvp4c готов к работе.
Последовательность команд
>>solinit = bvpinit([0 1 2 3 4],[1 0]); sol = bvp4c(@twoode,@twobc,solinit); решает двухточечную краевую задачу на интервале [0, 4].
В результате выполнения команды bvp4c возвращается структура sol
, содержащая поля x, y, yp
(вектор производных) и solver
(имя решателя, в данном случае – bvp4c
). График решения можно построить по команде plot(sol.x,sol.y)
, он показан на рис. 5.5 (верхняя кривая).
Для формирования начальной сетки по х мы использовали целочисленные точки [0 1 2 3 4], а начальный вид функции у и ее производной задается вектором [1 0], т.е. перед началом итераций функция y
была равна единице, а ее производная – нулю:
1
)
(

x
y
,
0
)
(


x
y
Изменим знак начальных значений функции на противоположный, приняв
1
)
(


x
y
,
0
)
(


x
y
, и вновь найдем решение:
>>solinit = bvpinit([0 1 2 3 4],[1 0]);
>>sol = bvp4c(@twoode,@twobc,solinit);
Оно будет существенно отличаться от первого (нижняя кривая на рис. 5.5). Этот пример показывает, что разные начальные приближения могут привести к различным решениям.
0 0.5 1
1.5 2
2.5 3
3.5 4
-2
-1 0
1 2
3
x y
A BVP with two solutions
y
0
(x))= -1
y
0
(x) = 1
Рис. 5.5
Заметим, что команда bvp4c при решении использует переменный шаг. В этом можно убедиться, набрав команду plot(sol.x)
. Для получения аппроксимации решения в конкретных точках можно использовать уже упоминавшуюся команду deval. Первый ее аргумент – структура sol
, а второй – вектор точек, в которых нужно вычислить аппроксимацию. Так, по командам xint = linspace(0,4); yint = deval(sol,xint);
вычисляется решение в 100 точках, равномерно расположенных на

96 интервале [0 4]. График первой компоненты решения строится с помощью команды plot(xint,yint(1,:)).
Оба эти решения можно получить аналитически, заменяя нелинейное уравнение
0



y
y
двумя линейными уравнениями











0
при
0 0
при
0
y
y
y
y
y
y
Их общие решения имеют вид
x
C
x
C
x
y
cos sin
)
(
2 1


и
t
t
e
C
e
C
x
y



4 3
)
(
, соответственно. Для получения решения при краевых условиях
2
)
4
(
,
0
)
0
(



y
y
необходимо рассмотреть два случая, в зависимости от знака производной
)
0
(
y

.
Случай 1.
0
)
0
(


y
, тогда
)
(
4 3
x
x
e
C
e
C
x
y



Коэффициенты
4 3
, C
C
находим из краевых условий
0
)
0
(

y
,
2
)
4
(


y
:
2
,
0 4
4 4
3 0
4 0
3







e
C
e
C
e
C
e
C
. Имеем
1 2
,
1 2
8 4
4 8
4 3





e
e
C
e
e
C
Отсюда
4
sh sh
2
)
(
x
x
y


. При любых
0

x
это решение отрицательно. Таким образом, нижняя кривая на рис. 5.5 – перевернутый гиперболический синус.
Случай 2. Пусть
0
)
0
(


y
, тогда
x
C
x
C
x
y
cos sin
)
(
2 1


, и из условия
0
)
0
(

y
имеем
x
C
x
y
C
sin
)
(
,
0 1
2


. Производная этого решения
x
C
x
y
cos
)
(
1


. Решение обращается в ноль при


x
. После этого оно становится отрицательным и принимает вид
x
x
e
C
e
C
x
y



4 3
)
(

. Из условий
2
)
4
(





y
и
0
)
0
(

y
находим константы
3 4
3
,
)
4
(
sh
1
C
C
C






. Производная этой функции имеет вид
)
4
(
sh
)
(








x
x
e
e
x
y
. Из равенства
)
(
)
0
(


y
y



находим
)
4
(
sh
2 1



C
Таким образом, верхняя кривая на рис. 5.5 состоит из двух частей:
4 0
при
)
4
(
sh
)
sh(x
2
)
(
;
0
при
)
4
(
sh sin
2
)
(















x
y
t
x
x
y
Обе части гладко сопрягаются в точке π на оси абсцисс, где у них совпадают не только значения, но и производные.