Sergey Moiseev

Sergey Moiseev

418 Reputation

14 Badges

17 years, 101 days
Sergey N. Moiseev received M.S., Ph.D. and Dr.Sc. degrees in radio physics from the Voronezh State University, Voronezh, Russia in 1986, 1993 and 2003, respectively. From 1984 to 2003 the topics of his research have included theory and methods of signal processing, nonlinear optimization, decision-making theory, time series prediction, statistical radio physics, ionosphere sporadic channel models. He is currently a principal scientist in the JSC Kodofon, Voronezh, Russia. His current research interests are wide spread in the area of the communications.

MaplePrimes Activity

These are answers submitted by Sergey Moiseev


I analysed the problem. The DirectSearch works correctly.

The SolveEquations command finds solutions by minimizing squared residuals.

The problem is that squared residuals have many minimums.

The plot below shows it.



# Generate the cardioid, and set up a rotation

  spX:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   cos(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <= 2*Pi, sin(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
  spY:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   sin(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <  2*Pi, cos(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
  p1:=plot( [ spX, spY, theta=0..2*Pi ],thickness=3,color=red,discont=true):
 # display( [ seq( rotate(p1,alpha),alpha=0..8*evalf(Pi), evalf(Pi/8))],insequence=true);

# Define a rotation matrix which allows the above
# plot to be rotated by any specified angle beta,
# and produce the "rotated" x and y-coordinates

  mRot:=Matrix(2,2, [ [ cos(beta), -sin(beta)],
                      [ sin(beta),  cos(beta)]  ]):
  rS:= mRot.< spX, spY >:

# Set the number of solutions to be used for all calculations.
  numSols:= 0..5:

Now we plots the squared residuals for one solution


eq:=eval( rS[1],beta=N*Pi/(op(2, numSols)/2)):
plot(eq_squared,theta=0..2*Pi,title="Squared residuals");




The plot shows that there are four minimums; two are exact and two are inconsistent (not exact).

Because the SolveEquations can solve inconsistent, overdetermined and underdetermined system of equations

the any of four solutions can be returned.

Solve this equation by fsolve and SolveEquations commands.


sol_fsolve := 3.612831552


[HFloat(0.45205120448170044), Vector[column](%id = 18446744074252680838), [theta = HFloat(5.883009651016983e-9)], 12]


The fsolve returns one exact solution and SolveEquations returns one inconsistent (not exact).

The large squared residual (first item of sol) shows it. But it is correct result for SolveEquations.

How can we make the SolveEquations returns all the exact solutions?

Assuming only two exact solutions we can use options AllSolutions, solutions=2, pointrange=[0..2*Pi]:






theta = HFloat(3.6128315516259537), theta = HFloat(2.6703537555499497)


The small squared residuals (first items of solutions) show that obtained solutions very likely are exact.



Now we can solve the problem by selecting one of two exact solution.

We select solution for which y-axis intercept is positive.

numSols:= 0..15:
for k from op(1, numSols) to op(2, numSols) do
  sol:=SolveEquations( eval( rS[1],beta=k*Pi/(op(2, numSols)/2)),
  if evalf(eval( rS[2],[beta=k*Pi/(op(2, numSols)/2),sol[1,3][1]]))>=0 then
  end if;
end do:
# ans;

  plt:= seq( display( [ pointplot( [ 0,evalf( eval
                    ( rS[2],[ beta=k*Pi/(op(2, numSols)/2),theta=ans[k+1]]))],
              rotate( p1, k*Pi/(op(2, numSols)/2))],scaling=constrained),
  display(plt, insequence=true);


The animation shows that all works correctly.





I do not detect any problem to solve this equation by SolveEquations command from DirectSearch package. Note that true solution lies within interval (828.1666667, 993.8000000] (see Kitonum answer).

eq:=m = floor(v1/x)+floor(v2/x);

             [0., [0.], [x = 950.984332050577], 20]

When the function is piecewise-flat it is more reliable to take the big initial step and global search strategy:

SolveEquations(eq, step=1000, strategy=globalsearch);

            [0., [0.], [x = 987.528685750666], 225]

The following command obtains approximately 100 correct solutions:

sol:=SolveEquations(eq, step=1000, AllSolutions);

The following command obtains more than 900 correct solutions:

sol:=SolveEquations(eq, step=1000, number=1000, AllSolutions);


The SolveEquations command from DirectSearch package solves these equations:


eq1:=x^2 + y^2 +z^2 =3;
eq2:=x^2 +y^2 +2*z^2 +x*y -y*z=3;

SolveEquations(eq1, assume=integer);

                            [0., [0.], [x = 1, y = 1, z = 1], 132]


                            [0., [0.], [x = 1, y = 0, z = -1], 154]

The DataFit command from the DirectSearch package finds the following solution for the model f

f := a+b*abs(x-d)^c;

subject to constraint  c>=1:

[a = 0.249966755570439e-1, b = 7.18511056172813*10^(-7), c = 1.68823162809910, d = 12.2879102948195]


The fit is good enough.

You can use DataFit command from DirectSearch package:


xy:=Matrix([[0.2e-1, .158], [0.2e-1, .159], [0.3e-1, .161], [0.3e-1, .164], [0.3e-1, .166], [0.4e-1, .169], [0.6e-1, .173], [0.8e-1, .178], [.1, .185], [.11, .187], [.14, .193], [.19, .2], [.28, .21], [.38, .223], [.44, .233], [.58, .244], [.82, .256], [1.4, .278], [1.71, .281], [1.78, .282], [1.78, .282], [1.81, .282]]);


f:=a+b*abs(x)^c;     # model function

sol:=DataFit(f, X, Y, x);

f1:=eval(f,sol[2]);     # model function with estimated parameters




Use fsolve() instead of solve():

fsolve({eq_1, eq_10, eq_11, eq_12, eq_2, eq_3, eq_4, eq_5, eq_6, eq_7, eq_8, eq_9}, {c1, c2, c3, c4, t1, t2, t4, t5, t6, t7, t8, vc});

y := A*x+B*(sum(exp(-n^2*(x-C))/n^2, n = 1 .. 10));
X := Vector([1, 2, 3, 4, 5], datatype = float);
Y := Vector([2, 2, 6, 6, 8], datatype = float);

sol:=DataFit(y, X, Y, x);

sol:=[0.64, [A=1.6, B=0.9, C=-62.1], 73]

By default the multimodal optimizer from DirectSearch package uses 100 random initial points in the range. It is too small for large number of roots. Therefore, we can increase number of initial points (by option number) in order to more reliably find all roots.


sol:=SolveEquations(eq, [x=-5..5], pointrange=[x=-5..5], AllSolutions, number=200);
sort([seq(rhs(sol[i][3][1]), i=1..25)]);

But randomness of initial points decrease reliability. It is more reliable to use regular set of initial points in range.

startpoints:=[seq([x=-5.+i/10.], i=1..100)];
sol:=SolveEquations(eq, [x=-5..5], AllSolutions, initialpoints=startpoints);
sort([seq(rhs(sol[i][3][1]), i=1..25)]);

The command

solve(w(x,y), F);

really means

solve(w(x,y)=0, F);

Therefore, Maple gave the correct solution for incorrect task. The correct task:

eq:=w(x,y)=4*F/Pi^4/E[2]/I/a/b*Sum(Sum(sin(m*Pi*xi/a)*sin(n*Pi*eta/b)*sin(m*Pi*x/a)*sin(n*Pi*y/b)/((m/a)^2+(n/b)^2)^2, n= 1..infinity), m= 1..infinity);
solve(eq, F);

SolveEquations([eqn1, eqn2, eqn3]);

sol:=DirectSearch[SolveEquations](sys, assume = positive);
Y:=rhs(sol[3][1]);     # first solution x1a
Y2:=rhs(sol[3][2]);   # second solution x1c, etc...

CDF(Z, 0.87)-CDF(Z, 0);
1-CDF(abs(Z), 1.39);

The DirectSearch cannot find feasible point. It seems that constraints are too tight. For such tight constraints DirectSearch has two option: penaltymethod=truefalse and feasibilitytolerance=nonnegative. The first option convert inequality constraints to part of objective function. This method allow inequality constraints to be violated. Try, for example,

constr:={a>=0,b>=0,c>=0,u>=0,t>=0,a^2+2*b*c*u >= b^2+c^2, a*t+b*u <= c, b^2*(t^2-u^2)/(t^2-1)+c^2 <= 2*b*c*u};
sol:=GlobalOptima(f, constr, penaltymethod);
eval(constr, sol[2]);

The option feasibilitytolerance=nonnegative sets the maximum absolute allowable inequality constraint violation. This option is more appropriate for your task. Try

sol:=GlobalOptima(f, constr, feasibilitytolerance=10^(-10), tolerances=10^(-14));
eval(constr, sol[2]);

I obtain the following feasible solution:

sol := [4.25363469533755, [a = 0.10665115417e-1, b = 0.439302870444997e-2, c = 0.100583437022062e-1, t = .911765508073957, u = 0.76088604888748e-1], 3302];

It seems that the solution can be improved.


CDF(Normal(64.3, 2.6), 56.5); # cumulative distribution function
1 2 3 4 Page 1 of 4