★ FEATURED ARTICLE ####################################################################### Neutral DDE 高维 Hopf 分岔 / 规范形计算脚本合并整理版## 说明# - 原文主程序块重复三次此处仅保留一份。# - 已修正的明显笔误# lanbda - lambda# linearAlgebra - LinearAlgebra# o - 0for 循环下界# for i_1 from 0 (1i) - from 0 to (1i)# £» - ;# 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 thenvarseq : seq(seq(xt[i,j], j 1..ndelay1), i 1..ndim), op(vpar)elif dtype neutral thenvarseq : seq(seq(xt[i,j], j 1..ndelay1), i 1..2*ndim), op(vpar)elseerror %1,type of DDE not recognized, ddetypeend if;varseqend 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 donterm : nops(expand(Eqs_temp[i]1))-1;for j to nterm doterm : op(j, expand(Eqs_temp[i]1));if degree(term, varlist) Norder thenEq_ex[i] : Eq_ex[i] termend ifend doend do;Eq_exend 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 dofor j to nops(expand(Eqs[i])1)-1 doterm : op(j, expand(Eqs[i])1);if 1 degree(term, varlist) thenEq_G[i] : Eq_G[i] termend ifend doend do;Eq_F : Eqs - Eq_G;Eq_F, Eq_Gend 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 dofor j to nops(expand(Eqs[i])1)-1 doterm : op(j, expand(Eqs[i])1);if 1 degree(term, varpar) thenif degree(term, svarlist) 1 thenEq_ex[i] : Eq_ex[i] termend ifend ifend doend do;Eq_ex, Eqs - Eq_exend proc:CoeffMulvar : proc(expression::algebraic,varlist::list, exponentslist::list)local nvar, eq, i;nvar : nops(varlist);eq : expression;for i to nvar doeq : coeff(eq, varlist[i], exponentslist[i])end do;eqend 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 (norder1)^nvar-1 dosummands : NULL;kk : j;for m to nvar dosummands : summands, modp(kk, norder1);kk : floor(kk/(norder1))end do;S[summands] : convert([summands], );if S[summands] norder thenT : combinat:-carprod([$([0, seq(i, i 2..ni)], nvar)]);while not T[finished] doidx : T[nextvalue]();if convert(idx, ) ninorder-neqsand sum(idx[nn]*summands[nn], nn 1..nvar) ninorder-neqs thenfor ii to rdim dofor cc to cdim doDEQs[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 doend doend ifend doend ifend do;DEQsend proc:SigSum : proc(eq, ndim::integer)local r_eq;if eq 0 thenr_eq : Vector[column](ndim)elser_eq : eqend if;r_eqend proc:# 3. 主计算流程 # ---------- 系统 Taylor 展开 ----------sys_rhs : map(mtaylor, sys_rhs,[seq(seq(op([xt[i,j] Eq_pos[i], xt[n_dimi,j] Eq_pos[i]]),i 1..n_dim), j 1..(n_delay1))],n_order1);n_ode : 2;# ---------- 变量序列 ----------var : VarSeq(n_dim, n_delay, var_par, dde_type);svar : seq(seq(xt[i,j], j 1..n_delay1), i 1..n_dim);dvar : seq(seq(xt[i,j], j 1..n_delay1), i (n_dim1)..2*n_dim);n_var : nops([var]);n_svar : nops([svar]);n_dvar : nops([dvar]);for i from 1 to n_delay1 doXT[i] : seq(xt[j,i], j 1..n_dim);end do;for i from 1 to n_delay1 doNXT[i] : seq(xt[j,i], j n_dim1..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 docat(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 dofor jj from 2 to n_order docat(h,jj,ii) : 0;end do;end do;# ---------- 齐次多项式生成 ----------S : table(): aa : 0;for j from 0 to (n_order1)^n_ode_var-1 dosummands : NULL;kk : j;for m from 1 to n_ode_var dosummands : summands, modp(kk, n_order1);kk : floor(kk/(n_order1));end do;S[summands] : convert([summands], Matrix);for k from 2 to n_order doif S[summands] k thenfor jj from 1 to n_dim doaa : 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 doE_subs : [op(E_subs), seq(xt[i,j] Eq_pos[i], j 1..n_delay)];end do;for i from 0 to n_delay doA[i] : convert(subs(E_subs, linalg[jacobian](L||0, XT[i1])), Matrix);B[i] : convert(subs(E_subs, linalg[jacobian](sys_D_lhs, XT[i1])), 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 doPhiyt[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 dofor j from 1 to n_dim dosubs_group : {op(subs_group), XT[i1][j] Phiyt[i][j]};subs_group_h : {op(subs_group_t), XT[i1][j] Phiyth[i][j]};end do;end do;subs_group_W : {};linalg[jacobian](L||0, XT[i1], Matrix);# ---------- 逐阶计算 ----------for i from 2 to n_order doG_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 doif dde_type neutral thent1 : 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_delay1), n_dim) SigSum(add(LinearAlgebra[Multiply](add(add(subs(subs_group_W, subs_group,DiffMulvar(VectorCalculus[Jacobian](Gx||(k1),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_delay1), n_dim) SigSum(add(add(LinearAlgebra[Multiply]((VectorCalculus[Jacobian](Gx||(I1),convert(NXT[jj], list)) add(add(subs(subs_group_W, subs_group,DiffMulvar(VectorCalculus[Jacobian](Gx||(k1),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-k1]),k 2..i-I-1), n_dim)))),jj 1..n_delay1), I 1..(i-2)), n_dim);K[i] : t1 t2 t3;elif dde_type retarded thent1 : 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[i1-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[i1-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 (1i)^n_ode_var-1 dosumindex : NULL;k : i_1;for j from 1 to n_ode_var dosumindex : sumindex, modp(k, i1)1;k : floor(k/(i1));end do;SS[sumindex] : convert([sumindex], ) - (n_ode_var);if SS[sumindex] i thenid : id 1;for m from 1 to n_ode doii : 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) thenResonant||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);elseUResonat||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 doindet_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 thenC_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_delay1), 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_delay1), ii 1..n_dim)], XT[1]);I_ode[i] : K[i] L_0 - DV;for i_2 from 0 to (1i)^n_ode_var-1 dosumindex_2 : NULL;k_2 : i_2;for j_2 from 1 to n_ode_var dosumindex_2 : sumindex_2, modp(k_2, i1)1;k_2 : floor(k_2/(i1));end do;SS2[sumindex_2] : convert([sumindex_2], ) - (n_ode_var);if SS2[sumindex_2] i thenfor mm_2 from 1 to n_dim doindet_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_delay1 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 dofor j_s from 1 to n_delay1 dosubs_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]);
阅读完成 · 觉得有帮助？
