restart:

# NEWTON'S 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)): 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.

LINSYS := proc(N,OK,A,Y) local K, I1, Z, IR, IA, J, C, L, JA:

K := N-1:

OK := TRUE:

I1 := 1:

while OK = TRUE and I1 <= K do

Z := abs(A[I1-1,I1-1]):

IR := I1:

IA := I1+1:

for J from IA to N do

if abs(A[J-1,I1-1]) > Z then

IR := J:

Z := abs(A[J-1,I1-1]):

fi:

od:

if Z <= 1.0e-20 then

OK := FALSE:

else

if IR <> I1 then

for J from I1 to N+1 do

C := A[I1-1,J-1]:

A[I1-1,J-1] := A[IR-1,J-1]:

A[IR-1,J-1] := C:

od:

fi:

for J from IA to N do

C :=A[J-1,I1-1]/A[I1-1,I1-1]:

if abs(C) <= 1.0e-20 then

C := 0:

fi:

for L from I1 to N+1 do

A[J-1,L-1] := A[J-1,L-1]-C*A[I1-1,L-1]:

od:

fi:

I1 := I1+1:

od:

if OK = TRUE

then if abs(A[N-1,N-1]) <= 1.0e-20 then

OK := FALSE:

else

Y[N-1] := A[N-1,N]/A[N-1,N-1]:

for I1 from 1 to K do

J := N-I1:

JA := J+1:

C := A[J-1,N]:

for L from JA to N do

C := C-A[J-1,L-1]*Y[L-1]:

od:

Y[J-1] := C/A[J-1,J-1]:

od:

fi:

if OK = FALSE then

print(`Linear system is singular`):

fi:

end:

print(`This is the Newton 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

P[I1,J] := diff(F[I1],evaln(x || J)):print(`partial derivative i,j `):print(I1,J):print(P[I1,J]):

P[I1,J] := unapply(P[I1,J],evaln(x || (1..N))):

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 the 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 do

print(`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(`Response is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive\134\134:name.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\134n`):

print(`1. Answer only\134n`):

print(`2. All intermeditate approximations\134n`):

print(`Enter 1 or 2\134n`):

FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):

fprintf(OUP, `NEWTONS METHOD FOR NONLINEAR SYSTEMS\134n\134n`):

if FLAG = 2 then

fprintf(OUP, `Iteration, Approximation, Error\134n`):

fi:

# Step 1

K := 1:

# Step 2

while OK = TRUE and K <= NN do

# Step 3

for I1 from 1 to N do

for J from 1 to N do

A[I1-1,J-1] := evalf(P[I1,J](seq(X[i],i=0..N-1))):

od:

A[I1-1,N] := evalf(-F[I1](seq(X[i],i=0..N-1))):

od:

# Step 4

LINSYS(N,OK,A,Y):

if OK = TRUE then

# Step 5

R := 0:

for I1 from 1 to N do

if abs(Y[I1-1]) > R then

R := abs(Y[I1-1]):

fi:

X[I1-1] := X[I1-1]+Y[I1-1]:

od:

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`, R):

fi:

# Step 6

if R < TOL then

OK := FALSE:

fprintf(OUP, `Iteration %d gives solution:\134n\134n`, K):

for I1 from 1 to N do

fprintf(OUP, ` %11.8f`, X[I1-1]):

od:

fprintf(OUP, `\134n\134nto within tolerance %.10e\134n`, TOL):

# Step 7

else

K := K+1:

fi:

od:

if K > NN then

# Step 8

fprintf(OUP, `Procedure does not converge in %d iterations\134n`, NN):

fi:

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created sucessfully`):

fi:

SVFUaGlzfmlzfnRoZX5OZXd0b25+TWV0aG9kfmZvcn5Ob25saW5lYXJ+U3lzdGVtcy5HNiI=

SUFJbnB1dH50aGV+bnVtYmVyfm5+b2Z+ZXF1YXRpb25zLkc2Ig==

SSVufj1+RzYi

IiIk

SWJvSW5wdXR+dGhlfmZ1bmN0aW9ufmZbaV0oeDEseDIsfi4uLn54bil+aW5+dGVybXN+b2Z+eDEsfngyLH4uLi5+Zm9yfml+PX5HNiI=

IiIi

STdGb3J+ZXhhbXBsZTp+eDEteDJeMisxRzYi

SS1mdW5jdGlvbn5pc35HNiI=

LChJI3gxRzYiIiIkLUkkY29zRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YkNiMqJkkjeDJHRiQiIiJJI3gzR0YkRi4hIiIkISImRjBGLg==

SWJvSW5wdXR+dGhlfmZ1bmN0aW9ufmZbaV0oeDEseDIsfi4uLn54bil+aW5+dGVybXN+b2Z+eDEsfngyLH4uLi5+Zm9yfml+PX5HNiI=

IiIj

STdGb3J+ZXhhbXBsZTp+eDEteDJeMisxRzYi

SS1mdW5jdGlvbn5pc35HNiI=

LCoqJEkjeDFHNiIiIiMiIiIqJCwmSSN4MkdGJUYnJEYnISIiRidGJiEjIiktSSRzaW5HNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2I0kjeDNHRiVGJyQiJDEiISIjRic=

SWJvSW5wdXR+dGhlfmZ1bmN0aW9ufmZbaV0oeDEseDIsfi4uLn54bil+aW5+dGVybXN+b2Z+eDEsfngyLH4uLi5+Zm9yfml+PX5HNiI=

IiIk

STdGb3J+ZXhhbXBsZTp+eDEteDJeMisxRzYi

SS1mdW5jdGlvbn5pc35HNiI=

LCgtSSRleHBHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHNiI2IywkKiZJI3gxR0YoIiIiSSN4MkdGKEYtISIiRi1JI3gzR0YoIiM/JCIrOGIoPloqISIqRi0=