maple 程序，编译，运行
来源：互联网 发布：淘宝签收后如何退款 编辑：程序博客网 时间：2024/06/05 11:34
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])); 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
3 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
4 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
5 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
6 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
7 DiffMulvar := proc (Eqs::{Matrix, Vector[column]}, varlist::list, ww, norder::integter, 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
8 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
9 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:=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(lanbda=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[Matrix](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:=[]; #IIndex||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 (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; 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 o 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 V1\\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]);
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:=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(lanbda=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[Matrix](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:=[];
#IIndex||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[jj1],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 (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;
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 o 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 V1\\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]);
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:=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(lanbda=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[Matrix](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:=[];
#IIndex||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[jj1],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 (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;
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 o 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 V1\\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]);
- maple 程序，编译，运行
- MAPLE
- java程序编译运行
- 程序编译运行相关
- Scala程序编译运行
- 编译运行java程序
- 编译和运行Java程序
- 手工编译运行java程序
- 动态编译、运行JAVA程序
- 动态编译，运行java程序
- java 命令行编译运行程序
- cuda程序的编译运行
- cuda程序的编译运行
- linux编译运行C程序
- 自己动手编译、运行Java程序
- sublime3编译运行scala程序
- ubuntu编译运行java程序
- 命令行编译运行Java程序
- hibernate 多对多的增删改操作
- 教室信息管理系统总结（结构体数组）
- 很全的mvn介绍
- C语言特性
- mysql int(3)与int(11)的区别
- maple 程序，编译，运行
- WCF学习附带收获
- 使用jackson转换xml为hashmap形式
- 存minicom log
- Python进阶With语法
- 犀牛——第3章类型、值和变量3.9变量声明
- easyui datagrid中动态改变某行的颜色
- shell oracle(sqlplus spool)
- 关于MemStore的刷新方式
