home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / RK4.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  1.2 KB  |  57 lines
  1. PROGRAM d15r1(input,output);
  2. (* driver for routine RK4 *)
  3. CONST
  4.    n=4;
  5. TYPE
  6.    glnarray = ARRAY [1..n] OF real;
  7. VAR
  8.    h,x : real;
  9.    i,j : integer;
  10.    y,dydx,yout : glnarray;
  11.  
  12. (*$I MODFILE.PAS *)
  13. (*$I BESSJ0.PAS *)
  14.  
  15. (*$I BESSJ1.PAS *)
  16.  
  17. (*$I BESSJ.PAS *)
  18.  
  19. PROCEDURE derivs(x : real; y : glnarray;
  20.        VAR dydx : glnarray);
  21. (* The arrays y and dydx must carry the dimension given to
  22. them in the calling routine. *)
  23. BEGIN
  24.    dydx[1] := -y[2];
  25.    dydx[2] := y[1]-(1.0/x)*y[2];
  26.    dydx[3] := y[2]-(2.0/x)*y[3];
  27.    dydx[4] := y[3]-(3.0/x)*y[4]
  28. END;
  29.  
  30. (*$I RK4.PAS *)
  31.  
  32. BEGIN
  33.    x := 1.0;
  34.    y[1] := bessj0(x);
  35.    y[2] := bessj1(x);
  36.    y[3] := bessj(2,x);
  37.    y[4] := bessj(3,x);
  38.    dydx[1] := -y[2];
  39.    dydx[2] := y[1]-y[2];
  40.    dydx[3] := y[2]-2.0*y[3];
  41.    dydx[4] := y[3]-3.0*y[4];
  42.    writeln;
  43.    writeln('Bessel function:':16,'j0':5,
  44.          'j1':12,'j3':12,'j4':12);
  45.    FOR i := 1 to 5 DO BEGIN
  46.       h := 0.2*i;
  47.       rk4(y,dydx,n,x,h,yout);
  48.       writeln;
  49.       writeln('for a step size of:',h:6:2);
  50.       write('rk4: ':11);
  51.       FOR j := 1 to 4 DO write(yout[j]:12:6);
  52.       writeln;
  53.       writeln('actual: ':11,bessj0(x+h):12:6,
  54.          bessj1(x+h):12:6,bessj(2,x+h):12:6,bessj(3,x+h):12:6)
  55.    END
  56. END.
  57.