Sergey Moiseev

## 418 Reputation

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.

## Solve by DirectSearch...

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,beta=N*Pi/(op(2, numSols)/2)): eq_squared:=eq^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:=fsolve(eq,theta=0..2*Pi); sol:=SolveEquations(eq,[theta=0..2*Pi]);  (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]); (2)
 > sol[1,3],sol[2,3]; (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,beta=k*Pi/(op(2, numSols)/2)),     [theta=0..2*Pi],AllSolutions,solutions=2,pointrange=[theta=0..2*Pi]);   if evalf(eval( rS,[beta=k*Pi/(op(2, numSols)/2),sol[1,3]]))>=0 then     ans:=[op(ans),rhs(sol[1,3])]:   else     ans:=[op(ans),rhs(sol[2,3])]:   end if; end do: # ans;
 > plt:= seq( display( [ pointplot( [ 0,evalf( eval                     ( rS,[ 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.

 >
 >

## SolveEquations from DirectSearch...

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

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

## try DirectSearch...

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

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

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);     # model function with estimated parameters

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

## solve...

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

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

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

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

## by DirectSearch package...

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]), 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]), i=1..25)]);

## correct solution...

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

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

## obtain each solution point from DirectSe...

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

## another way...

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

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

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

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

```with(Statistics):

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