Advanced Search

Sergey Moiseev

Sergey Moiseev

418 Reputation

14 Badges

18 years, 290 days
Contact Sergey Moiseev
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

Solve by DirectSearch...

Sergey Moiseev 418
April 17 2018
1 2

 

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.

 

> 

restart;with(plots):with(plottools):with(DirectSearch):

# 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

 

> 

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

1

  

 

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:=fsolve(eq,theta=0..2*Pi);
sol:=SolveEquations(eq,[theta=0..2*Pi]);
 

sol_fsolve := 3.612831552

  

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

 (1)

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]:

 

> 

sol:=SolveEquations(eq,[theta=0..2*Pi],
         AllSolutions,solutions=2,pointrange=[theta=0..2*Pi]);

_rtable[18446744074252703134]

 (2)
> 

sol[1,3][1],sol[2,3][1];

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

 (3)

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:
ans:=[]:
for k from op(1, numSols) to op(2, numSols) do
  sol:=SolveEquations( eval( rS[1],beta=k*Pi/(op(2, numSols)/2)),
    [theta=0..2*Pi],AllSolutions,solutions=2,pointrange=[theta=0..2*Pi]);
  if evalf(eval( rS[2],[beta=k*Pi/(op(2, numSols)/2),sol[1,3][1]]))>=0 then
    ans:=[op(ans),rhs(sol[1,3][1])]:
  else
    ans:=[op(ans),rhs(sol[2,3][1])]:
  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]]))],
                symbol=solidcircle,symbolsize=40,color=blue),
              rotate( p1, k*Pi/(op(2, numSols)/2))],scaling=constrained),
          k=numSols):
  display(plt, insequence=true);

 

The animation shows that all works correctly.

> 

 

> 

 


 

Download Answer_Archimedian.mw

SolveEquations from DirectSearch...

Sergey Moiseev 418
November 22 2013
1 2

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).

with(DirectSearch):
m:=5;v1:=4969;v2:=208;
eq:=m = floor(v1/x)+floor(v2/x);
SolveEquations(eq);

             [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);

 

try DirectSearch...

Sergey Moiseev 418
November 21 2013
2 0

The SolveEquations command from DirectSearch package solves these equations:

with(DirectSearch);

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]

SolveEquations(eq2,assume=integer);

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

DataFit from DirectSearch package...

Sergey Moiseev 418
November 14 2013
1 4

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]

f:=0.249966755570439e-1+7.185110562*10^(-7)*abs(x-12.2879102948195)^1.68823162809910;

The fit is good enough.

DataFit from DirectSearch package...

Sergey Moiseev 418
October 24 2013
2 1

You can use DataFit command from DirectSearch package:

with(DirectSearch):with(plots):

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]]);

X:=xy[..,1];
Y:=xy[..,2];

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

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

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

with(plots):
p1,p2:=pointplot(X,Y),plot(f1,x=0..2);
plots[display](p1,p2);

solve...

Sergey Moiseev 418
October 10 2013
1 0

eq:=x=(2012-x)^(1/2)*(2013-x)^(1/2)+(2013-x)^(1/2)*(2014-x)^(1/2)+(2012-x)^(1/2)*(2014-x)^(1/2):

solve(eq);

use fsolve...

Sergey Moiseev 418
June 11 2013
0 0

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});

solution by DataFit from DirectSearch pa...

Sergey Moiseev 418
May 24 2013
1 3

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);
ysol:=eval(y,sol[2]);

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

by DirectSearch package...

Sergey Moiseev 418
May 13 2013
0 0

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.

with(DirectSearch):
eq:=sin(Pi*x^2/2);


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)]);

correct solution...

Sergey Moiseev 418
May 08 2013
1 1

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);

solution by DirectSearch package...

Sergey Moiseev 418
April 10 2013
1 0

with(DirectSearch):
SolveEquations([eqn1, eqn2, eqn3]);

obtain each solution point from DirectSe...

Sergey Moiseev 418
March 27 2013
1 0

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

another way...

Sergey Moiseev 418
March 04 2013
2 0

with(Statistics):
Z:=RandomVariable(Normal(0,1));
                                     _R1
CDF(Z, 0.87)-CDF(Z, 0);
                                0.3078497979
1-CDF(abs(Z), 1.39);
Probability(abs(Z)>1.39);
                                0.1645288774
                                0.1645288774

too tight constraints...

Sergey Moiseev 418
October 26 2012
1 4

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,

with(DirectSearch);
f:=(3*a/sqrt(1-u^2)+b/sqrt(1-t^2))/c;
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.


shorter...

Sergey Moiseev 418
October 04 2012
1 0 
with(Statistics):

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