1. JJ-simplex
- Posted by "J. Janovec" <Jan.Janovec at UPCE.CZ> Oct 29, 1998
- 386 views
--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--