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_M  df_M       df_M |
             | ----  ----  ...  ---- |
             | dx_1  dx_2       dx_N |

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

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

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

Если вычисление не выполнено успешно, 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).

"jacobian type"

Строка, определяющая тип якобиана, используемого с методом интегрирования жёсткого обратного дифференцирования (BDF). Допустимые значения

"full"

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

"banded"

Приближаются или используются только диагональ и количество нижних и верхних поддиагоналей, указанных параметрами "lower jacobian subdiagonals" и "upper jacobian subdiagonals" соответственно, из функции пользовательского якобиана. Пользовательская функция якобиана может устанавливать все остальные частные производные в произвольные значения.

"diagonal"

Если функция якобиана предоставлена пользователем, данное значение не оказывает никакого влияния. Приближенный методом конечных разностей якобиан lsode ограничен диагональю, где каждая частная производная вычисляется путём применения конечного изменения ко всем элементам состояния вместе; если настоящий якобиан действительно всегда диагонален, это имеет тот же эффект, что и применение конечного изменения только к соответствующему элементу состояния, но более эффективно.

"lower jacobian subdiagonals"

Количество нижних поддиагоналей, используемых, если параметр "jacobian type" установлен на "banded". Значение по умолчанию — ноль.

"upper jacobian subdiagonals"

Количество верхних поддиагоналей, используемых, если параметр "jacobian type" установлен на "banded". Значение по умолчанию — ноль.

Вот пример решения системы из трёх дифференциальных уравнений с помощью 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)];

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