restart:

# CONTINUATION METHOD FOR SYSTEMS ALGORITHM 104

# 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 Runge-Kutta 4 iterations N.

# OUTPUT: Approximate solution X=(X(1),...,X(n)).

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

K := N-1:

OK := TRUE:

I2 := 1:

while OK = TRUE and I2 <= K do

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

IR := I2:

IA := I2+1:

for J from IA to N do

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

IR := J:

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

fi:

od:

if Z <= 1.0e-20 then

OK := FALSE:

else

if IR <> I2 then

for J from I2 to N+1 do

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

A[I2-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,I2-1]/A[I2-1,I2-1]:

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

C := 0:

fi:

for L from I2 to N+1 do

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

od:

fi:

I2 := I2+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 I2 from 1 to K do

J := N-I2:

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 Continuation 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.\134n`):

fi:

od:

for I2 from 1 to N do

print(`Input the function F_I in terms of x1 ... xn for I = `,I2):

F[I2] := scanf(`%a`)[1]: print(`Function is `): print(F[I2]):

od:

for I2 from 1 to N do

for J from 1 to N do

P[I2,J] := diff(F[I2],evaln(x || J)):

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

od:

for I2 from 1 to N do

F[I2] := unapply(F[I2],evaln(x || (1..N))):

od:

OK := FALSE:

while OK = FALSE do

print(`Input the number N for RK4.`):

NN := scanf(`%d`)[1]: print(`N = `):print(NN):

if NN > 0 then

OK := TRUE:

else

print(`Must be a positive integer.`):

fi:

od:

for I2 from 1 to N do

print(`Input initial approximation X(I) for I = `, I2):

X[I2-1] := scanf(`%f`)[1]: print(`Input is `): print(X[I2-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(`Input 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:

fprintf(OUP,`This is the Continuation Method using the Runge Kutta Method of Order 4 \134n\134n`):

fprintf(OUP,`N = %2d \134n\134n`,NN):

fprintf(OUP,`Iter# Approximation \134n\134n`):

K:=0:

fprintf(OUP, ` %2d`, K):

for I2 from 1 to N do

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

od:

fprintf(OUP,` \134n`):

# Step 1

H := 1/NN:

for I2 from 1 to N do

b[I2-1] := H*evalf(-F[I2](seq(X[i],i=0..N-1))):

od:

# Step 2

for K from 1 to NN do

# Steps 3 - 6

for I2 from 1 to N do

for J from 1 to N do

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

od:

for I2 from 1 to N do

A[I2-1,N] := b[I2-1]:

od:

LINSYS(N,OK,A,Y):

if OK = FALSE then

break:

fi:

for I2 from 1 to N do

K1[I2-1] := Y[I2-1]:

X1[I2-1] := X[I2-1]+0.5*K1[I2-1]:

od:

for I2 from 1 to N do

for J from 1 to N do

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

od:

for I2 from 1 to N do

A[I2-1,N] := b[I2-1]:

od:

LINSYS(N,OK,A,Y):

if OK = FALSE then

break:

fi:

for I2 from 1 to N do

K2[I2-1] := Y[I2-1]:

X2[I2-1] := X[I2-1]+0.5*K2[I2-1]:

od:

for I2 from 1 to N do

for J from 1 to N do

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

od:

for I2 from 1 to N do

A[I2-1,N] := b[I2-1]:

od:

LINSYS(N,OK,A,Y):

if OK = FALSE then

break:

fi:

for I2 from 1 to N do

K3[I2-1] := Y[I2-1]:

X3[I2-1] := X[I2-1]+K3[I2-1]:

od:

for I2 from 1 to N do

for J from 1 to N do

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

od:

for I2 from 1 to N do

A[I2-1,N] := b[I2-1]:

od:

LINSYS(N,OK,A,Y):

if OK = FALSE then

break:

fi:

# Step 7

for I2 from 1 to N do

K4[I2-1] := Y[I2-1]:

X4[I2-1] := X[I2-1]+(K1[I2-1]+2*K2[I2-1]+2*K3[I2-1]+K4[I2-1])/6:

X[I2-1] := X4[I2-1]:

od:

fprintf(OUP, ` %2d`, K):

for I2 from 1 to N do

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

od:

fprintf(OUP,` \134n`):

od:

# Step 8

if OUP <> default then

fclose(OUP):

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

fi:

SVdUaGlzfmlzfnRoZX5Db250aW51YXRpb25+TWV0aG9kfmZvcn5Ob25saW5lYXJ+U3lzdGVtcy5HNiI=

SUFJbnB1dH50aGV+bnVtYmVyfm5+b2Z+ZXF1YXRpb25zLkc2Ig==

SSVufj1+RzYi

IiIk

NiRJVklucHV0fnRoZX5mdW5jdGlvbn5GX0l+aW5+dGVybXN+b2Z+eDF+Li4ufnhufmZvcn5Jfj1+RzYiIiIi

SS1GdW5jdGlvbn5pc35HNiI=

LChJI3gxRzYiIiIkLUkkY29zRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YkNiMqJkkjeDJHRiQiIiJJI3gzR0YkRi4hIiIkISImRjBGLg==

NiRJVklucHV0fnRoZX5mdW5jdGlvbn5GX0l+aW5+dGVybXN+b2Z+eDF+Li4ufnhufmZvcn5Jfj1+RzYiIiIj

SS1GdW5jdGlvbn5pc35HNiI=

LCoqJEkjeDFHNiIiIiMiIiIqJCwmSSN4MkdGJUYnJEYnISIiRidGJiEjIiktSSRzaW5HNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2I0kjeDNHRiVGJyQiJDEiISIjRic=

NiRJVklucHV0fnRoZX5mdW5jdGlvbn5GX0l+aW5+dGVybXN+b2Z+eDF+Li4ufnhufmZvcn5Jfj1+RzYiIiIk

SS1GdW5jdGlvbn5pc35HNiI=

LCgtSSRleHBHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHNiI2IywkKiZJI3gxR0YoIiIiSSN4MkdGKEYtISIiRi1JI3gzR0YoIiM/JCIrOGIoPloqISIqRi0=

STxJbnB1dH50aGV+bnVtYmVyfk5+Zm9yflJLNC5HNiI=

SSVOfj1+RzYi

IiIl

NiRJSklucHV0fmluaXRpYWx+YXBwcm94aW1hdGlvbn5YKEkpfmZvcn5Jfj1+RzYiIiIi

SSpJbnB1dH5pc35HNiI=

JCIiIUYj

NiRJSklucHV0fmluaXRpYWx+YXBwcm94aW1hdGlvbn5YKEkpfmZvcn5Jfj1+RzYiIiIj

SSpJbnB1dH5pc35HNiI=

JCIiIUYj

NiRJSklucHV0fmluaXRpYWx+YXBwcm94aW1hdGlvbn5YKEkpfmZvcn5Jfj1+RzYiIiIk

SSpJbnB1dH5pc35HNiI=

JCIiIUYj

STpTZWxlY3R+b3V0cHV0fmRlc3RpbmF0aW9uRzYi

SSoxLn5TY3JlZW5HNiI=

SS0yLn5UZXh0fmZpbGVHNiI=

SS1FbnRlcn4xfm9yfjJHNiI=

SSpJbnB1dH5pc35HNiI=

IiIi

This is the Continuation Method using the Runge Kutta Method of Order 4 N = 4 Iter# Approximation 0 0.00000000 0.00000000 0.00000000 1 0.12499997 -0.00329005 -0.13092026 2 0.24999977 -0.00450740 -0.26185576 3 0.37499970 -0.00343035 -0.39276344 4 0.50000000 0.00000001 -0.52359878

JSFH