1. JJ-simplex
--Message-Boundary-23822
Content-type: text/plain; charset=US-ASCII
Content-description: Mail message body
Hi Arthur,
I have completed a short program for solving system of nonlinear
equations. It's of "first-try" quality, but it works. Look and try
simplex in nonlinear problem.
By Jan.
--Message-Boundary-23822
Content-type: text/plain; charset=US-ASCII
Content-description: Text from file 'nonlineq.EX'
-- NONLIN-SIMPLEX JJ/98 e-mail: jan.janovec at upce.cz
-- Calculation of the systems of (nonlinear) equations
-- Example: Equations X(1)*X(1) + 3*cos(X(2)) = 2
-- cos(X(1)) + 2*X(1)*X(2) = 4
-- Pragmatic function EF = E1*E1 + E2*E2 ... minimalised
-- where:
-- E1 = Z(1)*Z(1) + 3*cos(Z(2)) - 2
-- E2 = cos(Z(1)) + 2*Z(1)*Z(2) - 4
include get.e
include graphics.e
sequence C,
E,
P,
R,
X,
Z,
str
atom EF,
E1,
E2,
E3
integer N,
CL,
H,
L,
H2,
Poc,
FlgV,
FlgC,
key
Poc = 0
FlgV = 1
FlgC = 0
-- ***************************************************************
procedure Pragfun()
atom E1, E2 -- ...
E1 = Z[1]*Z[1] + 3*cos(Z[2]) - 2
E2 = cos(Z[1]) + 2*Z[1]*Z[2] - 4
EF = E1*E1 + E2*E2
end procedure
------------------------------------------------------------------
function inputs(sequence text)
object inps
inps=""
text_color(YELLOW)
puts (1,text)
inps = gets (0)
inps = inps[1..(length(inps)-1)]
text_color(WHITE)
return inps
end function
function input(sequence text)
sequence inp, cur
inp={1,0}
text_color(YELLOW)
puts (1,text)
while inp[1] do
cur=get_position()
inp = get (0)
position (cur[1],cur[2])
end while
text_color(WHITE)
return inp[2]
end function
function abs(atom A)
if A<0 then
A=-1*A
end if
return A
end function
-------------------------------------------------------------------
N = input("\n Number of variables (Ex.=>2) ? ")
-- Dimension
C = repeat(0,N)
E = repeat(0,N+1)
P = E
for I=1 to N+1 do
P[I] = repeat(0,N)
end for
R = C
X = C
Z = C
str = {}
procedure First_estimation()
puts(1,"\n")
puts(1,"First estimation:")
for I=1 to N do
puts(1,"\n")
str = " Z(" & sprintf("%2d", I) & ") = "
X[I] = input(str)
Z[I] = X[I]
end for
CL = 0
Pragfun()
E[1] = EF
E1 = EF
puts(1,"\n")
for I=1 to N do
puts(1,"Z(")
printf(1, "%2d)", I)
printf(1, "= %.4f\n", X[I])
end for
printf(1, "EF = %.8f\n", EF)
puts(1,"\n")
end procedure
procedure Start_val()
for J=1 to N do
P[1][J] = X[J]
end for
for I=2 to N+1 do
for J=1 to N do
P[I][J] = X[J]
end for
P[I][I-1] = 1.1 * X[I-1]
if abs(X[I-1]) < 1e-10 then
P[I][I-1] = 1e-4
end if
end for
end procedure
procedure PL_PH()
H = 1
L = 1
for I=1 to N+1 do
for J=1 to N do
X[J] = P[I][J]
Z[J] = X[J]
end for
Pragfun()
E[I] = EF
if E[I]<E[L] then
L = I
if E[I]>E[L] then
H = I
end if
end if
end for
end procedure
procedure PH2()
H2 = L
for I=1 to N+1 do
if E[I] >= E[H2] and I != H then
H2 = I
end if
end for
end procedure
procedure Centroid()
for J=1 to N do
C[J] = - P[H][J]
for I=1 to N+1 do
C[J] = C[J] + P[I][J]
end for
C[J] = C[J] / N
end for
end procedure
procedure Rebound()
for J=1 to N do
R[J] = 1.9985 * C[J] - 0.9985 * P[H][J]
Z[J] = R[J]
end for
Pragfun()
E2 = EF
end procedure
procedure Extension()
L = H
H = H2
for J=1 to N do
X[J] = 1.95 * R[J] - 0.95 * C[J]
Z[J] = X[J]
end for
Pragfun()
E3 = EF
end procedure
procedure Prnt()
puts(1,"\n")
for J=1 to N do
puts(1,"Z(")
printf(1, "%2d)", J)
printf(1, "= %.8f\n", P[L][J])
end for
printf(1, "EF = %.8f", E[L])
puts(1,"\n")
end procedure
procedure Truncation()
for J=1 to N do
R[J] = 0.5015 * C[J] + 0.4985 * P[H][J]
Z[J] = R[J]
end for
Pragfun()
E2 = EF
end procedure
procedure Reducing()
for I=1 to N+1 do
for J=1 to N do
P[I][J] = P[I][J] + 0.5 * ( P[L][J] - P[I][J] )
end for
end for
end procedure
--***************************************************
First_estimation()
Start_val()
while FlgV do
PL_PH()
while FlgV do
while FlgV do
PH2()
Centroid()
while FlgV do
Rebound()
if E2 < E[L] then
FlgC = 1
exit
end if
if E2 >= E[H] then
Truncation()
if E2 < E[L] then
FlgC = 1
exit
end if
if E2 >= E[H] then
Reducing()
end if
end if
for J=1 to N do
P[H][J] = R[J]
end for
E[H] = E2
if E2 <= E[H2] then
exit
end if
end while
if FlgC then
FlgC = 0
exit
end if
H = H2
end while
Extension()
if E3<E2 then
for J=1 to N do
P[L][J] = X[J]
end for
E[L] = E3
else
for J=1 to N do
P[L][J] = R[J]
end for
E[L] = E2
end if
Prnt()
if abs(E[L]-E1) < 1e-8 or E[L] < 1e-30 then
FlgV = 0
printf(1, "Number of iterations = %4d\n", Poc)
exit
end if
E1 = E[L]
Poc = Poc + 1
if integer(0.1*Poc) then
printf(1, "Number of iterations = %4d\n", Poc)
end if
CL = CL +1
if CL>1000 then
puts(1," Isn't convergency - stood !!!")
FlgV = 0
exit
end if
end while
end while
key = wait_key()
--Message-Boundary-23822--
