alg102.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> restart;
> # BROYDEN ALGORITHM 10.2
> #
> # To approximate the solution of the nonlinear system F(X) = 0
> # given an initial approximation X.
> #
> # INPUT:   Number n of equations and unknowns; initial
> #          approximation X = (X(1),...,X(n)); tolerance TOL;
> #          maximum number of iterations N.
> #
> # OUTPUT:  Approximate solution X = (X(1),...,X(n)) or a message
> #          that the number of iterations was exceeded.
> alg102 := proc() local OK, N, I, F, J, PD, TOL, NN, X, FLAG, NAME, OUP, A, V, B, I1, I2, C, K, SN, S, VV, Y, ZN, Z, P, U, KK;
> printf(`This is the Broyden Method for Nonlinear Systems.\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number n of equations.\n`);
> N := scanf(`%d`)[1];
> if N >= 2 then
> OK := TRUE;
> else
> printf(`N must be an integer greater than 1.\n`);
> fi;
> od;
> for I from 1 to N do
> printf(`Input the function F%d in terms of x1...x%d \n` ,I ,N);
> F[I] := scanf(`%a`)[1];
> od;
> for I from 1 to N do
> for J from 1 to N do
> PD[I,J] := diff(F[I],evaln(x || J));
> PD[I,J] := unapply(PD[I,J],evaln(x || (1..N)));
> od;
> od;
> for I from 1 to N do
> F[I] := unapply(F[I],evaln(x || (1..N)));
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input tolerance\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be positive.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Must be a positive integer.\n`);
> fi;
> od;
> for I from 1 to N do
> printf(`Input initial approximation X(%d).\n`, I);
> X[I-1] := scanf(`%f`)[1];
> od;
> if OK = TRUE then
> printf(`Select output destination\n`);
> printf(`1. Screen\n`);
> printf(`2. Text file\n`);
> printf(`Enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> printf(`Select amount of output\n`);
> printf(`1. Answer only\n`);
> printf(`2. All intermeditate approximations\n`);
> printf(`Enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> fprintf(OUP, `BROYDENS METHOD FOR NONLINEAR SYSTEMS\n\n`);
> if FLAG = 2 then
> fprintf(OUP, `Iteration, Approximation, Error\n`);
> fi;
> # Step 1
> # A will hold the Jacobian for the initial approximation
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := evalf(PD[I,J](seq(X[i],i=0..N-1)));
> od;
> # Compute V = F(X(0))
> V[I-1] := evalf(F[I](seq(X[i],i=0..N-1)));
> od;
> # Step 2
> for I from 1 to N do
> for J from 1 to N do
> B[I-1,J-1] := 0;
> od;
> B[I-1,I-1] := 1;
> od;
> I := 1;
> while I <= N and OK = TRUE do
> I1 := I+1;
> I2 := I;
> if I <> N then
> C := abs(A[I-1,I-1]);
> for J from I1 to N do
> if abs(A[J-1,I-1]) > C then
> I2 := J;
> C := abs(A[J-1,I-1]);
> fi;
> od;
> if C <= 1.0e-20 then
> OK := FALSE;
> else
> if I2 <> I then
> for J from 1 to N do
> C := A[I-1,J-1];
> A[I-1,J-1] := A[I2-1,J-1];
> A[I2-1,J-1] := C;
> C := B[I-1,J-1];
> B[I-1,J-1] := B[I2-1,J-1];
> B[I2-1,J-1] := C;
> od;
> fi;
> fi;
> else
> if abs(A[N-1,N-1]) <= 1.0e-20 then
> OK := FALSE;
> fi;
> fi;
> if OK = TRUE then
> for J from 1 to N do
> if J <> I then
> C := A[J-1,I-1]/A[I-1,I-1];
> for K from 1 to N do
> A[J-1,K-1] := A[J-1,K-1]-C*A[I-1,K-1];
> B[J-1,K-1] := B[J-1,K-1]-C*B[I-1,K-1];
> od;
> fi;
> od;
> fi;
> I := I+1;
> od;
> if OK = TRUE then
> for I from 1 to N do
> C := A[I-1,I-1];
> for J from 1 to N do
> A[I-1,J-1] := B[I-1,J-1]/C;
> od;
> od;
> else
> printf(`Jacobian has no inverse\n`);
> fi;
> if OK = TRUE then
> # Step 3
> K := 2;
> # Note: S = S(1)
> SN := 0;
> for I from 1 to N do
> S[I-1] := 0;
> for J from 1 to N do
> S[I-1] := S[I-1]-A[I-1,J-1]*V[J-1];
> od;
> SN := SN+S[I-1]^2;
> od;
> SN := sqrt(SN);
> for I from 1 to N do
> X[I-1] := X[I-1]+S[I-1];
> od;
> if FLAG = 2 then
> fprintf(OUP,` %d`,K-1);
> for I from 1 to N do
> fprintf(OUP,` %11.8f`,X[I-1]);
> od;
> fprintf(OUP,`\n %12.6e\n`,SN);
> fi;
> # Step 4
> while K < NN and OK = TRUE do
> # Step 5
> for I from 1 to N do
> VV := evalf(F[I](seq(X[i],i=0..N-1)));
> Y[I-1] := VV-V[I-1];
> V[I-1] := VV;
> od;
> # Note:  V = F(X(K)) and Y = Y(K)
> # Step 6
> ZN := 0;
> for I from 1 to N do
> Z[I-1] := 0;
> for J from 1 to N do
> Z[I-1] := Z[I-1]-A[I-1,J-1]*Y[J-1];
> od;
> ZN := ZN+Z[I-1]*Z[I-1];
> od;
> ZN := sqrt(ZN);
> # Note:  Z = -A(K-1)^(-1)*Y(K)
> # Step 7
> P := 0;
> # P will be S(K)^T*A(K-1)^(-1)*Y(K)
> for I from 1 to N do
> P := P-S[I-1]*Z[I-1];
> od;
> # Step 8
> for I from 1 to N do
> U[I-1] := 0;
> for J from 1 to N do
> U[I-1] := U[I-1]+S[J-1]*A[J-1,I-1];
> od;
> # Step 9
> od;
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := A[I-1,J-1]+(S[I-1]+Z[I-1])*U[J-1]/P;
> od;
> od;
> # Step 10
> SN := 0;
> for I from 1 to N do
> S[I-1] := 0;
> for J from 1 to N do
> S[I-1] := S[I-1]-A[I-1,J-1]*V[J-1];
> od;
> SN := SN+S[I-1]^2;
> od;
> SN := sqrt(SN);
> # Note: A = A(K)^(-1) and S = -A(K)^(-1)*F(X(K))
> # Step 11
> for I from 1 to N do
> X[I-1] := X[I-1]+S[I-1];;
> od;
> # Note:  X = X(K+1)
> KK := K+1;
> if FLAG = 2 then
> fprintf(OUP, ` %2d`, K);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n%12.6e\n`, SN);
> fi;
> if SN <= TOL then
> # Procedure completed successfully
> OK := FALSE;
> fprintf(OUP, `Iteration number %d`, K);
> fprintf(OUP, ` gives solution:\n\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n\nto within tolerance %.10e\n\n`, TOL);
> fprintf(OUP, `Process is complete\n`);
> else
> # Step 13
> K := KK;
> fi;
> od;
> if K >= NN then
> # Step 14
> fprintf(OUP, `Procedure does not converge in %d iterations\n`, NN);
> fi;
> fi;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created sucessfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg102();