Модель "хищник - жертва": конечно-разностная схема.

[byulent]

Цитата:
"Динамика изменения численности жертвы $n_1$ и хищника $n_2$ задаётся уравнениями:

$$\frac{\partial{n_1}}{\partial{t}}=b(n_1)n_1-n_1n_2+D_1\frac{\partial^2{n_1}}{\partial{x}^2},$$

$$\frac{\partial{n_2}}{\partial{t}}=-d(n_2)n_2+n_1n_2+D_2\frac{\partial^2{n_2}}{\partial{x}^2}$$

... Функцию роста $b(n)$ и смерти $d(n)$ часто аппроксимируют линейным или квадратичным способом. Для конкретности возьмём:

$$b(n_1)=a+bn_1+cn_1^2,$$

$$d(n_2)=d_0+d_1n_2$$

"
Так вот, пытаюсь составить конечно - разностную схему для всего этого. Получается какая-то дичь:

$$\frac{N_{i_1}^{j+1}-N_{i_1}^j}{\tau}=\frac{(aN_{i_1}^{j+1}+b(N_{i_1}^{j+1})^2-c(N_{i_1}^{j+1})^3-N_{i_1}^{j+1}N_{i_2}^{j+1})+(aN_{i_1}^{j}+b(N_{i_1}^{j})^2-c(N_{i_1}^{j})^3-N_{i_1}^{j}N_{i_2}^{j})}{2}+D_1\frac{N_{{i+1}_1}^{j}-2N_{i_1}^{j}+N_{{i-1}_1}^{j}}{h^2},$$

формула не влезла целиком

$$\frac{N_{i_2}^{j+1}-N_{i_2}^j}{\tau}=\frac{(d_0N_{i_2}^{j+1}+d_1(N_{i_2}^{j+1})^2+N_{i_1}^{j+1}N_{i_2}^{j+1})+(d_0N_{i_2}^{j}+d_1(N_{i_2}^{j})^2+N_{i_1}^{j}N_{i_2}^{j})}{2}+D_2\frac{N_{{i+1}_2}^{j}-2N_{i_2}^{j}+N_{{i-1}_2}^{j}}{h^2}$$

Как мне выразить $N_{i_1}^{j+1}$, $N_{i_2}^{j+1}$? И как будет зависеть $\tau$ от $D_1, D_2$ и $h$?

[Munin]

Чтобы формула влезала, её надо вручную разбить на несколько строчек. Например, так:

byulent в сообщении #1082012 писал(а):

$$\begin{gathered}\frac{N_{i_1}^{j+1}-N_{i_1}^j}{\tau}=\frac{(aN_{i_1}^{j+1}+b(N_{i_1}^{j+1})^2-c(N_{i_1}^{j+1})^3-N_{i_1}^{j+1}N_{i_2}^{j+1})+(aN_{i_1}^{j}+b(N_{i_1}^{j})^2-c(N_{i_1}^{j})^3-N_{i_1}^{j}N_{i_2}^{j})}{2}+{}\\{}+D_1\frac{N_{{i+1}_1}^{j}-2N_{i_1}^{j}+N_{{i-1}_1}^{j}}{h^2},\end{gathered}$$

$$\begin{gathered}\frac{N_{i_2}^{j+1}-N_{i_2}^j}{\tau}=\frac{(d_0N_{i_2}^{j+1}+d_1(N_{i_2}^{j+1})^2+N_{i_1}^{j+1}N_{i_2}^{j+1})+(d_0N_{i_2}^{j}+d_1(N_{i_2}^{j})^2+N_{i_1}^{j}N_{i_2}^{j})}{2}+{}\\{}+D_2\frac{N_{{i+1}_2}^{j}-2N_{i_2}^{j}+N_{{i-1}_2}^{j}}{h^2}\end{gathered}$$

[Pphantom]

В первом "диком" выражении перед $c$ должен быть плюс, а не минус. Потом не очень ясно, зачем Вы берете среднее арифметическое значений на шагах $j$ и $j+1$ - это иногда может оказаться полезным, но без определенной потребности это только усложнит дело.

В существующем варианте Вы записали схему неявного метода, соответственно, полученные нелинейные алгебраические уравнения на каждом шаге надо решать относительно $N^{j+1}_i$ (например, методом Ньютона или итерациями). Но лучше бы сначала убедиться, что явного метода (когда для членов без производных все выражается через $N^j_i$) недостаточно, чем городить огород с неявным.

Осталось разобраться с шагом по времени. Если метод явный, то возьмите для начала $\tau \lesssim \frac{h^2}{\sqrt{D_1 \, D_2}}$, а там посмотрим.

P.S. Для понимания последующего: тут исправлена описка в оценке шага по времени (сначала квадрата не было).

[Munin]

Что мне не нравится в постановке задачи с точки зрения её логики - это член

$$\dfrac{\partial n_2}{\partial t}=\ldots+n_1n_2+\ldots$$

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

$$\dfrac{\partial n_2}{\partial t}=\ldots+A\dfrac{n_1}{C+n_1}n_2+\ldots$$

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

[Pphantom]

Munin в сообщении #1082095 писал(а):

Что мне не нравится в постановке задачи с точки зрения её логики - это член

Совсем не факт, что ТС занимается моделированием биологических систем. Модель Лотки-Вольтерры довольно часто встречается и в других задачах, а, судя по предыдущим темам, ТС вообще интересуется ей с точки зрения прикладной математики.

[B@R5uk]

byulent в сообщении #1082012 писал(а):

Получается какая-то дичь

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

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

Если используете высокоуровневую среду программирования типа МАТЛАБ или ещё что-нибудь, то первую часть — интегрирующую процедуру — можно взять одну из множества готовых, поставляемых в комплекте.

[Munin]

Pphantom в сообщении #1082102 писал(а):

Совсем не факт, что ТС занимается моделированием биологических систем.

Не факт, конечно. Скорее всего, это просто задание на ДУЧП. Просто считаю нужным оговорить, что эти уравнения могут быть не очень хорошими в плане соответствия реальности. Такой reality check, такую оглядку полезно время от времени делать.

[Pphantom]

B@R5uk в сообщении #1082109 писал(а):

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

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

[B@R5uk]

Pphantom в сообщении #1082125 писал(а):

но, увы, пока не получается

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

[Pphantom]

B@R5uk в сообщении #1082145 писал(а):

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

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

[B@R5uk]

Pphantom в сообщении #1082152 писал(а):

но труднодостижимо даже в реальной жизни

Как говорил Сократ, кто хочет, найдёт способ, кто не хочет — причину.

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

[byulent]

Сегодня спросил у препода про вот это:

byulent в сообщении #1082012 писал(а):

$$b(n_1)=a+bn_1-cn_1^2,$$

$$d(n_2)=d_0+d_1n_2$$

(да, там минус, а не плюс, как я написал сначала).
Он сказал, что это относится только к элементу $N_{i_{1,2}}^j$. Задача немного упростилась:

$$\begin{gathered}\frac{N_{i_1}^{j+1}-N_{i_1}^j}{\tau}=aN_{i_1}^{j}+b(N_{i_1}^{j})^2-c(N_{i_1}^{j})^3-N_{i_1}^{j}N_{i_2}^{j}+D_1\frac{N_{{i+1}_1}^{j}-2N_{i_1}^{j}+N_{{i-1}_1}^{j}}{h^2},\end{gathered}$$

$$\begin{gathered}\frac{N_{i_2}^{j+1}-N_{i_2}^j}{\tau}=d_0N_{i_2}^{j}+d_1(N_{i_2}^{j})^2+N_{i_1}^{j}N_{i_2}^{j}+D_2\frac{N_{{i+1}_2}^{j}-2N_{i_2}^{j}+N_{{i-1}_2}^{j}}{h^2}\end{gathered}$$

Вопрос про шаг по времени до сих пор остаётся открытым.
Почему он должен выглядеть так:

Pphantom в сообщении #1082079 писал(а):

$\tau \lesssim \frac{h}{\sqrt{D_1 \, D_2}}$

, если конечно-разностная схема берётся из уравнения теплопроводности, а там это условие выглядит совсем по-другому: $\tau\leqslant\dfrac{h^2}{2D}$?

[Pphantom]

byulent в сообщении #1082199 писал(а):

(да, там минус, а не плюс, как я написал сначала).

Ок, тогда с разностной аппроксимацией тут все было в порядке.

byulent в сообщении #1082199 писал(а):

Он сказал, что это относится только к элементу $N_{i_{1,2}}^j$. Задача немного упростилась:

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

byulent в сообщении #1082199 писал(а):

Почему он должен выглядеть так:

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

[byulent]

Испробовал все варианты: $\tau=\dfrac{h^2}{2\max(D_1,D_2)}$, $\tau=\dfrac{h^2}{2(D_1+D_2)}$, ваш вариант - всё без толку. Значения улетают куда-то в бесконечность уже на, по-моему, 20-й итерации.
Попробовать, может, включить туда коэфф-ты $a$, $b$, $c$ и т. д.? Или это от них не зависит совершенно?

[Munin]

byulent в сообщении #1082199 писал(а):

Он сказал, что это относится только к элементу $N_{i_{1,2}}^j$.

Это и есть "явная схема". Как вам тут и сказали.

По условию: полагаю, Pphantom просто описался насчёт квадрата. Он же и из размерностей следует. Очевидно, он там должен быть.

-- 14.12.2015 22:59:18 --

byulent в сообщении #1082204 писал(а):

Испробовал все варианты: $\tau=\dfrac{h^2}{2\max(D_1,D_2)}$, $\tau=\dfrac{h^2}{2(D_1+D_2)}$, ваш вариант - всё без толку. Значения улетают куда-то в бесконечность уже на, по-моему, 20-й итерации.

В принципе, можно ещё коэффициенты типа $1/2,$ $1/10,$ $1/1000,$ попборовать. Но вообще, это скорее признак, что вы где-то в кодировании напортачили.