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

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

На заметку

где константа C определяется исходя из начальных условий. В частности, для

начального значения y(1) = 1 получаем решение y(x) = 1 .

x

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

>>f=@(x,y)(y.^2-2./x.^2);

>>[x,y]=ode45(f,1:0.1:10,1);

>>plot(x,y,'r-','LineWidth',2)

>>hold on

>>plot(x,1./x,'b:','LineWidth',2)

>>grid on

На рис. 6.6 приведен документ с командным кодом, представленным выше.

Рис. 6.6. Решение дифференциального уравнения с помощью функции ode45()

функции ode45() в команде [x,y]=ode45(f,1:0.1:10,1). Вторым аргументом функции ode45() передается список значений аргумента x , для которых вычисляются значения функции y . Третий аргумент функции ode45() определяет начальное значение функции.

На заметку

Точка начального значения по аргументу определяется по первому элементу списка – второго аргумента функции ode45().

В качестве результата функцией возвращается список узловых точек и вычисленных значений функции в этих точках. Таким образом, в результате выполнения команды [x,y]=ode45(f,1:0.1:10,1) в переменную x за-

246

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

писываются узловые точки (определяются вторым аргументом функции ode45()), а в переменную y записываются значения функции, вычисленные в узловых точках. Эти переменные можно использовать для создания графика, что собственно и делается с помощью команды plot(x,y,'r-', 'LineWidth',2). Кривая для вычисленной в результате решения дифференциального уравнения зависимости отображается сплошной линией красного цвета толщины 2. Для сравнения командой plot(x,1./x, 'b:','LineWidth',2) пунктирной кривой синего цвета отображается точная аналитическая зависимость для решения рассмотренной задачи. На рис. 6.7 показаны графики, построенные на основе точного решения дифференциального уравнения и вычисленного с помощью функции ode45().

Рис. 6.7. Результат решения дифференциального уравнения в числовом виде

Как видим, на интервале значений аргумента от 1 до 10 совпадение более чем приемлемое. Правда, в данном случае зависимость монотонная. Если изменить начальное значение для функции, можно получить не такое тривиальное решение. На рис. 6.8 представлен документ (продолжение предыдущего), в котором решается то же самое уравнение, но с начальным условием y(1) = −0.5 .

Все практически так же, как в предыдущем случае, только внесены соответствующие изменения в начальном условии и для сравнения создается

247

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

Рис. 6.8. Решение уравнения с новым начальным условием

3 2

график функции y(x) = (1 + x3)x − x , которая является решением соответствующей задачи. График приведен на рис. 6.9.

Рис. 6.9. Еще одно решение дифференциального уравнения

Как и в предыдущем случае, совпадение результатов хорошее.

На заметку

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

248

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

Решение системы дифференциальных уравнений

Ну что Вы теперь на это скажете, мой дорогой психолог?

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

Системы дифференциальных уравнений решаются практически так же, как и отдельные дифференциальные уравнения, с той лишь поправкой, что придется иметь дело с векторными функциями (списком функций). В качестве примера рассмотрим систему дифференциальных уравнений x′(t) = y(t) + 2 exp(t) и y′(t) = x(t) +t2 с начальными значениями

x(0) = −2 и y(0) = −1. Сразу отметим, что точное решение имеет вид x(t) = t exp(t) −t2 −2 и y(t) = (t −1)exp(t) −2t . Для решения этой системы предварительно определяем в редакторе m-файлов следующую векторную функцию:

function f=DESyst(t,y) f=[y(2)+2*exp(t);y(1)+t.^2]; end

код функции в окне редактора m-файлов представлен на рис. 6.10.

Рис. 6.10. Код векторной функции

Что касается определения функции DESyst(), то у нее два аргумента – независимая переменная t (значение переменной t ) и переменная y, которая, как предполагается, содержит два значения – для функции x(t) и для функции y(t) (в момент времени t ). Ссылка на значение функции x(t) выглядит как y(1), а ссылка на значение функции y(t) выглядит как y(2). Именно эти ссылки использованы при вычислении результата функции DESyst(). Возвращаемое значение является списком (вектор-столбец). Значения

249

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

элементов вычисляются в соответствии с уравнениями (правыми частями уравнений), формирующими систему.

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

>>[t,y]=ode45('DESyst',0:0.01:2,[-2;-1]);

>>plot(t,y,'LineWidth',2)

>>grid on

>>hold on

>>plot(t,t.*exp(t)-t.^2-2,'r:')

>>plot(t,(t-1).*exp(t)-2*t,'k--')

>>title('Решение системы дифференциальных уравнений')

Интерес здесь представляет команда [t,y]=ode45('DESy st',0:0.01:2,[-2;-1]). Она по форме мало отличается от тех, что использовались ранее для решения отдельного дифференциального уравнения. Но поскольку первый аргумент (ссылка 'DESyst') передает векторфункцию, то в переменную y также записывается несколько функций: первый столбец этой переменной определяет зависимость x(t), а второй столбец определяет зависимость y(t). Соответствующие узловые точки (значения независимой переменной t ) записываются в переменную t. Прочие команды предназначены для иллюстрации полученного решения в графическом виде. На рис. 6.11 показано командное окно с использованными инструкциями.

Рис. 6.11. Решение системы дифференциальных уравнений

На заметку

Командой plot(t,y,'LineWidth',2) отображаются сразу две кривые, для функций x(t) и y(t). Соответствующие данные "спрятаны" в переменной y. Кроме этого, командами plot(t,t.*exp(t)-t.^2-2,'r:') и plot(t,(t-1).*exp(t)-2*t,'k--') строятся кривые для точных решений. На рис. 6.12 числовое и аналитическое решения неразличимы.

250