24.1 Обыкновенные дифференциальные уравнения

Функцию lsode можно использовать для решения ОДУ вида

dx
-- = f (x, t)
dt

с помощью решателя ОДУ Хиндмарша LSODE.

[x, istate, msg] = lsode (fcn, x_0, t)

[x, istate, msg] = lsode (fcn, x_0, t, t_crit)

Решатель обыкновенных дифференциальных уравнений (ОДУ).

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

dx
-- = f (x, t)
dt

с

x(t_0) = x_0

Решение возвращается в матрице x, где каждая строка соответствует элементу вектора t. Первый элемент t должен быть t_0 и соответствовать начальному состоянию системы x_0, так что первая строка выходных данных — это x_0.

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

xdot = f (x, t)

в котором xdot и x — векторы, а t — скаляр.

Если fcn — это массив строк из двух элементов или массив ячеек из двух элементов строк, встроенных функций или обработчиков функций, первый элемент задаёт имя функции f, описанной выше, а второй элемент задаёт имя функции для вычисления якобиана f. Функция якобиана должна иметь вид

jac = j (x, t)

в котором jac — матрица частных производных

| df_1  df_1       df_1 |
             | ----  ----  ...  ---- |
             | dx_1  dx_2       dx_N |
             |                       |
             | df_2  df_2       df_2 |
             | ----  ----  ...  ---- |
      df_i   | dx_1  dx_2       dx_N |
jac = ---- = |                       |
      dx_j   |  .    .     .    .    |
             |  .    .      .   .    |
             |  .    .       .  .    |
             |                       |
             | df_N  df_N       df_N |
             | ----  ----  ...  ---- |
             | dx_1  dx_2       dx_N |

Второй аргумент задаёт начальное состояние системы x_0. Третий аргумент — вектор t, задающий значения времени, для которых требуется решение.

Четвёртый аргумент является необязательным и может быть использован для задания набора моментов времени, за которые решатель ОДУ не должен интегрировать. Это полезно для избежания трудностей с особенностями и точками, где существует разрыв в производной.

После успешного вычисления значение istate будет равно 2 (соответствует Fortran-версии LSODE).

Если вычисление не выполняется успешно, istate будет отличным от 2, а msg будет содержать дополнительную информацию.

Вы можете использовать функцию lsode_options для задания необязательных параметров для lsode.

См. также: daspk, dassl, dasrt.

lsode_options ()

val = lsode_options (opt)

lsode_options (opt, val)

Запрос или установка параметров для функции lsode.

При вызове без аргументов отображаются имена всех доступных параметров и их текущие значения.

При передаче одного аргумента возвращается значение параметра opt.

При вызове с двумя аргументами lsode_options устанавливает параметр opt в значение val.

Параметры включают

"absolute tolerance"

Абсолюльная погрешность. Может быть либо вектором, либо скаляром. Если это вектор, он должен соответствовать размерности вектора состояния.

"relative tolerance"

Параметр относительной погрешности. В отличие от абсолютной погрешности, этот параметр может быть только скаляром.

Тест локальной погрешности, применяемый на каждом шаге интегрирования, имеет вид

abs (local error in x(i)) <= ...
      rtol * abs (y(i)) + atol(i)

"integration method"

Строка, определяющая метод интегрирования для решения системы ОДУ. Допустимые значения

"adams"

"non-stiff"

Якобиан не используется (даже если он доступен).

"bdf"

"stiff"

Использование метода обратной дифференциации по времени (BDF) для жёстких систем. Если функция для вычисления якобиана не указана, lsode вычислит приближение матрицы якобиана методом конечных разностей.

"initial step size"

Шаг, который будет предпринят на первом шаге (по умолчанию определяется автоматически).

"maximum order"

Ограничение максимального порядка метода решения. При использовании метода Адамса этот параметр должен быть от 1 до 12. В противном случае — от 1 до 5 включительно.

"maximum step size"

Установка максимального шага позволит избежать прохождения очень больших областей (по умолчанию не указано).

"minimum step size"

Минимальный абсолютный размер шага (по умолчанию 0).

"step limit"

Максимальное количество шагов (по умолчанию 100000).

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

## oregonator differential equation
function xdot = f (x, t)

  xdot = zeros (3,1);

  xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ...
            - 8.375e-06*x(1)^2);
  xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
  xdot(3) = 0.161*(x(1) - x(3));

endfunction

и начальное условие x0 = [ 4; 1.1; 4 ], система уравнений может быть интегрирована с помощью команды

t = linspace (0, 500, 1000);

y = lsode ("f", x0, t);

Если вы попробуете это, вы увидите, что значение результата резко меняется между t = 0 и 5, а также примерно при t = 305. Более эффективным набором выходных точек может быть

t = [0, logspace(-1, log10(303), 150), ...
        logspace(log10(304), log10(500), 150)];

См. работу Алана Хиндмарша, ODEPACK, Систематизированный набор решателей ОДУ, в «Научном моделировании», Р. С. Степлеман, редактор, (1983), для получения дополнительной информации о внутренней работе lsode.

М-файл для дифференциального уравнения, используемого выше, включён в дистрибутив Octave в каталоге примеров под именем oregonator.m.