## 12515 Reputation

8 years, 209 days

## Solution...

 > restart

Minimize   y   for integers  x,y  satisfying:
0 < x < y,  154*x+69*y =0 mod  7^3 , 13*x+716*y =0 mod 13^3  , 23*x+3059*y =0 mod  23^3.

 >
 > restart;
 > eqs := 154*x+69*y, 13*x+716*y, 23*x+3059*y; m := 7^3, 13^3, 23^3;
 (1)

For an "ideal Maple" it should work:

 > # Optimization:-Minimize(y, {1<=x, x<=y, eqs =~ m*~(u,v,w)}, assume = integer, depthlimit=100000);
 > # But the problem it too large!
 > #################################################################################################
 > Eqs := 154*x+69*y, 13*x+716*y, x+133*y; M := 7^3, 13^3, 23^2;
 (2)
 > sys:=y =-154 * 1/69*x mod M[1] + M[1]*u,  y = -13/716*x mod M[2] + M[2]*v, y = -1/133*x mod M[3] + M[3]*w;
 (3)
 > isolve(sys[1]-sys[2]);
 (4)
 > eval(%, [_Z1=a,_Z2=b]);
 (5)
 > e1:=eliminate(%, {u,v});
 (6)
 > isolve(sys[1]-sys[3]);
 (7)
 > eval(%, [_Z1=c,_Z2=d]);
 (8)
 > e2:=eliminate(%, {u,w});
 (9)
 > ee:=e1[1] union e2[1], e1[2] union e2[2];
 (10)
 > EE:=eval(ee, [a=u, c=u]);
 (11)
 > X, X1 := solve~(EE[2],x)[];
 (12)
 > s:=eval(isolve(X=X1),[_Z1=p,_Z2=q,_Z3=r]);
 (13)
 > X0:=eval(X, s);
 (14)
 > Y0:=rhs(eval(eval(sys[1], s), x=X0));
 (15)
 > Optimization:-Minimize(Y0, { X0<=Y0, 1<=X0}, assume=integer, depthlimit=1000 );
 (16)

Unfortunately Minimize uses hardware integers and it results an overflow!
(This is actually a bug!)

Let's minimize the reduced X0, Y0 problem by brute force:

 > p:='p':q:='q':
 > (Xmin,Ymin,pmin,qmin):=(infinity\$4): for p from -1000 to 1000 do for q from -1000 to 1000 do   if X0<1 or X0 >= Y0 then next fi;   if Y0
 (17)
 >

## ExtremePoints for piecewise...

ExtremePoints is not reliable for piecewice defined functions.

 > restart;
 > EP:=Student:-Calculus1:-ExtremePoints:
 > piecewise(-1 <= x and x <= 2, x^2, undefined);
 (1)
 > EP(%); #OK
 (2)
 > piecewise(-1 < x and x < 2, x^2, undefined);
 (3)
 > EP(%); #OK
 (4)
 > piecewise(-1 <= x and x <= 2, x, undefined);
 (5)
 > EP(%); # expected [1,2]
 (6)
 > piecewise(x<1, undefined, x<=2, x*sin(1/x), undefined);
 (7)
 > EP(%, 1..2); # expected [1,2]
 (8)
 > EP(x*sin(1/x), 1..2); # OK
 (9)
 > piecewise(x<0, undefined, x=0, 0,  x<=1, x*sin(1/x), undefined);
 (10)
 > EP(%); # expected [1]
 (11)
 > piecewise(x<-1, undefined, x<0, x*sin(1/x), x=0, 0, undefined);
 (12)
 > EP(%); # expected [-1]
 (13)
 >

## Bug...

It is a bug. For solve too.

Both commands work if V is given as a set, or is omitted.
So, use:
PolynomialSystem(F, {V[]});  or  PolynomialSystem(F, {V[]}, explicit);
or PolynomialSystem(F);

 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 (4)
 > solutions := pdsolve({eq1, eq10, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9}, dependent_variables);
 (5)
 >
 (6)

```A219954list := nmax ->
[seq[scan=`+`](ifelse(n=1, 0, 3^numboccur(convert(n-1, 'base', 2), 1) - ifelse(2^ilog2(n)=n, n/2, 0)), n = 1..nmax)]:
```

A219954list(33);
[0, 2, 5, 12, 15, 24, 33, 56, 59, 68, 77, 104, 113, 140, 167,  240, 243, 252,
261, 288, 297, 324, 351, 432, 441, 468, 495,   576, 603, 684, 765, 992, 995]

## Do not use floats...

 >
 >
 (1)
 (2)
 > u_exact(t,X,Y);
 (3)
 >
 (4)
 > eval(PDE, u=u_exact)
 (5)
 > simplify(%);
 (6)
 >

## global...

Technically, you must use:

`U,S,Vt := SingularValues(A, output=[':-U',':-S',':-Vt']);`

(Of course, in your example, this does not make much sense because U, S, Vt are not needed/used.)

## A solution...

```restart;
a:=floor(33*sqrt(3)):
L2:=R2*16/6: # the edge^2
nr:=0: W:=Array(1..0):
for  x from 0 to a do
for  y from x to a do
for  z from y to a do
if x^2+y^2+z^2=R2 then nr++: W(nr):=[x,y,z] fi
od od od:
f:=u -> ([u[1],u[2],u[3]], [-u[1],u[2],u[3]],[u[1],-u[2],u[3]], [u[1],u[2],-u[3]],
[-u[1],-u[2],u[3]], [-u[1],u[2],-u[3]],[u[1],-u[2],-u[3]], [-u[1],-u[2],-u[3]]):
P:=f~([seq](W)): # the integer points on the sphere
n:=nops(P);
#                            n := 96

d2 :=(u,v) -> inner(u-v,u-v):
A:=Matrix( n, (i,j)-> `if`(d2(P[i],P[j])=L2, 1, 0) ):
with(GraphTheory):
G:=Graph(A):
K4 := CompleteGraph(4):
IsSubgraphIsomorphic(K4,G, isomorphism);
#             true, {1 = 95, 2 = 94, 3 = 93, 4 = 89}

T := P[[rhs~(%[2])[]]]; # A desired tetrahedron
# T := [ [33, 33, 33], [-33, -33, 33], [-33, 33, -33], [33, -33, -33]  ]

plots:-display(plottools:-tetrahedron(T), plottools:-sphere([0, 0, 0], sqrt(R2), transparency=0.8), color="Blue");
```

## DirectSearch...

The floor method cannot work because NLPSolve finds only local extrema and almost any point is a local max,
actually any point in ( (-5, oo) \ Z )^9.

The DirectSearch package (in Application Center) works, but the solution is of course not guaranteed.

```GlobalOptima(add(x[k],k=1..9),
[seq(x[k]::integer, k=1..9),seq(x[k]<=x[k+1], k=1..8), x[1]>=-5, x[9]<=50, add(x[k]^3,k=1..9)=0],
evaluationlimit=50000, maximize);
```

[12., [x[1] = -4, x[2] = 2, x[3] = 2, x[4] = 2, x[5] = 2, x[6] = 2, x[7] = 2, x[8] = 2, x[9] = 2], 2169]

## Hermite...

No, such polynomial would have been found by solve.
(You have a Hermite interpolation problem here!)

## u,v as exp...

```ex:=F(u^(1/n), v^(1/n))^n - F(u, v):
ex1:=eval(ex, [u=exp(-A), v=exp(-B)]):
expand(simplify(ex1)) assuming positive;  # 0
```

## Just use ReducedRowEchelon...

Just use ReducedRowEchelonForm and append the identity matrix to your square matrix.

 > restart;
 > with(LinearAlgebra):
 > n:=3:
 > A:=RandomMatrix(n);
 (1)
 > LinearAlgebra:-ReducedRowEchelonForm();
 (2)
 > B:=%[..,n+1..];
 (3)
 > A.B; #check

 (4)
 >

## You just need a pencil for this.c=5, a=1...

You just need a pencil for this.

c=5, a=12/2=6, b^2 = a^2 - c^2 = 11

So, x^2/36 + y^2/11 = 1.

## I see the QuantifierElimi...

I see the QuantifierElimination package (new in Maple 2023) as experimental.
Here is a much simpler problem and the solution obtained by  QuantifierElimination.

 > restart;
 > with( QuantifierElimination );
 (1)
 > f:=(x-1)/(x^2+x+1);
 (2)
 > fmax:=evala(maximize(f));
 (3)
 > forall(x, f <= K);
 (4)
 > sol:=QuantifierEliminate( % ); ### ?
 (5)
 > should_be = (K >= fmax);   # !!
 (6)
 > is(eval(sol, K=fmax));            # ckeck
 > is(eval(sol, K=fmax-10^(-20)));   #
 (7)

N.B. Why do you use such an unusual UserProfile?  The output of   2*x + 30  is:

## 3 half-lines...

The system is not very simple, but I agree that SolveTools:-SemiAlgebraic is very slow in general, and I had not the patience for a result.

Anyway, the result (for eq1) seems to be (and is confirmed by MMA):

```a := RootOf(_Z^3 - 4*_Z^2 - 27, 5.054 .. 5.063):
sol := {z=0, y = x, a < x}, {y=0, z = x, a < x}, {x=0, y = z, a < z};
```

plots:-spacecurve([[t,t,0], [t,0,t], [0,t,t]], t=a..3*a, color=[red,blue,green], thickness=5);

 1 2 3 4 5 6 7 Last Page 3 of 112
﻿