maple 程序，编译，运行

######################################################################
# Neutral DDE 高维 Hopf 分岔 / 规范形计算脚本（合并整理版）
#
# 说明：
# - 原文主程序块重复三次，此处仅保留一份。
# - 已修正的明显笔误：
# lanbda -> lambda
# linearAlgebra -> LinearAlgebra
# o -> 0（for 循环下界）
# for i_1 from 0 (1+i) -> from 0 to (1+i)
# £» -> ;
# V_Fn[i][M] / (V1||i)[M] / c_element[M] -> [m]
# ww.l -> ww
# tau[jj1] -> tau[jj-1]
# S 表构造循环 summands:=modp(...) -> summands:=summands,modp(...)
# linalg[jacobian]((L||0,...) 多余括号 / convert(...)) 多余括号
# cat(_C,...) -> cat(`_C`,...)
# - 仍需核对（原代码本身可能不完整/依赖外部定义）：
# * sys_rhs 中 delta*(xt[4,1], xt[4,2]) 的向量写法
# * L||0, F||i, Gx||i, V1_Fn 等符号在脚本内未定义
# * S[summands]:=convert([summands],Matrix) 用 Matrix 作表键
# * tem*tem 是否本意为 tem^2
# * P_si := Matrix([[psi[1]]],[[psi[2]]]) 的构造方式
######################################################################

# ===================== 1. 参数与系统定义 =====================
n_dim := 1;
n_delay := 1;
n_order := 3;
Eq_pos := `<,>`(0, 0);
var_par := [delta];
dde_type := "neutral";

sys_rhs := `<,>`(
xt[2,1],
beta[1]*xt[2,1]*beta[2]*(xt[1,1]-xt[1,1]^3)
+ delta*<xt[4,1], xt[4,2]> # 原为 delta*(xt[4,1], xt[4,2])
);

sys_D_lhs := `<,>`(
xt[1,1],
(1-beta[3*c])*xt[2,1] + beta[3*c]*xt[2,2]
);

# ===================== 2. 工具过程 =====================
VarSeq := proc(ndim::integer, ndelay::integer,
varpar::{set, list, Vector[column]},
ddetype::{name, string, symbol} := "retarded")
local dtype, vpar, i, j, varseq;
dtype := convert(ddetype, string);
vpar := convert(varpar, list);
if dtype = "restarded" then
varseq := seq(seq(xt[i,j], j = 1..ndelay+1), i = 1..ndim), op(vpar)
elif dtype = "neutral" then
varseq := seq(seq(xt[i,j], j = 1..ndelay+1), i = 1..2*ndim), op(vpar)
else
error "%1,type of DDE not recognized", ddetype
end if;
varseq
end proc:

DegreeSeparate := proc(Eqs::{Matrix, Vector[column]},
varlist::{set, list},
Norder::integer)
local Eqs_temp, ndim, Eq_ex, i, nterm, j, term;
Eqs_temp := convert(Eqs, Vector[column]);
ndim := LinearAlgebra:-Dimension(Eqs);
Eq_ex := Vector[column](ndim);
for i to ndim do
nterm := nops(expand(Eqs_temp[i]+1))-1;
for j to nterm do
term := op(j, expand(Eqs_temp[i]+1));
if degree(term, varlist) = Norder then
Eq_ex[i] := Eq_ex[i] + term
end if
end do
end do;
Eq_ex
end proc:

FGSeparate := proc(Eqs::(Vector[column]), varlist::{set, list})
local ndim, i, nvar, j, Eq_F, Eq_G, term;
ndim := LinearAlgebra:-Dimension(Eqs);
Eq_G := Vector[column](ndim);
for i to ndim do
for j to nops(expand(Eqs[i])+1)-1 do
term := op(j, expand(Eqs[i])+1);
if 1 <= degree(term, varlist) then
Eq_G[i] := Eq_G[i] + term
end if
end do
end do;
Eq_F := Eqs - Eq_G;
Eq_F, Eq_G
end proc:

FLSeparate := proc(Eqs::(Vector[column]),
svarlist::{set, list},
varpar::{set, list, Vector[column]})
local ndim, Eq_ex, i, j, term;
ndim := LinearAlgebra:-Dimension(Eqs);
Eq_ex := Vector[column](ndim);
for i to ndim do
for j to nops(expand(Eqs[i])+1)-1 do
term := op(j, expand(Eqs[i])+1);
if 1 <= degree(term, varpar) then
if degree(term, svarlist) = 1 then
Eq_ex[i] := Eq_ex[i] + term
end if
end if
end do
end do;
Eq_ex, Eqs - Eq_ex
end proc:

CoeffMulvar := proc(expression::algebraic,
varlist::list, exponentslist::list)
local nvar, eq, i;
nvar := nops(varlist);
eq := expression;
for i to nvar do
eq := coeff(eq, varlist[i], exponentslist[i])
end do;
eq
end proc:

DiffMulvar := proc(Eqs::{Matrix, Vector[column]},
varlist::list, ww,
norder::integer, neqs::integer, ni::integer)
local Eq_max, nvar, rdim, cdim, j, summands, kk, k, DEQs, S,
m, jj, idx, ik, mm, nn, T, ii, cc;
Eq_max := convert(Eqs, Matrix);
nvar := nops(varlist);
rdim := LinearAlgebra:-RowDimension(Eq_max);
cdim := LinearAlgebra:-ColumnDimension(Eq_max);
DEQs := Matrix(rdim, cdim);
for j from 0 to (norder+1)^nvar-1 do
summands := NULL;
kk := j;
for m to nvar do
summands := summands, modp(kk, norder+1);
kk := floor(kk/(norder+1))
end do;
S[summands] := convert([summands], "+");
if S[summands] = norder then
T := combinat:-carprod([`$`([0, seq(i, i = 2..ni)], nvar)]);
while not T[finished] do
idx := T[nextvalue]();
if convert(idx, "+") <= ni+norder-neqs
and sum(idx[nn]*summands[nn], nn = 1..nvar) = ni+norder-neqs then
for ii to rdim do
for cc to cdim do
DEQs[ii,cc] := DEQs[ii,cc]
+ diff(Eq_max[ii,cc],
[seq(`$`(varlist[i], summands[i]), i = 1..nvar)])
* mul(ww[i][idx[i]]^summands[i], i = 1..nvar)
* mul(1/factorial(summands[i]), i = 1..nvar)
end do
end do
end if
end do
end if
end do;
DEQs
end proc:

SigSum := proc(eq, ndim::integer)
local r_eq;
if eq = 0 then
r_eq := Vector[column](ndim)
else
r_eq := eq
end if;
r_eq
end proc:

# ===================== 3. 主计算流程 =====================

# ---------- 系统 Taylor 展开 ----------
sys_rhs := map(mtaylor, sys_rhs,
[seq(seq(op([xt[i,j] = Eq_pos[i], xt[n_dim+i,j] = Eq_pos[i]]),
i = 1..n_dim), j = 1..(n_delay+1))],
n_order+1);

n_ode := 2;

# ---------- 变量序列 ----------
var := VarSeq(n_dim, n_delay, var_par, dde_type);
svar := seq(seq(xt[i,j], j = 1..n_delay+1), i = 1..n_dim);
dvar := seq(seq(xt[i,j], j = 1..n_delay+1), i = (n_dim+1)..2*n_dim);
n_var := nops([var]);
n_svar := nops([svar]);
n_dvar := nops([dvar]);

for i from 1 to n_delay+1 do
XT[i] := <seq(xt[j,i], j = 1..n_dim)>;
end do;

for i from 1 to n_delay+1 do
NXT[i] := <seq(xt[j,i], j = n_dim+1..2*n_dim)>;
end do;

V := <seq(v[i], i = 1..n_dim)>;
W := LinearAlgebra[Transpose](<seq(w[i], i = 1..n_dim)>);
Y := <seq(y[i], i = 1..n_ode)>;
ylist := [seq(y[i], i = 1..n_ode)];

H := <seq(0*i, i = 1..n_dim)>;
for j from 2 to n_order do
cat(h,j) := <seq(cat(h,j,i), i = 1..n_dim)>;
H := H + eval(cat(h,j));
end do;

ode_var := seq(y[i], i = 1..n_ode), op(var_par);
s_ode_var := seq(y[i], i = 1..n_ode);
n_ode_var := nops([ode_var]);

for ii from 1 to n_dim do
for jj from 2 to n_order do
cat(h,jj,ii) := 0;
end do;
end do;

# ---------- 齐次多项式生成 ----------
S := table(): aa := 0;
for j from 0 to (n_order+1)^n_ode_var-1 do
summands := NULL;
kk := j;
for m from 1 to n_ode_var do
summands := summands, modp(kk, n_order+1);
kk := floor(kk/(n_order+1));
end do;
S[summands] := convert([summands], Matrix);
for k from 2 to n_order do
if S[summands] = k then
for jj from 1 to n_dim do
aa := h||k||jj;
aa := aa
+ (h[jj][summands])(theta)
* mul(ode_var[i]^(summands)[i], i = 1..n_ode_var);
h||k||jj := aa;
end do;
end if;
end do;
end do;

# ---------- 平衡点处的线性化 ----------
E_subs := [];
for i from 1 to n_dim do
E_subs := [op(E_subs), seq(xt[i,j] = Eq_pos[i], j = 1..n_delay)];
end do;

for i from 0 to n_delay do
A[i] := convert(subs(E_subs, linalg[jacobian](L||0, XT[i+1])), Matrix);
B[i] := convert(subs(E_subs, linalg[jacobian](sys_D_lhs, XT[i+1])), Matrix);
end do;

# ---------- 特征方程 ----------
ident := Matrix(n_dim, shape = identity);

Delta := lambda*sum('B[i]*exp(-lambda*tau[i])', 'i' = 0..n_delay)
- sum('A[i]*exp(-lambda*tau[i])', 'i' = 0..n_delay);

char_eq := combine(collect(LinearAlgebra[Determinant](Delta), lambda));
D_Delta := subs(lambda = I*omega, map(diff, Delta, lambda));

eq_im := evalc(subs(lambda = I*omega, char_eq));
eq_Re := coeff(expand(eq_im), I, 0);
eq_Im := coeff(expand(eq_im), I, 1);

SC_subs := solve([eq_Re, eq_Im],
{op(selectfun(eq_Re, cos)),
op(selectfun(eq_Re, sin)),
op(selectfun(eq_Im, sin)),
op(selectfun(eq_Im, cos))});

# ---------- 中心空间基 ----------
Dv := subs(lambda = I*omega, LinearAlgebra[Multiply](Delta, V));
V := subs(op(solve([seq(Dv[i], i = 1..n_dim-1)],
[seq(v[i], i = 2..n_dim)])),
v[1] = 1, V);

assume_variables := seq(indets(eq_im, name)[i]::real, i = 1..nops(indets(eq_im, name))),
theta::real, xi::real;

phi[1] := V*exp(I*omega*theta);
phi[2] := evalc(map(conjugate, V)*exp(-I*omega*theta)) assuming assume_variables;
phi := Matrix(phi[1], phi[2]);

wD := subs(lambda = I*omega, LinearAlgebra[Multiply](W, Delta));
W := subs(op(solve([seq(wD[i], i = 1..n_dim-1)],
[seq(w[i], i = 2..n_dim)])), W);

w1 := solve(LinearAlgebra[Multiply](LinearAlgebra[Multiply](W, D_Delta), V) = 1, w[1]);

psi[2] := map(conjugate, W)*exp(I*omega*xi) assuming assume_variables;

P_si := Matrix([[psi[1]]], [[psi[2]]]);
P_si_0 := subs(xi = 0, P_si);

B := Matrix([[I*omega, 0], [0, -I*omega]]);

Phiy := LinearAlgebra[Multiply](phi, Y);
for i from 0 to n_delay do
Phiyt[i] := subs(theta = -tau[i], Phiy);
Phiyth[i] := eval(subs(theta = -tau[i], Phiy + H));
end do;

subs_group := {}: subs_group_t := {};
for i from 0 to n_delay do
for j from 1 to n_dim do
subs_group := {op(subs_group), XT[i+1][j] = Phiyt[i][j]};
subs_group_h := {op(subs_group_t), XT[i+1][j] = Phiyth[i][j]};
end do;
end do;

subs_group_W := {};
linalg[jacobian](L||0, XT[i+1], "Matrix");

# ---------- 逐阶计算 ----------
for i from 2 to n_order do
G_element||i := [];
L_element||i := [];
coe_element||i := [];
Resonant||i := [];

j:'j';
g1[i] := <seq(j*0, j = 1..n_ode)>;
V_Fn[i] := <seq(j*0, j = 1..n_ode)>;
V1||i := <seq(j*0, j = 1..n_ode)>;

indet_h||i := {}; eq_h||i := {}; eq_dh||i := {}; indet_c||i := {};
end do;

for i from 2 to n_order do
if dde_type = "neutral" then
t1 := subs(subs_group, L||(i-1))
+ SigSum(add(subs(subs_group_L[k], L||(i-k)), k = 2..(i-1)), n_dim);

t2 := subs(subs_group, F||i)
+ SigSum(add(subs(subs_group_W, subs_group,
convert(DiffMulvar(F||k, [svar], ww, l, k, i),
Vector[column])),
k = 2..(i-1)), n_dim)
+ SigSum(add(add(subs(subs_group_W, subs_group,
convert(DiffMulvar(F||k, [svar], ww, j, k, i),
Vector[column])),
k = j..i-j), j = 2..floor(i/2)), n_dim);

t3 := SigSum(add(LinearAlgebra[Multiply](
subs([op(subs_group_W), op(subs_group)],
VectorCalculus[Jacobian](Gx||i, convert(NXT[jj], list))),
LinearAlgebra[Multiply](
subs(theta = -tau[jj-1], Phi),
LinearAlgebra[Multiply](B, Y))),
jj = 1..n_delay+1), n_dim)
+ SigSum(add(LinearAlgebra[Multiply](
add(add(subs(subs_group_W, subs_group,
DiffMulvar(
VectorCalculus[Jacobian](Gx||(k+1),
convert(NXT[jj], list)),
[svar], ww, j, k, i-1)),
k = j..(i-1-j)), j = 1..floor((i-1)/2)),
LinearAlgebra[Multiply](
subs(theta = -tau[jj-1], Phi),
LinearAlgebra[Multiply](B, Y))),
jj = 1..n_delay+1), n_dim)
+ SigSum(add(add(LinearAlgebra[Multiply](
(VectorCalculus[Jacobian](Gx||(I+1),
convert(NXT[jj], list))
+ add(add(subs(subs_group_W, subs_group,
DiffMulvar(
VectorCalculus[Jacobian](Gx||(k+1),
convert(NXT[jj], list)),
[svar], ww, j, k, I)),
k = j..I-j), j = 1..floor(I/2))),
subs(theta = -tau[jj-1],
(LinearAlgebra[Multiply](Phi, g1[i-1])
+ LinearAlgebra[Multiply](
LinearAlgebra[Multiply](Phi,
VectorCalculus[Jacobian](V1||(i-1), ylist))
+ VectorCalculus[Jacobian](V2||(i-1), ylist),
LinearAlgebra[Multiply](B, Y))
+ SigSum(add(LinearAlgebra[Multiply](
LinearAlgebra[Multiply](Phi,
VectorCalculus[Jacobian](V1||k, ylist))
+ VectorCalculus[Jacobian](V2||k, ylist),
g1[i-I-k+1]),
k = 2..i-I-1), n_dim)))),
jj = 1..n_delay+1), I = 1..(i-2)), n_dim);

K[i] := t1 + t2 + t3;

elif dde_type = "retarded" then
t1 := subs(subs_group, L||(i-1))
+ SigSum(add(subs(subs_group_L[k], L||(i-k)), k = 2..(i-1)), n_dim);
t2 := subs(subs_group, F||i)
+ SigSum(add(subs(subs_group_W, subs_group,
convert(DiffMulvar(F||k, [svar], ww, l, k, i),
Vector[column])),
k = 2..(i-1)), n_dim)
+ SigSum(add(add(subs(subs_group_W, subs_group,
DiffMulvar(F||k, [svar], j, k, i)),
k = j..i-j), j = 2..floor(i/2)), n_dim);
K[i] := t1 + t2;
end if;

f1[i] := LinearAlgebra[Multiply](P_si_0, K[i])
- SigSum(add(LinearAlgebra[Multiply](
VectorCalculus[jacobian](V1||k, ylist),
g1[i+1-k]), k = 2..(i-1)), n_ode);

V2||i := h||i;

ode := -LinearAlgebra[Multiply](Phi,
LinearAlgebra[Multiply](P_si_0, K[i]))
- SigSum(add(LinearAlgebra[Multiply](
VectorCalculus[Jacobian](V2||k, ylist),
g1[i+1-k]), k = 2..i-1), n_dim)
- LinearAlgebra[Multiply](VectorCalculus[Jacobian](V2||i, ylist),
LinearAlgebra[Multiply](B, Y))
+ map(diff, V2||i, theta);

id := 0;

# ---- 按度数分离各项 ----
for i_1 from 0 to (1+i)^n_ode_var-1 do
sumindex := NULL;
k := i_1;
for j from 1 to n_ode_var do
sumindex := sumindex, modp(k, i+1)+1;
k := floor(k/(i+1));
end do;
SS[sumindex] := convert([sumindex], "+") - (n_ode_var);

if SS[sumindex] = i then
id := id + 1;
for m from 1 to n_ode do
ii := "ii":
element[m] := <seq(0*ii, ii = 1..n_ode)>;
ii := "ii":
element[m][m] := product(ode_var[ii]^(sumindex[ii]-1), ii = 1..n_ode_var);
sus[m] := convert(linalg[jacobian](element[m], Y), "Matrix");
M_element[m] := LinearAlgebra[Multiply](
LinearAlgebra[Multiply](sus[m], B), Y)
- LinearAlgebra[Multiply](B, element[m]);
c_element[m] := M_element[m][m]/element[m][m];

if evalb(c_element[m] = 0) then
Resonant||i := [element[m], op(Resonant||i)];
tem := CoeffMulvar(f1[i][m], [ode_var],
map("-", [sumindex], 1));
g1[i][m] := g1[i][m]
+ tem*tem*product(ode_var[ij]^(sumindex[ij]-1),
ij = 1..n_ode_var);
else
UResonat||i := [element[m], op(UResonat||i)];
tem := CoeffMulvar(f1[i][m], [ode_var],
map("-", [sumindex], 1));
V_Fn[i][m] := V1_Fn[i][m]
+ tem*product(ode_var[ij]^(sumindex[ij]-1),
ij = 1..n_ode_var);
(V1||i)[m] := (V1||i)[m]
+ tem/c_element[m]
* product(ode_var[ij]^(sumindex[ij]-1),
ij = 1..n_ode_var);
end if;

G_element||i := [element[m], op(G_element||i)];
L_element||i := [M_element[m], op(L_element||i)];
coe_element||i := [c_element[m], op(coe_element||i)];
end do;

for mm from 1 to n_dim do
indet_h||i := {CoeffMulvar(ode[mm], [ode_var],
map("-", [sumindex], 1)),
op(indet_h||i)};
eq_h||i := {CoeffMulvar(V2||i[mm], [ode_var],
map("-", [sumindex], 1)),
op(eq_dh||i)};
eq_dh||i := {CoeffMulvar(map(diff, V2||i[mm], theta),
[ode_var], map("-", [sumindex], 1)),
op(eq_dh||i)};
end do;
end if;
end do;

# ---- 求解 h_i 中的常数项 ----
if i <> n_order then
C_constant := {seq(cat(`_C`, i_constant), i_constant = 1..id*n_dim)};
V2||i := subs(convert(dsolve(indet_h||i, eq_h||i), int),
map(eval, V2||i));

L_0 := subs([seq(seq(xt[ii,jj] = subs(theta = -tau[jj-1], V2||i[ii]),
jj = 1..n_delay+1), ii = 1..n_dim)], L||0);

DV := subs([seq(seq(xt[ii,jj] = subs(theta = -tau[jj-1],
subs(dsolve(indet_h||i, eq_h||i),
subs(solve(indet_h||i, eq_dh||i),
map(diff, V2||i, theta)))[ii]),
jj = 1..n_delay+1), ii = 1..n_dim)], XT[1]);

I_ode[i] := K[i] + L_0 - DV;

for i_2 from 0 to (1+i)^n_ode_var-1 do
sumindex_2 := NULL;
k_2 := i_2;
for j_2 from 1 to n_ode_var do
sumindex_2 := sumindex_2, modp(k_2, i+1)+1;
k_2 := floor(k_2/(i+1));
end do;
SS2[sumindex_2] := convert([sumindex_2], "+") - (n_ode_var);

if SS2[sumindex_2] = i then
for mm_2 from 1 to n_dim do
indet_c||i := {CoeffMulvar(I_ode[i][mm_2], [ode_var],
map("-", [sumindex_2], 1)),
op(indet_c||i)};
end do;
end if;
end do;

V2||i := subs(solve(indet_c||i, C_constant), map(eval, V2||i));
W||i := LinearAlgebra[Multiply](Phi, V1||i + V2||i);

for i_w from 1 to n_delay+1 do
(WW||i)[i_w] := subs(theta = -tau[i_w-1], W||i);
end do;

subs_group_L[i] := {};
ix := 1;
for i_s from 1 to n_dim do
for j_s from 1 to n_delay+1 do
subs_group_L[i] := {op(subs_group_L[i]),
XT[j_s][i_s] = (WW||i)[j_s][i_s]};
subs_group_W := {op(subs_group_W),
ww[ix][i] = (WW||i)[j_s][i_s]};
ix := ix + 1;
end do;
end do;
end if;

# ---- 输出 ----
print(G_element||i);
printf("The images of the basis of %d order under Lie bracket", i);
print(L_element||i);
printf("Resonant terms of %d order", i);
print(Resonant||i);
printf("Unresonant terms of %d order", i);
print(UResonat||i);
end do;

# ---------- 汇总 ----------
V_1 := add(V1||i, i = 2..n_order);
V_2 := add(V2||i, i = 2..n_order-1);
SYS1 := subs(V[1] = 1,
(LinearAlgebra[Multiply](B, Y) + sum(g1[z], z = 2..3)));

printf("*********************************\n");
printf("The critical values needed for the occurance of a Hopf bifurcation\n");
printf("*********************************\n");
print(SC_subs);

printf("*********************************\n");
printf("Bases of center space \n");
printf("*********************************\n");
print("Phi" = Phi); print("psi" = P_si);

printf("*********************************\n");
printf("nonlinear transformation V1\n");
printf("*********************************\n");
print(V_1);

printf("*********************************\n");
printf("nonlinear transformation V2\n");
printf("*********************************\n");
print(V_2);

printf("*********************************\n");
printf("The normal form of the required order \n");
printf("*********************************\n");
print(map(factor, map(simplify, SYS1[1], exp)));
print(map(factor, map(simplify, SYS1[2], exp)));

printf("where\n");
print([W[1] = W1]);