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

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

dx
-- = f (x, t)
dt

с помощью решателя ОДУ Hindmarsh 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 можно использовать функцию lsode_options.

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

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

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.