Изучение методов нахождения спектральных характеристик

17. Садовничий В.А. Теория операторов / Садовничий В. А. - М.: Высш. шк.,1999. - 368 с.

18. Смирнов В. И. Научное творчество Алексея Николаевича Крылова, УМН, Т 1. - вып. 3-4(13-14), 1946.- С. 3-12

19. Срочко В. А. Численные методы: Курс лекций. Иркутск: Иркут. ун-т, 2003.-168 с.

20. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры / Фаддеев Д.К., Фаддеева В.Н. - Л.: Наука, 1975. - 656 с.



Приложение

restart:

Digits:=8:

with(linalg):

assume(n,integer);

N:=5:

Hmat:=array(symmetric,sparse,1..N,1..N);

m:=1: w:=1: a:=2: b:=1/2: c:=5:

Hmat := array(symmetric, sparse, 1 .. 5, 1 .. 5, [])

V:=(zeta) -> m*w^2/2*zeta^2+a*zeta^4+b*zeta^6;

V12:=(x,y)->V(x)+V(y)+c*x*y;

H:=(arg,x,y)->-1/2*diff(arg,x$2)-1/2*diff(arg,y$2)+V12(x,y)*arg;

alpha:=2: psi[0]:=x->exp(-alpha*x^2);

psi[1]:=(x,y)->psi[0](x)*psi[0](y);

with(student,Doubleint):

NormXY:=proc();

value(Doubleint(args[1],x=-infinity..infinity,y=-infinity..infinity));

end:

N1:=sqrt(NormXY(psi[1](x,y)^2));

phi[1]:=unapply(psi[1](x,y)/(N1),x,y);

psi[2]:=unapply(simplify( H(phi[1](x,y),x,y)),x,y):

Hd:=NormXY( phi[1](x,y)*psi[2](x,y) ):

Hmat[1,1]:=evalf(Hd);

xi[2]:=unapply(simplify(psi[2](x,y)-Hd*phi[1](x,y)),x,y):

N2:=sqrt(NormXY( xi[2](x,y)^2 )):

phi[2]:=unapply(simplify(xi[2](x,y)/N2),x,y):

Hmat[1,2]:=evalf(N2);

for n from 2 to N-1 do

psi[n+1]:=unapply(simplify(H(phi[n](x,y),x,y)),x,y);

Hd:=NormXY( phi[n](x,y)*psi[n+1](x,y) );

Hmat[n,n]:=Hd;

xi[n+1]:=unapply(simplify(psi[n+1](x,y)-Hd*phi[n](x,y)-Hmat[n-1,n]*phi[n-1](x,y)),x,y);

N2:=sqrt(NormXY( xi[n+1](x,y)^2 ));

Hmat[n,n+1]:=N2;

phi[n+1]:=unapply(simplify(xi[n+1](x,y)/N2),x,y);

od:

psi[N+1]:=unapply(simplify(H(phi[N](x,y),x,y)),x,y):

Hd:=NormXY(phi[N](x,y)*psi[N+1](x,y)):

Hmat[N,N]:=Hd:

A:=matrix(N,N):

A:=evalm(map(evalf,Hmat));

A1:=matrix(N,N):

E:=matrix(N,N,[[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,0,1]]):

B:=matrix(N,N):

p:=array(1..N):

b:=matrix(N,N):

y0:=matrix(N,1,[1,0,0,0,0]):

y:=matrix(N,1):

xL:=matrix(N,N):

A1:=A:

for j from 1 to N do

P:=0;

for i from 1 to N do

P:=P+A1[i,i]:

od:

p[j]:=(1/j)*P:

B:=evalm(A1-p[j]*E):

for i from 1 to N do

b[i,j]:=B[i,1]:

od:

A1:=evalm(evalf(A&*B)):

od:

sys:=lambda^N:

for i from 1 to N do

sys:=sys-p[i]*lambda^(N-i);

od:

ch:=evalf(solve(sys=0,lambda)):

print(sys);

print(ch);

1.8030133, 6.1059782, 14.415994, 52.898308, 202.91829

for k from 1 to N do

y:=y0:

for j from 1 to N-1 do

y:=evalm(ch[k]*y):

for i from 1 to N do

y[i,1]:=y[i,1]+b[i,j];

od: od:

for i from 1 to N do

xL[i,k]:=evalm(y[i,1]): od:

od:

print(xL);



for j from 1 to N do

s:=0:

for i from 1 to N do

s:=s+(xL[i,j])^2:od:

nn:=1/sqrt(s):

for i from 1 to N do

xL[i,j]:=evalf(nn*xL[i,j]):

od: od:

print(xL):

q:=matrix(N,N):

xK:=matrix(N,N):

y:=matrix(N,1,[1,0,0,0,0]):

for i from 0 to N-1 do

for j from 1 to N do

A1[j,N-i]:=y[j,1]:

od:

y:=evalm(A&*y):

od:

y:=evalm((-1)*y):

p:=linsolve(A1,y):

sys:=lambda^N:

for i from 1 to N do

sys:=sys+p[i,1]*lambda^(N-i):

od:

ch:=evalf(solve(sys=0,lambda));

ch := 1.8023796, 6.1070927, 14.415473, 52.898342, 202.91828

for i from 1 to N do

q[1,i]:=1:

od:

for j from 2 to N do

for i from 1 to N do

q[j,i]:=ch[i]*q[j-1,i]+p[j-1,1]:

od: od:

for j from 1 to N do

for i from 1 to N do

s:=0:

for k from 1 to N do

s:=s+q[k,j]*A1[i,k]:

od:

xK[i,j]:=s:

od:od:

print(xK):

for j from 1 to N do

s:=0:

for i from 1 to N do

s:=s+(xK[i,j])^2:od:

nn:=1/sqrt(s):

for i from 1 to N do

xK[i,j]:=evalf(nn*xK[i,j]):

od: od:

print(xK):

B:=matrix(N,N):

y:=matrix(N,N):

x:=matrix(N,1):

xD:=matrix(N,N):

M:=proc(y,N,N0) local my,i;

my:=matrix(N0,N0,0):

for i from 1 to N0 do

my[i,i]:=1:

od:

for i from 1 to N0 do

my[N,i]:=-y[N+1,i]/y[N+1,N];

od:

my[N,N]:=1/y[N+1,N];

RETURN(my);

end:

evalm(M (A,N-1,N)):y:=A:

for i from 1 to N-1 do

B:=M(y,N-i,N);

y:=evalm(B^(-1)&*y&*B):

Mt(N-i):=B:

od:

sys:=lambda^N:

for i from 1 to N do

sys:=sys-y[1,i]*lambda^(N-i):

od:

print(sys);

ch:=solve(sys,lambda):

print(ch);

1.8024104, 6.1070423, 14.415537, 52.898303, 202.91830

for i from 1 to N do

for j from 1 to N do

x[j,1]:=evalf( ch[i]^(N-j) ):

od:

for i1 from 1 to N-1 do

x:=evalm(Mt(i1)&*x);

od:

for i1 from 1 to N do

xD[i1,i]:=x[i1,1];

od;

od;

evalm(xD);

for j from 1 to N do

s:=0:

for i from 1 to N do

s:=s+(xD[i,j])^2:od:

nn:=1/sqrt(s):

for i from 1 to N do

xD[i,j]:=evalf(nn*xD[i,j]):

od: od:

print(xD):

Размещено на Allbest.ru