alg104.m

Numerical Anaysis 8th Edition by Burden and Faires

% CONTINUATION METHOD FOR SYSTEMS ALGORITHM 10.1
%
% 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)); number of RK4 steps N.
%
% OUTPUT:  Approximate solution X=(X(1),...,X(n)) or a message
%          that the number of iterations was exceeded.
 syms('OK', 'N', 'I', 'J', 'P', 'NN', 'X','ZZ');
 syms('FLAG', 'NAME', 'OUP', 'K', 'A', 's', 'ss','mm');
 syms('KK1', 'I1', 'Z1', 'IR1', 'IA1', 'J1', 'C1', 'L1', 'JA1');
 syms('K1','K2','K3','K4','X1','KK','H','b');
 TRUE = 1;
 FALSE = 0;
 fprintf(1,'This is the Continuation 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;
 A = zeros(N,N+1);
 K1 = zeros(4,N);
 X = zeros(1,N);
 Y = zeros(1,N);
 X1 = zeros(1,N);
 b = zeros(1,N);
 mm = zeros(1,4);
 mm(1) = 0.5;
 mm(2) = 0.5;
 mm(3) = 1.0;
 mm(4) = 0;
% 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
% fprintf(1,'Input the function F_(%d) in terms of y1 ... y%d \n' ,I ,N);
% s(I) = input(' ','s');
% end;
% 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)';
 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 the number of RK4 steps.\n');
 NN = input(' ');
 if NN > 0
 OK = TRUE;
 else
 fprintf(1,'Must be a positive integer.\n');
 end;
 end;
 for I = 1 : N
 fprintf(1,'Input initial approximation X(%d).\n', I);
 X(I) = input(' ');
 end;
 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;
% STEP 1
 H = 1/NN;
 for I = 1 : N
 ZZ = -H*FN(I,N,X,s);
 b(I) = ZZ;
 end;
% STEP 2
 K = 1;
 while OK == TRUE & K <= NN
% STEPS 3 - 6
 for I = 1 : N
 X1(I) = X(I);
 end;
 KK = 1;
 while OK == TRUE & KK <= 4
 for I = 1 : N
 for J = 1 : N
 ZZ = JAC(I,J,N,X1,ss);
 A(I,J) = ZZ;
 end;
 A(I,N+1) = b(I);
 end;
 KK1 = N-1;
 OK = TRUE;
 I1 = 1;
 while OK == TRUE & I1 <= KK1
 Z1 = abs(A(I1,I1));
 IR1 = I1;
 IA1 = I1+1;
 for J1 = IA1 : N
 if abs(A(J1,I1)) > Z1
 IR1 = J1;
 Z1 = abs(A(J1,I1));
 end;
 end;
 if Z1 <= 1.0e-20
 OK = FALSE;
 else
 if IR1 ~= I1
 for J1 = I1 : N+1
 C1 = A(I1,J1);
 A(I1,J1) = A(IR1,J1);
 A(IR1,J1) = C1;
 end;
 end;
 for J1 = IA1 : N
 C1 = A(J1,I1)/A(I1,I1);
 if abs(C1) <= 1.0e-20
 C1 = 0;
 end;
 for L1 = I1 : N+1
 A(J1,L1) = A(J1,L1)-C1*A(I1,L1);
 end;
 end;
 end;
 I1 = I1+1;
 end;
 if OK == TRUE
 if abs(A(N,N)) <= 1.0e-20
 OK = FALSE;
 else
 Y(N) = A(N,N+1)/A(N,N);
 for I1 = 1 : KK1
 J1 = N-I1;
 JA1 = J1+1;
 C1 = A(J1,N+1);
 for L1 = JA1 : N
 C1 = C1-A(J1,L1)*Y(L1);
 end;
 Y(J1) = C1/A(J1,J1);
 end;
 end;
 end;
 if OK == FALSE
 fprintf(1,'Linear system is singular\n');
 end;
 if OK == TRUE
 for I = 1 : N
 K1(KK,I) = Y(I);
 X1(I) =  X(I) + mm(KK)*K1(KK,I);
 end;
 KK = KK + 1;
 end;
 end;
% STEP 7
 if OK == TRUE
 for I = 1 : N
 X(I) = X(I) + (K1(1,I)+2*K1(2,I)+2*K1(3,I)+K1(4,I))/6;
 end;
 fprintf(OUP, ' %2d',  K);
 for I = 1 : N
 fprintf(OUP, ' %11.8f', X(I));
 end;
 fprintf(OUP, '\n');
 end;
 K = K + 1;
 end;
% STEP 8
 if OUP ~= 1
 fclose(OUP);
 fprintf(1,'Output file %s created successfully \n',NAME);
 end;