Материал: Matlab. Практический подход. Самоучитель

Самоучитель Matlab

ответствующе решения на графике (а мы это увидим далее), то получим осциллирующие с затуханием функции.

Приведенное решение в числовом виде можно получить достаточно просто без использования специальных функций для решения систем дифферен-

и матрица коэффициентов правых частей системы дифференциальных

Затем воспользуемся формальной аналогией. Если бы у нас было не вектор-

ное уравнение, а скалярное df (t) = Af (t) с начальным условием f (0) = f0 , dt

то решение можно было бы записать практически сразу: f (t) = f0 exp(At) . Но у нас уравнение векторное, поэтому все операции в соответствующем выражении для решения следует понимать в смысле операций с матрицами.

решения в фиксированный момент времени t необходимо матрицу ˆ умно-

A

жить на скаляр t (выполняется поэлементное умножение). В результате по-

лучаем матрицу ˆ . Затем вычисляем матрицу ˆ (по правилу вычис-

At exp(At)

ления функции expm()) и умножаем ее (по правилу умножения матриц)

на вектор начальных значений f0 .

На заметку

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

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

function [x y]=ldes(A,t,z) N=length(t);

166

Глава 4. Элементы матричной алгебры

for i=1:N M=expm(A*t(i));

x(i)=z(1)*M(1,1)+z(2)*M(1,2);

y(i)=z(1)*M(2,1)+z(2)*M(2,2); end

end

Аргументов у функции три: матрица коэффициентов A, список значений аргумента t (узловые точки), в которых вычисляются значения искомых функций, а также вектор начальных значений функций z. Результат вычисления функции записывается в массивы x и y, каждый из которых представляет набор точек-значений соответственно первой x(t) и второй y(t) искомых функций в узловых точках.

Командой N=length(t) определяется количество узловых точек, для которых необходимо вычислить значения функций. Затем запускается оператор цикла, в котором индексная переменная i пробегает значения от 1 до N. На каждой итерации определяются значения искомых функций в соответствующей узловой точке. Так, при фиксированном значении индексной переменной командой M=expm(A*t(i)) вычисляется матричная экспонента (для данного момента времени t(i)). Затем командами x(i)=z(1)*M(1,1)+z(2)*M(1,2) и y(i)=z(1)*M(2,1)+z(2)*M(2,2)

вычисляются значения функций в узловой точке.

На заметку

При вычислении произведения матрицы M на вектор-столбец начальных условий z получаем вектор-столбец. Первый элемент этого вектора-столбца дает значение функции x(t) в узловой точке (элемент x(i)). Второй элемент вектора-столбца дает значение функции y(t) в узловой точке (элемент y(i)). Чтобы получить первый элемент вектора-столбца, необходимо перемножить и сложить соответствующие элементы первой строки матрицы M и элементы вектора начальных значений z. Первая строка матрицы M – это элементы M(1,1) и M(1,2). Первый элемент вектора начальных значений – это z(1), а второй элемент вектора начальных значений – это z(2). В результате получаем комбинацию z(1)*M(1,1)+z(2)*M(1,2). Аналогично вычисляется второй элемент вектора-столбца: перемножаются и суммируются элементы второй строки матрицы M и элементы вектора начальных значений. Желающие могут подумать, как тот же алгоритм можно было бы реализовать с использованием операции матричного произведения вместо явного выписывания сумм и произведений.

Окно редактора m-файлов с кодом функции ldes() показано на рис. 4.14.

После того, как функция ldes() создана, ее можно использовать для решения системы дифференциальных уравнений. В командном окне вводим следующие команды (жирным шрифтом выделен ввод пользователя):

167

Самоучитель Matlab

Рис. 4.14. Редактор m-файлов с кодом функции для вычисления решения дифференциального уравнения

>> A=[-2 3;-3 1]

A =

-2 3 -3 1

>>t=0:0.01:3*pi;

>>[x y]=ldes(A,t,[1 3]);

>>plot(t,x,t,y)

Здесь командой A=[-2 3;-3 1] создается матрица коэффициентов, командой t=0:0.01:3*pi формируется набор узловых точек, а затем командой [x y]=ldes(A,t,[1 3]) находим решение системы дифференциальных уравнений. Командное окно с кодом приведено на рис. 4.15.

Рис. 4.15. Командное окно с кодом для вычисления решения дифференциального уравнения

В результате выполнения команды plot(t,x,t,y) создается график для функциональных зависимостей x(t) и y(t), полученных в результате решения системы дифференциальных уравнений. Графики функций x(t) и y(t) представлены на рис. 4.16.

168

Глава 4. Элементы матричной алгебры

Рис. 4.16. Графическое представление для найденного решения

На заметку

Графики функций строятся на интервале от 0 до 3π. Наличие числа π здесь является чистой формальностью и никакой особой смысловой нагрузки не несет.

Преобразование матриц

Этого объяснить я Вам не могу, потому что сам толком ни черта не понимаю.

К/ф "Семнадцать мгновений весны"

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

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

169

Самоучитель Matlab

Рассмотрим следующий командный код более внимательно (жирным шрифтом выделен ввод пользователя):

>> E=eye(2)

E =

10

01

>>O=zeros(2)

O =

00

00

>>U=ones(2)

U =

11

ты нулевые). Функцией zeros() создается матрица с нулевыми элементами, а функцией ones() создается матрица, все элементы которой единичные. Размер матрицы указывается аргументом перечисленных функций. Таким образом, создается три матрицы размерами 2×2 каждая: единичная матрица E, нулевая матрица O и матрица из единиц U. Затем командой A=[E,2*U;3*U,O] создается блочная матрица A. Если посмотреть на эту команду формально, то речь идет о создании матрицы 2×2. Однако элементы матрицы сами являются матрицами. Причем матрица U при передаче элементом в блочную матрицу предварительно умножается на число. В результате получаем матрицу A, которая имеет размеры 4×4. Ее левый верхний блок – это матрица E. Верхний правый блок – матрица, состоящая из двоек (результат умножения матрицы U на 2). Левый нижний блок – матрица, состоящая из троек (результат умножения матрицы U на 3). Правый нижний блок – нулевая матрица O.

170