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