# Broyden-Algorithm-Non-Linear-System-Numerical-Analysis-MATLAB-Code, Notas de aula de Engenharia Aeronáutica e Programação de Computadores. Universidade Federal de Santa Catarina (UFSC)

PDF (121 KB)
6 páginas
907Número de visitas
Descrição
engenharia
20pontos
este documento
Pré-visualização3 páginas / 6

% 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,'2. All intermediate approximations\n');

fprintf(1,'Enter 1 or 2\n');

docsity.com

FLAG = input(' ');

fprintf(OUP, 'BROYDENS METHOD FOR NONLINEAR SYSTEMS\n\n');

if FLAG == 2

fprintf(OUP, 'Iteration, Approximation, Error\n');

end;

% STEP 1

% A will hold the Jacobian for the initial approximation.

for I = 1 : N

for J = 1 : N

ZZ = JAC(I,J,N,X,ss);

A(I,J) = ZZ;

end;

% Compute V = F(x(0))

V(I) = FN(I,N,X,s);

end;

% STEP 2

% Invert the Jacobian.

for I = 1 : N

for J = 1 : N

B(I,J) = 0;

end;

B(I,I) = 1;

end;

I = 1;

while I <= N & OK == TRUE

I1 = I+1;

I2 = I;

if I ~= N

C = abs(A(I,I));

for J = I1 : N

if abs(A(J,I)) > C

I2 = J;

C = abs(A(J,I));

end;

end;

if C <= 1.0e-20

OK = FALSE;

else

if I2 ~= I

for J = 1 : N

C = A(I,J);

A(I,J) = A(I2,J);

A(I2,J) = C;

C = B(I,J);

B(I,J) = B(I2,J);

B(I2,J) = C;

end;

end;

end;

else

if abs(A(N,N)) <= 1.0e-20

OK = FALSE;

end;

end;

docsity.com

if OK == TRUE

for J = 1 : N

if J ~= I

C = A(J,I)/A(I,I);

for K = 1 : N

A(J,K) = A(J,K)-C*A(I,K);

B(J,K) = B(J,K)-C*B(I,K);

end;

end;

end;

end;

I = I+1;

end;

if OK == TRUE

for I = 1 : N

C = A(I,I);

for J = 1 : N

A(I,J) = B(I,J)/C;

end;

end;

else

fprintf(1,'Jacobian has no inverse\n');

end;

if OK == TRUE

% STEP 3

K = 2;

% Note: S = S(1)

% Compute the product S = -Av and the L2 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);

for I = 1 : N

X(I) = X(I)+S(I);

end;

if FLAG == 2

fprintf(OUP,' %d',K-1);

for I = 1 : N

fprintf(OUP,' %11.8f',X(I));

end;

fprintf(OUP,'\n %12.6e\n',SN);

end;

% STEP 4

while K <= NN & OK == TRUE

% STEP 5

for I = 1 : N

VV = FN(I,N,X,s);

Y(I) = VV-V(I);

V(I) = VV;

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

comentários (0)
Até o momento nenhum comentário
Seja o primeiro a comentar!