В.Э. Адлер · Классические интегрируемые системы

Лекция 10 · 10 апреля 2023

Некоторые вычисления на Mathematica

Построение графиков

. Plot. Гармонический осциллятор

Определяем собственную функцию основного состояния и последующие при помощи оператора рождения. 
◼ присваивание
◼ определение функции
◼ Factor 
◼ дифференцирование D 
◼ тест “такой, что”

F[0]=Exp[-(x^2)/2];
F[n_]:=Factor[-D[F[n-1],x]+x F[n-1]]/;n>0

Можно убедиться, что результат выражается через многочлены Эрмита (они встроены в Mathematica, как и многие другие спецфункции).
◼ Table – создание списка 
◼ % – предыдущий результат  
◼ применение Expand “сзади”

F[0]
F[1]
F[2]
F[3]
F[4]
Table[F[j]-HermiteH[j,x]F[0],{j,0,5}]
%//Expand

Определяем нормировочные константы.
◼ Интегрирование

c[n_]:=((Integrate[((F[n])^2),{x,-∞,∞}])^(1/2))
Table[c[j],{j,0,5}]

Строим график.
◼ Plot и некоторые опции

osc=Plot[(* список выводимых функций *)
Evaluate[Append[Table[F[j]/c[j]+2j+1,{j,0,5}],(x^2)]], 
(* переменная *)
{x,-4.5,4.5}, 
(* опции *)
Filling->Table[j->2j-1,{j,0,6}], (* заливка *)
Ticks->{Automatic, (* метки по оси х *)
Table[If[OddQ[i],{i,i,{0,.015}},{i,,{0,.01}}],{i,0,13,1/2}]}, (* метки по оси y *)
PlotRange->{{-4.6,4.6},{-0.5,12.2}},  (* выводимый прямоугольник *)
ImageSize->400,(* выводимый прямоугольник *)
AxesLabel->{"x","u, λ"}, (* метки на осях *)
BaseStyle->{FontSize->16,FontFamily->"Times New Roman"} ] (* шрифт *)

Запись изображения на диск.
◼ <> конкатенация строк

Export[NotebookDirectory[]<>"osc.pdf",osc,"PDF"];

. ListPlot и Manipulate. Солитон цепочки Вольтерры

Определим функцию (Subscript[u, n])(t) по явным формулам и проверим, что она удовлетворяет цепочке (Subscript[u, n,t])=(Subscript[u, n])((Subscript[u, n+1])-(Subscript[u, n-1])).
◼ Clear – очистка определений 
◼ Together – приведение к общему знаменателю

Clear[a,c,t]
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])]

Определим функцию plot. Она выводит график при фиксированных значениях a,c,t. Для наглядности, два раза – в виде точек и в виде линии.
◼ Пример лямбда-функции:  {#,#}&[x]  даёт  {x,x}

pl[a0_,c0_,t0_]:=
ListPlot[
(* список точек (два раза) *)
{#,#}&[Table[{n,u[n]/.{t->t0,a->a0,c->c0}},{n,-30,30}]],
(* опции *)
PlotRange->{0.99,1.102},
Joined->{True,False},  (* точки из первого списка соединяем, из второго нет *)
PlotStyle->{Brown,{Blue,PointSize[0.012]}}, (* стиль, для каждого списка *)
BaseStyle->{FontSize->14,FontFamily->"Times New Roman"},
AxesLabel->{"n","u"},
ImageSize->400]

pl[1.5,1,0]

Запустим график в манипуляторе.
◼ Динамический контент желательно удалять после использования.

Manipulate[pl[a,c,t],
{{t,1},-5,5,Appearance->"Labeled"},
{{a,1.5},-1,5,Appearance->"Labeled"},
{{c,1},-3,3,Appearance->"Labeled"}]

. Plot3D. Солитон огибающей ((SuperPlus[NLS]))

Комплексные числа
◼  i (а также π,e) – зарезервированные обозначения
◼ ComplexExpand предполагает, что все буквы обозначат вещественные переменные

Clear[a,b,c,d]

u=((a Exp[I (b x+((b^2)-(a^2))t+c)])/(Cosh[a(x+2b t)+d]))
v=u/.I->-I

Simplify[-I D[u,t]+D[u,x,x]+2(u^2)v]   (* проверили, что это решение НУШ *)

Simplify[ComplexExpand[ReIm[u]]]  (* вещественная и мнимая части решения *)

plot[a_,b_,c_,d_,t_]:=Module[{p=((a Cos[b x+((b^2)-(a^2))t+c])/(Cosh[a(x+2b t)+d])),q=((a Sin[b x+((b^2)-(a^2))t+c])/(Cosh[a(x+2b t)+d]))},
Plot[{p,q,Sqrt[(p^2)+(q^2)]},{x,-7,7},
PlotStyle->{Blue,Green,Red},
PlotRange->{-1.15,1.15},
BaseStyle->{FontSize->16,FontFamily->"Times New Roman"},
ImageSize->400,
AxesLabel->{"x",None},
PlotLegends->LineLegend[{Blue,Green,Red},{"Re u","Im u","|u|"}]
]]

Manipulate[
plot[a,b,c,d,t],
{{a,1},0,2,Appearance->"Labeled"},
{{b,10},-15,15,Appearance->"Labeled"},
{{c,1},-3,3,Appearance->"Labeled"},
{{d,1},-3,3,Appearance->"Labeled"},
{{t,0},-1,1,Appearance->"Labeled"}]

Plot3D[((a Cos[b x+((b^2)-(a^2))t+c])/(Cosh[a(x+2b t)+d]))/.{a->1,b->3,c->0,d->0},
{x,-10,10},{t,-2,2},
PlotPoints->100,
Mesh->False,
PlotRange->{-1.2,1.2}]

. ParametricPlot, ParametricPlot3D. Уравнение Хопфа

(Subscript[u, t])=u (Subscript[u, x])   ⇔  u(x,t)=f(x+u(x,t)t)

Проверим, что это действительно решение.
◼ В данном случае переменные с индексом – просто обозначение, оно никак не интерпретируется машиной.  Можно было бы обозначить ux вместо (Subscript[u, x]). Использование индексов иногда приводит к недоразумениям.

Clear[x,y,t,u,f]
eq=u-f[x+u t];
eqx=D[eq,x]+D[eq,u](Subscript[u, x])
eqt=D[eq,t]+D[eq,u](Subscript[u, t])
Factor[eqt-u eqx]

Subscript[u,x]
D[%,x]

Зададим начальное условие u(x,0)=f(x). Построим кривые (x(s),u(s)), где параметр s=x+u t. Так как u=f(s), то x=s-t f(s).
◼ N[expr] переводит все точные числа в выражении в формат float
ToString[expr] возвращает строку из числа (или выражения)

f[x_]:=Exp[-0.05(((x-1))^2)](4Sin[x]+Cos[2x]+Cos[3x])

g1[t_]:=ParametricPlot[
(* две кривые: при t=0 и при заданном t *)
{{s,f[s]},{s-f[s]t,f[s]}},
(* интервал изменения параметра, совпадает с интервалом для x *)
{s,-4,7},
(* опции *)
PlotStyle->{Blue,Black},
PlotRange->{{-4.2,7.2},Automatic},
AspectRatio->Automatic, (* одинаковый масштаб по осям *)
AxesLabel->{"x","u"},
Epilog->{Text["\!\(\*SubscriptBox[\(t\), \(\\ \)]\)\!\(\*SubscriptBox[\(=\), \(\\ \)]\)"<>ToString[N[t]],{4.5,3.6},{-1,0}]}, (* вывод текста *)
BaseStyle->{FontSize->16,FontFamily->"Times New Roman"}]

Manipulate[g1[t],{{t,1},-1,1,Appearance -> "Labeled"}]

То же самое можно сделать в 3D.

(Subscript[u, t])=u (Subscript[u, x])+(u^2)(Subscript[u, y]) 	u(x,y,t)=f(x+u(x,y,t)t,y+(u^2)(x,y,t)t)

Проверка

Clear[f]
eq=u-f[x+u t,y+(u^2)t];
eqx=D[eq,x]+D[eq,u](Subscript[u, x])
eqy=D[eq,y]+D[eq,u](Subscript[u, y])
eqt=D[eq,t]+D[eq,u](Subscript[u, t])
Factor[eqt-u eqx-(u^2)eqy]

f[x_,y_]:=Sin[x]Cos[y]

Manipulate[ParametricPlot3D[{sx-f[sx,sy]t,sy-((f[sx,sy])^2)t,f[sx,sy]},{sx,-7,7},{sy,-7,7},
PlotRange->{{-7,7},{-7,7},{-1.1,1.1}},
AspectRatio->Automatic,
ImageSize->500],
{{t,0.7},-2,2,Appearance -> "Labeled"}]

Численное решение ОДУ

. NDSolve, ListPlot. Цепочка Вольтерры

(* задаём период и время счёта *)
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

(* функция для вывода. Для наглядности, массив выводится два раза: как ломаная линия 
и как набор жирных точек в ее вершинах *)
plotv1[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[plotv1[t],{{t,40},0,t1}]

Запись кадров в пдф-файлы, для последующей вставки в латех: (это довольно долго)

Do[Export[NotebookDirectory[]<>"VL_"<>ToString[j]<>".pdf",plotv1[j],"PDF"],{j,0,100}]

Также можно сохранить список кадров в виде гиф-файла:

Export[NotebookDirectory[]<>"VL.gif",Table[plotv1[j],{j,0,100}],"GIF"];

Если взять начальное условие с более широким горбиком, то мы придем примерно к той картине, что была в эксперименте Забуски-Краскала для КдФ. Напомним, что там решалось чисто разностное уравнение с начальными условиями в виде синусоидальной волны с большим периодом (порядка 100 узлов решетки) и не очень большой амплитудой. Такое начальное условие распадалось на несколько солитонов плюс “рябь” (число солитонов зависит от начальной амплитуды), которые двигались с разной скоростью и, благодаря периодичности, всё время нагоняли друг друга. Качественно, в цепочке Вольтерры наблюдается точно такая же картина.

M2=200;
t2=3000;

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]==1+0.07Cos[((2π)/M2) n]},{n,0,M2-1}];
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

plotv2[t_]:=ListPlot[{u2[t],u2[t]},
PlotRange->{{-0.1,M2+1},{0.9,1.3}},
PlotStyle->{{Thin,Black},{Blue,PointSize[0.007]}},
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.25},{-1,0}]}
]

Manipulate[plotv2[t],{{t,540},0,t2}]

. NDSolve, ListPlot3D. Коммутирующие векторные поля, 2-зонное решение КдФ

Задача. Пусть (Subscript[u, n])=(∂_x^n)(u). Рассмотрим уравнение Новикова

(1)	(Subscript[u, 4])-10u (Subscript[u, 2])-5(Subscript[u, 1]^2)+10(u^3)+(Subscript[c, 1]) ((Subscript[u, 2])-3(u^2))+(Subscript[c, 2])u+(Subscript[c, 3])=0,

где (Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]) произвольные постоянные (это один раз проинтегрированное стационарное уравнение для симметрии КдФ пятого порядка плюс младшие симметрии). Проверить, что (1) совместно с уравнением КдФ

(2)  	(Subscript[u, t])=(Subscript[u, 3])-6 u (Subscript[u, 1]).

Построить график любого решения u(x,t) удовлетворяющего обоим уравнениям.

Решение.

Первый способ символьной проверки. Определим эволюционное дифференцирование (Subscript[D, t]) на бесконечном наборе (Subscript[u, n]) и проверим выполнение тождества

(Subscript[D, t])(G)=((Subscript[D, x]^3)-6u (Subscript[D, x]))(G),

где G – левая часть (1). Тогда, если G=0, то и (Subscript[D, t])(G)=0, что и требуется.
◼ Block – объявляем переменную a локальной.

dif[f_]:=Block[{a},D[f/.u[n_]:>u[n]+a du[n],a]/.a->0]
dx[f_]:=dif[f]/.du[n_]:>u[n+1]
dx[f_,n_]:=Nest[dx,f,n]

ut=u[3]-6u[0]u[1];
dt[f_]:=dif[f]/.du[n_]:>dx[ut,n]

G=u[4]-10u[0]u[2]-5((u[1])^2)+10((u[0])^3)+c1(u[2]-3((u[0])^2))+c2 u[0]+c3;
Expand[dt[G]-dx[G,3]+6u[0]dx[G]]

Второй способ символьной проверки. Перепишем оба уравнения в виде конечномерных динамических систем. Так как четвертая производная выражается из ОДУ (1) через младшие, то динамическими переменными являются u=(Subscript[u, 0]),(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3]). Сразу можно написать векторное поле (Subscript[D, x]) (точнее – то, как оно действует на произвольную функцию от динамических переменных).
	После этого определяем (Subscript[D, t]). Уравнение (1) задаёт правило дифференцирования по t только для переменной (Subscript[u, 0]), а нам нужно ещё для (Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3]). Значит, нужно продолжить дифференцирование (Subscript[D, t]) на эти производные, пользуясь тем, что (Subscript[D, x]) уже определено и полагая, по определению, что (Subscript[D, t])((Subscript[u, k]))=(Subscript[D, x]^k)((Subscript[u, t])). Так мы и делали выше при определении (Subscript[D, t]). Сейчас лишь заменяем (Subscript[D, x]) на укороченное векторное поле, действующее на четырех переменных.

dx[f_]:=u[1]D[f,u[0]]+u[2]D[f,u[1]]+u[3]D[f,u[2]]+
(10u[0]u[2]+5((u[1])^2)-10((u[0])^3)-c1(u[2]-3((u[0])^2))-c2 u[0]-c3)D[f,u[3]]

dt[f_]:=ut D[f,u[0]]+dx[ut]D[f,u[1]]+dx[ut,2]D[f,u[2]]+dx[ut,3]D[f,u[3]]

Можно посмотреть на компоненты векторного поля (Subscript[D, t]):

dt[u[0]]
Collect[dt[u[1]],{u[3],u[2],u[1],u[0]},Factor]
Collect[dt[u[2]],{u[3],u[2],u[1],u[0]},Factor]
Collect[dt[u[3]],{u[3],u[2],u[1],u[0]},Factor]

Итак, построили два векторных поля. Совместность означает, что они коммутируют. Считаем коммутатор. Для этого нужно вычислить его значения на координатных функциях (= динамических переменных), то есть, [(Subscript[D, x]),(Subscript[D, t])]((Subscript[u, k])), k=0,1,2,3:

Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

Всё сошлось. На самом деле, ясно, что равенство нулю коммутатора на первых трёх координатах выполняется автоматически – в силу построения дифференцирования (Subscript[D, t]) по формуле (Subscript[D, t])((Subscript[u, k]))=(Subscript[D, x]^k)((Subscript[u, t])). Для примера, испортим что-то в уравнениях. Тогда коммутативность нарушится, но только на последней координате. Так что, можно было бы проверять только её.

(* портим ut *)
ut=u[3]-66u[0]u[1];
Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

(* откатываем назад *)
ut=u[3]-6u[0]u[1];
Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

Численное построение решения. (Более эффективный счёт получается при переходе к уравнениям Дубровина, которые дают также аналитический метод решения, см. лекцию 5. Но, здесь мы этого не касаемся.) Определяем решение задачи Коши по x:

nsolx[{u00_,u10_,u20_,u30_},{c1_,c2_,c3_},{x0_,x1_}]:=
NDSolve[
{u[0]'[x]==u[1][x],
u[1]'[x]==u[2][x],
u[2]'[x]==u[3][x],
u[3]'[x]==10u[0][x]u[2][x]+5((u[1][x])^2)-10((u[0][x])^3)-c1(u[2][x]-3((u[0][x])^2))-c2 u[0][x]-c3,
u[0][x0]==u00,
u[1][x0]==u10,
u[2][x0]==u20,
u[3][x0]==u30},
{u[0],u[1],u[2],u[3]},
{x,x0,x1}][[1]]

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

dt[u[0]]/.u[k_]:>u[k][t]
Expand[dt[u[1]]]/.u[k_]:>u[k][t]
Expand[dt[u[2]]]/.u[k_]:>u[k][t]
Expand[dt[u[3]]]/.u[k_]:>u[k][t]

nsolt[{u00_,u10_,u20_,u30_},{c1_,c2_,c3_},{t0_,t1_}]:=
NDSolve[
{u[0]'[t]==-6 u[0][t] u[1][t]+u[3][t],
u[1]'[t]==-c3-c2 u[0][t]+3 c1 ((u[0][t])^2)-10 ((u[0][t])^3)-((u[1][t])^2)-c1 u[2][t]+4 u[0][t] u[2][t],
u[2]'[t]==-c2 u[1][t]+6 c1 u[0][t] u[1][t]-30 ((u[0][t])^2) u[1][t]+2 u[1][t] u[2][t]-c1 u[3][t]+4 u[0][t] u[3][t],
u[3]'[t]==c1 c3+c1 c2 u[0][t]-4 c3 u[0][t]-3 (c1^2) ((u[0][t])^2)-4 c2 ((u[0][t])^2)+22 c1 ((u[0][t])^3)-40 ((u[0][t])^4)+c1 ((u[1][t])^2)-40 u[0][t] ((u[1][t])^2)+(c1^2) u[2][t]-c2 u[2][t]-8 c1 u[0][t] u[2][t]+10 ((u[0][t])^2) u[2][t]+2 ((u[2][t])^2)+6 u[1][t] u[3][t],
u[0][t0]==u00,
u[1][t0]==u10,
u[2][t0]==u20,
u[3][t0]==u30},
{u[0],u[1],u[2],u[3]},
{t,t0,t1}][[1]]

Решаем задачу Коши для (1) от x=(Subscript[x, 0]) до (Subscript[x, 1]). Значения вектора ((Subscript[u, 0]),(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])) в промежуточных точках, с некоторым шагом δ, принимаем в качестве начальных условий для (2). Решаем соответствующие задачи Коши от t=(Subscript[t, 0]) до (Subscript[t, 1]).
	Потом делаем то же самое в другом порядке, с шагом ϵ по t. Сравниваем получившиеся сетки. Если значения совпадают – значит функция u(x,t) определена корректно. Можно строить её 3D график по полученной сетке.
	Следует отметить, что в этой задаче многие решения имеют полюсные особенности. Ниже, значения параметров и начальных условий подобраны так, чтобы их не было.

x0=0;x1=55; (* интервал по x *)
Mx=200; (* число точек в сетке по x для 3D графика *)
mx=10; (* число точек для линейчатых графиков *)

t0=0;t1=25; (* интервал по t *)
Mt=200; (* число точек в сетке по t для 3D графика *)
mt=10; (* число точек для линейчатых графиков *)

cc={0.15,-0.7,0}; (* параметры *)
U0={-0.1,0.1,0.1,-0.125};  (* начальные данные при x0, t0 *)

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

(* решаем систему по t; 
решаем систему по x для последовательности значений t[0]= t0, ..., t[mt]=t1;
результат сохраняем в виде графического примитива Line *)

nst=nsolt[U0,cc,{t0,t1}];

nstx1=Table[
ns=nsolx[{u[0][t],u[1][t],u[2][t],u[3][t]}/.t->t0+j(t1-t0)/mt/.nst,cc,{x0,x1}];
Line[Table[{x0+i(x1-x0)/Mx,t0+j(t1-t0)/mt,u[0][x0+i(x1-x0)/Mx]/.ns},{i,0,Mx}]],
{j,0,mt}];

(* делаем то же самое, меняя ролями x и t *)
nsx=nsolx[U0,cc,{x0,x1}];

nsxt1=Table[
ns=nsolt[{u[0][x],u[1][x],u[2][x],u[3][x]}/.x->x0+i(x1-x0)/mx/.nsx,cc,{t0,t1}];
Line[Table[{x0+i(x1-x0)/mx,t0+j(t1-t0)/Mt,u[0][t0+j(t1-t0)/Mt]/.ns},{j,0,Mt}]],
{i,0,mx}];

gr0=Graphics3D[{Thick,Red,nstx1[[1]],Blue,nsxt1},
PlotRange->{{x0,x1},{t0,t1},{-1,2}},
Axes->True,
AxesLabel->Evaluate[Style[#,Italic,24]&/@{"x","t","u"}],
ImageSize->500,
BoxRatios->{1,1,0.2},
BaseStyle->{FontFamily->"Times",FontSize->16}]

gr1=Graphics3D[{Thick,Red,nstx1,Blue,nsxt1},
PlotRange->{{x0,x1},{t0,t1},{-1,2}},
Axes->True,
AxesLabel->Evaluate[Style[#,Italic,24]&/@{"x","t","u"}],
ImageSize->500,
BoxRatios->{1,1,0.2},
BaseStyle->{FontFamily->"Times",FontSize->16}]

Красные кривые (эволюция по x) пересекают синие (эволюция по t), что визуально подтверждает коммутативность потоков.

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

(* решаем систему по t; 
решаем систему по x для последовательности значений t[0]= t0, ..., t[Mt]=t1; 
сохраняем результат в виде 2D массива *)

nst=nsolt[U0,cc,{t0,t1}];
nstx2=Table[
ns=nsolx[{u[0][t],u[1][t],u[2][t],u[3][t]}/.t->t0+j(t1-t0)/Mt/.nst,cc,{x0,x1}];
Table[u[0][x0+i(x1-x0)/Mx]/.ns,{i,0,Mx}],{j,0,Mt}];

(* делаем то же самое, меняя ролями x и t *)
nsx=nsolx[U0,cc,{x0,x1}];
nsxt2=Table[
ns=nsolt[{u[0][x],u[1][x],u[2][x],u[3][x]}/.x->x0+i(x1-x0)/Mx/.nsx,cc,{t0,t1}];
Table[u[0][t0+j(t1-t0)/Mt]/.ns,{j,0,Mt}],{i,0,Mx}];

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

Max[Abs[nstx2-Transpose[nsxt2]]]

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

gr3=ListPlot3D[nstx2,
DataRange->{{x0,x1},{t0,t1}},
PlotRange->{-1,2},
Mesh->None,
AxesLabel->Evaluate[Style[#,Italic,24]&/@{"x","t","u"}],
ImageSize->500,
BaseStyle->{FontFamily->"Times",FontSize->16},
BoxRatios->{1,1,0.2},
PlotRangePadding->0]

gr2=Show[gr3,gr1]

Символьные вычисления

Определения из раздела 3.1 используются во всех дальнейших.

. (Subscript[D, x]),(SubStar[f])(g),E(f), интегрирование по частям

Определение (Subscript[D, x]), (SubStar[f])(g), E(f)

Некоторые конструкции функционального программирования:
◼ _u – пример паттерна: всё, что имеет голову u, например u[n], u[x,y], ... 
◼ Cases[f,...] – вхождения в выражение f
◼ Union – объединение элементов списка
◼ лямбда-функция: (выражение с #)&
◼ f/@xxx – применение f к каждому элементу списка xxx
◼ y@@xxx – замена головы выражения xxx  на y (в нашем случае List заменяется на Plus)
◼ Nest[f,x,n] – применяет f к x n раз

vars[f_]:=Union[Cases[f,_u,{0,∞}]]
dif[f_]:=Plus@@(D[f,#]dd[#]&/@vars[f])

dx[f_]:=D[f,x]+dif[f]/.dd[u[k_]]:>u[k+1]
dx[f_,k_]:=Nest[dx,f,k]

Теперь можем определить эволюционную производную

(SubStar[f])(g)=(Subscript[∇, g])(f)=(∂_0)(f)+... +(∂_n)(f)(Subscript[D, x]^n)(g),

и вариационную производную (оператор Эйлера)

E(f)=((δ f)/(δ u))=(∂_0)(f)-(Subscript[D, x])((∂_1)(f))+(Subscript[D, x]^2)((∂_2)(f))-...+(((-(Subscript[D, x])))^n)((∂_n)(f)),

где f=f(x,(Subscript[u, 0]),(Subscript[u, 1]),...,(Subscript[u, n])) и (∂_k) обозначает частную производную по (Subscript[u, k]).

star[g_,f_]:=dif[g]/.dd[u[k_]]:>dx[f,k]
vard[f_]:=Expand[Plus@@((-1)^#[[1]]dx[D[f,#],#[[1]]]&/@vars[f])]

Утверждение: выражение является полной производной по x тогда и только тогда, когда его вариационная производная равна 0:

E(f)=0   ⇔  ∃g:f=(Subscript[D, x])(g).

Алгоритм интегрирования по частям

Некоторые конструкции процедурного программирования:
◼ If[cond,t,f,u]  t,f,u – результаты, если cond=True,False, undefined 
◼ And[(Subscript[cond, 1]),(Subscript[cond, 2]),...]=(Subscript[cond, 1])&&(Subscript[cond, 2])&&...
◼ Module[{локальные переменные}, команды]
◼ While[cond, команды]

(* порядок по производным *)
ord[f_]:=If[#=={},-1,#[[-1]][[1]]]&[vars[f]]

int[f_]:=Module[{a=0,b=Expand[f],c,k}, 
While[k=ord[b];     (* вычисляем k - порядок b *)
And[k>0,                   (* выполняем цикл, пока выполняется условие, что k>0 и *)
c=D[b,u[k]];     (* b линейна по старшей переменной *)
Expand[D[c,u[k]]]==0
],
c=Integrate[c,u[k-1]];   (* интегрируем по переменной, следующей по старшинству *)
a=Expand[a+c];                        (* пересчет a и b *)
b=Expand[b-dx[c]]
];
If[ord[b]==-1,    (* если в результате остаток не зависит от u, то это функция от х *)
{Expand[a+Integrate[b,x]],0},
{a,b}  (* возвращаем пару (a,b)- результат и препятствие *)
]
]

Протестируем. Определим какое-нибудь не слишком сложное выражение (многочлен).

f=Expand[(((u[0]+((u[1])^2)+3x))^2)+2u[3]u[4]Exp[u[0]]]
vard[f]
int[f]

Как видим, f∉Im (Subscript[D, x]). Сделаем проверку для dx[f].

Expand[dx[f]]
vard[%]
int[%%]
%[[1]]-f

. Законы сохранения КдФ (обращение преобразования Миуры)

В лекции 2 было доказано, что ряд

v=-z+(Subscript[V, 0])+(Subscript[V, 1])/(2z)+(Subscript[V, 2])/(((2z))^2)+...,

удовлетворяющий уравнению

(Subscript[v, x])+(v^2)=(z^2)+u,

служит производящей функцией для плотностей законов сохранения уравнения КдФ

(1) 	(Subscript[u, t])=(Subscript[u, 3])-6u (Subscript[u, 1]).

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

(Subscript[V, 0])=0, (Subscript[V, 1])=-u,

(Subscript[V, n+1]) =(Subscript[V, 1])(Subscript[V, n-1])+...+(Subscript[V, n-1])(Subscript[V, 1])+(Subscript[V, n,x])=0,   n=1,2,3,... .

Вычислим несколько первых коэффициентов.

M=13;
V[1]=-u[0];
Do[V[n+1]=Factor[(dx[V[n]]+Sum[V[s]V[n-s],{s,1,n-1}])],{n,1,M-1}];

Чтобы было удобнее смотреть, заменим u[n] на (Subscript[u, n]):

Do[Print[(Subscript[V, n]),"= ",V[n]/.u[k_]:>(Subscript[u, k])],{n,1,M}]

Можно проверить, что построенный отрезок ряда действительно удовлетворяет уравнению (5), с точностью до членов (z^n) с n≥M:
◼ Series[f,{z,(Subscript[z, 0]),n}]  разложение в ряд Тейлора или Лорана

v=-z+Sum[V[n]/(((2z))^n),{n,1,M}];
Series[Expand[(v^2)+dx[v]-(z^2)-u[0]],{z,Infinity,M-1}]

Теория утверждает, что:

(i)  (Subscript[V, 2n])=(Subscript[D, x])((Subscript[a, n]))  (тривиальные плотности)

(ii)  (Subscript[D, t])((Subscript[V, 2n-1]))=(Subscript[D, x])((Subscript[b, n]))   и  (Subscript[V, 2n-1])∉Im (Subscript[D, x])

для некоторых (Subscript[a, n]),(Subscript[b, n])  – функций от u,(Subscript[u, 1]),(Subscript[u, 2]),… Проверяем:

Table[vard[V[n]],{n,2,M,2}]
Table[vard[V[n]]===0,{n,1,M,2}]
Table[vard[star[V[n],u[3]-6u[0]u[1]]],{n,1,M,2}]

Функции (Subscript[b, n]) (токи законов сохранения) можно вычислить, применяя интегрирование по частям. Например:

ρ=V[11]
Expand[star[%,u[3]-6u[0]u[1]]]
int[%]
σ=%[[1]];

Expand[star[ρ,u[3]-6u[0]u[1]]-dx[σ]]

. Построение и проверка высших симметрий КдФ. Оператор рекурсии

Напомним (см. лекцию 4), что уравнение КдФ (1) допускает бесконечную последовательность высших симметрий нечётного порядка по производным

(2)	(Subscript[u, (Subscript[t, k])])=(Subscript[f, k])(u,...,(Subscript[u, 2k-1])),   k=1,2... ,

то есть, таких уравнений, что [(Subscript[D, t]),(Subscript[D, (Subscript[t, k])])]=0. В частности,

(Subscript[u, (Subscript[t, 1])])=(Subscript[f, 1])=(Subscript[u, 1]),  
	(Subscript[u, (Subscript[t, 2])])=(Subscript[f, 2])=(Subscript[u, 3])-6u (Subscript[u, 1]),
	(Subscript[u, (Subscript[t, 3])])=(Subscript[f, 3])=(Subscript[u, 5])-10u (Subscript[u, 3])-20(Subscript[u, 1])(Subscript[u, 2])+30(u^2)(Subscript[u, 1]),  
	(Subscript[u, (Subscript[t, 4])])=(Subscript[f, 4])=(Subscript[u, 7])-14u (Subscript[u, 5])-42(Subscript[u, 1])(Subscript[u, 4])-70 (Subscript[u, 2])(Subscript[u, 3])+70(u^2)(Subscript[u, 3])+280u (Subscript[u, 1])(Subscript[u, 2])+70(Subscript[u, 1]^3)-140(u^3)(Subscript[u, 1]),

и так далее. Первые две симметрии классические, то есть, отвечают однопараметрическим группам точечных преобразований (сдвиги x->x+a и t->t+a). Следующие симметрии – высшие.

Есть несколько способов вычисления (Subscript[f, k]):

(i) исходя из определения симметрии;
	(ii) метод неопределённых коэффициентов, с использованием однородности w(u)=2, w((∂_x))=1;
	(iii) при помощи мастер-симметрии;
	(iv) оператор рекурсии (лекция 4);
	(v) квадратичное рекуррентное соотношение (лекция 4);
	(vi) расширение рекуррентного соотношения для законов сохранения.

Способы (iv)–(vi) примерно эквивалентны и основаны на рекуррентном соотношении

(3)	(Subscript[D, x])((Subscript[f, k+1]))=((Subscript[D, x]^3)-4u (Subscript[D, x])-6(Subscript[u, 1]))((Subscript[f, k])).

Применяя (Subscript[D, x]^(-1)), можно привести это к виду

(Subscript[f, k+1])=R((Subscript[f, k])),  R=(Subscript[D, x]^2)-4u-2(Subscript[u, 1])(Subscript[D, x]^(-1)).

Вычислим симметрии до (Subscript[f, 15]) за 1.5 сек. 
◼ Timing[expr] возвращает {затраченное время, результат}
◼ в нашем случае результат пуст, так как цикл Do ничего не возвращает, только что-то делает внутри себя
◼ перед обращением к Timing очищаем f и кэш, так как иначе при повторном запуске ячейки время окажется заниженным – Mathematica не будет считать заново, а воспользуется сохраненными в памяти результатами.

R[f_]:=Expand[dx[f,2]-4u[0]f-2u[1]int[f][[1]]]

M=15;
Clear[f]
ClearSystemCache[]

Timing[
f[1]=u[1];
Do[f[j]=R[f[j-1]],{j,2,M}]]

Посмотрим, что получилось.

Do[Print[f[j]/.u[n_]:>(Subscript[u, n])],{j,1,7}]

Length/@Table[f[j],{j,1,M}]  (* OEIS A182747 *)

Проверим, что это действительно симметрии.

ClearSystemCache[]
Timing[
Table[star[f[2],f[j]]-star[f[j],f[2]],{j,1,M}]//Expand]

Проверим, что и друг с другом они коммутируют

Table[star[f[i],f[j]]-star[f[j],f[i]],{i,1,5},{j,i+1,6}]//Expand

. Явное рекуррентное соотношение

Пусть  (Subscript[f, k])=(Subscript[D, x])((Subscript[a, k]))  и  a=1/2+(Subscript[a, 1])/(-4λ)+(Subscript[a, 2])/(((-4λ))^2)+...   – производящая функция для (Subscript[a, k]). Тогда рекуррентное соотношение (3) сводится к

2a (Subscript[a, x x])-(Subscript[a, x]^2)+4(λ-u)(a^2)=λ.

В свою очередь, это сводится (проверьте!) к

(Subscript[v, x])+(v^2)=(z^2)+u,    a(λ)(v(-z)-v(z))=z,   (z^2)=-λ.

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

Clear[a,V,ff]
ClearSystemCache[];

Timing[
V[1]=-u[0];Do[V[n+1]=Factor[(dx[V[n]]+Sum[V[s]V[n-s],{s,1,n-1}])],{n,1,2M-2}];
a[0]=-1/2;
Do[a[k]=Expand[-V[2k-1]+2Sum[a[k-s]V[2s-1],{s,1,k-1}]];
ff[k]=dx[a[k]],{k,1,M}]
]

Table[Expand[ff[k]-f[k]],{k,1,M}]