Хищник-жертва (Лотка-Вольтерра)

До сих пор популяция жила сама по себе. Но в природе всё связано: зайцы кормят рысей, рыси едят зайцев. Стало много зайцев — расплодились рыси — съели зайцев — рыси остались без еды и вымерли — зайцы снова расплодились. Эти качели и описывает классическая модель хищник-жертва.

Модель Лотки-Вольтерры — модель двух связанных популяций: жертвы размножаются и поедаются хищниками, а хищники растут за счёт съеденных жертв и вымирают без них. В результате обе численности колеблются циклически, причём пик хищников отстаёт от пика жертв.

Зачем эта модель? Она объясняет загадку, которую веками наблюдали натуралисты: численность зверей в лесу не постоянна и не растёт монотонно, а колеблется волнами. Классический пример — записи пушной компании Гудзонова залива: шкурки зайца и рыси то густо, то пусто, и пики рыси аккуратно отстают от пиков зайца на пару лет. Две простые «связи» порождают эти вечные качели.

Четыре стрелки взаимодействия

За шаг времени $dt$ жертвы (prey) и хищники (pred) меняются на четыре величины, по числу параметров. В непрерывном виде система выглядит так:

$$\begin{cases} \frac{d(\text{prey})}{dt} = a \cdot \text{prey} - b \cdot \text{prey} \cdot \text{pred} \\ \frac{d(\text{pred})}{dt} = c \cdot \text{prey} \cdot \text{pred} - d \cdot \text{pred} \end{cases}$$

$a$ — рождаемость жертв. Член $a \cdot \text{prey}$: будь жертвы одни, они росли бы экспоненциально (как в уроке 5.1).
$b$ — успешность охоты. Член $-b \cdot \text{prey} \cdot \text{pred}$: чем больше встреч «жертва $\times$ хищник», тем больше съедено. Это минус для жертв.
$c$ — эффективность размножения хищников. Член $+c \cdot \text{prey} \cdot \text{pred}$: те же встречи дают хищникам энергию для потомства. Это плюс для хищников.
$d$ — смертность хищников. Член $-d \cdot \text{pred}$: без еды хищники вымирают сами по себе.

Обратите внимание: член взаимодействия $\text{prey} \cdot \text{pred}$ один и тот же в обоих уравнениях, только у жертв он со знаком минус (их едят), а у хищников — с плюсом (они сыты). Это и есть «связь» двух популяций.

Запустим дискретную имитацию методом явного Эйлера: на каждом крошечном шаге $dt$ считаем изменения и прибавляем их.

def lotka_volterra(prey, pred, a, b, c, d, dt, steps):
    out = []
    for s in range(steps):
        dprey = (a * prey - b * prey * pred) * dt
        dpred = (c * prey * pred - d * pred) * dt
        prey += dprey
        pred += dpred
        out.append((s, prey, pred))
    return out

res = lotka_volterra(40, 9, 0.1, 0.01, 0.005, 0.1, 0.05, 400)
print(f"{'шаг':>4} {'жертвы':>8} {'хищники':>8}")
for s, pr, pd in res[::50]:
    print(f"{s:>4} {pr:>8.2f} {pd:>8.2f}")

Вывод:

 шаг   жертвы  хищники
   0    40.02     9.04
  50    39.76    11.62
 100    36.81    14.65
 150    31.59    17.53
 200    25.47    19.51
 250    19.85    20.17
 300    15.46    19.57
 350    12.37    18.12

Запаздывание фаз — главный сюжет

Прочитаем числа как историю. В начале жертв много (40), хищников мало (9). Хищникам есть что есть — они начинают размножаться: к шагу 100 их уже почти 15, к шагу 250 — больше 20. Но пока хищники набирали силу, они выедали жертв: жертвы сползают с 40 до 20 и ниже.

И вот ключевой момент: хищники продолжают расти даже тогда, когда жертв уже стало мало — по инерции, их ведь много и они ещё доедают остатки (пик хищников около шага 250–270). А затем, когда еды совсем не осталось, хищники начинают вымирать (к шагу 350 их уже 18 и падает). Это и есть запаздывание: пик хищников приходит позже пика жертв. Хищник всегда «реагирует с опозданием» на численность добычи.

жертвы:   /‾‾\__        пик раньше
                \__
хищники:    __/‾‾\      пик позже (с запаздыванием)
         __/
         ----------------> время

Если продолжить имитацию дальше, жертвы (оставшись почти без хищников) снова расплодятся, хищники с запозданием снова поднимутся — и цикл повторится. Получаются незатухающие колебания: классические циклы Лотки-Вольтерры.

Как работает под капотом

Метод явного Эйлера — самый простой способ имитировать непрерывную систему. Идея: «производная — это скорость изменения, так умножим скорость на маленький шаг времени $dt$ и прибавим». Мы дробим время на 400 крошечных кусочков по $dt = 0.05$ и на каждом делаем линейный шажок.

Почему $dt$ должен быть маленьким? Потому что между шагами мы считаем скорость постоянной, хотя на самом деле она всё время меняется. Чем мельче шаг, тем меньше ошибка. Если взять $dt$ слишком большим, явный Эйлер начинает врать — траектория искажается, и численная «популяция» уходит от верной динамики, хотя настоящая система устойчиво колеблется. Сравним аккуратный и грубый шаги на одинаковом промежутке времени:

def lv_last(dt, steps):
    prey, pred = 40.0, 9.0
    a, b, c, d = 0.1, 0.01, 0.005, 0.1
    for _ in range(steps):
        dprey = (a * prey - b * prey * pred) * dt
        dpred = (c * prey * pred - d * pred) * dt
        prey += dprey
        pred += dpred
    return prey, pred

p1 = lv_last(0.05, 400)   # мелкий шаг, время = 20
p2 = lv_last(2.0, 10)     # грубый шаг, то же время = 20
print(f"dt=0.05:  жертвы={p1[0]:.2f}  хищники={p1[1]:.2f}")
print(f"dt=2.00:  жертвы={p2[0]:.2f}  хищники={p2[1]:.2f}")

Вывод:

dt=0.05:  жертвы=10.36  хищники=16.29
dt=2.00:  жертвы=7.25  хищники=18.12

Оба прогона покрывают одинаковый промежуток «времени» (20 единиц), но дают заметно разные числа: грубый шаг искажает динамику. Вывод практический — в явном Эйлере $dt$ надо брать заведомо маленьким и при сомнении уменьшать его вдвое, проверяя, что результат почти не меняется.

Важно помнить: это имитация непрерывной системы. Настоящая модель Лотки-Вольтерры записывается дифференциальными уравнениями, и её подробный вывод, анализ устойчивости и фазовый портрет разбираются в курсе про дифференциальные уравнения. Здесь мы лишь «прокручиваем» её по шагам, чтобы своими глазами увидеть колебания и запаздывание.

Частые ошибки

Брать слишком большой $dt$. Явный Эйлер тогда искажает динамику и даёт нефизичный результат. $dt$ должен быть маленьким; проверяйте сходимость, уменьшая его вдвое.
Думать, что пики совпадают. Нет: пик хищников всегда отстаёт от пика жертв — в этом запаздывании вся суть модели.
Путать знаки членов взаимодействия. Один и тот же $\text{prey} \cdot \text{pred}$ у жертв идёт с минусом (их едят), а у хищников — с плюсом (они сыты).
Ждать, что система придёт к равновесию. В классической модели колебания не затухают — это вечные циклы, а не выход на постоянный уровень.
Забывать член $-d \cdot \text{pred}$. Без него хищники не вымирали бы без еды, и модель потеряла бы смысл «качелей».

Итоги

Хищник-жертва — две связанные популяции: жертвы плодятся ($a$) и поедаются ($b$), хищники растут за счёт еды ($c$) и вымирают без неё ($d$).
Член взаимодействия $\text{prey} \cdot \text{pred}$ общий для обоих уравнений: минус для жертв, плюс для хищников.
Численности колеблются циклически, причём пик хищников запаздывает относительно пика жертв.
Мы имитируем систему явным Эйлером — дробим время на шаги $dt$; слишком большой $dt$ искажает динамику.
Это дискретная имитация непрерывной модели; её строгий вывод — в курсе про дифференциальные уравнения.