DataMuseum.dk

top - metrics - download

    Length: 3840 (0xf00)
    Types: TextFile
    Names: »pminttest«

└─⟦2c55ea56f⟧ Bits:30001844 SW-save af projekt 1000, Alarm-system
    └─⟦6b41451d2⟧ 
        └─⟦this⟧ »pminttest« 

TextFile

p=algol list.yes bossline.yes
\f


begin
procedure intgr(t,c,f,n,h,t1,t2,eps,print); value t,n,h,t1,t2,eps;
real t,h,t1,t2,eps; array c; procedure f,print; integer n;
begin
  array d,d1,d2,c1,c2,b,b1,b2,a,a1,a2,p,p1,p2(1:n);
  real norm; integer i; boolean dbl;

  procedure l1(a,a1,a2,b,b1,b2,t); array a,a1,a2,b,b1,b2; real t;
  begin integer i; real x;
    for i:=1 step 1 until n do a(i):=b1(i)+b2(i)*h+b(i);
    f(a2,a,t+h);  norm:=0;
    for i:=1 step 1 until n do begin
      a(i):=x:=b1(i)+(b2(i)+a2(i))*h/2+b(i);
      a1(i):=b(i)-x+(b2(i)+a2(i))*h/2+b1(i);
      if abs(a2(i)-b2(i))/2 > norm then norm:=abs(a2(i)-b2(i))/2;
    end;
  end l1;

  procedure l2(a,a1,a2,b,b1,b2,c,c1,c2,t);
    array a,a1,a2,b,b1,b2,c,c1,c2; real t;
  begin integer i; real x;
    for i:=1 step 1 until n do a(i):=c1(i)+(9/4*b2(i)+3/4*c2(i))*h+c(i);
    f(a2,a,t+h); norm:=0;
    for i:=1 step 1 until n do begin
      x:=27/32*(b(i)-c(i)+b1(i)-c1(i))
         -(45/32*b2(i)+21/32*c2(i)-3/8*a2(i))*h;
      a(i):=a(i)+x;
      a1(i):=c(i)-a(i)+(9/4*b2(i)+3/4*c2(i))*h+x+c1(i);
      if abs x > norm then norm:=abs x;
    end;
    norm:=norm/h;
  end l2;

  if t >= t1 then print(t,c,t2,t2);
  if t+3*h > t2 then h:=(t2-t)*1.5 else h:=h*2;
  for i:=1 step 1 until n do c1(i):=0;
  f(c2,c,t);

  repeat
    repeat
      h:=h/2; l1(b,b1,b2,c,c1,c2,t);
    until norm <= eps;
    l1(a,a1,a2,b,b1,b2,t+h);
  until norm <= eps;
  dbl:=false;  t:=t+2*h;
  repeat
    l2(p,p1,p2,a,a1,a2,c,c1,c2,t);
    if norm > eps then begin
      for i:=1 step 1 until n do begin
        d(i):=c(i); d1(i):=c1(i); d2(i):=c2(i);
        c(i):=b(i); c1(i):=b1(i); c2(i):=b2(i);
      end;
      h:=h/2; l2(b,b1,b2,c,c1,c2,d,d1,d2,t-2*h); dbl:=false;
    end else begin
      t:=t+h;  if t >= t1 then print(t,p,t1,t2);
      if t < t2 then begin
        while t+h > t2 do begin
          for i:=1 step 1 until n do begin
            d(i):=b(i); d1(i):=b1(i); d2(i):=b2(i);
            b(i):=a(i); b1(i):=a1(i); b2(i):=a2(i);
            a(i):=c(i); a1(i):=c1(i); a2(i):=c2(i);
          end;
          h:=h/2;
          l2(c,c1,c2,d,d1,d2,a,a1,a2,t-4*h);
          l2(a,a1,a2,b,b1,b2,d,d1,d2,t-2*h);
          norm:=eps;
        end t+h > t2;
        if norm < eps/16 and t+2*h <= t2 and dbl then begin
          for i:=1 step 1 until n do begin
            a(i):=p(i); a1(i):=p1(i); a2(i):=p2(i);
            c(i):=d(i); c1(i):=d1(i); c2(i):=d2(i);
          end;
          h:=2*h;  dbl:=false;
        end else begin
          for i:=1 step 1 until n do begin
            d(i):=c(i); d1(i):=c1(i); d2(i):=c2(i);
            c(i):=b(i); c1(i):=b1(i); c2(i):=b2(i);
            b(i):=a(i); b1(i):=a1(i); b2(i):=a2(i);
            a(i):=p(i); a1(i):=p1(i); a2(i):=p2(i);
          end;
          dbl:=true;
        end eps/16 < norm < eps;
      end t<t2;
    end norm < eps;
  until t>= t2;
end intgr;

procedure f(x1,x,t); array x1,x; real t;
begin x1(1):=-x(2); x1(2):=x(1); end;

procedure print(t,x,t1,t2); real t,t1,t2; array x;
write(out,<:<10>:>,<<  -d.ddddd ddddd>,t,x(1),x(2));

array x(1:2);
x(1):=1; x(2):=0;
intgr(0,x,f,2,3,0,1.5707963,'-6,print);
end;

p
finis
▶EOF◀