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 (соответствует Fortran-версии LSODE).

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

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

См. Алана К. Хиндмарша, ODEPACK, Систематизированная коллекция решателей ОДУ, в Scientific Computing, Р. С. Степлеман, редактор, (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.