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

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

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

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

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

За шаг времени dt жертвы (prey) и хищники (pred) меняются на четыре величины, по числу параметров. Изменение жертв — dprey = (a * prey - b * prey * pred) * dt, изменение хищников — dpred = (c * prey * pred - d * pred) * dt.

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

Обратите внимание: член взаимодействия prey * 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 должен быть маленьким; проверяйте сходимость, уменьшая его вдвое.
Думать, что пики совпадают. Нет: пик хищников всегда отстаёт от пика жертв — в этом запаздывании вся суть модели.
Путать знаки членов взаимодействия. Один и тот же prey * pred у жертв идёт с минусом (их едят), а у хищников — с плюсом (они сыты).
Ждать, что система придёт к равновесию. В классической модели колебания не затухают — это вечные циклы, а не выход на постоянный уровень.
Забывать член -d * pred. Без него хищники не вымирали бы без еды, и модель потеряла бы смысл «качелей».

Итоги

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