张琪昌编的<分岔与混沌理论及应用>附录的后两个程序

octopussheng

看了论坛已经有了这本书附录中中心流形的计算程序，虽然有点问题，很多方程不能算，但是后面两个程序还没有经过大家的测试，先贴出来，免得大家敲程序，呵呵，很累的哦！

关于第一个中心流形计算程序，我把例4.4.1中korder取为3的时候，也会出来很有意思的结果，虽然是错误的，请有想法的朋友也试一下吧！

用别的系统，出不来结果！谁让咱的mathematica学的不咋样，刚开始接触！:@L :@L

言归正传  
附录2——5.6 计算高阶规范形的周期平均法
"The first three lines are the parameters and function determined by the system"
"The last line is the results which is saved with the filename lf in the current
directory. The first part of olf is , and the second part of olf is , wich is the
high order averaging equations."
k=The order of the high order averaging equation
w=Natural frequency of the nonlinear oscillator
f={0,The nonlinear function}
y={y1,y2}
eta={{Cos[w t],Sin[w t]},{-Sin[w t],Cos[w t]}}
x=eta.y
x1=x[[1]]
x2=x[[2]]
enta=Simplify[Inverse[eta]]
g=enta.f>>og
ff=Table[0,{k},{2}]
gg=Table[0,{k},{2}]
hh=Table[0,{k},{2}]
dh=Table[0,{k},{2},{2}]
dhf={0,0}
lf={0,0}
gg[[1]]=g/.ep->0>>og1
tt=2 Pi/w
fit={{Cos[ct],Sin[ct]},{-Sin[ct]/r,Cos[ct]/r}}
ff[[1]]=1/tt Integrate[gg[[1]],{t,0,tt}];
lf=Simplify[Expand[fit.ff[[1]]/.{y1->r Cos[ct],y2->r Sin[ct]}]];
Do[
   yy1=y1;
   yy2=y2;
   hh[[i-1]]=1/tt Integrate[tao(gg[[i-1]]-ff[[i-1]])/.t->(t+tao),{tao,0,tt}];
   dh[[i-1]]=Outer[D,hh[[i-1]],{y1,y2}];
   Do[
      yy1=yy1+ep^j hh[[j,1]]
   ,{j,1,i-1}];
   Do[
      yy2=yy2+ep^j hh[[j,2]]
   ,{j,1,i-1}];
   Do[
      dhf=dhf+dh[[i-j]].ff[[j]]
   ,{j,1,i-1}];
   gg[]=(1/(i-1)! D[g/.{y1->yy1,y2->yy2},{ep,i-1}]/.ep->0)-dhf;
   ff[]=1/tt Integrate[gg[],{t,0,tt}];
   lf=lf+ep^(i-1) Simplify[Expand[fit.ff[]/.{y1->rCos[ct],y2->rSin[ct]}]]
,{i,2,k}]
ep lf>>olf

[ 本帖最后由 无水1324 于 2007-7-17 11:42 编辑 ]

octopussheng

5.9 计算半单系统规范形的通用程序方法


Print["***********************************************************"];
Print["               Read the input functions                    "];
Print["***********************************************************"];
SetDirectory["path"];
ReadList["file"];
Print["***********************************************************"];
Print["Transform the real functions form to complex functions form"];
Print["***********************************************************"];
epsilon=\[Epsilon];
FA=Table[0,{i,n}];
F=Table[f,{i,1,n}];
X=Table[x,{i,1,n}];
U=Table[u,{i,1,n}];
Jacobi[funs_List,vars_List]:=Outer[D,funs,vars]
A=Jacobi[F,X];;
Do[
   A=A/.{x[k]->0};
,{k,1,n}]
B=Eigenvalues[A];
Table[bb[k]=0,{k,1,n}];
Do[
   If[Im[B[[k]]]!=0,
        bb[k]=1/2;,
        bb[k]=1;];
,{k,1,n}];
BC=DiagonalMatrix[Table[bb[k],{k,1,n}]];
EV=Eigenvectors[A];
EV=Transpose[EV].BC;
IEV=Inverse[EV];
X=EV.U;
Do[
   x=X[]
,{i,n}]
F=IEV.F;
J0=DiagonalMatrix;
F=F-J0.U;
Do[
   FFI=CoefficientList[F[],\[Epsilon]];
   le=Length[FFI];
   Do[
      fci[i,k1]=Part[FFI,k1];
      Do[
         fci[i,k1]=fci[i,k1]/.{u[k]->v[k]};
      ,{k,1,n}];
   ,{k1,1,le}];
,{i,1,n}];
Print["***********************************************************"];
Print["   Compute the normal form and nonlinear transformation    "];
Print["***********************************************************"];
Co3_:=D[b[k,nn],{u,nk}]/nk!
ne[k]_:=j-Sum[nk,{i,1,k-1}];
Table[ne=0,{i,n,8}];
Table[nk=0,{i,1,n}];
Do[TP[j]=Table[0,{i,1,n}],{j,0,norder-1}];
Do[TNF[j]=Table[0,{i,1,n}],{j,0,norder-1}];
Do[
   If[j>2,
      Do[
         Do[
            fc[k,k1]=fci[k,k1];
            If[k1<j,
               Do[vv[k]=u[k]+Sum[\[Epsilon]^j1* TP[j1][[k]],{j1,1,j-k1}],{k,1,n}];
               Do[fc[k,k1]=fc[k,k1]/.{v->vv};,{i,1,n}];
               Do[fc[k,k1]=fc[k,k1]/.{v->vv};,{i,1,n}];
              ];
         ,{k1,2,le[k]}];
         FF[k]=Table[fc[k,k1],{k1,2,le[k]}];
         EE=Table[If[i<j,\[Epsilon]^i,0],{i,1,le[k]-1}];
         f[k]=FF[k].EE;
      ,{k,1,n}];
      F=Table[f,{i,1,n}];
      Do[FF=CofficientList[F[],\[Epsilon]],{i,1,n}];
      Do[FF=FFI,{i,1,n}];
   ];
   Do[le1[k]=Length[FF[k]],{k,1,n}];
   Do[If[j>le1[k],g[k]=0;,g[k]=Part[FF[k],j]];,{k,1,n}];
   FA=Sum[Jacobi[TP[j-1-j1],U].TNF[j1],{j1,1,j-2}];
   If[FA==0,FA=Table[0,{i,n}];];
   Table[p[k]=0,{k,n}];
   Table[nf[k]=0,{k,n}];
   Do[
      Do[
         Do[
            Do[
               Do[
                  Do[
                     Do[
                        Do[
                           nk[n]=j-Sum[nk,{i,1,n-1}];
                           Do[
                              nn=Sum[nk*10^(n-i),{i,1,n}]+10^n;
                              b[k,nn]=g[k]-FA[[k]];
                              Do[b[k,nn]=Co3,{i,1,n}];
                              \[Lambda]=Sum[B[]*nk,{i,1,n}]-B[[k]];
                              If[
                                 \[Lambda]==0,
                                 un=b[k,nn];
                                 hp=0;
                                 hp=b[k,nn]/\[Lambda];
                                 un=0;
                                 ];
                              xs=Product[u^nk,{i,1,n}];
                              p[k]=p[k]+hp*xs;
                              nf[k]=nf[k]+un*xs;
                           ,{k,1,n}];
                        ,{nk[8],0,ne[8]}];
                     ,{nk[7],0,ne[7]}];
                  ,{nk[6],0,ne[6]}];
               ,{nk[5],0,ne[5]}];
            ,{nk[4],0,ne[4]}];
         ,{nk[3],0,ne[3]}];
      ,{nk[2],0,ne[2]}];
   ,{nk[1],0,j}];
   TP[j-1]=Table[p[k],{k,1,n}];
   TNF[j-1]=Table[nf[k],{k,1,n}];
,{j,2,norder}];

"************************************************************"
"             The normal form in complex form               ";
"**********************************************************" ;
"****************************************************">>d.txt;
"         The normal form in complex form           ">>>d.txt;
"***************************************************">>>d.txt;
Do[
   tf[k]=Sum[\[Epsilon]^jTNF[j],{j,1,norder-1}][[k]]
,{k,1,n}];
Do[
   "D[u["<>ToString<>"],t]=">>>d.txt;
   tf=B[]*u+tf>>>d.txt;
,{i,1,n}];

"***************************************************">>>d.txt;
"      Transform back to system in real form        ">>>d.txt;
"***************************************************">>>d.txt;
Clear[F,f,X,x,Y,y,v];
X=Table[x,{i,1,n}];
F=Simplify[Table[tf,{i,1,n}]];
Do[F=F/.{u[k]->y[k]};,{k,1,n}];
Y=IEV.X;
Do[y[k]=Y[[k]],{k,1,n}];
"D[Y,t]=">>>d.txt
F1=Simplify[EV.F]>>>d.txt

"***************************************************">>>d.txt;
"      Transform back to system in Polar form       ">>>d.txt;
"***************************************************">>>d.txt;
conjug=Complex[a_,b]_:>Complex[a,-b];
nz=0;ni=0;cnz=0;cni=0;cm=0;
Do[
   If[Re[B[[k]]]! =0,
      If[Im[B[[k]]]==0,cnz=cnz+1;];,
      cm=cm+1;
      If[Im[B[[k]]]==0,nz=nz+1;];
   ];
,{k,n}];
sm=n-cm;
ni=(cm-nz)/2;
cni=(cm-cnz)/2;
k2=0;k1=0;
Do[
   If[Im[B[[k]]]>0,
       k2=k1+1;
       v[k]=r[k2] Exp[I \[Theta][k2]];,
       If[Im[B[[k]]]==0,
           v[k]=u[k];,
           k1=k1+1;
           v[k]=r[k1] Exp[-I \[Theta][k1]];
       ];
   ];
,{k,n}];
k1=0;k2=0;nzn=nz+cnz;
Do[
   If[Im[B[[k]]]==0,
       ga=tf[k];
       Do[ga=ga/.{u->v};,{i,1,n}];
       k2=k2+1;
       "D[u["<>ToString[k2]<>"],t]=">>>d.txt;
       g[k2]=Expand[TrigReduce[ComplexExpand[ga]]]>>>d.txt;
       If[Im[B[[k]]]<0,
           ga=tf[k];
           Do[ga=ga/.{u->v};,{i,1,n}];
           k1=k1+1;
           gc=ga Exp[I \[Theta][k1]];
           gd=gc/.conjug;
           "D[r["<>ToString[k1]<>"],t]=">>>d.txt;
           g[nzn+2 k1-1]=Expand[TrigReduce[ComplexExpand[Simplify[(gc+gd)/2]]]]>>>d.txt;
           "D[\[Theta]["<>ToString[k1]<>"],t]=">>>d.txt;
           g[nzn+2 k1]=Expand[TrigReduce[ComplexExpand[Simplify[-(gc-gd)/(2*I*r[k1])]]]]>>>d.txt;
       ];
   ];
,{k,n}];

[ 本帖最后由 无水1324 于 2007-7-17 11:43 编辑 ]

无水1324

与另外的帖子重复没有？

octopussheng

没有，我看那个只有第一个中心流形的计算程序，不放心无水可以检查一下，呵呵！

无水1324

原帖由 octopussheng 于 2007-7-17 13:43 发表

                               
登录/注册后可看大图

没有，我看那个只有第一个中心流形的计算程序，不放心无水可以检查一下，呵呵！

一直没有时间做这个问题，最近在用奇异性理论分析，遇到很多困难，看书还是模糊的

octopussheng

奇异性理论我就看过一遍，因为暂时没有用上，所以就没有太用心去看，呵呵！

中心流形、规范形理论倒是看的仔细点！

不过现在规范形还是没有用好！还在学习！