Baixe Broyden-Algorithm-Non-Linear-System-Numerical-Analysis-MATLAB-Code e outras Notas de aula em PDF para Engenharia Aeronáutica, somente na Docsity! % 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. syms('OK', 'N', 'I', 'ZZ', 'J', 's', 'ss', 'TOL', 'NN', 'X'); syms('FLAG', 'NAME', 'OUP', 'A', 'V', 'B', 'I1', 'I2'); syms('C', 'K', 'SN', 'S', 'VV', 'Y', 'ZN', 'Z', 'P', 'U', 'KK'); TRUE = 1; FALSE = 0; fprintf(1,'This is the Broyden Method for Nonlinear Systems.\n'); fprintf(1,'The functions could be input or defined in code.\n'); fprintf(1,'This code assumes input of functions - see \n'); fprintf(1,'comments in code for alternate version.\n'); fprintf(1,'This program also uses M-files JAC.M and FN.M \n'); fprintf(1,'If the number of equations exceeds 7 then JAC.M\n'); fprintf(1,'and FN.M must be changed.\n'); OK = FALSE; while OK == FALSE fprintf(1,'Input the number n of equations.\n'); N = input(' '); if N >= 2 OK = TRUE; else fprintf(1,'N must be an integer greater than 1.\n'); end; end; for I = 1 : N fprintf(1,'Input the function F_(%d) in terms of y1...y%d \n',I,N); s(I) = input(' ','s'); end; % Define components of F as follows: % s(1) = '3*y1-cos(y2*y3)-0.5'; % s(2) = 'y1^2-81*(y2+0.1)^2+sin(y3)+1.06'; % s(3) = 'exp(-y1*y2)+20*y3+(10*pi-3)/3'; for I = 1 : N for J = 1 : N fprintf(1,'Input the partial of F_(%d) with respect to x_%d \n',I,J); fprintf(1,'in terms of y1, ..., y%d \n',N); ss((I-1)*N+J) = input(' ','s'); end; end; % Define the entries of the Jacobian in row major ordering. % ss(1) = '3'; % ss(2) = 'y3*sin(y2*y3)'; % ss(3) = 'y2*sin(y2*y3)'; % ss(4) = '2*y1'; % ss(5) = '-162*(y2+0.1)'; docsity.com % ss(6) = 'cos(y3)'; % ss(7) = '-y2*exp(-y1*y2)'; % ss(8) = '-y1*exp(-y1*y2)'; % ss(9) = '20'; OK = FALSE; while OK == FALSE fprintf(1,'Input tolerance\n'); TOL = input(' '); if TOL > 0 OK = TRUE; else fprintf(1,'Tolerance must be positive.\n'); end; end; OK = FALSE; while OK == FALSE fprintf(1,'Input the maximum number of iterations.\n'); NN = input(' '); if NN > 0 OK = TRUE; else fprintf(1,'Must be a positive integer.\n'); end; end; X = zeros(1,N); A = zeros(N,N); B = zeros(N,N); V = zeros(1,N); S = zeros(1,N); Y = zeros(1,N); U = zeros(1,N); Z = zeros(1,N); for I = 1 : N fprintf(1,'Input initial approximation X(%d).\n', I); X(I) = input(' '); end; if OK == TRUE fprintf(1,'Select output destination\n'); fprintf(1,'1. Screen\n'); fprintf(1,'2. Text file\n'); fprintf(1,'Enter 1 or 2\n'); FLAG = input(' '); if FLAG == 2 fprintf(1,'Input the file name in the form - drive:\\name.ext\n'); fprintf(1,'for example A:\\OUTPUT.DTA\n'); NAME = input(' ','s'); OUP = fopen(NAME,'wt'); else OUP = 1; end; fprintf(1,'Select amount of output\n'); fprintf(1,'1. Answer only\n'); fprintf(1,'2. All intermediate approximations\n'); fprintf(1,'Enter 1 or 2\n'); docsity.com end; % Note: V = F(X(K)) and Y = Y(K) % STEP 6 % Form Z = -Ay and norm ZN of Z. ZN = 0; for I = 1 : N Z(I) = 0; for J = 1 : N Z(I) = Z(I)-A(I,J)*Y(J); end; ZN = ZN+Z(I)*Z(I); end; ZN = sqrt(ZN); % Note = Z = -A(K-1)^(-1)*Y(K) % STEP 7 P = 0; % P will be S(K)^T*A(K)^(-1)*Y(K) for I = 1 : N P = P-S(I)*Z(I); end; % STEP 8 for I = 1 : N U(I) = 0; for J = 1 : N U(I) = U(I)+S(J)*A(J,I); end; end; % STEP 9 for I = 1 : N for J = 1 : N A(I,J) = A(I,J)+(S(I)+Z(I))*U(J)/P; end; end; % STEP 10 % Form S = -Av and norm SN of S. SN = 0; for I = 1 : N S(I) = 0; for J = 1 : N S(I) = S(I)-A(I,J)*V(J); end; SN = SN+S(I)^2; end; SN = sqrt(SN); % Note = A = A(K)^(-1) and S = -A(K)^(-1)*F(X(K)) % STEP 11 for I = 1 : N X(I) = X(I)+S(I);; end; % Note: X = X(K+1) KK = K+1; if FLAG == 2 fprintf(OUP, ' %2d', K); for I = 1 : N docsity.com fprintf(OUP, ' %11.8f', X(I)); end; fprintf(OUP, '\n%12.6e\n', SN); end; if SN <= TOL % procedure completed successfully OK = FALSE; fprintf(OUP, 'Iteration number %d', K); fprintf(OUP, ' gives solution:\n\n'); for I = 1 : N fprintf(OUP, ' %11.8f', X(I)); end; fprintf(OUP, '\n\nto within tolerance %.10e\n\n', TOL); fprintf(OUP, 'Process is complete\n'); else % STEP 13 K = KK; end; end; if K >= NN % STEP 14 fprintf(OUP, 'Procedure does not converge in %d iterations\n', NN); end; end; end; if OUP ~= 1 fclose(OUP); fprintf(1,'Output file %s created successfully \n',NAME); end; docsity.com