DataMuseum.dk

top - metrics - download

    Length: 5376 (0x1500)
    Types: TextFile
    Names: »alggaussbar«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦7e928b248⟧ »algbib« 
            └─⟦this⟧ 

TextFile

;kemlab5 1
lookup gosdiag
if ok.no
(gosdiag=set 6
gosdiag= algol extgosdiag index.no)
gaussbar=set 47
gaussbar=algol index.no
\f




GAUSS-BARRIER af 13/10-72.

begin
comment
     Program til beregning af egenværdier og egenvektorer
for operatoren

   1/2 P**2 + A X**2 + B X**4 + C X**6 + D exp(-E X**2) 

som basis benyttes egenfunktioner til operatoren 

   1/2 P**2 +   1/2 X**2                         ;

zone res(128,1,stderror);
real A,B,C,D,E,V,ny,Q2,Q4,Q6,eps,a,b,c,d,sq1,sq2,sq3;
integer N,m1,m2,outmatrix,outvector,i,j,k,k1,k2,m,n;
array head(1:12);
boolean even,sp,nl,closeres;
system(2,n,head); even:= true; sp:= false add 32; nl:= false add 10;
readhead(in,head,1);
read(in,N,A,B,C,D,E,ny,m1,m2,outmatrix,outvector);
closeres:= outmedium(res);  a:= m2-m1+6;
n:= entier (-a + sqrt((a-6)*a-100+n//2));
if N>n then N:= n;
m:= if ny>0 then (m2+m2-1) else (m1+m1-2);
begin
comment oprettelse af H-matrix.;
array h(1:N*(N+1)//2),ev(m1:m2),e(m1+m1-2:m),x(m1:m2,1:N);
matrix:
if D=0 then begin
   for i:= N*(N+1)//2 step -1 until 1 do h(i):= 0
end else begin
h(1):= D/sqrt(1+E);  if -, even then h(1):= h(1)/(1+E);
i:= j:= 1;  k:= m:= if even then 0 else 1;
for m:= m+2 while j<N do begin
   n:= k; i:= i+1; a:= sqrt(m*(m-1));
   h(i):= -E*(m+n-1)*h(i-j)/((1+E)*a);
   j:= j+1;  a:= a*(1-E);
   for n:= n+2 while n<=m do begin
      i:= i+1;
      h(i):= (h(i-j)*a - E*(m+n-1)*h(i-1))/((1+E)*sqrt(n*(n-1)))
   end n
end m
end exponentdel;
sq1:= sq2:= sq3:= 0;  a:= A+0.5; b:= (A-0.5)/2; 
n:= if even then -2 else -1; k:= 0;
for i:=1 step 1 until N do begin 
   n:=n+2; k:= k+i;
   h(k):= h(k) + (((n+1)*n+1.5)*C*2.5+a)*(n+0.5)
               + B*(n*(n+1)+0.5)*1.5;
   if i>1 then
   h(k-1):= h(k-1) + (((n-1)*n+1)*C*1.875+(n-0.5)*B + b)*sq1;
   if i>2 then
   h(k-2):= h(k-2) + ((n-1.5)*C*0.75+B*0.25)*sq1*sq2;
   if i>3 then
   h(k-3):= h(k-3) + C*0.125*sq1*sq2*sq3;
   sq3:= sq2; sq2:= sq1;  sq1:= sqrt((n+1)*(n+2))
end i;
i:=1;
if even or outvector>0 or (m2-m1)>5 then
write(res,<:<12>:>,nl,2,string head(increase (i)),nl,1)
else write(res,nl,3);
if even then write(res,<:Even eigenfunctions and -values:>) else
write(res,<:Odd eigenfunctions and -values:>);
write(res,nl,1,<:Matrix dimension:>,<<  ddd>,N,nl,1,
          <:A,B,C,D,E:>,<< -d.ddddd'-d>,A,B,C,D,E);
if (D>0 or A<0) and even then begin
   a:= 0;  eps:= D*E;  b:= abs(A-eps)/(B+B+E*eps) + 1;
   for c:= ((3*C*a+B+B)*a+A-eps)/((3*C*a+B)*2+E*eps)
   while abs c < b and (c-a)*E < 1000 do begin
      a:= a-c; b:= abs c; eps:= exp(-E*a)*E*D
   end;
   if a<0 then a:= b:= 0 else b:= sqrt(a);
   V:= D - ((C*a+B)*a+A)*a - D*exp(-E*a);
   write(res,nl,1,<:V, Xmin  :>,<< -d.ddddd'-d>,V,b)
end else V:= 0;
if outmatrix>0 then begin
for k1:=1,k1+6 while k1<N do begin write(res,<:<10>:>);
k2:=k1-1;
for i:=k1 step 1 until N do begin write(res,<:<10>:>);
k:=i*(i-1)//2; if k2<k1+5 then k2:=k2+1;
for j:=k1 step 1 until k2 do write(res,<<  -b.dddd'-d>,h(k+j));
end; end; end;
if m1=0 then goto ud;
eps:= h((N+1)*N//2);  gosdiag(N,m1,m2,h,ev,x);
if outvector>0 then begin
k1:= m1-4;  n:= if even then -2 else -1;
rep:  write(res,nl,2);
k1:= k1+4; k2:= k1+3; if k2>m2 then k2:= m2;
for i:=k1 step 1 until k2 do
   write(res,sp,6,<:n = :>,<<dd>,n+i+i,sp,6);
write(res,nl,2);
for i:=k1 step 1 until k2 do
   write(res,<<  -d.dddddddddd'-d>,ev(i));
for j:=1 step 1 until N do  begin
   write(res,nl,1);
   for i:=k1 step 1 until k2 do
      write(res,<<    -d.dddddddd>,x(i,j),sp,3)
end;
if k2<m2 then goto rep
end else begin
for i:= m1 step 1 until m2 do begin
write(res,nl,2,<<-d.dddddddddd'-d>,ev(i),
          <<  -d.d'-dd>,eps*x(i,N)**2);
Q2:= 0.5; Q4:= 0.75; Q6:= 1.875;
a:= b:= c:= sq1:= sq2:= sq3:= 0;
n:= if even then -2 else -1; if ny>0 then e(n+i+i):= ev(i);
for j:=1 step 1 until N do begin
   d:= c; c:= b; b:= a; a:= x(i,j); n:= n+2;
   Q2:= Q2+n*a*a+sq1*a*b;
   Q4:= Q4+1.5*(n+1)*n*a*a+(n+n-1)*sq1*a*b+0.5*sq1*sq2*a*c;
   Q6:= Q6+((n+1.5)*n*2.5+5)*n*a*a+((n-1)*n+1)*3.75*sq1*a*b
          +(n-1.5)*1.5*sq1*sq2*a*c+0.25*sq1*sq2*sq3*a*d;
   sq3:= sq2; sq2:= sq1; sq1:= sqrt((n+1)*(n+2))
end;
write(res,sp,4,<<      d.dddddddd'-d>,Q2,nl,1,sp,30,
          Q4,nl,1,sp,30,Q6);
end end;
ud:  if even then begin even:= false ;goto matrix end;
if outvector=0 and ny>0 then begin
   N:= m2+m2-1; n:= m1+m1-2; a:= ny/(e(n+1)-e(n)); e(n):= e(n)*a;
   write(res,nl,3,<:V = :>,<<ddddddd.dd>,V*a,
             nl,2,<:n = :>,<<dd>,n,sp,4,<<ddddddd.dd>,e(n));
   for i:=n+1 step 1 until N do begin
      e(i):= e(i)*a;
      write(res,nl,1,sp,22,<<ddddd.dd>,e(i)-e(i-1),
                nl,1,<:n = :>,<<dd>,i,sp,4,<<ddddddd.dd>,e(i))
end end;
write(res,<:<25>:>); close(res,closeres);
end end;
permanent gaussbar.15
▶EOF◀