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₀ и должен соответствовать начальному состоянию системы x₀, так что первая строка результата — x₀.

Первый аргумент 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₀. Третий аргумент — вектор t, задающий значения времени, для которых требуется решение.

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

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

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

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

См. Алана К. Хиндмарша, ODEPACK, Систематизированный набор решателей ОДУ, в «Вычислительной науке», Р. С. Стэплман, редактор, (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.