97920

Распределенная модель «хищник -жертва»

Дипломная

Экономическая теория и математическое моделирование

Развитие математико-экологических моделей можно проследить по эволюции тех научных и прикладных вопросов, для ответа на которые эти модели создавались. Вопросы эти усложнялись по мере развития экологии и совершенствования методики моделирования. Если вначале сами вопросы и результаты математического моделирования представляли отвлеченный теоретический интерес...

Русский

2015-10-26

4.16 MB

3 чел.

Федеральное агентство по образованию

Государственное образовательное учреждение

высшего профессионального образования

«Нижегородский государственный университет

им. Н. И. Лобачевского»

Механико-математический факультет

Кафедра численного моделирования физико-механических процессов

Дипломная работа

на тему

Распределенная модель “хищник -жертва”

 

Исполнитель: студентка 655 группы

Цупко Д. А.

Руководитель: доц. каф. ЧМФМП,

к.ф.-м.н.   Сабаева Т.А.

                                                    Нижний Новгород, 2010

Содержание

Содержание……………………………………………………………………..2

Введение………………………………………………………………………...3

  1.Постановка задачи…………………………………………………………..5

  2. Простейшая модель "хищник-жертва”………………………………........6

  3. Усложненная модель "хищник-жертва”……………………………….......8

  4. Релаксационные колебания в системе “хищник жертва”.

   Переход к моделям теории катастроф………………………………….…11

5. Катастрофа типа сборки в системе "хищник — жертва"……………...…14

6. Распределенные модели……………………………………………………19

7. Различные типы граничных условий и их экологическая

      интерпретация. Возможные динамические режимы….…………………..21

  8. Система "хищник—жертва" с диффузионной неустойчивостью………..23

  9. Численный  метод решения……………………………………..….…….....28

Заключение………………………………………………………………....….....38

Список литературы………………………………….……………………….......39

Приложение……………………………………………………………….….…..40

Введение

Развитие математико-экологических моделей можно проследить по эволюции тех научных и прикладных вопросов, для ответа на которые эти модели создавались. Вопросы эти усложнялись по мере развития экологии и совершенствования методики моделирования. Если вначале сами вопросы и результаты математического моделирования представляли отвлеченный теоретический интерес, то в дальнейшем они стали носить конкретный практический характер. Значительная часть работ по моделированию природных экосистем имеет прикладной характер. Эти работы ставят перед собой практические задачи - построение прогнозов поведения во времени реальных биологических систем.

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

Модели развития экосистем можно разделить на два больших класса:

  •  дискретные (на основе систем алгебраических уравнений, которые в большинстве случаев можно получить путем дискретизации дифференциальных уравнений; в пространственно распределенном случае это характерные модели, называемые клеточными автоматами).
  •  непрерывные (на основе обыкновенных дифференциальных уравнений для сосредоточенных моделей и дифференциальных уравнений в частных производных для пространственно-распределенных систем).

Перед применением как непрерывных, так и дискретных  моделей необходимо доказать их адекватность, т.е. их соответствие реальному объекту, куда входят  способы оценки адекватности (например, путем сравнения с  наблюдениями).  Здесь же надо выяснить вопросы  непрерывной зависимости от начальных данных, параметров системы, допустимые области их задания и другие вопросы.

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

Взаимодействия между хищниками и их жертвами (так называемые отношения «хищник – жертва») приводят к тому, что эволюция хищников и жертв происходит сопряжено -  в процессе её хищники совершенствуют способы нападения, а жертвы – способы защиты. Следствием этих отношений являются сопряженные изменения численности популяций хищников и жертв.

Математическая модель совместного существования двух биологических видов (популяций) типа "хищник - жертва" впервые была получена А.Лоткой (1925 г.), который использовал ее для описания динамики взаимодействующих биологических популяций. Чуть позже и независимо от Лотки аналогичные (и более сложные) модели были разработаны итальянским математиком В. Вольтерра (1926 г.), глубокие исследования которого в области экологических проблем заложили фундамент математической теории биологических сообществ или так называемой математической экологии.[1],[2]

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

1.Постановка задачи

В работе производится численное решение системы, описывающей частный случай распределенной модели с возможностью диффузионной неустойчивости, строится вариант схемы аппроксимации задачи, рассматриваются условия ее применимости. Исследуется связь между сосредоточенной и пространственно распределенной моделью “хищник-жертва. Постановка задачи выполняется на двумерном ареале распространения особей, численности которых изменяются не только во времени, но и неравномерны в пространстве. И в связи с этим  рассматриваются варианты так называемой “диффузионной” неустойчивости, условия ее возникновения, а так же случаи, когда пространственно распределенная модель сводится к сосредоточенной.

      Результаты, полученные ранее для сосредоточенной модели,  используются в качестве

      теста к итоговой пространственно распределенной постановке задачи.

2.Простейшая модель "хищник-жертва

В данной главе кратко рассматриваются итоги исследования предыдущих курсовых работ, так как некоторые полученные в них результаты будут  использоваться при решении поставленной задачи. Ранее исследовалась простейшая модель типа “хищник – жертва”, которая не учитывает внутривидовую конкуренцию и предполагает неограниченное количество пищи для жертв. Эта модель состоит из двух дифференциальных уравнений:

                                            (1)       

Система имеет четыре параметра и два начальных условия.

Здесь X— численность  жертв. Y — численность хищников, k1 ..., k4  —  положительные параметры. Вообще говоря, выделение коэффициентов рождаемости и смертности представляет собой проблему. Существующие модели (в частности, рассматриваемая модель хищника и жертвы) записываются без явного разделения правых частей уравнений на коэффициенты рождаемости и смертности. Вариант, предложенный Лоткой, состоит в том, чтобы правую часть уравнений по формальному алгебраическому признаку разделить на две: члены со знаком «+» отнести к коэффициенту рождаемости, а члены со знаком « - » отнести к коэффициенту смертности. Подход в высшей степени произвольный, вряд ли может быть обоснован с экологически содержательной позиции, но он имеет право на существование.  В отсутствие хищников    численность жертв неограниченно возрастает: dX/dt>0, если Y=0. В отсутствие жертв хищники вымирают: dY/dt<0 при всех t, если . Соответственно, число вариантов решения системы очень велико. Приведем график решения, показывающего характерную периодическую особенность поведения численности популяций хищников и жертв. Вначале идет размножение жертв. Их число увеличивается - хищников еще мало. Некоторое время спустя начинается очень быстрый рост числа хищников. Быстро поедая жертв, и, лишаясь пищи, хищники начинают умирать, их численность падает до минимального значения. При минимуме хищников вновь начинается постепенный рост популяции "жертв". Затем все повторяется. Решение уравнений модели носит колебательный характер. (На рис.1: горизонтальная ось – время, вертикальная – численности жертв(верхний график) и хищников(нижний график).

Рис. 1.Численное решение простейшей модели хищник – жертва при k1 =0.1, k2 =1.0,

k3 =0.5,k4 =2.0

В данном случае существуют два положения равновесия. Линеаризация системы (1) в окрестности каждого из них и последующее составление характеристических уравнений  дают: Р1(0, 0) — седло; P2 (k3/k4 , k1/k2) — центр. Точка Р1 определяет положение равновесия, которое характеризуется полным истреблением жертв (X= 0) и вымиранием хищников (Y = 0). Точка  Р2 отображает стационарный режим  сосуществования  хищников и их жертв  с некоторыми ненулевыми численностями.[1]

Колебания численностей хищников и их жертв весьма часто встречаются в природе, и модель отображает этот факт. Однако рассмотренная система как и всякая консервативная система, является негрубой: при малых изменениях ее правых частей происходят качественные изменения в ее фазовом портрете. По-видимому, модель реальной биологической системы должна быть грубой, а колебания должны определяться не начальными условиями, а внутренними свойствами системы, т.е. это должны быть автоколебания. Данное соображение послужило стимулом для разработки новых, автоколебательных, моделей систем типа “хищник — жертва”. [2],[3]

3. Усложненная модель хищник – жертва

Ранее была описана простейшая модель совместного существования двух биологических видов, один из которых питается особями другого. Эта простейшая модель не учитывает многих факторов, реально существующих в природе.

Отличие данной модели от простейшей заключается в использовании более сложных гипотез о динамике роста жертвы:

- емкость среды ограничена некоторой величиной, и безграничный рост жертвы в отсутствие хищника невозможен;

- существует нижняя критическая численность  жертвы, и, если число особей падает по каким-либо причинам ниже этой величины, популяция вымирает.

- существует периодическое (сезонное) изменение емкости среды, которая определяется как максимальное число особей, которые могут существовать в данной среде.

Рис.2. Численное решение усложненной модели с простой зависимостью численностей от времени.

Рис.3. Фазовый портрет для усложненной модели с простой зависимостью численностей от времени.

Пример простого регулярно повторяющегося решения. Вначале, когда хищников мало, численность жертв имеет сезонные колебания. Постепенно число хищников растет. Оно также зависит от сезона. Число хищников нарастает с увеличением жертв, уменьшая их количество. Колебания численности имеют простую зависимость от времени (на рис.2 – верхний и нижний графики – соответственно зависимости численности жертв и хищников в данной модели).

Наконец, пример сложного нерегулярного, неповторяющегося решения (рис.4), (рис.5). Колебания численности имеют сложную зависимость от времени.

  Рис. 4. Численное решение усложненной модели со сложной зависимостью численностей от времени.

Рис.5. Фазовый портрет для усложненной модели со сложной зависимостью численностей от времени.

Если амплитуда внешнего воздействия возрастает, колебательный цикл в системе “хищник-жертва” разрушается, возникает нерегулярный, хаотический режим. Таким, образом, сильные периодические колебания условий среды могут приводить к псевдослучайным, нерегулярным колебаниям динамики численности популяций.

4. Релаксационные колебания в системе “хищник жертва”.

Переход к моделям теории катастроф

Пусть х(t), у(t) — численности популяций "жертвы" и "хищника" соответственно. Наиболее популярная модель системы "хищник — жертва" записывается в виде

                                (2)

                 

Здесь α(х) — мальтузианская функция "жертвы", V(х) — трофическая функция, т.е. скорость потребления жертвы одним хищником, k  — к.п.д. переработки биомассы жертвы в новую биомассу хищника,  m — коэффициент его естественной смертности.

Рассмотрим взаимодействие двух факторов в системе "хищник —  жертва": саморегулирования в популяции жертв и S-образность трофической функции хищника. Зададим максимально простую параметризацию мальтузианской функции жертв:  - логистический рост, и трофической функции хищника: —  S-образная функция. Делая замену переменных:

, ,

и обозначая m/kA=b, вместо системы (2) получаем

           (3)

   

Предположим, что жертва размножается со скоростью гораздо большей, чем хищник, т.е. α0 >> 1. Тогда 1/ α0< < 1 можно рассматривать как малый параметр, переменную X как "быструю", а Y — "медленную". Динамика системы в значительной степени будет определяться видом равновесного многообразия:

,      

которое распадается на два:

;   .  

Второе многообразие зависит только от одного параметра Xkемкости среды для популяции жертв. На рис.6. в плоскости (X, Y) изображены два вида этой кривой при различных значениях Xk. Пусть Хk>27, и равновесное многообразие имеет вид, изображенный на рис. 7. Так как поведение системы зависит еще от параметра b (значение ), то рассмотрим три варианта поведения при различных b = m/kА.

1.  (см. рис. 7, на котором для удобства дальнейшего рассмотрения равновесное многообразие изображено в координатах X = X(Y)),

где, , . Экологически это означает, что хищник хорошо адаптирован к среде, "давление" хищника на жертву велико, так что равновесная численность жертвы мала (b и X* связаны монотонной зависимостью) . Предположим, что сначала система находится в точке А. С ростом переменной Y переменная X убывает, пока не достигнет точки В, где переменная Y пересекает значение С2 , и система совершает "катастрофический" скачок на нижнюю ветвь многообразия в точку В/. Дальнейшее возрастание переменной Y уводит систему  к равновесной точке E1.

Рис. 6. Равновесное многообразие                    Рис. 7.  Типы динамик    

Y = F (X)для системы (3);                                    системы "хищник - жертва" при

;                                     различных значениях параметра b = m/kА.

2. Пусть теперь , где , , т.е. хищник плохо адаптирован к среде, его "давление" на жертву мало, и равновесная численность жертвы велика. Предположим, что сначала система находилась в точке D. Так как , то переменная Y будет убывать, и система следует по нижней ветви многообразия к точке G. В этой точке переменная Y пересекает значение C1, и система "катастрофически быстро" перескакивает в точку G/ на верхнюю ветвь многообразия, после чего спокойно движется по этой ветви к точке Е2 - положению равновесия.

3. Пусть теперь хищник имеет промежуточную адаптацию, т.е. ; . Вначале система находилась в точке А. Затем она начала двигаться к точке В, в этой точке происходит "катастрофа", система перескакивает на нижнюю ветвь, но, в отличие от ситуации, описанной в пункте 1, система начинает двигаться от точки В/ не вправо, а влево, к точке С. В точке G новая "катастрофа", система перескакивает в точку G/ на верхнюю ветвь, и снова начинает двигаться к точке В.  Цикл замкнулся, возник так называемый "релаксационный" цикл. Описываемые им колебания, для которых характерны чередования достаточно гладких режимов и "катастроф", называются обычно "релаксационными". Существующее здесь равновесие (точка Е3) неустойчиво, но сам релаксационный цикл устойчив.[6]

Здесь можно говорить, что в точках Y = С1,2 отображение Х(Y) прямой на прямую имеет особенность. Эта особенность является складкой и представляет собой простейший тип катастрофы. Но, эта особенность устойчива, т.е. она не исчезает при малых шевелениях самого отображения. Устойчивость такого типа называется еще структурной.[1]

                  5. Катастрофа типа сборки в системе "хищник — жертва"

До сих пор мы рассматривали равновесное многообразие Х = Х(Y) как отображение прямой на прямую при фиксированном Xk. Будем теперь менять параметр Xk, тогда многообразие Х=X(Y, Хk) представляет собой поверхность в трехмерном пространстве   (X, Y, Xk), задающую отображение двумерного многообразия в двумерное. Если нарисовать поверхность Х = Х(Y, Xk) в координатах (X, Y, Xk), то мы получим картину, изображенную на рис.8. В точке С это отображение имеет особенность — сборку, и соответственно, здесь мы получаем более сложную "катастрофу" — типа сборки. Эта особенность устойчива.

Рис.8 Поверхность X = Х(Y, Xk) и ее проекция на плоскость {Y, Xk}. На кривых С1L1 и С1K1, отображение особенности типа складки, в точке С1 – сборка.

Заметим, что направление движения по этому многообразию зависит от величины параметра b - показателя адаптации хищника к среде, который, в свою очередь, тесно связан со значением равновесной численности жертвы X*.

До сих пор наши рассуждения были тесно связаны с исходной динамической системой. Сделаем шаг в направлении более общей ситуации. Так как второе уравнение системы (3) дает нам лишь информацию о направлении движения по равновесному многообразию, то откажемся от него, а из первого уравнения сохраним связь между численностями хищника и жертвы и "емкостью" среды, обеспечивающую нам стационарную численность жертвы. По сути дела, мы оставляем только уравнение баланса жертвы при равенстве "прихода" и "расхода":

α(Х) X (прирост за счет размножения) = Y(Х) Y (выедание хищниками).

Если не рассматривать пока малоинтересное состояние с X*= 0, то при принятой параметризации связь между X, Y и Xk запишется в виде

 

или, после замены

,  ,

в виде

     

.

Это в точности поверхность катастрофы сборки. Тем самым мы подтвердили сделанный выше на основе лишь геометрических соображений вывод.

Вычислим координаты точки сборки С. В этой точке , откуда сразу получаем

,  , ,

или в старых координатах

, , .

На первый взгляд кажется, что теория катастроф не дает никакой новой информации по сравнению с проведенным нами обычным исследованием, но это не так. При нашем исследовании мы получили особенность типа сборки, параметризовав мальтузианскую и трофическую функции. Но в любой реальной ситуации мы никогда не знаем точного их вида, известны лишь качественные особенности (например, то, что мальтузианская функция уменьшается с ростом численности и достигает нуля при ее конечных значениях, а трофическая функция — это S-образная функция с насыщением). И неизвестно, сохранится ли эта особенность для функций, несколько отличающихся от используемых нами при анализе, не возникнут ли другие особенности? А теория катастроф позволяет утверждать, что эта особенность сохранится, а другие особенности неустойчивы и разрушатся при вариациях мальтузианской и трофической функции, или же просто при возмущениях, действующих на систему. Другими словами, какими бы мы ни выбрали трофическую и мальтузианскую функции (важно только, чтобы они были качественно похожими), поведение системы "хищник — жертва" качественно будет аналогичным описанному выше.

Дадим наглядную интерпретацию этим утверждениям. Будем характеризовать систему "хищник — жертва" тремя параметрами:

X, Y — численностями жертвы и хищника соответственно, и Xk — "емкостью" среды для жертвы, т.е. устойчивой предельной численностью жертвы в отсутствие хищника. Естественно, что между этими параметрами существует зависимость, описываемая поверхностью в трехмерном пространстве с координатами (X, Y, Xk). Спроектируем эту поверхность вдоль оси X на плоскость (Y, Xk). Выше было показано для частного случая, что сборка существует. Можно утверждать, что сборка, расположенная так, как это изображено на рис. 8, удовлетворительно описывает динамику системы "хищник — жертва" в случае общего положения, когда не известны точно мальтузианская и трофическая функции.[1]

Рассмотрим прикладное применение  нашей теории. Самым перспективным методом в борьбе с сельскохозяйственными и лесными вредителями в настоящее время считается биологический метод. Суть его состоит в том, что для вредителя — "жертвы" подбирается "хищник", который может регулировать численность популяции жертвы, снижая ее численность. Другими словами, искусственно создается система "хищник — жертва". Предположим, что подобран такой хищник, который хорошо адаптирован к среде, так что параметр m/kА мал. Соответственно будет низкой и равновесная численность жертвы. Будем считать, что его трофическая функция S-образна.

Посмотрим, как будет меняться численность вредителя в зависимости от его "емкости" среды и численности хищника. Если при малой "емкости", когда и в отсутствие хищника вредитель не достигал большой численности, увеличивать численность хищника, то численность вредителя монотонно и довольно медленно уменьшается (кривая АВ на рис. 8). Однако если "емкость" велика (а по-видимому, реально нас интересует именно это ситуация), то возникают качественно новые явления. Из рис. 9 видно, что численность вредителя будет меняться скачком (кривая A/DD/B/), поэтому при определенной численности хищника достаточно малого ее увеличения, чтобы получить высокую эффективность метода. Но в этой "катастрофичности" заключена и серьезная опасность. Дело в том, что до наступления "катастрофы" мы не видим особого эффекта от применения биологического метода и поневоле можем сделать вывод о неэффективности самого метода в целом.

Рассмотрим теперь ситуацию, когда резко увеличивается "емкость" среды для жертвы. Для листогрызущих насекомых это может, например, произойти, когда несколько дней подряд была теплая и влажная погода, и появилось большое число свежих и сочных листьев. Пусть мы наблюдаем за системами "хищник - жертва", которые в начальный момент имели одинаковые численности жертв и почти одинаковые численности хищников (точки 1 и 2 на рис. 8). Нетрудно видеть, что из-за наличия особенности при увеличении  (достаточно быстром, чтобы переменная Y практически не менялась) одна система перейдет на верхний лист многообразия (из точки 1 в точку 3), а вторая -на нижний (из точки 2 в точку 4). Этот пример иллюстрирует еще одну особенность систем с катастрофами: при одном и том же изменении параметра даже из близких начальных состояний система может прийти в весьма далекие друг от друга конечные состояния.

Рассмотрим  еще один пример системы с "катастрофическим" поведением. Предположим, что мы имеем пруд, в котором фитопланктон служит пищей для рыб. Тогда фитопланктон — это жертва, а рыбы — хищник. Обычно размножение фитопланктона лимитируется количеством растворенных биогенных веществ — азота и фосфора. Одним из способов повышения продуктивности пруда — увеличения биомассы рыбы — служит искусственная минеральная подкормка, когда в пруд вводят минеральные удобрения, надеясь, что это приведет к увеличению продуктивности фитопланктона и, как следствие, к увеличению биомассы рыбы. С нашей точки зрения это воздействие эквивалентно скачкообразному увеличению "емкости" среды для фитопланктона. Обычно при этом об эффективности подкормки судят по увеличению биомассы (или точнее плотности) фитопланктона. Рассмотрим возможные случаи.

Пусть сначала система находилась в точке А (см. рис. 8). После внесения подкормки она скачком перешла в точку А' и попала в область катастрофических режимов, где ее эволюция будет происходить вдоль кривой A/DD/B/. В точке D система проходит через катастрофу, в результате которой резко уменьшается плотность фитопланктона. И если в этот момент мы эту плотность измерим, то придем к выводу, что наше управление было неэффективным, хотя на самом деле биомасса рыб Y продолжает увеличиваться.

Похожую ошибку можно совершить, если из начального состояния, задаваемого точкой 2 на рис. 8, при увеличении Xk попасть на нижний лист многообразия, на котором плотность фитопланктона мала. Избежать этих ошибок (как, впрочем, и ошибки при экспериментах по определению эффективности биологического метода борьбы с вредителями) можно, если избегать значений Y, близких к , или, в более реальных обозначениях, близких к . Мальтузианская скорость роста фитопланктона обычно известна, так же как известны и параметры трофической функции рыб при выедании ими фитопланктона.[2]

Нужно отметить, что во всех приведенных выше качественных рассуждениях содержится довольно много неточностей, что, впрочем, характерно для подавляющего числа работ по приложениям теории катастроф. Однако качественные эффекты в динамике системы "хищник – жертва" здесь отражены верно. Можно точно сказать, что при малой "емкости" среды для жертвы, при сильном саморегулировании ее популяции в системе отсутствуют "катастрофические" режимы (скачки, релаксационные циклы). При увеличении "емкости", когда саморегуляция в популяции жертвы ослабевает, в системе появляются "катастрофы" (циклы, скачки, явления гистерезиса, когда поведение системы определяется ее предысторией; дивергенция, когда эволюция, начавшись в близких состояниях, может приводить к далеким конечным состояниям.

  1.  Распределенные модели

Любая особь популяции  в процессе своей жизнедеятельности перемещается по ареалу – имеет свой радиус индивидуальной активности. Если радиус больше, чем характерный размер ареала, то популяцию можно считать хорошо перемешанной и при ее описании можно не учитывать пространственных эффектов. Если же мы имеем обратную ситуацию, то рассмотрение пространственного распределения популяции становится необходимым.

Пусть в промежутке между рождением и размножением каждая особь перемещается в случайном направлении( все направления равновероятны) на некоторое расстояние. Тогда, если p(r) dr – вероятность перемещения на расстояние, лежащее

между r и  r+dr, за единицу времени, измеряемого в поколениях, тогда

- есть среднеквадратичное перемещение за поколение или радиус индивидуальной активности.[1]

Предположим, что процесс миграции и процессы размножения и гибели не связаны между собой. Пусть динамика плотности локально( в точке) описывается уравнением

,

Тогда приращение плотности за время  в точке (х,у)за счет миграции и локального роста равно

, где

Раскладывая N(x',y',t) в окрестности (х,у) и предполагая, что третий момент достаточно мал,<< можно ограничиться в этом разложении членами второго порядка пои .Тогда, переходя к пределу при , получаем уравнение вида  - уравнение нелинейной диффузии, где

, а D=  - коэффициент диффузии.[1]

Ранее  была рассмотрена динамика неравновесной экологической структуры “хищник-жертва”, распространяющейся по бесконечному ареалу. Эти структура есть не что иное, как переходный  процесс между стационарными состояниями, каждое из которых однородно по пространству.

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

Введем новое обозначение для распределенных моделей. Пусть (x,y,t ) — это плотность особей i-го вида в точке (х,y) G, где область G — ареал обитания рассматриваемого сообщества.

Формально диссипативные структуры - это ограниченные устойчивые пространственно неоднородные  решения   системы нелинейных диффузионных уравнений вида

   (4)

i= 1,…,n   

с соответствующими граничными и начальными условиями.

Очевидно, что N* ={ * }, где F(N*) = 0  (i=1,..,n) – однородное по

пространству стационарное решение системы . Тогда одним из возможных пусковых механизмов возникновения пространственно-неоднородной устойчивой структуры может служить потеря устойчивости этого решения под действием пространственных возмущений. [1],[2].

7. Различные типы граничных условий и их экологическая интерпретация. Возможные динамические режимы

В предыдущей главе употреблялся термин "соответствующие граничные условия". Рассмотрим этот вопрос. Мы считаем, что G — ограниченная область, тогда характер решений будет сильно зависеть от граничных условий. Наиболее простым (и наиболее естественным с экологической точки зрения) будет условие непроницаемости границы области, т.е.

   , где п — вектор нормали к границе G.

Это условие непроницаемости; оно описывает полную изоляцию ареала, абсолютный изоляционный барьер.

Но вполне возможна, например, и другая ситуация, когда на границе ареала возникают участки с оптимальными условиями обитания для некоторых видов, оптимальные биотопы или заповедники. Ее можно описать граничными условиями типа

  ,  где за   удобно принять равновесные значения плотностей.

Конечно, гипотезы, лежащие в основе модели, которая описывается приведенными выше уравнениями (4) с рассмотренными граничными условиями, с экологической точки зрения слишком упрощены. В самом деле, мы предполагаем, что одни и те же типы взаимодействий существуют в различных точках пространства и в разные моменты времени, что перемещение особей в пространстве происходит самым простым образом, что миграция случайна и т.п. Все это с позиций конкретной экологии можно критиковать, но тем не менее даже такая упрощенная модель позволяет получить ряд интересных режимов пространственной динамики популяций. Все многообразие этих режимов даже в таких простых моделях в настоящее время полностью не исследовано. Однако ясно, что, кроме простейших стационарных однородных по пространству решений, здесь могут существовать неоднородные стационарные решения, сложные колебательные режимы и динамическая стохастичность.

Наиболее изучены динамические режимы, возникающие в таких системах тогда, когда стационарное однородное решение теряет устойчивость. Другими словами, когда спектр системы, линеаризованной в окрестности равновесия, находится в открытой левой полуплоскости, за исключением нескольких собственных значений. Именно эта ситуация и называется в широком смысле слова диффузионной неустойчивостью. Если множество таких неустойчивых собственных значений состоит больше чем из одного элемента, то исследовать характер динамического поведения системы аналитически не удается — остаются только численные методы. Если же спектр состоит из объединения одного простого собственного значения, равного нулю, и части, лежащей в открытой левой полуплоскости, то в этом случае возможно аналитически исследовать устойчивость возникающих здесь решений и найти их вид. Эти неоднородные стационарные решения называют "мягкими" диссипативными структурами. [1].


8. Диффузионная неустойчивость в сообществе типа "хищник — жертва"

Пусть сообщество, описываемое моделью (4), имеет двумерный ареал, так что

Ni= Ni (x,у,t). Для исследования устойчивости решения N = const=N* линеаризуем (4) в его окрестности :

 ,      i = 1,…n      (5)

Здесь Ui = Ni - , aij = () N* . Будем искать решение (5) в виде

(x,y,t)=соs1х + k2y + )уе (6)

Такой вид позволяет удовлетворить достаточно широкому классу начальных условий (все симметричные функции, разложимые в ряд Фурье). Подставляя (6) в (5), мы получим, что должны быть собственными значениями матрицы  [- ], где к2 =  , - символ Кронекера. Ясно, что если даже все собственные значения матрицы ||aij || лежат в левой полуплоскости (т.е. равновесие N* в отсутствие диффузии устойчиво), то при Di 0 вполне вероятно, что при определенных волновых числах к некоторые из  будут лежать правее мнимой оси и амплитуда всех возмущений с этими волновыми числами будут возрастать, т.е. возникает типичное явление неустойчивости. Неустойчивость такого типа называется диффузионной.

Рассмотрим возникновение диффузионной неустойчивости в системе "хищник - жертва", но сначала выпишем условия, обеспечивающие возникновение диффузионной неустойчивости в системе  при п =2.  Ясно, что равновесие (* ) локально (т.е. в отсутствие диффузии) устойчиво, если

<0,      > 0 (7) 

Диффузионная неустойчивость возникает в случае, если хотя бы одно из условий

-  <0,                         

 > 0  (8)

нарушается. Очевидно, что если имеет место первое из неравенств (7), то будет выполнено и первое из (8). Следовательно, нужно, чтобы при некоторых к2 многочлен

был отрицательным. Условие того, что всегда найдется хотя бы одно , при котором  < 0, т.е. условие возникновения диффузионной неустойчивости запишется в виде

>  (9)

Вернемся теперь к системе "хищник —жертва". Предположим, что в ней имеется нетривиальное равновесие. Из чисто биологических соображений следует, что

<0,   

>0,      

  

Мы считаем, что  N1 — это численность жертвы, a N2 — численность хищника. Последнее неравенство указывает, что вид "хищник" либо самолимитирован (строгое неравенство), либо в его популяции внутривидовая конкуренция отсутствует, и его динамика определяется только скоростью поглощения трофического ресурса и естественной смертностью, пропорциональной численности (равенство нулю). Если и относительно жертвы сделать аналогичные предположения (т.е. либо самолимитирование, либо отсутствие конкуренции), то    < 0

Итак, если  < 0, то неравенство (9) при любых D1 и D2 не будет иметь место и диффузионная неустойчивость не может возникнуть. Для ее возникновения нужно, по крайней мере, чтобы

> 0. Но, с другой стороны, если > 0, то для локальной устойчивости необходимо, чтобы < 0.

Чтобы дать содержательные интерпретации условиям возникновения  диффузионной   неустойчивости, мы  конкретизируем вид

функций F1 и F2 и рассмотрим несколько вариантов системы "хищник—жертва''.

Пусть

                        

Если = const, m = const, V(N) — трофическая функция одного из двух типов, то = 0. Уже отсюда видно, что здесь не будет эффектов, связанных с диффузионной неустойчивостью. Пусть теперь по-прежнему = const, но m = m(N), m' (N) > 0. Тогда

= - m'(N) N < 0.

Рассмотрим выражение для

= -

Ясно, что если трофическая функция первого типа, то  >0 для любых N и N, а если второго - то для тех N, для которых < 0. Следовательно, здесь уже можно ожидать эффектов, связанных с возникновением диффузионной неустойчивости.

Если же трофическая функция вольтерровского типа, то при = const, = 0. При = ( N1),  =  '(N) N  и  > 0 только при '(N) > 0, т.е. диффузионная неустойчивость возникает лишь тогда, когда популяция жертвы — это популяция типа Олли. Для таких популяций характерно, что существуют численности, при которых они растут быстрее, чем экспонента. И окончательно, в исследуемой системе “хищник — жертва” возникновение диффузионной неустойчивости (при локальной устойчивости равновесия) возможно лишь в том случае, когда естественная смертность хищника возрастает с ростом его численности быстрее, чем линейная функция, и  трофическая функция отличается от вольтеровской, либо когда популяция жертвы — это популяция типа Олли. [1]

Рассмотрим конкретный пример системы "хищник — жертва", в которой при определенных значениях параметров может возникнуть диффузионная неустойчивость. Пусть             ()=, V()=, m()=т0     

Рис.9. Области диффузионной неустойчивости (заштрихованы) в плоскости параметров (а) ( = D2/D1, 0 ),при  т0 = 1   и (б) ( , т0),при а0 = 1 .

Очевидно, что при  , () 0 , т.е. при достаточно больших численностях, популяция жертвы (в отсутствие хищника) возрастает экспоненциально. При малых же   ()  и   .

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

Далее, мы видим, что трофическая функция — вольтерровская, но коэффициент смертности хищника не постоянен (как в вольтерровской модели), а линейно возрастает с ростом численности. Это вполне обычная ситуация, если среди хищников существует сильная внутривидовая конкуренция не за трофический, а за какой-либо другой ресурс.

Не нарушая общности, можно положить К = 1,  = 1, = 1,  

Легко показать, что при 0m0 > 1 и  0m0 > 1 локальное равновесие всегда существует и устойчиво. А если выполняется неравенство

то в системе возникает диффузионная неустойчивость. Очевидно, что при ( а остальные параметры порядка единицы) последнее неравенство будет выполняться наверняка. Интересно, сколь велики будут области диффузиционной неустойчивости в пространстве параметров = D2/ D1, 0 и m0.

Пусть сначала т0 = 1, а затем 0 =1. На рис.9  изображены эти области. Из них видно, что при выборе одинаковых масштабов параметров , 0  и m0 и при  ~ 1 ширина областей диффузионной неустойчивости по 0  и т0 очень мала, области узкие, и вероятность попасть туда низка. Однако с ростом , т.е. когда различие в подвижностях видов возрастает, ширина областей также растет. [1]

Таким образом, распределенная модель является обобщением точечной модели, с тем отличием, что численности хищников и жертв становятся функциями не только времени, но и координат.

На основе этого вывода, преобразуем сосредоточенную модель (3) в пространственно – распределенную вида(4). Введем обозначения, принятые для пространственно распределенных моделей, т.е. численность жертв описывается функцией N(x,y,t), численность хищников -  M(x,y,t), где x,у – пространственные координаты, t – время.

Таким образом, получаем:

    (10)

             Данная пространственно распределенная задача принимает вид (3), если численности особей (N,M) являются константами или имеют линейную зависимость от координат ( это приведет к обращению в ноль вторых производных по координатам).

Таким образом, находя численное решение (10), все выводы, полученные для сосредоточенной модели(3), мы будем использовать в качестве теста для проверки этого решения, выбирая функции численностей определенным образом, что и является одной из поставленных целей.

9.Численный метод решения

Практическая часть курсовой работы состоит в:

1) численном решении системы вида

,

описывающей  поведение системы   ''хищник-жертва '' , распределенной на ограниченном  двумерном ареале.

Здесь во избежание сложной индексации   приняты обозначения:  N(x,y,t) – численность популяции жертв , M (x,y,t) –численность популяции хищников. Ареал распространения особей двумерный.

- константы, характеризующие различные параметры популяции.

Исследование вариантов развития сообщества производится при :

1.Граничных условиях Неймана  ,что  обеспечивает полную изоляцию ареала распространения особей.

2. Начальных условиях  φ1(x,y)

                                         φ2(x,y)

2)построении функций N(x,y), M(x,y) на отдельных временных слоях.

Отдельно проиллюстрируем случай, когда вторые частные производные по х, y функций N(x,y,t), M(x,y,t) равны 0, то есть наша модель  сведется  к нераспределенной, варианты которой были рассмотрены в  предыдущей курсовой работе. Этот вариант и будет тестировать полученное решение.

Численное решение

Алгоритм решения приведен для первого уравнения системы, описывающего        численность жертв  - N (x,y,t).

Для M(x,y,t) вся схема абсолютно аналогична.

1) Сначала строим конечно-разностную схему нашей задачи.

Будем решать двумерную систему дифференциальных уравнений в частных производных одной из разностных схем  метода переменных направлений, называемую продольно-поперечной разностной схемой (схемой Писмена-Рэчфорда).

Эта схема сочетает лучшие качества явной схемы — экономичность  и неявной — устойчивость.[4],[5].

Введем сетку по времени

={t=k, k=0,1…K-1,K=T}

и пространственную сетку

={u=(x,y), x=ih, y=jh},

i=0,1,…n-1, j=0,1,…m-1, nh=a, mh=b.

Будем обозначать N= N (u, t)(все принципы рассуждений и обозначений аналогичны и для функции численности популяции хищника M).

В методе переменных направлений переход от n к n+1 слою осуществляется в два этапа.  На первом этапе определяются промежуточные значения N из системы уравнений

(N- N)/0.5= D(+) +F(N,M,)

а на втором этапе, пользуясь найденными значениями N, находим N

из системы

(N- N)/0.5= D(+)+ F(N,M,)

для внутренних  точек сетки ( когда i=1,…n-2, j=1,…m-2). Обозначим их .

Система уравнений для первого этапа неявна только по переменному х, а для второго – только по у, поэтому данные системы можно решить последовательным применением одномерных прогонок сначала по направлению х, потом по направлению у.

Делим ареал на (n*m) кусков, тогда шаги по пространству h=a/n,

h=b/m. Аппроксимируем вторые производные  центральными разностями:

где i=1,…n-2

                                   где j=1,…m-2

в точках i=0, i=n-1,а j – любое; j=0,j=m-1,а  i - любое (т.е. на границе ареала) аппроксимация будет другой, учитывая граничные условия. Построение выполняется абсолютно аналогично одномерному случаю, который был полностью рассмотрен ранее. Например:

Теперь аппроксимируем производную по времени разностью вперед:       

Поскольку система уравнений для первого этапа неявна только по х, а второго– только по у, то их можно решать последовательным применением одномерных прогонок сначала по направлению х (то есть при каждом фиксированном j по переменному i), потом по направлению у(при каждом фиксированном i по переменному j).

Таким образом, получаем две системы линейных уравнений для нахождения матриц  N, M       - значений функций численности популяций жертв и хищников  на полушаге по времени.    Запишем систему первого этапа для N, для M отличие только в нелинейном члене:

(+)+

…………………………………………

(+)+

+

………………………………………….

(+)+

Это матрица для j=0; последовательно меняя вторые индексы, будем получать системы для всех фиксированных j. Таким образом получим  N, M.

Повторяя эти рассуждения, из второго этапа получим N,M. Первый шаг сделан  - найдены функции численностей на первом временном слое. Далее все повторяется снова – для второго слоя и т.д.

Схема переменных направлений абсолютно устойчива. Каждое из уравнений  первого и второго этапа имеет первый порядок точности по времени и второй по пространству, т.е. 0( +h). [5]

Система сведена к виду [A1]{ N}={r1}

Итак, при использовании продольно – поперечной схемы получаем 3-х диагональную матрицу A1, систему решаем методом прогонки.

Обозначим элементы матрицы A1=b, i=0..n-1.

A=a, i=1…n-1, a1=0

A=c, i=1…n-1, cn=0.

i-ый элемент вектор – столбца r1 соответственно обозначим d.

Прогоночные коэффициенты вычисляются по формулам

i=0:

A1= - c/b, B1=d/b

i= 1..n-2:

Ai= - c/(b+aA), Bi=(d-aB)/(b+aA)

i=n-1

An-1=0, Bn-1=(d-aB)/(b+aA), то есть N= Bn-1.

После вычисления всех прогоночных коэффициентов, учитывая формулу N=AiN+1+Bi, i=0…n-1, производим обратный пересчет. Полная аналогия для второй системы.

Для данных систем коэффициенты подбираются таким образом, чтоб  выполнялось:

Корректность и устойчивость метода прогонки (достаточное условие):

Пусть коэффициенты удовлетворяют:

  1.  | a|0, | b|>0, | c |0, i=0..n-1.
  2.  | b|| a| + | c|, i=1...n-2.
  3.  | b, || c |, | b|| a|

причем хотя бы в одном из неравенств 2) или 3) выполняется строгое неравенство, т.е. матрица системы имеет диагональное преобладание.

Тогда гарантируется корректность и устойчивость метода.

Графический результат решения, полученный с помощью пакета Matlab:

1.N=const,M=const (тестовая задача)

При таких начальных функциях  обращаются в ноль, и задача сводится к нераспределенному случаю.

Значения параметров:

D=0.001;         

D=0.001;

τ=0.001;

h=0.1;

h=0.1;

=100

 k=2;

 A=2; m/kA=b=0.25( хорошая адаптация хищника)=30;

 

 m=1;

=30;

300 10

Приход к точке равновесия Е1-увеличение численности хищников,уменьшение численности жертв

N=230;M=23

N=140;M=24

N=135;M=38

Приход к точке равновесия Е1-увеличение численности хищников,уменьшение численности жертв

2.(N сonst, M const одновременно)

Высокий радиус индивидуальной активности хищников

При таких начальных функциях  обращаются в ноль, и задача сводится к нераспределенному случаю.

Значения параметров:

D=0.001;         

D=2;

τ=0.001;

h=0.1;

h=0.1;

=100

k=2;

 A=2;

 m=4;

=30;

Начальные функции:

Численность жертв                             Численность хищников

3.Диффузионная неустойчивость. (N сonst, M const одновременно)

Низкий  радиус индивидуальной активности хищников

Значения параметров:

D=0.001;       

D=0.001;

τ=0.001;

h=0.1;

h=0.1;

=35

 k=5

 A=4

 b=3

=45

Начальные функции:           

Численность жертв                                          Численность хищников

Динамика развития системы:

Вымирание хищников

Заключение

В работе исследованы возможности возникновения диссипативных структур в различных распределенных вариантах модели “хищник – жертва”, для одной из них построена численная схема решения, приведены  графические итоги этого решения, а результаты для сосредоточенной модели, полученные в предыдущей курсовой работе, использованы в качестве тестовой задачи.

Полученные результаты позволяют еще раз подтвердить гипотезу о том, что в явлениях живой природы, которые так или иначе можно трактовать как диссипативные структуры, наряду с локальными взаимодействиями, такую же, а может быть и более важную роль играет форма границы. Именно она порождает диссипативные структуры. А любая пространственно распределенная  модель сводится к сосредоточенной, если функции численностей либо не зависят от координат, либо зависимость носит линейный характер.

Таким образом, цель работы выполнена.

Список литературы

1.Свирежев, Ю.М. Нелинейные волны, диссипативные структуры и катастрофы в экологии / Ю.М. Свирежев. – М.: Наука. Главная редакция физико-математической литературы,1987. – 144-228.

2.Вольтерра, В. Математическая теория борьбы за существование  / Пер. с французского под ред. Ю.М.Свирежева. – М.: Наука, 1976. –286-288.

3.Марчук, Г.И. Математическое моделирование в проблеме окружающей среды/ Г.И. Марчук. - М.: Наука, 1982. – 320-326.

4.Березин, И.С. Методы вычислений/ И.С. Березин, Н.П. Жидков. – М.:Государственное издательство физико-математической литературы,1959. – 516-521.

5.Самарский, А.А. Численные методы/А.А. Самарский, А.В. Гулин. – М.: Наука. Главная редакция физико-математической литературы,1989. – 369-378.

6. Свирежев, Ю.М.  Устойчивость биологических сообществ/ Ю.М. Свирежев, Д.О.Логофет - М.: Наука. Главная редакция физико-математической литературы,1978. – 352-354.

Приложение

Код программы

r=15;//размер ареала распространения особей

N1=[200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200; 200,200,200,200,200,200,200,200,200,200,200,200,200,200,200;];

//массив 15*15.Элемент N1[i,j]-количественный показатель концентрации особей

жертв в точке рассматриваемого ареала с координатами (i*h1;j*h2)

M1=[10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10; 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10;];

//массив 15*15.Элемент M1[i,j]-количественный показатель концентрации особей

хищников в точке рассматриваемого ареала с координатами (i*h1;j*h2)

d1=0.001;//коэффициент диффузии жертв

d2=0.001;//коэффициент диффузии хищников

dt=0.001;// шаг по времени

h1=0.1;//пространственный шаг по координате х

h2=0.1;//пространственный шаг по координате у

K=10;//количество рассматриваемых временных слоев

m0=2;//константный параметр популяции

al0=1;//константный параметр популяции

for k=1:K//цикл по времени

for i=1:r //цикл по координате х

   for j=1:r//внутренний цикл по У

       N2(i,j)=0;//изначально задаем в массивах нулевые значения

       M2(i,j)=0;

   end

end

 

for j=1:r//Ищем численность жертв на первом полушаге по времени. Одномерная прогонка по направлению х(т.е. по i при каждом фиксированном j)

  A(1)=(d1/(h1*h1))/(1/dt+d1/(h1*h1)); //подсчет первой группы  прогоночных коэффициентов

   B(1)=(N1(1,j)*(al0*N1(1,j)/(1+N1(1,j))-M1(1,j))+N1(1,j)/dt+(N1(j,2)-N1(j,1))/(h2*h2))/(1/dt+d1/(h1*h1));

   for i=2:r-1//подсчет прогоночных коэффициентов со второй группы по предпоследнюю

       A(i)=(d1/(h1*h1))/(1/dt+2*d1/(h1*h1)-(d1/(h1*h1))*A(i-1));

       B(i)=(N1(i,j)*(al0*N1(i,j)/(1+N1(i,j))-M1(i,j))+N1(i,j)/dt+(N1(j,i+1)-2*N1(j,i)+N1(j,i-1))/(h2*h2)+(d1/(h1*h1))*B(i-1))/(1/dt+2*d1/(h1*h1)-(d1/(h1*h1))*A(i-1));

   end

   A(r)=0; //подсчет последней группы  прогоночных коэффициентов

  B(r)=(N1(r,j)*(al0*N1(r,j)/(1+N1(r,j))-M1(r,j))+N1(r,j)/dt+(N1(j,r)-N1(j,r-1))/(h2*h2)+(d1/(h1*h1))*B(r-1))/(1/dt+d1/(h1*h1)-(d1/(h1*h1))*A(r-1));

   N2(r,j)=B(r);

   for i=r-1:-1:1//обратный пересчет

       N2(i,j)=A(i)*N2(i+1,j)+B(i);//значения матрицы N-численности жертв на полушаге по времени

   end

end    

 

for j=1:r//аналогично предыдущему для матрицы М-численности хищников. Ищем численность хищников на полушаге по времени. Одномерная прогонка по направлению х(т.е. по i при каждом фиксированном j).Первый полушаг

   A(1)=(d2/(h2*h2))/(1/dt+d2/(h2*h2)); //подсчет первой группы  прогоночных коэффициентов

   B(1)=(M1(1,j)*(N1(1,j)-m0*M1(1,j))+M1(1,j)/dt+(M1(j,2)-M1(j,1))/(h2*h2))/(1/dt+d2/(h2*h2));

   for i=2:r-1//подсчет прогоночных коэффициентов со второй группы по предпоследнюю

       A(i)=(d2/(h2*h2))/(1/dt+2*d2/(h2*h2)-(d2/(h2*h2))*A(i-1));

       B(i)=(M1(i,j)*(N1(i,j)-m0*M1(i,j))+M1(i,j)/dt+((M1(j,i+1)-2*M1(j,i)+M1(j,i-1))/(h2*h2)+(d2/(h2*h2))*B(i-1)))/(1/dt+(h2*h2)-(d2/(h2*h2))*A(i-1));

   end

   A(r)=0; //подсчет последней группы  прогоночных коэффициентов

   B(r)=(M1(r,j)*(N1(r,j)-m0*M1(r,j))+M1(r,j)/dt+((M1(j,r-1)-M1(j,r))/(h2*h2)+(d2/(h2*h2))*B(r-1)))/(1/dt+(h2*h2)-(d2/(h2*h2))*A(r-1));

   M2(r,j)=B(r);

   for i=r-1:-1:1//обратный пересчет

       M2(i,j)=A(i)*M2(i+1,j)+B(i);

   end

end

 

 

 

   N1=N2;M1=M2;//заполняем матрицы полученными значениями. Для второго полушага эти значения будут исходными

   

   

 

for j=1:r//Ищем численность жертв на втором полушаге по времени. Одномерная прогонка по направлению у(т.е. по j при каждом фиксированном i)    A(1)=(d2/(h2*h2))/(1/dt+d2/(h2*h2)); //подсчет первой группы  прогоночных коэффициентов

   B(1)=(N1(j,1)*(al0*N1(j,1)/(1+N1(j,1))-M1(j,1))+N1(j,1)/dt+(N1(2,j)-N1(1,j))/(h1*h1))/(1/dt+d2/(h2*h2));

   for i=2:r-1//подсчет прогоночных коэффициентов со второй группы по предпоследнюю

       A(i)=(d2/(h2*h2))/(1/dt+2*d2/(h2*h2)-(d2/(h2*h2))*A(i-1));

       B(i)=(N1(j,i)*(al0*N1(j,i)/(1+N1(j,i))-M1(j,i))+N1(j,i)/dt+(N1(i+1,j)-2*N1(i,j)+N1(i-1,j))/(h1*h1)+(d2/(h2*h2))*B(i-1))/(1/dt+2*d2/(h2*h2)-(d2/(h2*h2))*A(i-1));

   end

   A(r)=0; //подсчет последней группы  прогоночных коэффициентов

   B(r)=(N1(j,r)*(al0*N1(j,r)/(1+N1(j,r))-M1(j,r))+N1(j,r)/dt+(N1(r,j)-N1(r-1,j))/(h1*h1)+(d2/(h2*h2))*B(r-1))/(1/dt+d2/(h2*h2)-(d2/(h2*h2))*A(r-1));

   N2(j,r)=B(r);

   for i=r-1:-1:1//обратный пересчет

       N2(j,i)=A(i)*N2(j,i+1)+B(i);

   end

end

 

for j=1:r //аналогично предыдущему для матрицы М-численности хищников. Ищем численность хищников на полушаге по времени. Одномерная прогонка по направлению y(т.е. по j при каждом фиксированном i).Второй полушаг

 A(1)=(d1/(h1*h1))/(1/dt+d1/(h1*h1)); //подсчет первой группы  прогоночных коэффициентов

   B(1)=(M1(j,1)*(N1(j,1)-m0*M1(j,1))+M1(j,1)/dt+(M1(2,j)-M1(1,j))/(h1*h1))/(1/dt+d1/(h1*h1));

   for i=2:r-1//подсчет прогоночных коэффициентов со второй группы по предпоследнюю

       A(i)=(d1/(h1*h1))/(1/dt+2*d1/(h1*h1)-(d1/(h1*h1))*A(i-1));

       B(i)=(M1(j,i)*(N1(j,i)-m0*M1(j,i))+M1(j,i)/dt+((M1(i+1,j)-2*M1(i,j)+M1(i-1,j))/(h1*h1)+(d1/(h1*h1))*B(i-1)))/(1/dt+(h1*h1)-(d1/(h1*h1))*A(i-1));

   end

   A(r)=0; //подсчет последней группы  прогоночных коэффициентов

   B(r)=(M1(j,r)*(N1(j,r)-m0*M1(j,r))+M1(j,r)/dt+((M1(r-1,j)-M1(r,j))/(h1*h1)+(d1/(h1*h1))*B(r-1)))/(1/dt+(h1*h1)-(d1/(h1*h1))*A(r-1));

   M2(j,r)=B(r);

   for i=r-1:-1:1//обратный пересчет

       M2(j,i)=A(i)*M2(j,i+1)+B(i);

   end

end

   

   N1=N2;M1=M2;//значения численностей жертв и хищников на следующем шаге по времени

   

end

 

[x,y]=meshgrid(0:h1:(r-1)*h1,0:h2:(r-1)*h2);//рисуем сетку

mesh(x,y,N1)//ставим значения численности жертв N1 в точке (х,у)

figure;//соединяем по точкам

mesh(x,y,M1) //ставим значения численности жертв хищников М1 в точке (х,у)



 

А также другие работы, которые могут Вас заинтересовать

82392. Типичная русская еда 49.5 KB
  Задачи направленные на достижение метапредметных результатов: регулятивные УУД формирование умения принимать и сохранять учебные задачи; формирование умения планировать контролировать и оценивать свою деятельность деятельность группы; познавательные УУД формирование умения использовать различные...
82393. Односоставные предложения 87.5 KB
  Обобщить изученный материал об односоставных предложениях. В чём же преимущество увиденной рекламы Чем она больше всего запоминается Простотой изложения лаконичностью легкостью ненавязчивостью Какие предложения использованы в этом рекламном ролике Односоставные В этой рекламе прозвучали стихи...
82394. Н.В.Гоголь «Мертвые души». История замысла, жанр и композиция поэмы, логика последовательности ее глав 38.91 KB
  Цели: помочь учащимся разобраться в замысле, композиции, смысле названия, жанре произведения Н.В.Гоголя; создать условия, способствующие развитию культуры монологической речи, умения обобщать, сопоставлять и делать выводы; воспитывать духовно-нравственную личность, адаптированную к современным...
82395. Неравенства. Конспект урока математики в 8 классе 684 KB
  Цель урока: развивать умение обобщать правильно отбирать способы решения неравенств; расширять общий кругозор. Задачи урока: Образовательные: Расширить обобщить и систематизировать знания о линейных неравенствах; Повторить понятие неравенства алгоритм решения неравенства с одной переменной...
82396. Герб. Флаг. Гимн 32.05 KB
  Цель: формирование представлений о происхождении и истории символики государства России; формирование устойчивой учебно-познавательной мотивации к предмету. Оборудование: компьютер физическая карта России мультимедиапроектор энциклопедия книги о символике.
82397. Волны в океане 57.5 KB
  Цели и задачи: объяснить образование ветровых волн цунами приливов; Описание строения ветровых волн; Практическое использование приливной волны; Совершенствовать умение работы с картами атласа. Приходилось ли вам наблюдать волны на поверхности реки моря Что вы можете о них вспомнить?
82398. А.И. Куприн «Барбос и Жулька» 44.55 KB
  Так о каком же животном пойдёт речь сегодня Какие породы собак вам известны И тема нашего урока А. На свете множество собак И на цепи и просто так: Собак служебных – пограничных Дворовых шариков обычных И молодых пугливых шавок Что тявкать любят из подлавок И тех изнеженных болонок Чей нос курнос...
82399. Царь и кузнец (притча) 100.7 KB
  Это воистину волшебные создания С давних времен бабочки ассоциировалось с легкостью и беззаботностью. Среди скандинавов и славян бабочки олицетворяли душу человека а также служили символом всех влюбленных. С древних времен бабочки вызывали восхищение и являлись символом просветления...
82400. ПРИЕМЫ И ПРИНЦИПЫ АНАЛИЗА ХУДОЖЕСТВЕННОГО ПРОИЗВЕДЕНИЯ 74 KB
  Теория искусства при характеристике идейного содержания произведения в первую очередь выявляет авторское понимание объяснение жизни выраженное в произведении и авторский приговор оценку ее не противопоставляя идейную и художественную стороны.