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

Глава 6. Интегрирование и дифференциальные уравнения

Рис. 6.12. Графическое представление системы уравнений

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

Уравнения высоких порядков

Слушайте, что это Вы со мной загадками разговариваете? Честное слово, мне ничего не ясно!

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

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

251

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

следующую задачу: y′′(x) + 2y′(x) + 2y(x) = 0 с начальными условиями y(0) = 2 и y′(0) = 1 .

На заметку

У этой задачи есть точное решение y(x) = exp(−x)(3 sin(x) + 2 cos(x)) . Его и попытаемся найти.

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

необходимо дополнить начальными условиями. Они в новых обозначениях запишутся так: y1(0) = 2 (следствие условия y(0) = 2 ) и y2(0) = 2 (следствие условия y′(0) = 1 ). На этом математические формальности заканчиваются, и можно приступить непосредственно к поиску решения в Matlab.

Как и в предыдущем случае, начинаем с определения векторной функции, которая определяет систему уравнений. Вот этот код:

function f=LDE(t,y) f=[y(2);-2*(y(1)+y(2))]; end

Здесь все достаточно просто. Результатом функции LDE() возвращается вектор-столбец из двух элементов, каждый из которых задает правую часть соответствующего уравнения. При этом независимая переменная в явном виде не используется. Окно редактора m-файлов с данным программным кодом представлено на рис. 6.13.

Рис. 6.13. Программный код функции, определяющей систему уравнений

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

252

Глава 6. Интегрирование и дифференциальные уравнения

>>[x,y]=ode45('LDE',0:0.01:7,[2;1]);

>>plot(x,y)

>>grid on

>>hold on

>>plot(x,exp(-x).*(3*sin(x)+2*cos(x)),'r:','LineWidth',3)

>>title('Решение уравнения второго порядка')

>>legend('функция y(x)','производная y''(x)','точное решение')

Здесь командой [x,y]=ode45('LDE',0:0.01:7,[2;1]) решается система уравнений, а результат записывается в переменную y. При этом элемент y(1) содержит значение функции y(x) (или, что то же самое, функции y1(x) ), а элемент y(2) содержит значение производной y′(x) (функции y2(x) ). Кроме этой команды, инструкцией plot(x,y) строится график функции y(x) и ее производной y′(x) . Командой plot(x,exp(-x).* (3*sin(x)+2*cos(x)),'r:','LineWidth',3) строится график функции y(x) на основе точного (аналитического) выражения для решения дифференциального уравнения. Используется пунктирная кривая красного цвета толщины 3. Командное окно с соответствующими командами показано на рис. 6.14.

Рис. 6.14. Команды для решения системы уравнений

Графики для полученного числового решения (функция y(x), а также ее производная y′(x) ) и точного решения представлены на рис. 6.15.

Для удобства на графике отображается также и легенда. Это позволяет заметить, что найденное числовое решение практически не отличается (на использованном интервале значений аргумента) от аналитического решения.

На заметку

Как уже отмечалось, функции для решения дифференциальных уравнений, называние которых начинается с аббревиатуры ode (сокращение от ordinary differential equation, что означает обыкновенное дифференциальное уравнение), используются аналогично тому, как это делается для функции ode45().

253

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

Рис. 6.15. Графическое представление числового решения и точного результата

Кроме того, следует помнить, что, помимо описанного способа передачи аргументов встроенным функциям для решения дифференциальных уравнений, могут передаваться и дополнительные параметры – например, допускается в точном виде задавать точность вычислений. Что касается использования той ли иной функции, то обычно выбирают ту, которая позволяет получить приемлемой точности результат за оптимальное время. К сожалению, нередко единственным способом проверить эффективность разных функций является метод "проб и ошибок".

Снова об интегралах

Я просто не предполагал, что моя догадка настолько попадет в цель.

К/ф "Приключения Шерлока Холмса и доктора Ватсона. Кровавая надпись"

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

Для решения поставленной задачи можно использовать разные подходы. Здесь мы создадим специальную функцию, в которой при вычислении

254

Глава 6. Интегрирование и дифференциальные уравнения

интеграла вызывается одна из встроенных функций для решения дифференциальных уравнений (а именно, уже знакомая нам функция ode45()). Окно редактора m-файлов с кодом этой функции показано на рис. 6.16.

Рис. 6.16. Код функции для вычисления интегралов

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

function res=Int(F,a,b) f=@(x,y)F(x); [x,y]=ode45(f,a:(b-a)/100:b,0); res=y(length(y)); plot(x,F(x),'b-','LineWidth',2); hold on; plot(x,y,'r:','LineWidth',3); plot(b,res,'rs','LineWidth',4); plot([a,b],[res,res],'r-'); plot([b,b],[0,res],'r-');

grid on;

title('Вычисление интеграла','FontWeight','Bold','FontSize',12); legend('функция','первообразная','результат',2); text(a,res,num2str(res),'VerticalAlignment','bottom',...

'EdgeColor','red','BackgroundColor','yellow','FontWeight','BOLD'); hold off;

end

255