Задача об эпидемии – решение системы дифференциальных уравнений
В пункте 5 рис. 5.2 столбцы матрицы Z (время, число больных и число здоровых) разнесены по отдельным векторам и отображены графически, что позволяет проследить динамику развития эпидемии.
На рис. 5.2 получены несколько иные результаты, чем на рис. 5.1, хотя характер кривых сохранился: максимум больных (3119) наблюдается не на 9-й, а на 7-й день, на 13-й день мы имеем не 105, а 209 больных. Это объясняется и различными значениями точности расчетов (на рис. 5.1 делалось 13 шагов интегрирования, а на рис 5.2 ¾ 500) и различными примененными методиками (Эйлер против Рунге и Кутта[5]).
В функцию rkfixed заложен широко распространенный метод решения дифференциальных уравнений – метод Рунге ¾ Кутта[6]. Несмотря на то что это не самый быстрый метод, функция rkfixed почти всегда справляется с поставленной задачей. Однако есть случаи, когда лучше использовать более сложные методы. Эти случаи попадают под три широкие категории: система может быть жесткой[7]
(Stiffb, Stiffr), функции системы могут быть гладкими
(Bulstoer) или плавными (Rkadapt). Нередко приходится пробовать на одном дифференциальном уравнении (одной системе) несколько методов, чтобы определить, какой метод лучше (быстрее, точнее). Примерно так мы сравнивали в этюде 3 разные способы поиска оптимумов функции.
Когда известно, что решение гладкое, используется функция Bulstoer, куда заложен метод Булирша ¾ Штёра, а не Рунге ¾ Кутты, используемый функцией rkfixed. В этом случае решение будет точнее. Список аргументов и матрица, получаемая при работе с функцией Bulstoer, такие же, как и при работе с rkfixed.
Можно решить задачу более точно (более быстро), если уменьшать шаг (у нас это Dt) там, где производная меняется быстро, и увеличивать шаг там, где она ведет себя более спокойно. Для этого предусмотрена функция Rkadapt (adaption –адаптация). Но, несмотря на то что при решении дифференциального уравнения функция Rkadapt использует непостоянный шаг, она тем не менее представит ответ для точек, находящихся на одинаковом расстоянии, заданном пользователем. Аргументы и матрица, возвращаемая функцией Rkadapt, такие же, как при rkfixed.
При решении жестких систем следует использовать одну из двух встроенных функций, разработанных специально для таких случаев: Stiffb и Stiffr.
Они используют метод Булирша ¾ Штёра (b) или Розенброка (r). Форма матрицы-решения, полученной с помощью этих функций, идентична матрице, полученной через rkfixed. Однако Stiffb и Stiffr требуют дополнительного аргумента J:
Stiffb(x, tнач, tкон, n, f, J)
Stiffr(x, tнач, tкон, n, f, J),
где
x – вектор n начальных значений;
tнач, tкон – конечные точки интервала, в котором должно быть найдено решение дифференциальных уравнений. Начальные значения x определяются в точке tнач;
n – количество точек за начальной точкой, в которых должно быть определено решение. Это определяет число рядов (1 + n) матрицы, которую генерируют функции Stiffb и Stiffr;
f(t, x) – функция-вектор правых частей системы;
J(t, x) – матрица-функция размерности n×(n+1), в которой содержится матрица Якоби правых частей дифференциальных уравнений.
В функциях для решения дифференциальных уравнений, описанных ранее, предполагалось, что необходимо получить таблицу значений x(t) относительно значений t, расположенных на одинаковом расстоянии друг от друга на отрезке интегрирования, ограниченном точками tнач и tкон. Но часто требуется найти решение только в конечной точке tкон. Хотя описанные функции определяют значение x(tкон), они проделывают при этом много ненужной работы, высчитывая промежуточные значения x(t).
Если необходимо узнать только значение x(tкон), используются следующие функции:
bulstoer( x, tнач, tкон, acc, f, kmax, save)
rkadapt( x, tнач, tкон, acc, f, kmax, save)
stiffb( x, tнач, tкон, acc, f, J, kmax, save)
stiffr ( x, tнач, tкон, acc, f, J, kmax, save),
где
x – см. выше;
tнач, tкон – конечные точки интервала, на котором должно быть найдено решение дифференциальных уравнений. Начальные значения t определяются в точке tнач;
acc – контролирует точность решения. При небольших значениях acc делаются более мелкие шаги вдоль траектории, что увеличивает точность решения;
f(t, x) – см. выше;
J(t, x) – см. выше;
kmax – максимальное количество промежуточных точек, в которых будет найдено решение. Значение kmax устанавливает верхнюю границу количества рядов матрицы, получаемой этими функциями;
save – наименьшее допустимое расстояние между величинами, в которых должно быть найдено решение. Значение save устанавливает нижнюю границу разницы между любыми двумя числами в первой колонке матрицы решения.
Вернемся к нашему примеру с эпидемией. Если известно начальное число больных (N=50 – см. пункт 1 на рис. 5.1), то это значит, что их сразу пересчитали. А если это так, то их знают поименно. Но в этом случае больные должны быть изолированы, и никакой эпидемии не будет вообще (Пр=0). В реальной задаче мы можем знать число жителей в городе (число здоровых на день начала эпидемии) и число больных в какой-то последующий день, когда становится ясно, что разразилась эпидемия. Другими словами, нам что-то известно о параметрах на краях отрезка, охватывающего некий динамический процесс.
На рис. 5.3. реализован метод последовательных приближений для решения по разностной схеме такой краевой задачи: в начале эпидемии в городе 20 000 здоровых жителей, а в конце эпидемии – 100 больных. Спрашивается, сколько больных было в начале эпидемии.
