Как решать дифференциальные уравнения в Matlab

Matlab позволяет решать обыкновенные дифференциальные уравнения различного порядка с начальными условиями (задачи Коши).
Решатели Matlab реализуют различные методы решения систем дифференциальных уравнений (см. таблицу 1).

Таблица 1.

Для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode15s, ode23s, ode23t, ode23tb.

Все решатели могут решать системы уравнений явного вида у' = F(t, y). Решатели ode15s и ode23t способны найти корни дифференциально-алгебраических уравнений M(t)y' = F(t, у), где М называется матрицей массы. Решатели ode15s, ode23s, ode23t и ode23tb могут решать уравнения неявного вида M(t,y) у' = F(t, у).
Особенности решателей следующие:
- ode23tb, ode23s служат для решения жестких дифференциальных уравнений;
· ode15s - жестких дифференциальных и дифференциально-алгебраических уравнений;
· ode23t - умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

Последовательность решения обыкновенных дифференциальных уравнений следующая:
1. Приведение дифференциального уравнения к системе дифференциальных уравнений первого порядка. Для этого вводится столько дополнительных функций, каков порядок уравнения.
2. Создание специальной файл-функции для системы уравнений. Файл-функция содержит два входных аргумента: переменную t, по которой производится дифференцирование, даже если она входит в уравнение неявно, и вектор, размер которого равен числу неизвестных функций системы.
3. Вызов нужного решателя (встроенной функции). Входными аргументами решателя являются имя файл-функции, вектор с начальным и конечным значениями переменной

Пример. Решить дифференциальное уравнение вида

при начальных условиях x0=0, y(x0)=1 с шагом интегрирования h=0,1 на интервале [0,1]. Заметим, что аналитическое решение этого уравнения запишется так:

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

· options - аргумент, создаваемый функцией odeset - позволяет вывести параметры, установленные по умолчанию;

· tspan - вектор, определяющий интервал интегрирования [t0 tfinal]. Для получения решений в конкретные моменты времени t0, tl,..., tfinal (расположенные в порядке уменьшения или увеличения) нужно использовать tspan = [t0 tl ... tfinal];

· у0 - вектор начальных условий;

· Т, Y - матрица решений Y, где каждая строка соответствует времени, возвращенном в векторе-столбце Т.

Решение:

1. Создадим М-файл с функцией для вычисления правой части дифференциального уравнения, для этого выполним команду главного окна Matlab File > New> M-file.

2. В окне редактирования кода запишем код (рис. 1):

function dy=difur(x,y)
dy=2*x^2+2*y;

Рис. 1.

3. Сохраним М-файл с именем difur.

4. В командном окне введем инструкции (рис.2.).

>> % Интервал интегрирования
>> tspan = [0, 1];
>> y0 = [0; 1]; % начальные условия
>> [x,y] = ode45(@difur, tspan, y0, []);
>> % График
>> plot(x,y(:,1))

Рис. 2.

В результате выполнения инструкций будет построен график (рис. 3.)

Рис. 3.

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

Заметим, что дифференциальное уравнение можно решить в MatLab и символьно. Для этого нужно применить стандартную процедуру dsolve.
Запишем в командном окне инструкцию (рис. 4)
>> dsolve('Dy=-2*x ^2-2*y','y(0)=1')

В итоге, после ее выполнения получим аналитическую запись решения дифференциального уравнения:
ans =-x^2+exp(-2*t)*(x^2+1)

Рис. 4.