home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / POLYROOT.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-06-26  |  4.7 KB  |  188 lines
  1. program polyroot(input,output);
  2. (* Program uses the Newton-Raphson technique to find a simple real root *)
  3. (* of a polynomial and removes it by synthetic division.  The procedure *)
  4. (* may be continued until all such roots are found or until the reduced *)
  5. (* polynomial is of degree 2.  If the polynomial can be reduced to one *)
  6. (* of degree 2, the quadratic formula is applied to yield all roots. *)
  7.  
  8. label 10,20;
  9. type vector = array[0..10] of real;
  10. var coef: vector;
  11.     x0,error: real;
  12.     degree,delta,i,maxiter,counter: integer;
  13.     q1,q2,q3,q4: char;
  14.  
  15. procedure readin(var a1: vector; n1: integer);
  16. (* Reads in polynomial coefficients from the keyboard. *)
  17.  
  18. var i: integer;
  19.  
  20. begin
  21.    for i := 0 to n1 do
  22.       begin
  23.          write('Input the A(',i,') coefficent: ');
  24.          readln(a1[i]);
  25.       end;
  26.    writeln('The polynomial coefficients are: ');
  27.  
  28.    writeln;
  29.    for i := 0 to n1 do
  30.       begin
  31.          write(a1[i]:9,'  ');
  32.       end;
  33.    writeln;
  34. end;
  35.  
  36.  
  37.  
  38. procedure newtraf(coef: vector; error1: real; var root: real; var k1:integer;
  39.                   n1,max: integer);
  40. (* Searches for a real root of a polynomial by the Newton-Raphson technique. *)
  41.  
  42. var yd,fcn,deriv,a1: real;
  43.     i: integer;
  44.  
  45. begin
  46.    a1 := 1.0;
  47.    fcn := a1*coef[0];
  48.    for i := 1 to n1 do
  49.       begin
  50.          a1 := root*a1;  (* finds root^i *)
  51.          fcn := fcn + a1*coef[i];   (* value of trunc polynomial *)
  52.       end;
  53.    a1 := 1.0;
  54.    deriv := a1*coef[1];
  55.    for i := 1 to (n1-1) do
  56.       begin
  57.          a1 := root*a1;
  58.          deriv := deriv + (i+1)*a1*coef[i+1]; (* value of derivative *)
  59.       end;
  60.    yd := deriv;
  61.    if deriv = 0.0 then
  62.      begin
  63.         yd := 0.0000001;  (* prevent divide by zero *)
  64.         writeln('Try a new guess!');
  65.      end;
  66.    root := root - fcn/yd;
  67.    k1 := k1 + 1;
  68.    writeln('Value of fcn: ',fcn:9);
  69.    if (k1 < max) and (abs(fcn/yd) > error1) then
  70.       newtraf(coef,error1,root,k1,n1,max);
  71. end;
  72.  
  73.  
  74.  
  75. procedure syndiv(var c:vector; root: real; var n1: integer);
  76. (* Carries out synthetic division of a real polynomial by a first-degree *)
  77. (* real polynomial. *)
  78.  
  79. label 10;
  80. type vect2 = array[0..1] of real;
  81. var a: vector;
  82.     b: vect2;
  83. var i,j: integer;
  84.  
  85. begin
  86.    b[1] := 1.0;
  87.    b[0] := -root;
  88.    for i := n1 downto 1 do
  89.       begin
  90.          a[i-1] := c[i]/b[1];
  91.          if i = 1 then
  92.             goto 10;
  93.          for j := 0 to 1 do
  94.             begin
  95.                c[i-j] := c[i-j] - a[i-1]*b[1-j];
  96.             end;
  97.          10:writeln;
  98.       end;
  99.    n1 := n1 - 1;
  100.    for i := 0 to n1 do
  101.       begin
  102.          c[i] := a[i];
  103.       end;
  104. end;
  105.  
  106. procedure quadroot(c: vector); 
  107. (* Applies quadratic formula to a real second-degree polynomial, *)
  108. (* in the form  c[2]*x*x + c[1]*x + c[0] = 0 *)
  109. var y,y1,x1,x2,z1,z2: real;
  110.   begin
  111.     y := c[1]*c[1]/4-c[2]*c[0];
  112.     y1 := sqrt(abs(y));
  113.     z1 := -c[1]/(2*c[2]);
  114.     z2 := y1/c[2];
  115.     if y > 0 then
  116.       begin
  117.         x1 := z1 + z2;
  118.         x2 := z1 - z2;
  119.         writeln('The roots are: ',x1:10,' and ',x2:10);
  120.       end
  121.     else
  122.       begin
  123.         write('The roots are: ',z1:9,' + j',z2:9);
  124.         writeln(' and  ',z1:9,' - j',z2:9);
  125.      end;
  126.   end;
  127.  
  128. begin  (* main *)
  129.    write('Input degree of polynomial: ');
  130.    readln(degree);
  131.    writeln;
  132.    readin(coef,degree);
  133.    write('Input convergence factor: ');
  134.    readln(error);
  135.    writeln;
  136.    write('Input maximum number of iterations: ');
  137.    readln(maxiter);
  138.    writeln;
  139. 10:counter := 0;
  140.    write('Input the initial guess: ');
  141.    readln(x0);
  142.    writeln;
  143. 20:newtraf(coef,error,x0,counter,degree,maxiter);
  144.    writeln;
  145.    writeln('The calculated root is: ',x0:10);
  146.    writeln('The number of iterations is: ',counter);
  147.    writeln;
  148.    writeln('Try another initial guess (y,n)? ');
  149.    readln(q1);
  150.    if q1 = 'y' then
  151.       goto 10;
  152.    write('More iterations (y,n)? ');
  153.    readln(q2);
  154.    if q2 = 'y' then
  155.       begin
  156.          writeln('How many more iterations: ');
  157.          readln(delta);
  158.          maxiter := maxiter + delta;
  159.          writeln;
  160.          write('New convergence factor (y,n)?');
  161.          readln(q3);
  162.          if q3 = 'n' then
  163.             goto 20;
  164.          write('Input new convergence factor: ');
  165.          readln(error);
  166.          goto 20;
  167.       end;
  168.    syndiv(coef,x0,degree);
  169.    writeln('The coefficients of the resulting polynomial are: ');
  170.    writeln;
  171.    for i := 0 to degree do
  172.       begin
  173.          write(coef[i]:8:3,'  ');
  174.       end;
  175.    writeln;
  176.    if degree > 2 then
  177.       begin
  178.          write('Continue to search for roots (y,n)? ');
  179.          readln(q4);
  180.          if q4 = 'y' then
  181.             goto 10;
  182.       end;
  183.    if degree = 2 then
  184.       writeln;
  185.       quadroot(coef);
  186. end.
  187.  
  188.