DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »tchstiffp«

└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
    └─⟦b2ec5d50f⟧ 
        └─⟦09b4e9619⟧ »thcømat« 
            └─⟦this⟧ 

TextFile

decompose=algol
external
boolean procedure decompose(A,P,i);
array A; integer array P; integer i;
; end

solve=algol
external
procedure solve(A,P,i,B);
array A,B; integer array P; integer i;
; end

chstiffp=algol
external integer procedure chstiffp(n,fjac,xin,xout);
value n;
integer n;
procedure fjac,xin,xout;
begin
     integer i,j,sg;
     real r,q,h,t,tb,hb,dt,te,eps;
     boolean test;
     integer array P,OSG(1:n);
     real array X,XH,X2H,XH2,K1,K2(1:n),A(1:n,1:n);

procedure chstep(X,h,Y);
value h;
real h;
real array X,Y;
begin fjac(true,X,t,K1,A);
     for i:= 1 step 1 until n do
     for j:= 1 step 1 until n do
        A(i,j):= -0.788675134*A(i,j);
     for i:= 1 step 1 until n do A(i,i):= 1/h+A(i,i);
     if -,decompose(A,P,0) then
     begin chstiffp:= -n-1; goto chexit end;
     solve(A,P,0,K1);
     for i:= 1 step 1 until n do 
        Y(i):= X(i)-1.15470054*K1(i);
     fjac(false,Y,t,K2,A);
     solve(A,P,0,K2);
     for i:= 1 step 1 until n do
        Y(i):= X(i)+0.75*K1(i)+0.25*K2(i);
end chstep;
     chstiffp:= 0; test:= false;
     xin(X,tb,dt,te,hb,eps); t:= tb; h:= hb;
     if te<=t then goto chexit;
     if h>0.5*dt then h:= 0.5*dt; 
     for i:=1 step 1 until n do OSG(i):= 0;
rep: chstep(X,h,XH);
     chstep(XH,h,XH2);
     chstep(X,2*h,X2H);
     r:= 0;
     for i:= 1 step 1 until n do
     if X(i)<>0 then
     begin q:= abs(X2H(i)-XH2(i))/X(i);
        if q>r then r:= q;
     end;
     if r<>0 then
     begin q:= h*(14*eps/r)**0.25;
        if test or q<'-10*hb then xout(X,t,-h,0);
        if q<0.8*h then
        begin if q<0.01*h then
           h:= 0.01*h else  h:= q; goto rep end;
     end else q:= 100*h;
     if t+h>tb+dt then
     begin h:= tb+dt-t; chstep(X,h,XH); end;
     for i:= 1 step 1 until n do
     begin sg:= sign(XH(i)-X(i));
        if sg<>OSG(i) then
        begin xout(X,t,h,-sg*i); OSG(i):= sg end;
     end;
     t:= t+h;
     for i:= 1 step 1 until n do X(i):= XH(i);
     if t=tb+dt then
     begin xout(X,t,h,0); tb:= tb+dt; end;
     if q>hb then q:= hb;
     if t+q<tb+dt then h:= q else h:= tb+dt-t;
     if t<te then goto rep;
chexit:
end;
end;
▶EOF◀