home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / TRAPZ2.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1985-06-26  |  3.6 KB  |  127 lines
  1. program trapz2(input,output);
  2. (* Program solves second-order state model  x' = Ax + f using the method *)
  3. (* of trapezoidal integration. *)
  4.  
  5. type mat = array[1..2,1..2] of real;
  6. var x10,x20,t0,tr,h: real;
  7.     np: integer;
  8.     a: mat;
  9.  
  10. function fn1(x: real): real;  (* Forcing function # 1 *)
  11.    begin
  12.       fn1 := 20.0;
  13.    end;
  14.  
  15. function fn2(x: real): real;  (* Forcing function # 2 *)
  16.    begin
  17.       fn2 := 20*exp(-5*x)-100.0;
  18.    end;
  19.  
  20. procedure integ2(b,c: mat; h1,t1,tr1,y10,y20: real; np1: integer); forward;
  21. (* Uses trapezoidal integration to solve a second-order state model with *)
  22. (* initial conditions y10 and y20 at t1.  Integration increment: h1 *)
  23. (* Run time: tr1.  Number of points between print: np1. *)
  24.  
  25. var h2,y1,y2,y11,y12,y13,y21,y22,y23,t: real;
  26.     k,i: integer;
  27.  
  28. procedure constants(a: mat; h1: real; np1: integer);
  29. (* Evaluates two second-order matrices, and calculates inverse of one. *)
  30.  
  31. var b,c,d: mat;
  32.     h2,dtc,t1 : real;
  33. begin
  34.    h2 := h1/2;
  35.    c[1,1] := 1 - a[1,1]*h2;   (*  find matrix C = (U - A*h/2)  *)
  36.    c[1,2] := -a[1,2]*h2;
  37.    c[2,1] := -a[2,1]*h2;
  38.    c[2,2] := 1 - a[2,2]*h2;
  39.    d[1,1] := 1 + a[1,1]*h2;   (*  find matrix D = (U + A*h/2)  *)
  40.    d[1,2] := a[1,2]*h2;
  41.    d[2,1] := a[2,1]*h2;
  42.    d[2,2] := 1 + a[2,2]*h2;
  43.    dtc := c[1,1]*c[2,2] - c[1,2]*c[2,1];   (*  find det C  *)
  44.    if dtc = 0 then writeln('Divide by zero.  Change increment.');
  45.    b[1,1] := c[2,2]/dtc;      (*  find the inverse of matrix C  *)
  46.    b[1,2] := -c[1,2]/dtc;
  47.    b[2,1] := -c[2,1]/dtc;
  48.    b[2,2] := c[1,1]/dtc;
  49.    c[1,1] := b[1,1]*d[1,1] + b[1,2]*d[2,1]; (* Find C-inverse x D *)
  50.    c[1,2] := b[1,1]*d[1,2] + b[1,2]*d[2,2];
  51.    c[2,1] := b[2,1]*d[1,1] + b[2,2]*d[2,1];
  52.    c[2,2] := b[2,1]*d[1,2] + b[2,2]*d[2,2];
  53.    integ2(b,c,h1,0,tr,x10,x20,np1);
  54. end;
  55.  
  56. procedure change(var h1: real; var np1: integer); forward;
  57. (* Halves the increment of integration, if desired. *)
  58.  
  59. var b: char;
  60.  
  61. procedure integ2;  (* Comment above *)
  62.  
  63. begin
  64.    y1 := y10;
  65.    y2 := y20;
  66.    t := t1;
  67.    h2 := h1/2;
  68.    writeln(t:6:3,'':5,y1:9,'  ',y2:9);
  69.    while t <= t1 + tr1 + h1*np1 do
  70.       begin
  71.          for i := 1 to np1 do
  72.             begin   (*  find Y(k+1) = C*Y(k) + B*(F(k) + F(k+1))*h/2  *)
  73.                y11 := c[1,1]*y1 + c[1,2]*y2;
  74.                y12 := b[1,1]*(fn1(t) + fn1(t+h1))*h2;
  75.                y13 := b[1,2]*(fn2(t) + fn2(t+h1))*h2;
  76.                y21 := c[2,1]*y1 + c[2,2]*y2;
  77.                y22 := b[2,1]*(fn1(t) + fn1(t+h1))*h2;
  78.                y23 := b[2,2]*(fn2(t) + fn2(t+h1))*h2;
  79.                y1 := y11 + y12 + y13;
  80.                y2 := y21 + y22 + y23;
  81.                t := t + h1;
  82.             end;
  83.          writeln(t:6:3,'':5,y1:9,'  ',y2:9);
  84.       end;
  85.    change(h1,np1);
  86. end;
  87.  
  88.  
  89. procedure change;  (* Comment above *)
  90.  
  91. begin
  92.    write('Try a smaller increment? ');
  93.    readln(b);
  94.    writeln;
  95.    if b = 'y' then
  96.       begin
  97.          h1 := h1/2;
  98.          np1 := 2*np1;
  99.          constants(a,h1,np1);
  100.       end;
  101. end;
  102.  
  103.  
  104.  
  105. begin          (*  main  *)
  106.    writeln('Forcing functions must be in function subprograms.');
  107.    writeln;
  108.    write('Input coefficients a(1,1) and a(1,2) (with space between): ');
  109.    readln(a[1,1],a[1,2]);
  110.    writeln;
  111.    write('Input coefficients a(2,1) and a(2,2): ');
  112.    readln(a[2,1],a[2,2]);
  113.    writeln;
  114.    write('Input initial values x1(0) and x2(0): ');
  115.    readln(x10,x20);
  116.    writeln;
  117.    write('Input run time: ');
  118.    readln(tr);
  119.    writeln;
  120.    write('Input time increment and no. of increments between print: ');
  121.    readln(h,np);
  122.    writeln;
  123.    writeln('   Time':10,'   x1':10,'   x2':10);
  124.    writeln;
  125.    constants(a,h,np);
  126. end.
  127.