home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / TRAPZ1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-06-26  |  2.0 KB  |  84 lines
  1. program trapz1(input,output);
  2. (* Program solves first-order differential equation  x' = ax + f *)
  3. (* by trapezoidal integration. *)
  4.  
  5. var a,x0,t0,h,tr: real;
  6.     np: integer;
  7.  
  8. function fcn(z: real): real;  (* Define forcing function *)
  9.    begin
  10.       fcn := 10;
  11.    end;
  12.  
  13. procedure change(var h1: real; var np1: integer); forward;
  14. (* Halves the integration increment, and calls integration procedure. *)
  15.  
  16. var b: char;
  17.     a1,t1,tr1,x1: real;
  18.  
  19. procedure integ(x1,a1,h1,t1,tr1: real; np1: integer);
  20. (* Uses trapezoidal integration with increment h1 to solve x' = ax + f. *)
  21. (* Initial value is x1 at t1, run time is tr1, and np1 is the number of *)
  22. (* points between print. *)
  23.  
  24. var i: integer;
  25. x,k1,k2,x2,t: real;
  26. begin
  27.    t := t1;
  28.    x2 := x1;
  29.    writeln(t0:5,'':5,x0:10);
  30.    while t <= t1 + tr1 + h1*np1 do
  31.       begin
  32.          k1 := 1 + a1*h1/2;
  33.          k2 := 1 - a1*h1/2;
  34.          for i := 1 to np1 do
  35.             begin
  36.                x := (k1*x2 +(fcn(t) + fcn(t+h1))*h1/2)/k2;
  37.                t := t + h1;
  38.                x2 := x;
  39.             end;
  40.          writeln(t:5,'':5,x:10);
  41.       end;
  42.    change(h1,np1);
  43. end;
  44.  
  45. procedure change;  (* Comment above *)
  46.  
  47. begin
  48.    writeln('Try a smaller increment (y,n)?');
  49.    readln(b);
  50.    writeln;
  51.    if b = 'y' then
  52.       begin
  53.          h1 := h1/2;
  54.          np1 := 2*np1;
  55.          integ(x0,a,h1,t1,tr,np1);
  56.       end;
  57. end;
  58.  
  59. begin (* main *)
  60.    writeln('Forcing function must be in function subroutine.');
  61.    writeln;
  62.    write('Input the coefficient a (with sign): ');
  63.    readln(a);
  64.    writeln;
  65.    write('Input the initial value of x: ');
  66.    readln(x0);
  67.    writeln;
  68.    write('Input initial time: ');
  69.    readln(t0);
  70.    write('Input run time: ');
  71.    readln(tr);
  72.    writeln;
  73.    write('Input time increment: ');
  74.    readln(h);
  75.    writeln;
  76.    write('Input number of increments between print: ');
  77.    readln(np);
  78.    writeln;
  79.    writeln('Time   ':10,'Solution  ':10);
  80.    writeln;
  81.  
  82.    integ(x0,a,h,t0,tr,np);
  83. end.
  84.