Моделирование динамики популяции на Python

В современном биологическом и экологическом исследовании моделирование динамики популяции является незаменимым инструментом. С помощью математических моделей можно оценить‚ как изменяется численность вида под воздействием факторов среды‚ конкуренции‚ хищничества и генетических процессов. В этой статье мы подробно разберём‚ какие модели используют учёные‚ какие требования предъявляются к программному обеспечению и как реализовать простую‚ но гибкую программу предсказания популяции на языке Python.

Теоретическая база

Существует несколько классических моделей‚ которые часто применяются в задачах предсказания численности:

  • Экспоненциальный рост – модель без ограничения ресурсов (модель Мalthus). Формула: N(t) = N₀·e^{rt}.
  • Логистическая модель (модель Рикера) – учитывает ограниченность среды. Уравнение: dN/dt = r·N·(1 ― N/K)‚ где K – ёмкость среды.
  • Модель хищник‑жертва (Лотка‑Вольтерра) – система двух дифференциальных уравнений‚ описывающих взаимодействие популяций.
  • Возрастно‑половая структура – более сложные модели‚ учитывающие различия в выживаемости и плодовитости разных возрастных групп.

Выбор модели зависит от целей исследования и доступных данных. В большинстве практических задач начинают с логистической модели‚ а при необходимости добавляют дополнительные параметры.

Что должна уметь программа?

  1. Ввод исходных параметров: начальная численность N₀‚ темп роста r‚ ёмкость среды K‚ длительность прогноза (кол-во шагов времени).
  2. Выбор модели: экспоненциальная‚ логистическая‚ Лотка‑Вольтерра и т.д.
  3. Численное решение уравнений: использование методов Эйлера‚ Рунге‑Кутты 4‑го порядка или готовых ODE‑солверов из scipy.integrate.
  4. Визуализация результатов: графики изменения численности во времени (библиотека matplotlib).
  5. Экспорт данных: возможность сохранить таблицу предсказаний в CSV.

Алгоритм работы программы

Ниже представлена общая схема:

Запросить у пользователя тип модели и параметры.
Сформировать функцию правой части ОДУ в зависимости от выбранной модели.
Выбрать численный метод (по умолчанию – RK4 из scipy.integrate.solve_ivp).
Запустить интегрирование на заданном интервале времени.
Построить график N(t) и вывести таблицу значений.
При необходимости сохранить результаты в файл.

Реализация на Python

Для примера реализуем две модели: экспоненциальный рост и логистический рост. Код написан так‚ чтобы легко добавить новые модели.


#!/usr/bin/env python3

import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import csv

def exponential(t‚ N‚ r):
 """dN/dt = r * N"""
 return r * N

def logistic(t‚ N‚ r‚ K):
 """dN/dt = r * N * (1 ─ N / K)"""
 return r * N * (1 ─ N / K)

def get_float(prompt‚ default=None):
 while True:
 try:
 line = input(f"{prompt} " + (f"[{default}] " if default is not None else ""))
 if not line and default is not None:
 return float(default)
 return float(line)
 except ValueError:
 print("Ошибка: введите числовое значение.")

def main:
 print("=== Программа предсказания размера популяции ===")
 print("Выберите модель:")
 print("1 – Экспоненциальный рост")
 print("2 – Логистический рост")

 choice = input("Введите номер модели (1/2): ").strip

 if choice not in ('1'‚ '2'):
 print("Неверный выбор.")
 sys;exit(1)

 # Общие параметры
 N0 = get_float("Начальная численность N0:"‚ 10)
 r = get_float("Темп роста r (в 1/ед.времени):"‚ 0.1)
 t_max = get_float("Длительность прогноза (в единицах времени):"‚ 50)
 steps = int(get_float("Количество точек вывода:"‚ 500))

 t_eval = np.linspace(0‚ t_max‚ steps)

 if choice == '1':
 # Экспоненциальный рост
 sol = solve_ivp(fun=lambda t‚ N: exponential(t‚ N‚ r)‚
 t_span=(0‚ t_max)‚
 y0=[N0]‚
 t_eval=t_eval‚
 method='RK45')
 model_name = "Экспоненциальный рост"
 else:
 # Логистический рост
 K = get_float("Ёмкость среды K:"‚ 100)
 sol = solve_ivp(fun=lambda t‚ N: logistic(t‚ N‚ r‚ K)‚
 t_span=(0‚ t_max)‚
 y0=[N0]‚
 t_eval=t_eval‚
 method='RK45')
 model_name = "Логистический рост"
 if not sol.success:
 print("Ошибка интегрирования:"‚ sol.message)
 sys.exit(1)

 N = sol.y[0]

 plt.figure(figsize=(8‚5))
 plt.plot(sol.t‚ N‚ label='N(t)'‚ color='green')
 plt.title(f'Прогноз популяции ({model_name})')
 plt.xlabel('Время')
 plt.ylabel('Численность')
 plt.grid(True)
 plt.legend
 plt.tight_layout
 plt.show

 # Сохранение в CSV
 save = input("Сохранить результаты в CSV? (y/n): ").strip.lower
 if save == 'y':
 filename = input("Введите имя файла (по умолчанию result.csv): ");strip
 if not filename:
 filename = "result.csv"
 with open(filename‚ mode='w'‚ newline='') as f:
 writer = csv.writer(f)
 writer.writerow(['time'‚ 'population'])
 writer.writerows(zip(sol.t‚ N))
 print(f"Результаты записаны в {filename}")

if __name__ == "__main__":
 main

Краткое пояснение к коду:

  • Для решения ОДУ используется scipy.integrate.solve_ivp с методом Рунге‑Кутты 4‑го порядка (по умолчанию RK45).
  • Параметры запрашиваются у пользователя‚ но можно легко заменить ввод на чтение из конфигурационного файла.
  • График строится с помощью matplotlib‚ что позволяет быстро оценить динамику.
  • Опция экспорта в CSV упрощает дальнейший анализ в Excel или R.

Тестирование и проверка корректности

Для уверенности в правильности работы программы рекомендуется:

  1. Сравнить результаты логистической модели с аналитическим решением: N(t) = K / (1 + ((K ─ N₀)/N₀)·e^{-rt}). При небольшом шаге интегрирования отклонения должны быть < 0.1%.
  2. Провести чувствительный анализ: изменить r и K в диапазоне ±10% и наблюдать‚ как меняется кривая.
  3. Для модели хищник‑жертва добавить вторую ОДУ и проверить наличие осцилляций‚ характерных для системы Лотка‑Вольтерра.

Расширения и дальнейшее развитие

После создания базовой версии программы можно добавить следующие возможности:

  • Возрастно‑половая структура – использовать матрицы Левенштейна (Leslie matrix) для расчёта переходов между возрастными классами.
  • Стохастические модели – включить случайные флуктуации (модель Брауна‚ процесс Пуассона) для более реалистичного предсказания небольших популяций.
  • Оптимизация параметров – реализовать подгонку модели к реальным наблюдениям с помощью scipy.optimize.curve_fit или генетических алгоритмов.
  • Графический пользовательский интерфейс (GUI) – построить простое окно на tkinter или PyQt5 для удобного ввода параметров.
  • Веб‑интерфейс – разместить приложение на Flask/Django‚ чтобы делиться результатами через браузер.

Программное предсказание размера популяции – это сочетание биологической интуиции и математической точности. Начав с простых моделей (экспоненциальный и логистический рост)‚ вы получаете быстрый инструмент для оценки влияния ключевых параметров среды. Дальнейшее усложнение модели (включение взаимодействий‚ возрастных структур‚ стохастики) позволяет адаптировать программу под конкретные исследовательские задачи.

Приведённый выше код полностью готов к запуску на любой современной системе с установленным Python 3‚ numpy‚ scipy и matplotlib. Вы можете сразу использовать его для учебных целей‚ а затем постепенно расширять‚ опираясь на теоретические основы‚ описанные в литературе (см. «Моделирование динамики популяции»‚ «Популяционная динамика» Г. Ю. Ризниченко и др.).


Дата генерации статьи: