Численное решение ОДУ и их систем

Для решения обыкновенных дифференциальных уравнений (ОДУ или ODE), а также систем ОДУ в MATLAB реализованы различные методы (аналитические или численные). Реализации численных методов названы решателямиОДУ.

SOLVER (решатель) означает один из возможных численных методов решения ОДУ: ode45, ode23, ode113, ode23tb, ode15s, ode23s, ode23t , bvp4c или pdepe. Решатели реализуют следующие методы решения систем дифференциальных уравнений, причем для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode 15s, ode23s, ode23t. ode23tb:

· ode45 – одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка. Это классический метод, рекомендуемый для начальной пробы решения. Во многих случаях он дает хорошие результаты;

· ode23 – одношаговые явные методы Рунге-Кутта 2-го и 4-го порядка. При умеренной жесткости системы ОДУ и низких требованиях к точности этот метод может дать выигрыш в скорости решения;

· ode113 – многошаговый метод Адамса-Башворта-Мултона переменного порядка. Это адаптивный метод, который может обеспечить высокую точность решения;

· ode23tb – неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем;

· ode15s – многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения;

· ode23s – одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности решения жесткой системы дифференциальных уравнений;

· ode23t – метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих колебательные системы с почти гармоническим выходным сигналом;

· bvp4c – служит для решения проблемы граничных значений систем дифференциальных уравнений вида , (краевая задача);

· pdepe – нужен для решения систем параболических и эллиптических дифференциальных уравнений в частных производных, введен в ядро системы для поддержки новых графических функций Open GL, пакет расширения Partial Differential Equations Toolbox содержит более мощные средства.

Все решатели могут решать системы уравнений явного вида . Решатели ode15s и ode23t способны найти корни дифференциально-алгебраических уравнений , где М называется матрицей массы. Решатели ode15s, ode23s, ode23t и ode23tb могут решать уравнения неявного вида . Решатели ode23tb, ode23s служат для решения жестких дифференциальных уравнений, ode15s – для жестких дифференциальных и дифференциально-алгебраических уравнений, ode23t – для умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

Для решения ОДУ необходимо в окне команд MATLAB набрать одну из следующих команд, причем вместо solver необходимо подставить имя конкретного решателя (например ode45):

[X,Y] = solver(@F, xspan, у0)

[X,Y] = solver(@F, xspan, y0, options)

[X,Y] = solver(@F, xspan, y0, options, pl,p2...)

[X,Y,XE,YE,IE] = solver(@F, xspan, y0, options)

[X,Y,S] = solver(…)

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

· @F – дескриптор OДУ (ссылка на имя M-файла, совпадающая с именем функции F – дифференциального уравнения);

· xspan – вектор, определяющий интервал интегрирования [x0 xfinal]. Для получения решений в конкретные моменты времени x0, xl, ... , xfinal (расположенные в порядке уменьшения или увеличения) нужно использовать xspan = [x0 xl ... xfinal] или с шагом h [x0: h: xfinal];

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

· options – свойства (опции), создаваемые функцией odeset (еще одна функция odeget или bvpget (только для bvp4c) позволяет вывести параметры, установленные по умолчанию или с помощью функции odeset /bvpset);

· p1, р2,…, – произвольные параметры, передаваемые в функцию F;

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

· [X,Y] = solver(@F, xspan, у0) – интегрирует систему дифференциальных уравнений вида на интервале xspan с начальными условиями у0. @F – дескриптор OДУ-функции. Каждая строка в массиве решений Y соответствует значению аргумента, возвращаемому в векторе-столбце X;

· [X,Y] = solver(@F, xspan, y0, options) – дает решение, подобное описанному выше, но с параметрами, определяемыми значениями аргумента options, созданного функцией odeset. Обычно используемые параметры включают допустимое значение относительной погрешности RelTol (по умолчанию le-3) и вектор допустимых значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны 1е-6);

· [X,Y] = solver(@F, xspan, y0, options, pl,p2...) – дает решение, подобное описанному выше, передавая дополнительные параметры p1, р2,... в m-файл F всякий раз, когда он вызывается (вычисление функций правых частей). Используйте options=[], если никакие параметры не задаются (т.е. опции используются по умолчанию);

· [X,Y,XE,YE,IE] = solver(@F, xspan, y0, options) – при условии, что опция Events установлена в структуре options и позволяет зафиксировать момент, когда какая-либо из указанных переменных принимает значение 0; XE – вектор столбец моментов времени, когда фиксировалось событие, YE – перечень соответствующих переменных, IE – индексы событий, указывающие, в каком направлении изменялось решение в момент пересечения нуля (-1 убывало, +1 возрастало, 0 только фиксировалось);

· [X,Y,S] = solver(…) – позволяет получить статистические характеристики решателя; при этом вектор-столбец S имеет следующие компоненты (количество успешных шагов, количество неудачных попыток, количество вызовов функции F(X,Y), количество вызовов функции Якоби dF/dx, количество LU-разложений, количество решений вспомогательных систем уравнений);

Команда ODESET – задает, изменяет или выводит структуру свойств (опций) для решателя SOLVER обыкновенного дифференциального уравнения или их систем:

options = odeset('name1',value1,'name2',value2,...)

options = odeset(oldopts,'name1',value1,...)

options = odeset(oldopts,newopts)

odeset

options = odeset('name1',value1,'name2',value2,...) – формирует массив записей (структуру) options, в котором свойствам (опциям) с именами 'name1', 'name2', … приписываются некоторые значения value1, value2,…(достаточно указать только первый символ); неприсвоенному значению приписывается пустой массив [];

options = odeset(oldopts,'name1',value1,...) – заменяет значения отдельных свойств (опций);

options = odeset(oldopts,newopts) – заменяется набор старых свойств набором новых, но если при этом для какого-либо из новых свойств указано [], то сохраняется прежнее значение этого свойства;

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

Границы погрешностей интегрирования

Выходные параметры решателя

Характеристики якобиана (для функций ode15s и ode23)

Фиксация события

Характеристики матрицы при производных (для функций ode15s и ode23)

Шаг интегрирования

Параметры функции ode15s

Свойство	Значение	Описание
MaxOrder	1/2/3/4/{5}	Максимальный порядок формулы интегрирования
BDF	on | {off}	Определяет, должны ли использоваться формулы обратного дифференцирования вместо используемых по умолчанию формул прямого дифференцирования

Например, команда options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); – устанавливает свойству (опции) 'RelTol' (допустимая относительная погрешность для всех (так как скаляр) компонентов вектора решения) значение 1e-4, а свойству 'AbsTol' (абсолютная погрешность для соответствующих компонентов вектора) значения 1e-4, 1e-4, 1e-5.

Команда PLOT предназначена для построения двумерных графиков в линейном маштабе:

plot(Y)

plot(X1,Y1,...)

plot(X1,Y1,LineSpec1,...)

plot(...,'PropertyName',PropertyValue,...)

h = plot(...)

plot(Y) – строит график элементов одномерного массива Y в зависимости от номера элемента, если Y – двумерный массив, то строится график для столбцов;

plot(X1,Y1,...) – соответствует построению графика функции, когда одномерный массив Х1 соответствует значениям аргумента, а одномерный массив Y1 – значениям функции; если Х1 – одномерный массив, а Y1 двумерный, то строятся графики для столбцов массива Y1 в зависимости от значений Х1; если Х1 двумерный, Y1 одномерный, то строятся графики для столбцов Х1 в зависимости от Y1; если оба двумерные массивы, то строятся зависимости значений столбцов Y1 от значений столбцов Х1;

plot(X1,Y1,LineSpec1,...) – позволяет построить график функции, указав способ отображения линий, способ отображения маркера точек, цвет линий и маркера с помощью строковой переменной LineSpec, которая может включать до трех символов из следующей таблицы:

Если Х1 – одномерный массив, а Y1 – двухмерный, то для отдельной установки свойств каждому графику (столбцу) Y1 необходимо указать номер столбца Y1, (например для plot(Х1,Y1(:,1),’-’,Х1,Y1(:,2),’ –.’) первому столбцу Y1 назначается непрерывный тип линии, а второму столбцу – штрихпунктирный;

plot(X1,Y1,LineSpec1,X1,Y1,LineSpec1,...) – позволяет построить на одном графике несколько графических объектов Line для функций y1(x1), y2(x2),…, определив для каждой из них свой способ отображения.

plot(...,'PropertyName',PropertyValue,...) – позволяет задать значения свойств графического объекта Line, соответствующего построенному графику.

h = plot(...) – возвращает вектор дескрипторов для всех графических объектов Line текущего объекта.

Кроме команды PLOT могут быть использованы другие команды элементарной графики: FPLOT (построение графиков функций), LOGLOG (график в логарифмическом масштабе), SEMILOGX, SEMILOGY (график в полулогарифмическом масштабе), POLAR (график в полярных координатах), PIOTYY (график с двумя осями ординат), PLOT3 (построение линий и точек в трехмерном пространстве), MESHGRID (формирование прямоугольной сетки), MESH, MESHC, MESHZ (трехмерная сетчатая поверхность), SURF, SURFC (трехмерная сплошная поверхность), SHADING (затенение поверхности), SURFL (сплошная поверхность с подсветкой).

Покажем применение решателя ОДУ на примере решения задачи Коши

y¢=x² y³ sin³ (x+y), y(0)=1 на отрезке [0,3] методом Рунге-Кутты с постоянным шагом h=0,3.

Перед решением запишем дифференциальное уравнение в виде ode-функции. Для этого в главном меню MATLAB выберем File àNew à M-File (Файл àНовый à M-Файл) и введем

Сохраним m-файл-функцию в необходимом каталоге (он при помощи списка на стандартной панели инструментов MATLAB должен быть выбран как «Текущий каталог») под именем, совпадающим с именем функции (LR1.M). Так как необходимо решить одно ОДУ первого порядка, то матрица производных dydx вырождается в вектор столбец из одного элемента dydx=zeros(1,1), для которого и задается вид дифференциального уравнения dydx(1)=x^2*y(1)^3*(sin(x+y(1)))^3.

Тогда в окне команд для получения решения (решателем ode45) необходимо выполнить следующую команду:

»[X,Y]=ode45(@lrl,[0:0.3:3],[1]);

где @lr1 – значение дескриптора (@F) OДУ-функции; [0:0.3:3] – есть значения вектора (xspan), определяющего интервал интегрирования от 0 с шагом 0.3 до 3; [1] – значения вектора (у0) начальных условий y(0)=1.

Для получения сопровождающего решение графика следует выполнить команду PLOT, например в виде:

»plot(X,Y);

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

Тогда для решения задачи Коши для системы из трех ОДУ

y₁¢ = -y₂₊ sin (x y₃), y₁(0)=1

y₂¢ = y₁² y₂(0)=0

y₃¢ = -y₃ - y₁ y₃(0)=1

на отрезке [0,3] методом Рунге-Кутты с постоянным шагом h=0,3 получаем М-файл в виде:

В поле команд MATLAB набираем

»[X,Y]=ode45(@lr2,[0:0.3:3],[1 0 1]);

»plot(X,Y(:,1),'-o',X,Y(:,2),'--ms',X,Y(:,3),'--d');

Получаем следующий график решений, где например, для первого столбца Y назначается непрерывная линия с типом маркера кружок (Y(:,1), '-o'):