Разностная схема системы реакция-диффузия на диске

@ge4r · Регистрация: 04.12.2008

Студворк — интернет-сервис помощи студентам

Добрый день,
с численными методами знаком плохо, поэтому прошу Вашей помощи.
Есть задача такого вида:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{\partial{u}}{\partial{t}}=\frac{\partial^2{u}}{\partial{r^2}}+\frac{1}{r}\frac{\partial{u}}{\partial{r}}+\frac{1}{r^2}\frac{\partial^2{u}}{\partial{\theta^2}}+u(1-u-cv);\\ \\\frac{\partial{v}}{\partial{t}}=\mu(\frac{\partial^2{v}}{\partial{r^2}}+\frac{1}{r}\frac{\partial{v}}{\partial{r}}+\frac{1}{r^2}\frac{\partial^2{v}}{\partial{\theta^2}})+av(1-cu-bv);\\ \\\\$

С граничными условиями
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\left. \frac{\partial{u(r,\theta,t)}}{\partial{r}} \right |_{r=1}=0 \\\left. \frac{\partial{v(r,\theta,t)}}{\partial{r}} \right |_{r=1}=0 \\\left. \frac{\partial{u(r,\theta,t)}}{\partial{r}} \right |_{r=0}=0 \\\left. \frac{\partial{v(r,\theta,t)}}{\partial{r}} \right |_{r=0}=0 \\$

И начальными кусочно-заданными условиями:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u(r,\theta,0)=\begin{cases}\\0,&\text{if $r \leq r_0$,$0<\theta \leq \Pi/6$}\\1,&\text{if $r \leq r_0$,$\Pi/6<\theta \leq 2\Pi$}\\0,&\text{if $r > r_0$,$0<\theta \leq 2\Pi$}\\\end{cases}\\v(r,\theta,0)=\begin{cases}1,&\text{if $r \leq r_0$,$0<\theta \leq \Pi/6$}\\0,&\text{if $r \leq r_0$,$\Pi/6<\theta \leq 2\Pi$}\\0,&\text{if $r > r_0$,$0<\theta \leq 2\Pi$}\\\end{cases}$

Построил для этой задачи разностную схему, хотелось бы узнать - верно или нет, особенно граничные условия вызывают сомнения.
Итак, сначала задал сетку:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\text{Grid:}\\r_{min}=0;r_{max}=1;n=100; h_r=\frac{r_{max}-r_{min}}{n}\\\theta_{min}=0;\theta_{max}=2\Pi;m=600; h_\theta=\frac{\theta_{max}-\theta_{min}}{m}\\t_{min}=0;t_{max}=1;l=200; h_t=\frac{t_{max}-t_{min}}{l}\\$

Затем сами уравнения:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{u_{i,j}^{t+1}-u_{i,j}^{t}}{h_t}=\frac{u_{i+1,j}^{t}-2u_{i,j}^{t}+u_{i-1,j}^{t}}{h^2_r}+\frac{1}{i \cdot h_r}\frac{u_{i+1,j}^{t}-u_{i,j}^{t}}{h_r}+\\+\frac{1}{(i \cdot h_r)^2} \frac{u_{i,j+1}^{t}-2u_{i,j}^{t}+u_{i,j-1}^{t}}{h_\theta}+u_{i,j}^{t}(1-u_{i,j}^{t}-c \cdot v_{i,j}^{t})\\$ ;

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{v_{i,j}^{t+1}-v_{i,j}^{t}}{h_t}=\mu \cdot(\frac{v_{i+1,j}^{t}-2v_{i,j}^{t}+v_{i-1,j}^{t}}{h^2_r}+\frac{1}{i \cdot h_r}\frac{v_{i+1,j}^{t}-v_{i,j}^{t}}{h_r}+\\+\frac{1}{(i \cdot h_r)^2} \frac{v_{i,j+1}^{t}-2v_{i,j}^{t}+v_{i,j-1}^{t}}{h_\theta})+a \cdot v_{i,j}^{t}(1-c \cdot u_{i,j}^{t}-b \cdot v_{i,j}^{t})\\$ ;

Граничные условия:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{u_{2,j}^{t}-u_{1,j}^{t}}{h_r}=0 \Rightarrow u_{2,j}^{t}=u_{1,j}^{t} \\\frac{u_{n,j}^{t}-u_{n-1,j}^{t}}{h_r}=0 \Rightarrow u_{n,j}^{t}=u_{n-1,j}^{t}\\\\\frac{v_{2,j}^{t}-v_{1,j}^{t}}{h_r}=0 \Rightarrow v_{2,j}^{t}=v_{1,j}^{t}\\\frac{v_{n,j}^{t}-v_{n-1,j}^{t}}{h_r}=0 \Rightarrow v_{n,j}^{t}=v_{n-1,j}^{t}\\\\u_{i,1}^{t}=u_{i,m}^{t}\\v_{i,1}^{t}=v_{i,m}^{t}\\$

Ну и наконец начальные функции:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u(i \cdot h_r,j \cdot h_\theta,0)=\begin{cases}0,&\text {if $i \cdot h_r<r_0,0 \leq j \cdot h_\theta \leq \Pi/6$}\\1,&\text {if $i \cdot h_r<r_0,\Pi/6 < j \cdot h_\theta < 2\Pi$}\\0,&\text {if $i \cdot h_r>r_0,0 \leq j \cdot h_\theta \leq 2\Pi$}\\\end{cases}$ ;
$https://www.cyberforum.ru/cgi-bin/latex.cgi?v(i \cdot h_r,j \cdot h_\theta,0)=\begin{cases}1,&\text {if $i \cdot h_r<r_0,0 \leq j \cdot h_\theta \leq \Pi/6$}\\0,&\text {if $i \cdot h_r<r_0,\Pi/6 < j \cdot h_\theta < 2\Pi$}\\0,&\text {if $i \cdot h_r>r_0,0 \leq j \cdot h_\theta \leq 2\Pi$}\\\end{cases}$ ;

И вот тут возникает вопрос о дальнейшем задании этой схемы в каком-либо матпакете.
Правильно ли я понимаю ход действий:
1. Используя начальные условия, заполняем первый слой при t=0;
2. Используя граничные условия для производных получаем, что внешнии слои по r (первый и последний) равны пред-внешним (второй и предпоследний), но какие там будут значения лежать?
3. Используя граничные условия периодичности получим что крайние слои по тета тоже равны, но опять не понятно какие там значения.
4. Ну и наконец в системе выражаем u[i,j,t+1]=..., v[i,j,t+1]=... и запускам цикл по 3 переменным. И вот тут еще один вопрос. В каком порядке запускать циклы? Полагаю, что самый внешний будет по t, а внутренние по i и j.Но в каком порядке? Сначала по i, потом по j или наоборот?

Заранее спасибо за ответы.

Добавлено через 22 часа 22 минуты
Составил алгоритм расчета задачи. Верно ли организовал циклы, не подскажете?

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
clc;
clear all;
 
%Задание сетки
.
.
.
 
%Заполняем первый временной слой
for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           u0=0; 
        elseif ((r(i)<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           u0=1;
        elseif ((r(i)>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           u0=0; 
        end
        u(i,j,1)=u0;
    end
end
 
for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           v0=1; 
        elseif ((r(i)<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           v0=0;
        elseif ((r(i)>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           v0=0; 
        end
        v(i,j,1)=v0;
    end
end
 
%Сама схема
for t=1:nt-1 
    for i=2:nr-1
        for j=2:nth-1
        u(i,j,t+1)=...;
        v(i,j,t+1)=...;
        %Граничные условия Неймана по r
        u(2,j,t)=u(1,j,t);
        u(nr,j,t)=u(nr-1,j,t);
        v(2,j,t)=v(1,j,t);
        v(nr,j,t)=v(nr-1,j,t);
        end
    %Граничные условия по theta
    u(i,1,t)=u(i,nth,t);
    v(i,1,t)=v(i,nth,t);
    end
end

Добавлено через 4 часа 21 минуту
Нашел момент, в котором возникает проблема. Условия по тета имеют вид:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u(r,0,t)=u(r,2\Pi,t)\\v(r,0,t)=v(r,2\Pi,t)\\$
Если задаю их в разностной схеме как
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u_{i,1}^{t}=u_{i,n}^{t}\\v_{i,1}^{t}=v_{i,m}^{t}\\$
то получается, что они никак не опираются на другие подсчитанные значения и остаются неопределенными. Как быть в такой ситуации?

@ge4r · 08.11.2013, 18:44 **[ТС]**

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

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
clc;
clear all;
 
%constants
r0=0.1;
a=1.01;
b=1;
mu=1;
c=100;
 
%step-size for r
n=100;
r_max=1;
hr=r_max/n;
r=0:hr:r_max;
nr=max(size(r));
%step-size for theta
m=200;
th_max=2*pi;
hth=th_max/m;
th=0:hth:th_max;
nth=max(size(th));
%step-size for t
l=100;
t_max=5;
ht=t_max/l;
time=0:ht:t_max;
nt=max(size(time));
 
u=zeros(nr,nth,nt);
v=zeros(nr,nth,nt);
 
%Initialization
for i=1:nr
    for j=1:nth
        if ((i*hr<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           u0=0; 
        elseif ((i*hr<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           u0=1;
        elseif ((i*hr>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           u0=0; 
        end
        u(i,j,1)=u0;
    end
end
 
for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           v0=1; 
        elseif ((r(i)<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           v0=0;
        elseif ((r(i)>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           v0=0; 
        end
        v(i,j,1)=v0;
    end
end
 
for t=1:nt-1
    
   for i=2:nr-1
       u(i,1,t+1)=u(i,1,t)+ht*((u(i+1,1,t)-2*u(i,1,t)+u(i-1,1,t))/hr^2+(1/(i*hr)))*((u(i+1,1,t)-u(i,1,t))/hr)+(1/(i*hr)^2)*((u(i,2,t)-2*u(i,1,t)+u(i,nth,t)/hth)+u(i,1,t)*(1-u(i,1,t)-c*v(i,1,t)));
       v(i,1,t+1)=v(i,1,t)+ht*((v(i+1,1,t)-2*v(i,1,t)+v(i-1,1,t))/hr^2+(1/(i*hr)))*((v(i+1,1,t)-v(i,1,t))/hr)+(1/(i*hr)^2)*((v(i,2,t)-2*v(i,1,t)+v(i,nth,t)/hth)+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));
       u(i,nth,t+1)=u(i,nth,t)+ht*((u(i+1,nth,t)-2*u(i,nth,t)+u(i-1,nth,t))/hr^2+(1/(i*hr)))*((u(i+1,nth,t)-u(i,nth,t))/hr)+(1/(i*hr)^2)*((u(i,1,t)-2*u(i,nth,t)+u(i,nth-1,t)/hth)+u(i,nth,t)*(1-u(i,nth,t)-c*v(i,nth,t)));
       v(i,nth,t+1)=v(i,nth,t)+ht*((v(i+1,nth,t)-2*v(i,nth,t)+v(i-1,nth,t))/hr^2+(1/(i*hr)))*((v(i+1,nth,t)-v(i,nth,t))/hr)+(1/(i*hr)^2)*((v(i,1,t)-2*v(i,nth,t)+v(i,nth-1,t)/hth)+a*v(i,nth,t)*(1-c*u(i,nth,t)-b*v(i,nth,t)));
       for j=2:nth-1
        u(i,j,t+1)=u(i,j,t)+ht*((u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/hr^2+(1/(i*hr)))*((u(i+1,j,t)-u(i,j,t))/hr)+(1/(i*hr)^2)*((u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t)/hth)+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));
        v(i,j,t+1)=v(i,j,t)+ht*((v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/hr^2+(1/(i*hr)))*((v(i+1,j,t)-v(i,j,t))/hr)+(1/(i*hr)^2)*((v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t)/hth)+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));
       end
   end
   %Boundary
   for j=1:nth
        u(1,j,t)=u(2,j,t);
        u(nr,j,t)=u(nr-1,j,t);
        v(1,j,t)=v(2,j,t);
        v(nr,j,t)=v(nr-1,j,t);
   end  
end

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11680&d=1772460536 Одним из. . .	Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .
SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .	SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .	Конвертировать закладки radiotray-ng в m3u-плейлист damix 19.02.2026 Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .	Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .