home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / ODEINT.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  1.4 KB  |  63 lines
  1. PROGRAM d15r4(input,output);
  2. (* driver for ODEINT *)
  3. CONST
  4.    n=4;
  5. TYPE
  6.    glnarray = ARRAY [1..n] OF real;
  7.    glarray = glnarray;
  8. VAR
  9.    dxsav,eps,h1,hmin,x1,x2 : real;
  10.    i,kmax,kount,nbad,nok : integer;
  11.    ystart : glnarray;
  12.    xp : ARRAY [1..200] OF real;
  13.    yp : ARRAY [1..10,1..200] OF real;
  14.  
  15. (*$I MODFILE.PAS *)
  16. (*$I BESSJ0.PAS *)
  17.  
  18. (*$I BESSJ1.PAS *)
  19.  
  20. (*$I BESSJ.PAS *)
  21.  
  22. PROCEDURE derivs(x: real; y: glnarray; VAR dydx: glnarray);
  23. (* Programs using routine DERIVS must define the type
  24. TYPE
  25.    glnarray = ARRAY [1..4] OF real;
  26. in the calling routine. *)
  27. BEGIN
  28.    dydx[1] := -y[2];
  29.    dydx[2] := y[1]-(1.0/x)*y[2];
  30.    dydx[3] := y[2]-(2.0/x)*y[3];
  31.    dydx[4] := y[3]-(3.0/x)*y[4]
  32. END;
  33.  
  34. (*$I RK4.PAS *)
  35.  
  36. (*$I RKQC.PAS *)
  37.  
  38. (*$I ODEINT.PAS *)
  39.  
  40. BEGIN
  41.    x1 := 1.0;
  42.    x2 := 10.0;
  43.    ystart[1] := bessj0(x1);
  44.    ystart[2] := bessj1(x1);
  45.    ystart[3] := bessj(2,x1);
  46.    ystart[4] := bessj(3,x1);
  47.    eps := 1.0e-4;
  48.    h1 := 0.1;
  49.    hmin := 0.0;
  50.    kmax := 100;
  51.    dxsav := (x2-x1)/20.0;
  52.    odeint(ystart,n,x1,x2,eps,h1,hmin,nok,nbad);
  53.    writeln;
  54.    writeln('successful steps:',' ':13,nok:3);
  55.    writeln('bad steps:',' ':20,nbad:3);
  56.    writeln('stored intermediate values:','   ',kount:3);
  57.    writeln;
  58.    writeln('x':8,'integral':18,'bessj(3,x)':15);
  59.    FOR i := 1 to kount DO BEGIN
  60.       writeln(xp[i]:10:4,yp[4,i]:16:6,bessj(3,xp[i]):14:6)
  61.    END
  62. END.
  63.