Цепочка Вольтерры

К лекции 9 (2024)

VL:	(Subscript[u, n,t])=(Subscript[u, n]) ((Subscript[u, n+1])-(Subscript[u, n-1]))

Законы сохранения, симметрии.
Несколько простых численных экспериментов.

Численный счёт, решение солитонного типа

Clear[u]

(* задаём период и время счёта *)
M1=100;
t1=1000;

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

(* задаём начальные условия. Попробуйте что-нибудь поменять *)
ic=Table[{u[n][0]==1+0.5Exp[-0.01(((n-M1/2))^2)]},{n,0,M1-1}];

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

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

gr1[t_]:=ListPlot[{u1[t],u1[t]},
PlotRange->{{-0.1,M1+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.8M1,1.95},{-1,0}]}
]

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

Распад пичка

Что получится, если возмутить постоянное решение в одном-единственном узле?

Clear[u]
M2=500;
t2=1000;

eqs=Table[u[n]'[t]==u[n][t](u[Mod[n+1,M2]][t]-u[Mod[n-1,M2]][t]),{n,0,M2-1}];

ic=Table[{u[n][0]==If[n==Quotient[M2,2],1.5,1]},{n,0,M2-1}];

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

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

gr2[t_]:=ListPlot[{u2[t],u2[t]},
PlotRange->{{-0.1,M2+1},{0.7,1.6}},
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.8M2,1.95},{-1,0}]}
]

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

Игра с чётностью

В цепочке Вольтерра есть стационарное решение (Subscript[u, 2n])=a,  (Subscript[u, 2n+1])=b с разными постоянными a и b. Если немного возмутить его, взяв начальное условие в виде двух плавных кривых, отдельно для чётных и нечётных узлов, то эти узлы так и останутся разделёнными. Конечно, в уравнении КдФ такое невозможно. Для решений такого типа тоже можно написать непрерывный предел, но уже к системе на две полевые переменные.

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

M3=500;
t3=10000;

eqs=Table[u[n]'[t]==u[n][t](u[Mod[n+1,M3]][t]-u[Mod[n-1,M3]][t]),{n,0,M3-1}];

ic=Table[{u[n][0]==
If[OddQ[n],
1+0.1Exp[-0.005(((n-M3/3))^2)]+0.1Exp[-0.005(((n-2M3/3))^2)],   (* формула для нечётных узлов *) 
1.2-0.05Exp[-0.005(((n-M3/2))^2)]]},            (* формула для чётных узлов *)
{n,0,M3-1}];

sol3=NDSolve[Join[eqs,ic],Table[u[n],{n,0,M3-1}],{t,0,t3},MaxSteps->∞][[1]];

u3[t_]:=Table[{n,u[n][t]},{n,0,M3-1}]/.sol3
u3o[t_]:=Table[{n,u[n][t]},{n,1,M3-1,2}]/.sol3
u3e[t_]:=Table[{n,u[n][t]},{n,0,M3-1,2}]/.sol3

gr3[t_]:=ListPlot[{u3[t],u3o[t],u3e[t],u3o[t],u3e[t]},
PlotRange->{{-0.1,M3+1},{0.85,1.55}},
PlotStyle->{{Thin,Black},{Thin,Black},{Thin,Black},
{Blue,PointSize[0.01]},{Red,PointSize[0.01]}},
Joined->{True,True,True,False,False},
AxesLabel->{" \!\(\*SubscriptBox[\(n\), \(\\ \)]\)","(Subscript[u, n])"},
BaseStyle->{FontSize->14,FontFamily->"Times New Roman"},
ImageSize->400,
Epilog->{Text["\!\(\*SubscriptBox[\(t\), \(\\ \)]\)= "<>ToString[t],{0.8M3,1.50},{-1,0}]}
]

Manipulate[gr3[t],{{t,0},0,t3,Appearance->{"Labeled","Open"},AnimationRate->10}]

Преобразование типа Миуры

Определим оператор линеаризации и разностную вариационную производную.

vars[f_]:=Union[Flatten[Cases[f,Blank[#],{0,∞}]&/@{u}]]
df[f_]:=Plus@@(D[f,#]dif[#]&/@vars[f])
star[f_,g_]:=df[f]/.dif[u[k_]]:>(g/.n->k)

vard[f_]:=Expand[Plus@@(T[D[f,#],n-#[[1]]]&/@vars[f])]

T[f_,k_]:=f/.n->n+k

Вариационная производная уничтожает полную разность

vard[u[n]u[n+1]]

vard[u[n]-u[n+1]]

vard[
(T[#,1]-#)&[G[u[n-1],u[n],u[n+1],u[n+2]]]
]

Определим mVL (пишем u вместо f, U вместо u)

ut=u[n](u[n]-a)(u[n+1]-u[n-1])

По решению (Subscript[f, n]) этой цепочки определим переменную (Subscript[u, n])=(Subscript[f, n])((Subscript[f, n+1])-a) и проверим, что она удовлетворяет VL

U[n_]:=u[n](u[n+1]-a)
Factor[star[U[n],ut]-U[n](U[n+1]-U[n-1])]

Точно так же для переменной (Subscript[u, n])=((Subscript[v, n])-a)(Subscript[v, n+1])

U[n_]:=(u[n]-a)u[n+1]
Factor[star[U[n],ut]-U[n](U[n+1]-U[n-1])]

Представление нулевой кривизны

L[n_]:={{0,1},{-u[n],λ}}

U[n_]:={{-((λ )^2)/2+u[n-1]-1,λ},
{-λ u[n],((λ )^2)/2+u[n]-1}}

Factor[star[L[n],u[n](u[n+1]-u[n-1])]-U[n+1].L[n]+L[n].U[n]]

Плотности законов сохранения

Решаем уравнение (дискретное преобразование Миуры)

(Subscript[u, n])=(Subscript[f, n]) ((Subscript[f, n+1])-z)

в виде ряда

(Subscript[f, n])=-((Subscript[u, n])/z) (1+(Subscript[F, n]^((1)))/(z^2)+...)

Плотности определяются, как коэффициенты разложения

log (Subscript[f, n])=-log(- z)+log (Subscript[u, n])+log(1+((Subscript[F, n]^((1)))/(z^2))+((Subscript[F, n]^((2)))/(z^4))+((Subscript[F, n]^((3)))/(z^6)))=-log(- z)+log (Subscript[u, n])-(Subscript[r, n]^((1)))/(z^2)+(Subscript[r, n]^((2)))/(z^4)-(Subscript[r, n]^((3)))/(z^6)+...

F[0,n_]:=1
F[k_,n_]:=-u[n+1]Expand[Sum[F[s,n]F[k-1-s,n+1],{s,0,k-1}]]/;k>0

F[0,n]
F[1,n]
F[2,n]
F[3,n]

В принципе, для r[k] можно выписать рекуррентные формулы, но мы просто воспользуемся функцией Series.

M=6;

FF=Sum[F[s,n](z^(-2s)),{s,0,M}];

R=Series[Log[FF],{z,∞,2(M-1)}];

r[0]=Log[u[n]];

Do[r[k]=T[Expand[(((-1))^k)Coefficient[R,z,-2k]],-1],{k,1,2(M-1)}]

В результате, получаем такие плотности (если хочется больше, увеличим M):

r[0]
r[1]
r[2]
r[3]
r[4]
r[5]

Проверяем, что это плотности: дифференцируем в силу VL, должны получиться полные разности. Следовательно, если применить вариационную производную, должны получиться нули.

star[Table[r[k],{k,0,M-1}],u[n](u[n+1]-u[n-1])];
vard/@%

Высшие симметрии

Высшие симметрии определяются как

(Subscript[u, n,(Subscript[t, k])])=(Subscript[u, n])(T-(T^(-1)))((Subscript[u, n])δ (Subscript[r, n]^((k)))/δ (Subscript[u, n]))

uT[k_]:=Factor[u[n](T[#,1]-T[#,-1])&[u[n]vard[r[k]]]]

uT[0]
uT[1]
uT[2]
uT[3]

Проверка коммутативности

Expand[star[uT[1],uT[2]]-star[uT[2],uT[1]]]
Expand[star[uT[1],uT[3]]-star[uT[3],uT[1]]]
Expand[star[uT[1],uT[4]]-star[uT[4],uT[1]]]

Expand[star[uT[2],uT[3]]-star[uT[3],uT[2]]]
Expand[star[uT[2],uT[4]]-star[uT[4],uT[2]]]
Expand[star[uT[3],uT[4]]-star[uT[4],uT[3]]]

Солитон

Проверка связи

Clear[u]
ut=u[n](u[n+1]-u[n-1]);

F=u[n](u[n+1]+u[n]+u[n-1])+A u[n]+B+(((-1))^n)C;

Factor[star[F,ut]-u[n](T[F,1]-T[F,-1])]

Пара легких упражнений.

1) Определим (Subscript[u, n])(t) по явным формулам

(Subscript[u, n])(t)=(((Subscript[w, n+1])(t)(Subscript[w, n-2])(t))/((Subscript[w, n])(t)(Subscript[w, n-1])(t))),   (Subscript[w, n])(t)=1+c (a^n)exp((a-(a^(-1)))x)

где a и c произвольные постоянные, и проверим, что она удовлетворяет VL и связи.

Clear[a,c,u,w]
w[n_]:=1+c (a^n)Exp[(a-(a^(-1)))t];
u[n_]:=((w[n+1]w[n-2])/(w[n]w[n-1]))
Together[D[u[n],t]-u[n](u[n+1]-u[n-1])]

Together[F/.{A->-B-3,C->0}]

Together[%/.B->a+(a^(-1))+1]

Определим функцию для построения графика решения. Для наглядности, массив выведем два раза: как ломаную и как набор точек. Построим график в манипуляторе, чтобы можно было менять t и параметры.

pl[aa_,cc_,tt_,minu_,maxu_]:=
ListPlot[
Evaluate[{#,#}&[Table[{n,u[n]/.{t->tt,a->aa,c->cc}},{n,-30,30}]]],
PlotRange->{minu,maxu},
Joined->{True,False},
PlotStyle->{Black,{Blue,PointSize[0.015]}},
BaseStyle->{FontSize->14,FontFamily->"Times New Roman"},
Ticks->{Automatic,Range[0.9,1.1,0.05]},
AxesLabel->{"n","u"},
ImageSize->400]

Manipulate[pl[aa,cc,tt,minu,maxu],
{{tt,0,"t"},-15,15,Appearance->"Labeled"},
{{aa,1.5,"a"},-5,5,Appearance->"Labeled"},
{{cc,1,"c"},-3,3,Appearance->"Labeled"},
{{minu,.99},-10,1},
{{maxu,1.11},1,20}
]

2) То же самое для решения

(Subscript[u, n])=-((b(n+d)(e^(2b t)))/(a + (e^(2b t)))) при чётных n и (Subscript[u, n])=((b (a n + c))/(a + (e^(2b t)))) при нечётных n.

Это решение растет по n, то есть, оно не относится к солитонному типу. Но, тем не менее, это некоторое явное решение.

Clear[a,b,c,d,v]
u[n_]:=If[EvenQ[n],-((b(n+d)Exp[2b t])/(a+Exp[2b t])),((b(a n+c))/(a+Exp[2b t]))]

Table[Together[D[u[n],t]-u[n](u[n+1]-u[n-1])],{n,-10,10}]

Manipulate[
ListPlot[
Evaluate[{#,#}&[Table[{n,u[n]/.{t->tt,a->aa,b->bb,c->cc,d->dd}},{n,-30,30}]]],
PlotRange->{-13,13},
Joined->{True,False},
PlotStyle->{Black,{Blue,PointSize[0.015]}},
BaseStyle->{FontSize->14,FontFamily->"Times New Roman"},
Ticks->{Automatic,Range[0.9,1.1,0.05]},
AxesLabel->{"n","u"},
ImageSize->400],
{{tt,0},-5,5,Appearance->"Labeled"},
{{aa,1.5},1,5,Appearance->"Labeled"},
{{bb,1.5},1,5,Appearance->"Labeled"},
{{cc,1.5},1,5,Appearance->"Labeled"},
{{dd,1},-3,3,Appearance->"Labeled"}
]