DataMuseum.dk

top - metrics - download

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

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

TextFile

rkfifth=algol
procedure for integrating a system of first
order differential equations
10 11 72 13 30 00
external
procedure rkfifth(f,x,y,b,eps,n,fi);
value b,n,fi;
procedure f; array y;  real x,b,eps; integer n; boolean fi;
begin
  integer i; real h,xl,int,fh,fhm,mu,d0; own real oldh;
  array d, yl,k0,k1,k2,k3,k4,k5(1:n);
  procedure stp(xincr,yincr,k);
  real xincr,yincr; array k;
  begin
    x:=xincr*h+xl;
    for i:=1 step 1 until n do y(i):=yl(i)+yincr*h;
    f(x,y,k)
  end step;
  xl:=x; int:=b-xl;
  if fi then oldh:=int;
  for i:=1 step 1 until n do
  begin
    yl(i):=y(i); d(i):=0
  end;
  d0:=0; h:=abs(oldh);
  if int<0 then
  begin
    int:=-int;
    h:=-h
  end;
  int:=eps*int;
  fi:=true;
  for mu:=abs(h) while fi do
  begin comment integrationloop;
    if mu<int then
    begin
      eps:=2*eps;
      int:=2*int;
      h:=(if h<0 then -int else int)*5
    end adjustmentafeps;
    if (h-(b-xl))*h>=0 then
    begin
      oldh:=h;
      h:=b-xl;
      fi:=false
    end setlaststep;
    f(xl,yl,k0);
    stp(0.2222222222,k0(i)*2/9,k1);
    stp(0.3333333333,(k0(i)+k1(i)*3)/12,k2);
    stp(0.5,(k0(i)+k2(i)*3)/8,k3);
    stp(0.8,(53*k0(i)-135*k1(i)+126*k2(i)+56*k3(i))/125,k4);
    stp(1.0,(133*k0(i)-378*k1(i)+276*k2(i)+112*k3(i)+25*k4(i))/168,k5);
    fhm:=0;
    for i:=1 step 1 until n do
    begin
      fh:=abs(21*k0(i)-162*k2(i)+224*k3(i)-125*k4(i)
        +42*k5(i))/14/(abs(k0(i))+1)/eps;
      if fh>fhm then fhm:=fh
    end;
    mu:=1/(1+fhm)+0.5;
    if fhm<2 then
    begin comment stepaccepted;
      stp(1.0,(-63*k0(i)+189*k1(i)-36*k2(i)-112*k3(i)
        +50*k4(i))/28,k5);
      for i:=1 step 1 until n do
      begin comment finalincrement;
        fh:=(35*k0(i)+162*k2(i)+125*k4(i)+14*k5(i))/336*h+d(i);
        fhm:=y(i):=yl(i)+fh;
        d(i):=fh-(fhm-yl(i))
      end;
      fh:=h+d0; x:=xl+fh; d0:=fh-(x-xl); xl:=x;
      for i:=1 step 1 until n do yl(i):=y(i)
    end stepaccepted
    else fi:=true;
    h:=mu*h
  end integrationloop
end rkfifth;
end;
▶EOF◀