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 (согласуется с фортранской версией LSODE).

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

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

См. Alan C. Hindmarsh, ODEPACK, Систематизированное собрание решателей ОДУ, в Scientific Computing, R. S. Stepleman, редактор, (1983) или https://computing.llnl.gov/projects/odepack для получения дополнительной информации о внутренней работе lsode.

Пример: Решение уравнения Ван дер Поля

fvdp = @(y,t) [y(2); (1 - y(1)^2) * y(2) - y(1)];
t = linspace (0, 20, 100);
y = lsode (fvdp, [2; 0], t);

См. также: 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)];

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