##############################################################################
# ## IFB_COMP ## ## v0.1 ## Sep 07 2004 ## Martin Rasmussen, University of Augsburg, Germany ## martin.rasmussen@math.uni-augsburg.de ## ##############################################################################necessary libraries:with(linalg): with(plots):Dimensions of the system:Dim[0]:=1: # dimension of the pseudo-stable subspace
Dim[1]:=1: # dimension of the pseudo-unstable subspace
Linear parts:A:=k->Matrix(1,1,[[-0.5]]): # linearization (pseudo-stable)
B:=k->Matrix(1,1,[[-1.0]]): # linearization (pseudo-unstable)Nonlinear part:F1:=(x1,x2,y1,y2,k)->y1^2: # pseudo-stable nonlinear part (1st component)
F2:=(x1,x2,y1,y2,k)->.0: # pseudo-stable nonlinear part (2nd component) (if needed)
G1:=(x1,x2,y1,y2,k)->x1*y1: # pseudo-unstable nonlinear part (1st component)
G2:=(x1,x2,y1,y2,k)->.0: # pseudo-unstable nonlinear part (2nd component) (if needed)Coefficient sequences (if needed):cc[1]:=k->0.92+0.7*(2/Pi)*arctan(k); # coefficient sequences
cc[2]:=k->0.9+0.2*(2/Pi)*arctan(k);
cc[3]:=k->0.178+0.1*(2/Pi)*arctan(k);The following series of procedures has to be executed...graphs:=proc()
global c,d,e,graphf,graphg,graphh;
graphf:=(wb,k,x,y)->1/2*c(2,0,wb,k)*x^2+c(1,1,wb,k)*x*y+1/2*c(0,2,wb,k)*y^2+1/6*c(3,0,wb,k)*x^3+1/2*c(2,1,wb,k)*x^2*y+1/2*c(1,2,wb,k)*x*y^2+1/6*c(0,3,wb,k)*y^3+1/24*c(4,0,wb,k)*x^4+1/6*c(3,1,wb,k)*x^3*y+1/4*c(2,2,wb,k)*x^2*y^2+1/6*c(1,3,wb,k)*x*y^3+1/24*c(0,4,wb,k)*y^4+1/120*c(5,0,wb,k)*x^5+1/24*c(4,1,wb,k)*x^4*y+1/12*c(3,2,wb,k)*y^2*x^3+1/12*c(2,3,wb,k)*y^3*x^2+1/24*c(1,4,wb,k)*y^4*x+1/120*c(0,5,wb,k)*y^5+1/720*c(6,0,wb,k)*x^6+1/120*c(5,1,wb,k)*x^5*y+1/48*c(4,2,wb,k)*y^2*x^4+1/36*c(3,3,wb,k)*y^3*x^3+1/48*c(2,4,wb,k)*y^4*x^2+1/120*c(1,5,wb,k)*y^5*x+1/720*c(0,6,wb,k)*y^6+1/5040*c(7,0,wb,k)*x^7+1/720*c(6,1,wb,k)*x^6*y+1/240*c(5,2,wb,k)*x^5*y^2+1/144*c(4,3,wb,k)*y^3*x^4+1/144*c(3,4,wb,k)*y^4*x^3+1/240*c(2,5,wb,k)*y^5*x^2+1/720*c(1,6,wb,k)*y^6*x+1/5040*c(0,7,wb,k)*y^7;
graphg:=(wb,k,x)->1/2*d(2,wb,k)*x^2+1/6*d(3,wb,k)*x^3+1/24*d(4,wb,k)*x^4+1/120*d(5,wb,k)*x^5+1/720*d(6,wb,k)*x^6+1/(5040)*d(7,wb,k)*x^7;
graphh:=(wb,k,x)->1/2*e(2,wb,k)*x^2+1/6*e(3,wb,k)*x^3+1/24*e(4,wb,k)*x^4+1/120*e(5,wb,k)*x^5+1/720*e(6,wb,k)*x^6+1/(5040)*e(7,wb,k)*x^7;
end proc:CalcStart:=proc(km,gen,ord,wb)::integer;
local tmp;
if wb=0 then
return km;
else
tmp:=km-gen*(ord-1);
return tmp;
end if;
end proc:CalcStopp:=proc(kp,gen,ord,wb)::integer;
local tmp;
if wb=0 then
tmp:=kp+gen*(ord-1);
return tmp;
else
tmp:=kp+ord-2;
return tmp;
end if;
end proc:count:=proc(N,num,n)
local i,anz;
anz:=0;
for i from 1 by 1 to n do
if N[i]=num then
anz:=anz+1;
end if;
end do;
return anz;
end proc:valid1:=proc(N,ord)
local count::integer;
if (N[2]>0 and Dim[0]=1) or (N[4]>0 and Dim[1]=1) then #0 und 1 vertauscht
return 0;
end if;
count := N[1]+N[2]+N[3]+N[4];
if (count > ord) or (count < 2) then
return 0;
else
return 1;
end if;
end proc:valid2:=proc(N,l,n)
local i,cur;
cur:=l+1;
for i from n by -1 to 1 do
if N[i]<cur-1 then
return 0;
elif N[i]=cur-1 then
cur:=cur-1;
end if;
end do;
if cur>1 then
return 0;
else
return 1;
end if;
end proc:Differ:=proc(kmx::integer,kpx::integer,ord::integer)
local i,k,N,st;
global DF,DG;
for i from 1 by 1 to 4 do
N[i]:=0;
end do;
st:=4;
while st>0 do
st:=4;
while N[st]=ord do
N[st]:=0;
st:=st-1;
end do;
N[st]:=N[st]+1;
if valid1(N,ord)=1 and st>0 then
for k from kmx by 1 to kpx do
DF(N[1],N[2],N[3],N[4],0,k):=evalf(subs(x1=0,x2=0,y1=0,y2=0,diff(F1(x1,x2,y1,y2,k),x1$N[1],x2$N[2],y1$N[3],y2$N[4])));
DG(N[1],N[2],N[3],N[4],0,k):=evalf(subs(x1=0,x2=0,y1=0,y2=0,diff(G1(x1,x2,y1,y2,k),x1$N[1],x2$N[2],y1$N[3],y2$N[4])));
if Dim[0]=2 then
DF(N[1],N[2],N[3],N[4],1,k):=evalf(subs(x1=0,x2=0,y1=0,y2=0,diff(F2(x1,x2,y1,y2,k),x1$N[1],x2$N[2],y1$N[3],y2$N[4])));
end if;
if Dim[1]=2 then
DG(N[1],N[2],N[3],N[4],1,k):=evalf(subs(x1=0,x2=0,y1=0,y2=0,diff(G2(x1,x2,y1,y2,k),x1$N[1],x2$N[2],y1$N[3],y2$N[4])));
end if;
end do;
end if;
end do;
end proc:Init:=proc(kmx,kpx,wb)
local k;
global Spm,gpm;
for k from kmx by 1 to kpx do #Initialisierung Spm #Ordnung,Basis,Koordinate,Zeit
Spm(1,0,0,k):=0;
if Dim[wb] = 2 then
Spm(1,1,0,k):=0;
elif Dim[1-wb] = 2 then
Spm(1,0,1,k):=0;
end if;
end do;
for k from kmx by 1 to kpx do #Initialisierung gpm #Ordnung,Basis,Koordinate,Zeit
if wb=0 then
if Dim[wb]=1 then
gpm(1,0,0,k):=A(k)[1,1];
else
gpm(1,0,0,k):=A(k)[1,1];
gpm(1,1,0,k):=A(k)[1,2];
gpm(1,0,1,k):=A(k)[2,1];
gpm(1,1,1,k):=A(k)[2,2];
end if;
else
if Dim[wb]=1 then
gpm(1,0,0,k):=B(k)[1,1];
else
gpm(1,0,0,k):=B(k)[1,1];
gpm(1,1,0,k):=B(k)[1,2];
gpm(1,0,1,k):=B(k)[2,1];
gpm(1,1,1,k):=B(k)[2,2];
end if;
end if;
end do;
end proc:EvalMultilinear:=proc(MlForm,Laenge,Vektoren,Dimension)
local gbas,i,M,st,tmp;
tmp:=0;
for i from 1 by 1 to Laenge do
M[i]:=1;
end do;
M[0]:=0;
st:=Laenge;
while st>0 do
gbas:=0;
for i from 1 by 1 to Laenge do
if M[i]=2 then
gbas:=gbas+1;
end if;
end do;
tmp:=tmp+MlForm(gbas)*mul(Vektoren(i,M[i]),i=1..Laenge);
st:=Laenge;
while M[st]=Dimension do
M[st]:=1;
st:=st-1;
end do;
M[st]:=M[st]+1;
end do;
return tmp;
end proc:CalcG:=proc(ord,bas,koor,time,flag,wb) #flag, 0: summe bis ord, 1: summe bis ord-1 plus zus. Term
local ausg::float,anz,Dimges,Dx1,Dx2,Dy1,Dy2,gbas,i,j,l,N,M,MlForm,st,st2,tmp::float,x,y;
ausg:=0;
Dimges:=Dim[0]+Dim[1];
for l from 2 by 1 to (ord-flag) do
for i from 1 by 1 to ord do
N[i]:=1;
end do;
st:=ord;
while st>0 do
st:=ord;
while N[st]=l do
N[st]:=1;
st:=st-1;
end do;
N[st]:=N[st]+1;
if valid2(N,l,ord)=1 and st>0 then
for i from 1 to l do
anz[i]:=count(N,i,ord);
M[i]:=0; #Anzahl zweiter Basisvektoren
for j from 1 by 1 to ord do
if N[j]=i then
if bas>=j then
M[i]:=M[i]+1;
end if;
end if;
end do;
end do;
for i from 1 to l do
y(i,1):=gpm(anz[i],M[i],0,time);
if Dim[wb]=2 then
y(i,2):=gpm(anz[i],M[i],1,time);
end if;
if wb=0 then
if Dim[0]=1 then
if anz[i]=1 then
x(i,1):=1;
else
x(i,1):=0;
end if;
x(i,2):=Spm(anz[i],M[i],0,time);
if Dim[1]=2 then
x(i,3):=Spm(anz[i],M[i],1,time);
end if;
else
if anz[i]=1 then
x(i,1):=1-M[i];
x(i,2):=M[i];
else
x(i,1):=0;
x(i,2):=0;
end if;
x(i,3):=Spm(anz[i],M[i],0,time);
end if;
else
if Dim[0]=1 then
if anz[i]=1 then
if Dim[1]=1 then
x(i,2):=1;
else
x(i,2):=1-M[i];
x(i,3):=M[i];
end if;
else
x(i,2):=0;
x(i,3):=0;
end if;
x(i,1):=Spm(anz[i],M[i],0,time);
else
x(i,1):=Spm(anz[i],M[i],0,time);
x(i,2):=Spm(anz[i],M[i],1,time);
if anz[i]=1 then
x(i,3):=1;
else
x(i,3):=0;
end if;
end if;
end if;
end do;
for j from 1 by 1 to l do
M[j]:=1;
end do;
M[0]:=0;
st2:=l;
tmp:=0;
while st2>0 do
Dx1:=0; Dx2:=0; Dy1:=0; Dy2:=0;
for j from 1 by 1 to l do
if M[j]=1 then
Dx1:=Dx1+1;
elif M[j]=2 then
if Dim[0]=2 then
Dx2:=Dx2+1;
else
Dy1:=Dy1+1;
end if;
elif M[j]=3 then
if Dim[0]=2 then
Dy1:=Dy1+1;
else
Dy2:=Dy2+1;
end if;
end if;
end do;
if (flag=0 and wb=0) or (flag=1 and wb=1) then
tmp:=tmp + DF(Dx1,Dx2,Dy1,Dy2,koor,time) *mul(x(j,M[j]),j=1..l);
else
tmp:=tmp + DG(Dx1,Dx2,Dy1,Dy2,koor,time)*mul(x(j,M[j]),j=1..l);
end if;
st2:=l;
while M[st2]=Dimges do
M[st2]:=1;
st2:=st2-1;
end do;
M[st2]:=M[st2]+1;
end do;
if flag = 1 then
for j from 0 to l*(Dim[wb]-1) do
MlForm(j):=Spm(l,j,koor,time+1)
end do;
tmp:=tmp-EvalMultilinear(MlForm,l,y,Dim[wb]);
end if;
ausg:=ausg+tmp;
end if;
end do;
end do;
return ausg;
end proc:
CalcH:=proc(ord,bas,koor,time,wb)
local ausg;
if wb=0 then
ausg:=DG(ord-bas,bas,0,0,koor,time);
else
ausg:=DF(0,0,ord-bas,bas,koor,time);
end if;
ausg:=ausg + CalcG(ord,bas,koor,time,1,wb);
return ausg;
end proc:
Main:=proc(km,kp,gen,ord,wb)
#km: Startzeit
#kp: Stoppzeit
#gen: Genauigkeit
#ord: Ordnung
#wb: 0 -> pseudo-stabil, 1 -> pseudo-instabil
local Hpm,i,j,kmx::integer,kpx::integer,k,l,m,n,Matr1,Matr2,MlForm,Vektoren;
if wb=0 then
printf("Calculate pseudo-stable bundle up to order %d, init time: %d, final time: %d, precision: %d\n", ord,km,kp,gen);
else
printf("Calculate pseudo-unstable bundle of to order %d, init time: %d, final time: %d, precision: %d\n", ord,km,kp,gen);
end if;
kmx := CalcStart(km,gen,ord,wb);
kpx := CalcStopp(kp,gen,ord,wb);
graphs();
Differ(kmx,kpx,ord);
Init(kmx,kpx,wb);
for n from 2 by 1 to ord do
printf("Calculate order %d: ",n);
if n>2 then
printf("gpm,");
for k from CalcStart(km,gen,ord-n+2,wb) by 1 to CalcStopp(kp,gen,ord-n+2,wb) do
for i from 0 by 1 to n*(Dim[wb]-1) do
gpm(n-1,i,0,k):=CalcG(n-1,i,0,k,0,wb);
if Dim[wb]=2 then
gpm(n-1,i,1,k):=CalcG(n-1,i,1,k,0,wb);
end if;
end do;
end do;
end if;
printf("Hpm,");
for k from CalcStart(km,gen,ord-n+2,wb) by 1 to CalcStopp(kp,gen,ord-n+2,wb) do
for i from 0 by 1 to n*(Dim[wb]-1) do
Hpm(n,i,0,k):=CalcH(n,i,0,k,wb);
if Dim[1-wb]=2 then
Hpm(n,i,1,k):=CalcH(n,i,1,k,wb);
end if;
end do;
end do;
#Berechnung von Spm
printf("Spm,");
for k from CalcStart(km,gen,ord-n+2,wb) by 1 to CalcStopp(kp,gen,ord-n+2,wb) do
for i from 0 by 1 to n*(Dim[wb]-1) do
Spm(n,i,0,k):=0;
Spm(n,i,1,k):=0;
if wb=0 then
for j from k+gen by -1 to k do
if Dim[0]=1 then Matr1:=Matrix(1,1,[[1]]); else Matr1:=Matrix(2,2,[[1,0],[0,1]]); end if;
Matr2:=inverse(B(j));
for l from j-1 by -1 to k do
Matr1:=multiply(Matr1,A(l));
Matr2:=multiply(inverse(B(l)),Matr2);
end do;
#Initialisierung der Vektoren
for l from 1 by 1 to n do
if l<=i then #nur falls Dim[wb]=Dim[0]=2
Vektoren(l,1):=Matr1[1,2];
Vektoren(l,2):=Matr1[2,2];
else
Vektoren(l,1):=Matr1[1,1];
if Dim[wb]=2 then
Vektoren(l,2):=Matr1[2,1];
end if;
end if;
end do;
if Dim[1-wb]=1 then
Matr1:=Matrix(1,1,[[0]]);
else
Matr1:=Matrix(2,1,[[0],[0]]);
end if;
for l from 0 by 1 to Dim[1-wb]-1 do #l bezeichnet Koordinate
#Initialisierung der Multilinearform
for m from 0 by 1 to n*(Dim[wb]-1) do
MlForm(m):=Hpm(n,m,l,k);
end do;
Matr1[l+1,1]:=EvalMultilinear(MlForm,n,Vektoren,Dim[wb]);
end do;
Matr2:=multiply(Matr2,Matr1);
Spm(n,i,0,k):=evalf(Spm(n,i,0,k)-Matr2[1,1]);
if Dim[1-wb]=2 then
Spm(n,i,1,k):=evalf(Spm(n,i,1,k)-Matr2[2,1]);
end if;
end do;
else
for j from k-gen by 1 to k-1 do
if Dim[0]=1 then Matr1:=Matrix(1,1,[[1]]); else Matr1:=Matrix(2,2,[[1,0],[0,1]]); end if;
Matr2:=inverse(B(k-1));
for l from k-1 by -1 to j+1 do
Matr1:=multiply(Matr1,A(l));
Matr2:=multiply(inverse(B(l-1)),Matr2);
end do;
#Initialisierung der Vektoren
for l from 1 by 1 to n do
if l<=i then #nur falls Dim[wb]=Dim[1]=2
Vektoren(l,1):=Matr2[1,2];
Vektoren(l,2):=Matr2[2,2];
else
Vektoren(l,1):=Matr2[1,1];
if Dim[wb]=2 then
Vektoren(l,2):=Matr2[2,1];
end if;
end if;
end do;
if Dim[1-wb]=1 then
Matr2:=Matrix(1,1,[[0]]);
else
Matr2:=Matrix(2,1,[[0],[0]]);
end if;
for l from 0 by 1 to Dim[1-wb]-1 do #l bezeichnet Koordinate
#Initialisierung der Multilinearform
for m from 0 by 1 to n*(Dim[wb]-1) do
MlForm(m):=Hpm(n,m,l,k);
end do;
Matr2[l+1,1]:=EvalMultilinear(MlForm,n,Vektoren,Dim[wb]);
end do;
Matr1:=multiply(Matr1,Matr2);
Spm(n,i,0,k):=evalf(Spm(n,i,0,k)+Matr1[1,1]);
if Dim[1-wb]=2 then
Spm(n,i,1,k):=evalf(Spm(n,i,1,k)+Matr1[2,1]);
end if;
end do;
end if;
end do;
end do;
printf("\n");
end do;
for k from km by 1 to kp do
if Dim[wb]=2 then
for i from 2 by 1 to 7 do ##7=maximale Ordnung
for j from 0 by 1 to i do
if i>ord then
c(i-j,j,wb,k):=0;
else
c(i-j,j,wb,k):=evalf(Spm(i,j,0,k));
end if;
end do;
end do;
else
for i from 2 by 1 to 7 do ##7=maximale Ordnung
if i>ord then
d(i,wb,k):=0;
e(i,wb,k):=0;
else
d(i,wb,k):=evalf(Spm(i,0,0,k));
if Dim[1-wb]=2 then
e(i,wb,k):=evalf(Spm(i,0,1,k));
end if;
end if;
end do;
end if;
end do;
end proc:
Output:=proc(wb,time,koorm,koorp,gen)
local i,j,points;
if Dim[wb]=2 then
print(graphf(wb,time,x,y));
plot3d(graphf(wb,time,x,y),x=koorm..koorp,y=koorm..koorp, axes=BOXED);#,view=[koorm..koorp,koorm..koorp,-1..1]);
else
print(graphg(wb,time,x));
if Dim[1-wb]=1 then
plot(graphg(wb,time,x), x=koorm..koorp);
else
print(graphh(wb,time,x));
points:={seq([koorm+(koorp-koorm)*x/gen, graphg(wb,time,koorm+(koorp-koorm)*x/gen), graphh(wb,time,koorm+(koorp-koorm)*x/gen)], x=0..gen)};
pointplot3d(points,axes=NORMAL);
end if;
end if;
end proc:Calculate the pseudo-stable bundle:Main(0,0,10,6,0); # arguments: km,kp,prec,ord,wbOutput the pseudo-stable bundleOutput(0,0,-0.25,0.25,1000); #Argumente: wb,time,koorm,koorp,precCalculate the pseudo-unstable bundle:Main(0,0,10,6,1); # arguments: km,kp,prec,ord,wbOutput the pseudo-unstable bundleOutput(1,0,-0.25,0.25,1000); #Argumente: wb,time,koorm,koorp,gen