Примеры уравнений с дискретными переменными

В.Э. Адлер · Классические интегрируемые системы · Лекция 10 · 12 апреля 2021

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

Разделы открываются/закрываются щелчком по скобке справа.
Исполняемые ячейки выделены серым, они активируются нажатием [Shift-Enter].
Генерируемые ячейки с анимацией лучше удалять после просмотра, чтобы они не мешали следующим вычислениям.

Цепочка Вольтерры (Volterra lattice)

(Subscript[u, n,t])=(Subscript[u, n]) ((Subscript[u, n+1])-(Subscript[u, n-1])),  n∈ℤ, t∈ℝ

Эволюция по t

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

Мы выберем начальное условие в виде плавного горбика и посмотрим, как он разваливается на несколько солитонов. Картина очень напоминает солитоны КдФ, что не удивительно, так как цепочка Вольтерры переходит в КдФ в непрерывном пределе.

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

(* задаём период и время счёта; можно менять *)
M=300;
t1=1000;

(* уравнения для n=0,...,M-1, с учётом периодического замыкания *)
eqs=Table[u[n]'[t]==u[n][t](u[Mod[n+1,M]][t]-u[Mod[n-1,M]][t]),{n,0,M-1}];

(* задаём начальные условия; можно менять *)

ic=Table[{u[n][0]==1+0.5Exp[-0.01(((n-M/2))^2)]},{n,0,M-1}];

(* строим решение при помощи NDSolve *)
sol=NDSolve[Join[eqs,ic],Table[u[n],{n,0,M-1}],{t,0,t1},MaxSteps->∞][[1]];

(* массив точек *)
u1[t_]:=Table[{n,u[n][t]},{n,0,M-1}]/.sol

(* функция для вывода. Для наглядности, массив выводится два раза: как ломаная линия 
и как набор жирных точек в ее вершинах *)
gr1[t_]:=ListPlot[{u1[t],u1[t]},
PlotRange->{{-0.1,M+1},{0.95,2.05}},
PlotStyle->{{Thin,Black},{Blue,PointSize[0.01]}},
Joined->{True,False},
AxesLabel->{" \!\(\*SubscriptBox[\(n\), \(\\ \)]\)","(Subscript[u, n])"},
BaseStyle->{FontSize->14,FontFamily->"Times New Roman"},
ImageSize->400,
Epilog->{Text["\!\(\*SubscriptBox[\(t\), \(\\ \)]\)= "<>ToString[t],{0.8M,1.95},{-1,0}]}
]

(* анимация *)
Manipulate[gr1[t],{{t,0},0,t1,Appearance->"Open"}]

Эволюция по n

Можно трактовать цепочку иначе: как отображение ((Subscript[u, n-1]),(Subscript[u, n]))↦((Subscript[u, n]),(Subscript[u, n+1])), определяемое формулой

(1)	(Subscript[u, n+1])=(Subscript[u, n,t])/(Subscript[u, n])+(Subscript[u, n-1]).

Это такое рекуррентное соотношение, которое позволяет вычислять (Subscript[u, n]) по заданным начальным данным (Subscript[u, 0]),(Subscript[u, 1]).
У этого отображения есть интересное свойство: если пара (u,v)=((Subscript[u, n-1]),(Subscript[u, n])) удовлетворяет уравнению в частных производных

(2)	(Subscript[u, s])=-(Subscript[u, t t])+(Subscript[((u^2)+2u v), t]),   (Subscript[v, s])=(Subscript[v, t t])+(Subscript[((v^2)+2u v), t])

(здесь s дополнительная переменная, а t теперь играет роль пространственной переменной в (2)), то это же верно и для пары  (u,v)=((Subscript[u, n]),(Subscript[u, n+1])). Правда, не гарантируется, что это решение будет регулярным, даже если исходное было.

Задача. Докажите это свойство в общем виде (прямым вычислением).

Пример. Пусть u=0, тогда v удовлетворяет уравнению Бюргерса, решение которого можно получить из решения уравнения теплопроводности (Subscript[w, s])=(Subscript[w, t t]) по формуле v=(Subscript[w, t])/w. Например, годится w= 6 s t+2a s+(t^3)+a (t^2)+b t+c.  Строим по этим начальным условиям решение цепочки:

w=c+2 a s+b t+6 s t+a (t^2)+(t^3)
D[w,s]-D[w,t,t]

u[0]=0;
u[1]=((D[w,t])/w);

Table[u[n+1]=Together[D[u[n],t]/u[n]+u[n-1]],{n,1,8}]

u[7]

Хм... Неожиданно, через несколько шагов мы пришли опять в 0 и дальше считать не можем! Это, конечно, благодаря специальному выбору начального условия (можете попробовать взять другое). Проверяем выполнение уравнения (2) для того, что получилось:

Table[Together[{
D[u[n],s]+D[u[n],t,t]-D[((u[n])^2)+2u[n]u[n+1],t],
D[-u[n+1],s]+D[u[n+1],t,t]+D[((u[n+1])^2)+2u[n]u[n+1],t]}],
{n,0,6}]

Можно обобщить наше случайное наблюдение. Похоже, что если взять (Subscript[u, 0])=0 и (Subscript[u, 1])=(Subscript[w, t])/w, где w – многочлен по t степени k, то вычисления по формуле (1) приводят через несколько шагов к (Subscript[u, 2k+1])=0.

u[0]=0;

w=Sum[a[i](t^i),{i,0,4}]

u[1]=((D[w,t])/w);

Table[u[n+1]=Together[D[u[n],t]/u[n]+u[n-1]],{n,1,8}]

u[9]

u[0]=0;

w=(t^5)+1;

u[1]=((D[w,t])/w);

Table[u[n+1]=Together[D[u[n],t]/u[n]+u[n-1]],{n,1,10}]

u[11]

Задача (сложная). Попробуйте доказать это.

Разностная схема Забуски–Краскала (Zabusky-Kruskal scheme)

Термин “солитон” был введен в статье Забуски-Краскала 1965, в которой уравнение КдФ  (Subscript[u, t])=(Subscript[u, 3])+6u (Subscript[u, 1]) исследовалось численно при помощи следующей аппроксимации:

(Subscript[u, n])(1)-(Subscript[u, n])(-1)=a ((Subscript[u, n+2])-2(Subscript[u, n+1])+2(Subscript[u, n-1])-(Subscript[u, n-2]))+2a (h^2)((Subscript[u, n+1])+(Subscript[u, n])+(Subscript[u, n-1]))((Subscript[u, n+1])-(Subscript[u, n-1])). 	 ( ({{ ,  , ●,  ,  }, {●, ●, ●, ●, ●}, { ,  , ●,  ,  }}) )

Это чисто дискретное уравнение на переменные (Subscript[u, n,m]), n,m ∈ℤ  (для краткости, обозначено (Subscript[u, n])=(Subscript[u, n,m]),(Subscript[u, n])(±1)=(Subscript[u, n,m±1])). Само это уравнение не интегрируемо, но оно переходит в КдФ в пределе (Subscript[u, n,m])=u(x,t), x=n h, t=a m (h^3).

Требуется построить периодическое по n решение разностной задачи по следующим данным:

– начальное условие u(x,(Subscript[t, 0]))=f(x) на отрезке [(Subscript[x, min]),(Subscript[x, max])] такое, что f((Subscript[x, min]))=f((Subscript[x, max]));  
	– число N точек (Subscript[x, n]) сетки по x; считаем, что (Subscript[x, 0])=(Subscript[x, min]), (Subscript[x, N])=(Subscript[x, max]), 
	   при этом достаточно работать с n=0,...,N-1;
	– параметр a; 
	– временной интервал [(Subscript[t, min]),(Subscript[t, max])] (включающий (Subscript[t, 0])) и число M точек (Subscript[t, m]), 
	   в которых сохраняется решение (считая с (Subscript[t, 0])).

Шаг по x равен h=((Subscript[x, max])-(Subscript[x, min]))/N. Шаг по t равен a (h^3). Следовательно, чем больше a, тем выше скорость счета, но сильно увеличить этот параметр не удается, схема ломается.

Функции:

ZK[f,{(Subscript[x, min]),(Subscript[x, max]),N},a,{(Subscript[t, min]),(Subscript[t, 0]),(Subscript[t, max]),M}] возвращает решение в виде списка вида {{(Subscript[x, min]),h,a,(Subscript[t, 0]),(Subscript[M, m]),(Subscript[M, 0])},массив M×N}. Здесь (Subscript[M, 0]) – число пропускаемых шагов по времени, то есть, полное число шагов равно (Subscript[M, tot])=(M-1) (Subscript[M, 0])+1; число шагов для положительных и отрицательных t равны (Subscript[M, p])(Subscript[M, 0]) и (Subscript[M, m])(Subscript[M, 0]), соответственно.  Основные вычисления делаются в компилированной функции ZK0.

gr[sol,j,{(Subscript[u, min]),(Subscript[u, max])}] возвращает график j-го кадра в решении; также выводится соответствующее значение t=(Subscript[t, 0])+a (h^3)(j-1-(Subscript[M, m])) (Subscript[M, 0]).

gr3[sol,{(Subscript[u, min]),(Subscript[u, max])},options] возвращает трехмерный график решения.

ZK[f_,{xmin_,xmax_,Nx_},a_,{tmin_,t0_,tmax_,M_}]:=
Module[{n,
h=N[(xmax-xmin)/Nx],
Mtot,Mp,Mm,M0},
Mtot=Round[(tmax-tmin)/(a (h^3))]+1;
M0=Round[(Mtot-1)/(M-1)];
Mp=Round[(tmax-t0)/(a (h^3)M0)];
Mm=M-Mp-1;
{{xmin,a,h,t0,Mm,M0},
ZK0[a,h,Mp,Mm,M0,Table[f/.{x->xmin+n h},{n,0,Nx-1}]]}]

ZK0=Compile[
{{a,_Real},{h,_Real},
{Mp,_Integer},{Mm,_Integer},{M0,_Integer},
{F,_Real,1}},
Module[{U0=F,U=F,sol={F},U1,m0=0},
Do[U1=U0+2a (h^2)(RotateLeft[U,1]+U+RotateLeft[U,-1])(RotateLeft[U,1]-RotateLeft[U,-1])+a(RotateLeft[U,2]-2RotateLeft[U,1]+2RotateLeft[U,-1]-RotateLeft[U,-2]);
U0=U;
U=U1;
m0=m0+1;
If[m0==M0,m0=0;AppendTo[sol,U]],
{Mp*M0}];
U0=sol[[2]];
U=sol[[1]];
m0=0;
Do[U1=U0-2a (h^2)(RotateLeft[U,1]+U+RotateLeft[U,-1])(RotateLeft[U,1]-RotateLeft[U,-1])-a(RotateLeft[U,2]-2RotateLeft[U,1]+2RotateLeft[U,-1]-RotateLeft[U,-2]);
U0=U;
U=U1;
m0=m0+1;
If[m0==M0,m0=0;PrependTo[sol,U]],
{Mm*M0}];
sol]];

gr[sol_,j_,{umin_,umax_},opt___]:=
Module[{xmin,a,h,t0,Mm,M0,U=sol[[2,j]],Nx,n},
{xmin,a,h,t0,Mm,M0}=sol[[1]];
Nx=Length[U];
U=Transpose[{Table[xmin+n h,{n,0,Nx-1}],U}];
ListPlot[U,
Joined->True,
PlotRange->{{xmin,xmin+Nx h},{umin,umax}},
AxesLabel->Evaluate[Style[#,FontFamily->"Times",Italic,16]&/@{"x","u"}],
Epilog->{Text["\!\(\*SubscriptBox[\(t\), \(\\ \)]\)="<>ToString[Round[t0+a (h^3)(j-1-Mm) M0,0.01]],{xmin+0.85Nx h,umin+0.9(umax-umin)},{-1,0}]},
opt]]

gr3[sol_,{umin_,umax_},opt___]:=
Module[{xmin,a,h,t0,Mm,M0},{xmin,a,h,t0,Mm,M0}=sol[[1]];
Show[
{ListPlot3D[sol[[2]],
DataRange->{{xmin,xmin+h Length[sol[[2,1]]]},{t0-a (h^3)Mm M0,t0+a (h^3)(Length[sol[[2]]]-Mm-1) M0}},
PlotRange->{umin,umax},
Mesh->None,
AxesLabel->Evaluate[Style[#,Italic,24]&/@{"x","t","u"}],
ImageSize->400,
BaseStyle->{FontFamily->"Times",FontSize->16},
opt],
Graphics3D[{Thick,Red,Line[Table[{xmin+h n,t0,sol[[2,Mm+1,n]]},{n,1,Length[sol[[2,1]]]}]]}]
}]]

. Солитон

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

Clear[u]
u[x_,t_,d_,k_]:=((2(k^2))/((Cosh[k x+4(k^3)t+d])^2))
Simplify[D[u[x,t,d,k],t]-D[u[x,t,d,k],x,x,x]-6u[x,t,d,k] D[u[x,t,d,k],x]]

n=600;
a=0.01;
t1=20;
Round[t1/(a (((100/n))^3))]
Timing[sol1=ZK[u[x,0,-40,1],{-50,50,n},a,{0,0,t1,100}];]

gg[j_]:=Show[{
Plot[u[x,sol1[[1,4]]+(j-1) sol1[[1,2]]((sol1[[1,3]])^3)sol1[[1,6]],-40,1],{x,-50,50},
PlotStyle->Red,PlotRange->{-0.15,2.15}],
gr[sol1,j,{-0.15,2.15}]}]

Наложим графики численного решения и точного (красное).

Manipulate[gg[j],
{{j,sol1[[1,5]]+1},1,Length[sol1[[2]]],1,Appearance->"Labeled"}]

. Начальные данные Забуски–Краскала

Попытка воспроизвести результат из статьи 1965 года. До некоторого t это удается, но дальше что-то ломается. Возможно, нужно выбирать a меньше. Из статьи это трудно понять, схема там описана довольно подробно, но конкретные числа не указаны. По сравнению со статьей сделано некоторое масштабирование.

c=0.15;
n=400;
a=0.005;
t1=20;
Round[t1/(a (((2π/c/n))^3))]

Timing[sol2=ZK[Cos[c x],{0,((2π)/c),n},a,{0,0,12,800}];]
Length[sol2[[2]]]

gr3[sol2,{-2,3},ViewPoint->{0,-3,4}]

Manipulate[gr[sol2,j,{-1.55,2.85},
PlotStyle->Blue,
ImageSize->400,
BaseStyle->{FontFamily->"Times",FontSize->14}
],
{{j,sol2[[1,5]]+1},1,Length[sol2[[2]]],1,Appearance->{"Open","Labeled"}}]

Система шаров-ящиков (Box-ball system)

Имеется бесконечная лента, разбитая на клетки-ящики. В начальный момент в конечном числе клеток размещено по одному шару. Переход к следующему моменту времени заключается в том, что каждый шар по очереди, начиная с крайнего левого и кончая крайним правым, переносится в ближайшую свободную клетку справа от него, при этом каждый шар передвигается только один раз за ход [Takahashi, Satsuma 1990].

Группа из n шаров движется со скоростью n клеток за ход. Допустим, что справа от неё находится группа с m шарами. Понятно, что если m≥n, то обе группы будут двигаться независимо, а если m<n, то в какой-то момент случится столкновение. Что при этом произойдет? Эти две группы могли бы слипнуться, развалиться на другие группы или отдельные шары. Но нет – оказывается, что после некоторого обмена шарами, опять образуется две группы с тем же самым числом шаров, но более быстрая расположена уже справа. В результате взаимодействия происходит фазовый сдвиг. То же самое будет, если взять несколько групп: число шаров в них восстанавливается после взаимодействия. Таким образом, группы шаров ведут себя вполне аналогично солитонам, причем число шаров в группе играет роль амплитуды. Возможно, это самая простая (по формулировке) интегрируемая модель, которую можно вообразить.

Пусть (Subscript[x, n]^t)∈{0,1} обозначает состояние n-й клетки в момент t (0 означает, что клетка пуста) и (∑_(n=-∞)^∞)(Subscript[x, n]^t)<∞ (число шаров конечно). Тогда

(Subscript[x, n]^t)={({{1, если (Subscript[x, n]^t)=0  и  (∑_(k=-∞)^(n-1))((Subscript[x, k]^t)-(Subscript[x, k]^(t+1)))>0, }, {0, иначе.}})

move:=Module[{x1=x,M=Length[x]},For[n=1,n<M+1,n++,x1[[n]]=If[x[[n]]==0&&Sum[x[[k]]-x1[[k]],{k,1,n-1}]>0,1,0]];x=x1]

BB[n_]:=Join[{x},Table[move,{k,1,n}]]

gr[b_,j_,h_]:=Graphics[
Flatten[
Table[If[b[[t,n]]==0,{},Disk[{n,h-t},2/5]],
{t,1,j},
{n,1,Length[b[[1]]]}]],
PlotRange->{{0,Length[b[[1]]]+1},{-1,h}},
ImageSize->600
]

x={1,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
b=BB[23];
MatrixForm[b/.{0->".",1->"●"}]

x={1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
b=BB[80];

Manipulate[gr[b,j,Length[b]],{j,1,Length[b],Appearance->"Labeled"}]