alg028.pas

Numerical Anaysis 8th Edition, by Burden and Faires (Pascal Source)

      for I := 1 to N do
         begin
            writeln(OUP,'Coefficient of x^',I-1,' is ',A[I]:12:8);
         end;
      writeln(OUP,'The output is i, approximation x(i), f(x(i))');
      writeln(OUP);
      writeln(OUP,'Three columns means the results are real numbers,');
      writeln(OUP,'Five columns means the results are complex ');
      writeln(OUP,'numbers with real and imaginary parts of x(i)');
      writeln(OUP,'followed by real and imaginary parts f(x(i)).');
      writeln(OUP);
   end;
begin
   INPUT;
   if OK then
     begin
        OUTPUT;
{  Evaluate F using Horner's Method and save in the vector F           }
        for I := 1 to 3 do
           begin
              F[I] := A[N];
              for J := 2 to N do
                 begin
                    K := N-J+1;
                    F[I] := F[I]*X[I]+A[K]
                 end
           end;
{  Variable ISW is used to note a switch to complex arithmetic
    ISW=0 means real arithmetic, and ISW=1 means complex arithmetic    }
        ISW := 0;
{       STEP 1                                                         }
        H[1] := X[2]-X[1];
        H[2] := X[3]-X[2];
        DEL1[1] := (F[2]-F[1])/H[1];
        DEL1[2] := (F[3]-F[2])/H[2];
        DEL := (DEL1[2]-DEL1[1])/(H[2]+H[1]);
        I := 3;
        OK := true;
{       STEP 2                                                         }
        while ((I <= M) and (OK)) do
           begin
{        STEPS 3-7 for real arithmetic                                 }
              if (ISW = 0) then
                 begin
{              STEP 3                                                  }
                    B := DEL1[2]+H[2]*DEL;
                    D := B*B-4.0*F[3]*DEL;
{              test to see if need complex arithmetic                  }
                    if (D >= 0.0) then
                       begin
{                   real arithmetic - test to see if straight line     }
                          if (abs(DEL) <= ZERO) then
                             begin
{                          straight line - test if horizontal line     }
                                   if (abs(DEL1[2]) <= ZERO) then
                                      begin
                                         writeln('Horizontal Line');
                                         OK := false
                                      end
                                   else
                                      begin
{                                straight line but not horizontal      }
                                   writeln('line - not horizontal');
                                         X[4] := (F[3]-DEL1[2]*X[3])
                                              /DEL1[2];
                                         H[3] := X[4]-X[3]
                                      end
                             end
                          else
                             begin
{                          not straight line                           }
                                D := sqrt(D);
{                          STEP 4                                      }
                                E := B+D;
                                if (abs(B-D) > abs(E)) then E := B-D;
{                          STEP 5                                      }
                                H[3] := -2.0*F[3]/E;
                                X[4] := X[3]+H[3]
                             end;
                          if (OK) then
                             begin
{                          evaluate f(x(I))=F[4] by Horner's method    }
                                F[4] := A[N];
                                for J := 2 to N do
                                   begin
                                      K := N-J+1;
                                      F[4] := F[4]*X[4]+A[K]
                                   end;
                                writeln(OUP,I,X[4],F[4]);
{                          STEP 6                                      }
                                if (abs(H[3]) < TOL) then
                                   begin
                                      writeln(OUP);
                                      writeln(OUP,'Method Succeeds');
                                      write(OUP,'Approximation is within ');
                                      writeln(OUP,TOL);
                                      writeln(OUP,'in ',I,' iterations');
                                      OK := false
                                   end
                                else
{                             STEP 7                                   }
                                   begin
                                      for J := 1 to 2 do
                                         begin
                                            H[J] := H[J+1];
                                            X[J] := X[J+1];
                                            F[J] := F[J+1]
                                         end;
                                      X[3] := X[4];
                                      F[3] := F[4];
                                      DEL1[1] := DEL1[2];
                                      DEL1[2] := (F[3]-F[2])/H[2];
                                      DEL := (DEL1[2]-DEL1[1])
                                           /(H[2]+H[1])
                                   end
                             end
                       end
                    else
                       begin
{                    switch to complex arithmetic                      }
                          ISW := 1;
                          for J := 1 to 3 do
                             begin
                                ZR[J] := X[J]; ZI[J] := 0.0;
                                GR[J] := F[J]; GI[J] := 0.0
                             end;
                          for J := 1 to 2 do
                             begin
                                CDEL1R[J] := DEL1[J]; CDEL1I[J] := 0.0;
                                CH1R[J] := H[J]; CH1I[J] := 0.0
                             end;
                          CDELR := DEL; CDELI := 0.0
                       end
                 end;
              if ((ISW = 1) and (OK)) then
{               STEPS 3-7 complex arithmetic                           }
                 begin
{              test if straight line                                   }
                    if (CABS(CDELR,CDELI) <= ZERO) then
                       begin
{                   straight line - test if horizontal line            }
                          if (CABS(CDEL1R[1],CDEL1I[1]) <= ZERO) then
                             begin
                                writeln('horizontal line - complex');
                                OK := false
                             end
                          else
{                       straight line but not horizontal               }
                             begin
                                writeln('line - not horizontal');
                                CMULT(CDEL1R[2],CDEL1I[2],
                                      ZR[3],ZI[3],ER,EI);
                                CSUB(GR[3],GI[3],ER,EI,FR,FI);
                                CDIV(FR,FI,CDEL1R[2],
                                     CDEL1I[2],ZR[4],ZI[4]);
                                CSUB(ZR[4],ZI[4],ZR[3],ZI[3],
                                     CH1R[3],CH1I[3])
                             end
                       end
                    else
{                not straight line                                     }
                       begin
{                   STEP 3                                             }
                          CMULT(CH1R[2],CH1I[2],CDELR,CDELI,ER,EI);
                          CADD(CDEL1R[2],CDEL1I[2],ER,EI,CBR,CBI);
                          CMULT(GR[3],GI[3],CDELR,CDELI,ER,EI);
                          CMULT(ER,EI,4.0,0.0,FR,FI);
                          QR := CBR; QI := CBI;
                          CMULT(CBR,CBI,QR,QI,ER,EI);
                          CSUB(ER,EI,FR,FI,CDR,CDI);
                          CSQRT(CDR,CDI,FR,FI);
{                    STEP 4                                            }
                          CDR := FR; CDI := FI;
                          CADD(CBR,CBI,CDR,CDI,CER,CEI);
                          CSUB(CBR,CBI,CDR,CDI,FR,FI);
                          if (CABS(FR,FI) > CABS(CER,CEI)) then
                              CSUB(CBR,CBI,CDR,CDI,CER,CEI);
{                    STEP 5                                            }
                          CDIV(GR[3],GI[3],CER,CEI,ER,EI);
                          CMULT(ER,EI,-2.0,0.0,CH1R[3],CH1I[3]);
                          CADD(ZR[3],ZI[3],CH1R[3],CH1I[3],ZR[4],ZI[4])
                       end;
                    if (OK) then
                       begin
{                    evaluate f(x(i))=f(4) by Horner's method          }
                          GR[4]:= A[N]; GI[4] := 0.0;
                          for J := 2 to N do
                             begin
                                K := N-J+1;
                                CMULT(GR[4],GI[4],ZR[4],ZI[4],ER,EI);
                                GR[4] := ER+A[K];
                                GI[4] := EI
                             end;
                          writeln(OUP,i:4,' ',ZR[4]:14:8,' ',ZI[4]:14:8,' ',
                                  GR[4]:14:8,' ',GI[4]:14:8);
{                       STEP 6                                         }
                          if (CABS(CH1R[3],CH1I[3]) <= TOL) then
                             begin
                                writeln(OUP);
                                writeln(OUP,'Method Succeeds');
                                writeln(OUP,'Approximation is within ',TOL);
                                writeln(OUP,'in ',I,' iterations');
                                OK := false
                             end
                          else
{                       STEP 7                                         }
                             begin
                                for J := 1 to 2 do
                                   begin
                                      CH1R[J] := CH1R[J+1];
                                      CH1I[J] := CH1I[J+1];
                                      ZR[J] := ZR[J+1];
                                      ZI[J] := ZI[J+1];
                                      GR[J] := GR[J+1];
                                      GI[J] := GI[J+1]
                                   end;
                                ZR[3] := ZR[4];
                                ZI[3] := ZI[4];
                                GR[3] := GR[4];
                                GI[3] := GI[4];
                                CDEL1R[1] := CDEL1R[2];
                                CDEL1I[1] := CDEL1I[2];
                                CSUB(GR[3],GI[3],GR[2],GI[2],ER,EI);
                                CDIV(ER,EI,CH1R[2],CH1I[2],
                                     CDEL1R[2],CDEL1I[2]);
                                CSUB(CDEL1R[2],CDEL1I[2],
                                     CDEL1R[1],CDEL1I[1],ER,EI);
                                CADD(CH1R[2],CH1I[2],CH1R[1],
                                     CH1I[1],FR,FI);
                                CDIV(ER,EI,FR,FI,CDELR,CDELI)
                             end
                       end
                 end;
{         STEP 7 CONTINUED                                             }
              I := I + 1
           end;
{    STEP 8                                                            }
        if ( I > M ) and OK then
          writeln(OUP,'Method failed to give accurate approximation.');
        close(OUP)
     end
end.