Solve Systems of Non-Linear Equations

[rawi]

The program solves systems of non-linear equations. Equations are used in a form where the result should be zero. Nelder Mead Optimization is used to find a minimum of the sum of squared results. If the minimum is zero, the solution is found.

Example (taken from a YouTube video):
A quadrilateral is given with the following characteristics:
A, B, C, and D are located on a circle line. The length of the sides is as follows:
AB: 4, BC: 5, CD: 6, DA: 3.
What is the radius R of the circle?

The center of the circle is set to (0,0); A is set to (-R, 0). Then we get the following system of equations of 7 equations with 7 unknowns ( (X1,Y1) coordinates of A, (X2,Y2) of B etc.):
1) (-R-X2)^2+Y2^2=16
2) (X2-X3)^2+(Y2-Y3)^2=25
3) (X3-X4)^2+(Y3-Y4)^2=36
4) (X4+R)^2+Y4^2=9
5) X2^2+Y2^2=R^2
6) X3^2+Y3^2=R^2
7) X4^2+Y4^2=R^2

Application:
Parameters R, X2, X3, X4, Y2, Y3, Y4 are stored in List L1. From a rough drawing, we get the following initial estimates: R=3, X2=1; X3=3; X4=-3; Y2=-1; Y3=4; Y4=3.
We set: {3,1,3,-3,-1,4,3) ▶ L1
Further, the function FUNES is needed with the equations. Equations are changed, so that the results are 0. Results are stored in List L0.

Code:

EXPORT FUNES()
BEGIN
LOCAL Q1,I,S1;
// PUT HERE IN EQUATIONS f(L1(1),...L1(n))=0
(-L1(1)-L1(2))^2+L1(5)^2-15▶L0(1);
(L1(2)-L1(3))^2+(L1(5)-L1(6))^2-25▶L0(2);
(L1(3)-L1(4))^2+(L1(6)-L1(7))^2-36▶L0(3);
(L1(4)+L1(1))^2+L1(7)^2-9▶L0(4);
L1(2)^2+L1(5)^2-L1(1)^2▶L0(5);
L1(3)^2+L1(6)^2-L1(1)^2▶L0(6);
L1(4)^2+L1(7)^2-L1(1)^2▶L0(7);
// End of system of equations
SIZE(L1)▶S1;
L0(1)^2▶Q1;
FOR I FROM 2 TO S1 DO
Q1+L0(I)^2▶Q1;
END;
√Q1▶Q1;
RETURN Q1;
END;

Start the program NLES. You are asked for the following parameters:
S0: Initial step width of Nelder Mead Optimization. We start with 1.
E0: Termination criterion. Optimization ends if the square root of the sum of squared results of the equations is less than E0, the step length of the Nelder Mead Procedure is less than E0/10, or the maximum number of iterations is reached. We put in 0.0001.
ITM: Maximum number of iterations. We put in 500.

The interim result of the minimum is shown every 25 iterations. After a few seconds, you get the following information on the screen:

--- DONE -----------
Iter: 478
Eps: 8,4281055292E-6
Fmin: 1.70417577211E-4
Parameter: L1
Parameter L2

Interpretation:
Optimization stopped after 478 iterations because of the limit of step length of the Nelder Mead procedure, not because zero was reached (Fmin>1E-4).

So it might make sense to start a new one with a lower step length.
We take S0=0.25, E0=0.0001, ITM=400. After 46 iterations we get Fmin=1.385E-4, so our criterion is nearly reached and we stop.

Parameters are in list L1.
L1(1)=3.2704; this is the estimate for the radius of the circle.
Inspection of list L0 shows that the biggest deviation in the equations is 8.4E-5.

Code:

EXPORT NLES(S0,E0,ITM)
BEGIN
// Function with Equations: FUNES()
// S0 Starting value of simplex size
// E0 Termination criterion of simplex
// ITM Maximum number of iterations
// Parameters of function: L1
// Function values: L0
LOCAL I1,I2,I3,I4,IP,J1,J2,D0;
LOCAL fmin,fmax,fdiff;
LOCAL imin,imax,flag,ZAK;
Size(L1)▶D0;
1▶ZAK;
0▶J1;
0▶I4;
MAKEMAT(0,D0,D0+1)▶M1;
MAKELIST(0.,J1,1,D0)▶L2;
MAKELIST(0.,J1,1,D0+1)▶L3;
MAKELIST(0.,J1,1,D0)▶L0;
For I1 FROM 1 TO D0 DO //Simplex
J1+I1▶J1;
1▶J2;
FOR I2 FROM I1+1 TO D0+1 DO
√(J1)/I1/J2*S0▶M1(I1,I2);
I1+1▶J2;
END;
FOR I2 FROM 1 TO D0+1 DO
M1(I1,I2)+L1(I1)▶M1(I1,I2);
END;
END;
0▶flag;
REPEAT // Main loop
If flag=0 THEN 
FOR I1 FROM 1 TO D0+1 DO
FOR I2 FROM 1 TO D0 DO
M1(I2,I1)▶L1(I2);
END;
FUNES()▶L3(I1);
END;
END;
0▶flag;
1ᴇ99▶fmin;
−1ᴇ99▶fmax; 
FOR I1 FROM 1 TO D0+1 DO
IF L3(I1)<fmin Then
I1▶imin;
L3(I1)▶fmin;
END;
IF L3(I1)>fmax Then
I1▶imax;
L3(I1)▶fmax;
END;
END;
FOR I1 FROM 1 TO D0 DO // Reflection center
0▶L2(I1);
FOR I2 FROM 1 TO D0+1 DO
L2(I1)+M1(I1,I2)▶L2(I1);
END;
(L2(I1)-M1(I1,imax))/D0▶L2(I1);
2.*L2(I1)-M1(I1,imax)▶L1(I1); // New point
END;
FUNES()▶Y;
IF Y < fmin THEN // New point is min 
FOR I1 FROM 1 TO D0 DO
3.*L2(I1)-2.*M1(I1,imax)▶L1(I1);
END;
Y▶fmin;
FUNES()▶Y;
IF Y > fmin THEN
FOR I1 FROM 1 TO D0 DO
L1(I1)-L2(I1)+M1(I1,imax)▶L1(I1);
END;
ELSE
Y▶fmin;
END;
END;
IF Y > fmax THEN // New point is max
FOR I1 FROM 1 TO D0 DO
1.5*L3(I1)-0.5*M1(I1,imax)▶L1(I1);
END;
FUNES()▶Y;
If Y > fmax Then
FOR I1 FROM 1 TO D0 DO
FOR I2 FROM 1 TO D0+1 DO
(M1(I1,I2)+M1(I1,imin))*0.5▶M1(I1,I2);
END;
END;
1▶flag; 
END; 
END;
IF Y < fmax THEN // Updating simplex
FOR I1 FROM 1 TO D0 DO
L1(I1)▶M1(I1,imax);
END;
Y▶L3(imax);
END;
0▶ZAK; // Termination criterion
FOR I1 FROM 1 TO D0 DO
FOR I2 FROM 1 TO D0+1 DO
ZAK+(M1(I1,I2)-M1(I1,imin))^2▶ZAK;
END;
END;
√(ZAK/D0)▶ZAK;
I4+1▶I4;
PRINT(I4);
PRINT(fmin);
(I4>=ITM)+(ZAK<E0)+(fmin<E0)▶IP;
UNTIL IP>0;
FOR I1 FROM 1 TO D0 DO // Final work
M1(I1,imin)▶L1(I1);
END;
fmin▶Y;
fmax-fmin▶fdiff;
Print("--- DONE -----------");
Print("Iter: "+I4);
Print("Eps: "+ZAK);
Print("Fmin: "+fmin);
Print("Paramater: L1");
Print("Fct. val: L0");
END;