DataMuseum.dk

top - metrics - download

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

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

TextFile

chstiffpp=set 5
scope day chstiffpp
chstiffpp=algol list.no
external integer procedure chstiffpp(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,epsr,epsa;
     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 chstiffpp:= -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;

     chstiffpp:= 0; test:= false;
     xin(X,tb,dt,te,hb,epsr,epsa); 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);
     if test then xout(XH,t,-h,1);
     chstep(XH,h,XH2);
     if test then xout(XH2,t,-h,2);
     chstep(X,2*h,X2H);
     if test then xout(X2H,t,-h,3);
     r:= 0;
     for i:= 1 step 1 until n do
     begin
        if X(i)<>0 and X(i)>epsa then q:= abs(X2H(i)-XH2(i))/X(i)
        else if XH(i)<>0 and abs XH(i)>epsa then
           q:= abs((X2H(i)-XH2(i))/XH(i))
        else q:= 0;
        if q>r then r:= q;
     end;
     if r<>0 then
     begin q:= h*(14*epsr/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):= if XH(i)>0 then XH(i) else 0;
     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◀