vv

14242 Reputation

20 Badges

10 years, 345 days

MaplePrimes Activity


These are answers submitted by vv

Export the document as Maple code (.mpl).
Then open the .mpl. In Maple 2023, it will appear in a new Code Edit window.
Edit the code, i.e. replace `–` with -

There are not other errors (see the Diagnostics tab), so you can copy & paste the code in the worksheet and execute.

P.S. I use 2D input strictly for presentation purposes (and almost never the Document mode).

ex:=x*y*z^2*w^3:
degree(ex);
 #                              7
factors(ex)[2];
 #               [[y, 1], [z, 2], [w, 3], [x, 1]]

Miny:=proc(eqs::list, m::list(posint), R:=1000)
description "minimize y, for integers  x,y  satisfying 0<x<y, eqs[1] = 0 mod m[1], ...";
local p,q,X,Y,u,k, ir, Xmin,Ymin,pmin,qmin;
if nops(eqs)<>nops(m) then error "Lists must have same size" fi;
if nops(indets(eqs)) <> 2 then error "Must have 2 variables" fi; 
eval( isolve({(eqs =~ m*~[seq(u[k],k=1..nops(m))])[]}), [_Z1=p,_Z2=q] ):
(X,Y):= (eval([indets(eqs)[]],%)[]);
ir:=IntegerRelations:-LLL([[coeff(X,p),coeff(Y,p)], [coeff(X,q),coeff(Y,q)]]):
(X,Y):= (ir[1,1]*p+ir[2,1]*q, ir[1,2]*p+ir[2,2]*q);
(Xmin,Ymin,pmin,qmin):=(infinity$4):
for p from -R to R do
for q from -R to R do
  if eval(X<1 or X >= Y) then next fi;
  if eval(Y<Ymin) then (Xmin,Ymin,pmin,qmin):=eval([X,Y,p,q])[] fi
od od:
['X','Y','p','q'] =~ eval([Xmin,Ymin,pmin,qmin]);
end proc:

eqs := [154*x+69*y, 13*x+716*y, 23*x+3059*y]:
m := [7^3, 13^3, 23^3]:
Miny(eqs,m);

              [X = 15893, Y = 18837, p = 0, q = 1]

eqs := [154*x+69*y, 13*x+716*y, 23*x+3059*y, 2295*x + 4522*y, 6479*x + 5396*y]:
m := [7^3, 13^3, 23^3, 17^3,  29^3]:
Miny(eqs,m);

        [X = 186957261, Y = 190687133, p = -12, q = -1]

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;

154*x+69*y, 13*x+716*y, 23*x+3059*y

 

343, 2197, 12167

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

154*x+69*y, 13*x+716*y, x+133*y

 

343, 2197, 529

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

y = 301*x+343*u, y = 494*x+2197*v, y = 175*x+529*w

(3)

isolve(sys[1]-sys[2]);

{u = _Z1, v = 122*_Z1+193*_Z2, x = -1387*_Z1-2197*_Z2}

(4)

eval(%, [_Z1=a,_Z2=b]);

{u = a, v = 122*a+193*b, x = -1387*a-2197*b}

(5)

e1:=eliminate(%, {u,v});

[{u = a, v = 122*a+193*b}, {x+1387*a+2197*b}]

(6)

isolve(sys[1]-sys[3]);

{u = _Z1, w = 49*_Z1+126*_Z2, x = 203*_Z1+529*_Z2}

(7)

eval(%, [_Z1=c,_Z2=d]);

{u = c, w = 49*c+126*d, x = 203*c+529*d}

(8)

e2:=eliminate(%, {u,w});

[{u = c, w = 49*c+126*d}, {-x+203*c+529*d}]

(9)

ee:=e1[1] union e2[1], e1[2] union e2[2];

{u = a, u = c, v = 122*a+193*b, w = 49*c+126*d}, {-x+203*c+529*d, x+1387*a+2197*b}

(10)

EE:=eval(ee, [a=u, c=u]);

{u = u, v = 122*u+193*b, w = 49*u+126*d}, {-x+203*u+529*d, x+1387*u+2197*b}

(11)

X, X1 := solve~(EE[2],x)[];

-1387*u-2197*b, 203*u+529*d

(12)

s:=eval(isolve(X=X1),[_Z1=p,_Z2=q,_Z3=r]);

{b = p, d = -1513*p-1590*q, u = 502*p+529*q}

(13)

X0:=eval(X, s);

-698471*p-733723*q

(14)

Y0:=rhs(eval(eval(sys[1], s), x=X0));

-210067585*p-220669176*q

(15)

Optimization:-Minimize(Y0, { X0<=Y0, 1<=X0}, assume=integer, depthlimit=1000 );

[0, [p = 0, q = 0]]

(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<Ymin then (Xmin,Ymin,pmin,qmin):=(X0,Y0,p,q) fi
od od:
(x,y,'p','q') =~ (Xmin,Ymin,pmin,qmin);

x = 1103, y = 26390, p = 562, q = -535

(17)

 


 

Download diophant.mw

 

ExtremePoints is not reliable for piecewice defined functions.

 

restart;

EP:=Student:-Calculus1:-ExtremePoints:

piecewise(-1 <= x and x <= 2, x^2, undefined);

piecewise(-1 <= x and x <= 2, x^2, undefined)

(1)

EP(%); #OK

[-1, 0, 2]

(2)

piecewise(-1 < x and x < 2, x^2, undefined);

piecewise(-1 < x and x < 2, x^2, undefined)

(3)

EP(%); #OK

[0]

(4)

piecewise(-1 <= x and x <= 2, x, undefined);

piecewise(-1 <= x and x <= 2, x, undefined)

(5)

EP(%); # expected [1,2]

[2]

(6)

piecewise(x<1, undefined, x<=2, x*sin(1/x), undefined);

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

[1, 2]

(9)

piecewise(x<0, undefined, x=0, 0,  x<=1, x*sin(1/x), undefined);

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

piecewise(x < -1, undefined, x < 0, x*sin(1/x), x = 0, 0, undefined)

(12)

EP(%); # expected [-1]

[]

(13)

 

 

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

restart

variables := [t, x, y, z]; MyPsi := proc (t, x, y, z) options operator, arrow; psi(t, x, y, z) end proc; zeta_0 := proc (t, x, y, z) options operator, arrow; Zeta_0(t, x, y, z) end proc; zeta_1 := proc (t, x, y, z) options operator, arrow; Zeta_1(t, x, y, z) end proc; zeta_2 := proc (t, x, y, z) options operator, arrow; Zeta_2(t, x, y, z) end proc; zeta_3 := proc (t, x, y, z) options operator, arrow; Zeta_3(t, x, y, z) end proc

[t, x, y, z]

 

proc (t, x, y, z) options operator, arrow; psi(t, x, y, z) end proc

 

proc (t, x, y, z) options operator, arrow; Zeta_0(t, x, y, z) end proc

 

proc (t, x, y, z) options operator, arrow; Zeta_1(t, x, y, z) end proc

 

proc (t, x, y, z) options operator, arrow; Zeta_2(t, x, y, z) end proc

 

proc (t, x, y, z) options operator, arrow; Zeta_3(t, x, y, z) end proc

(1)

dzeta_0_dt := diff(zeta_0(t, x, y, z), t); dzeta_1_dt := diff(zeta_1(t, x, y, z), t); dzeta_2_dt := diff(zeta_2(t, x, y, z), t); dzeta_3_dt := diff(zeta_3(t, x, y, z), t); dzeta_0_dx := diff(zeta_0(t, x, y, z), x); dzeta_1_dx := diff(zeta_1(t, x, y, z), x); dzeta_2_dx := diff(zeta_2(t, x, y, z), x); dzeta_3_dx := diff(zeta_3(t, x, y, z), x); dzeta_0_dy := diff(zeta_0(t, x, y, z), y); dzeta_1_dy := diff(zeta_1(t, x, y, z), y); dzeta_2_dy := diff(zeta_2(t, x, y, z), y); dzeta_3_dy := diff(zeta_3(t, x, y, z), y); dzeta_0_dz := diff(zeta_0(t, x, y, z), z); dzeta_1_dz := diff(zeta_1(t, x, y, z), z); dzeta_2_dz := diff(zeta_2(t, x, y, z), z); dzeta_3_dz := diff(zeta_3(t, x, y, z), z)

diff(Zeta_0(t, x, y, z), t)

 

diff(Zeta_1(t, x, y, z), t)

 

diff(Zeta_2(t, x, y, z), t)

 

diff(Zeta_3(t, x, y, z), t)

 

diff(Zeta_0(t, x, y, z), x)

 

diff(Zeta_1(t, x, y, z), x)

 

diff(Zeta_2(t, x, y, z), x)

 

diff(Zeta_3(t, x, y, z), x)

 

diff(Zeta_0(t, x, y, z), y)

 

diff(Zeta_1(t, x, y, z), y)

 

diff(Zeta_2(t, x, y, z), y)

 

diff(Zeta_3(t, x, y, z), y)

 

diff(Zeta_0(t, x, y, z), z)

 

diff(Zeta_1(t, x, y, z), z)

 

diff(Zeta_2(t, x, y, z), z)

 

diff(Zeta_3(t, x, y, z), z)

(2)

eq1 := dzeta_0_dt = MyPsi(t, x, y, z); eq2 := -A^2*dzeta_1_dt+dzeta_0_dx = 0; eq3 := -B^2*dzeta_0_dy+dzeta_0_dy = 0; eq4 := -C^2*dzeta_3_dt+dzeta_0_dz = 0; eq5 := dzeta_1_dx+(diff(A, t))*zeta_0/A = MyPsi(t, x, y, z); eq6 := A^2*dzeta_1_dy+B^2*dzeta_2_dx = 0; eq7 := A^2*dzeta_1_dz+C^2*dzeta_3_dx = 0; eq8 := dzeta_2_dy+(diff(B, t))*zeta_0/B = MyPsi; eq9 := B^2*dzeta_2_dz+C^2*dzeta_3_dy = 0; eq10 := dzeta_3_dz+(diff(C, t))*zeta_0/C = MyPsi(t, x, y, z)

diff(Zeta_0(t, x, y, z), t) = psi(t, x, y, z)

 

-A^2*(diff(Zeta_1(t, x, y, z), t))+diff(Zeta_0(t, x, y, z), x) = 0

 

-B^2*(diff(Zeta_0(t, x, y, z), y))+diff(Zeta_0(t, x, y, z), y) = 0

 

-C^2*(diff(Zeta_3(t, x, y, z), t))+diff(Zeta_0(t, x, y, z), z) = 0

 

diff(Zeta_1(t, x, y, z), x) = psi(t, x, y, z)

 

A^2*(diff(Zeta_1(t, x, y, z), y))+B^2*(diff(Zeta_2(t, x, y, z), x)) = 0

 

A^2*(diff(Zeta_1(t, x, y, z), z))+C^2*(diff(Zeta_3(t, x, y, z), x)) = 0

 

diff(Zeta_2(t, x, y, z), y) = MyPsi

 

B^2*(diff(Zeta_2(t, x, y, z), z))+C^2*(diff(Zeta_3(t, x, y, z), y)) = 0

 

diff(Zeta_3(t, x, y, z), z) = psi(t, x, y, z)

(3)

dependent_variables := [Zeta_0(t, x, y, z), Zeta_1(t, x, y, z), Zeta_2(t, x, y, z), Zeta_3(t, x, y, z), A(t), B(t), C(t), psi(t, x, y, z)]; independent_variables := [t, x, y, z]

[Zeta_0(t, x, y, z), Zeta_1(t, x, y, z), Zeta_2(t, x, y, z), Zeta_3(t, x, y, z), A(t), B(t), C(t), psi(t, x, y, z)]

 

[t, x, y, z]

(4)

solutions := pdsolve({eq1, eq10, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9}, dependent_variables);

{Zeta_0(t, x, y, z) = (1/2)*_C1*t^2+(1/2)*(2*_C2*x+2*_C3*z+2*_C4)*t+(1/2)*(A^2*x^2+C^2*z^2)*_C1+A^2*_C10*x+C^2*_C6*z+_C12, Zeta_1(t, x, y, z) = (1/2)*((_C2*x^2+(2*_C1*t+2*_C3*z+2*_C4)*x+2*_C10*t+2*_C11)*A^2+(-C^2*z^2+t^2)*_C2-2*C^2*_C7*z-2*B^2*_C9*y)/A^2, Zeta_2(t, x, y, z) = MyPsi*y-C^2*_C5*z/B^2+_C9*x+_F2(t), Zeta_3(t, x, y, z) = (1/2)*((_C3*z^2+(2*_C1*t+2*_C2*x+2*_C4)*z+2*_C5*y+2*_C6*t+2*_C7*x+2*_C8)*C^2+(-A^2*x^2+t^2)*_C3)/C^2, psi(t, x, y, z) = _C1*t+_C2*x+_C3*z+_C4}

(5)

for i to nops(solutions) do print(solutions[i]) end do

Zeta_0(t, x, y, z) = (1/2)*_C1*t^2+(1/2)*(2*_C2*x+2*_C3*z+2*_C4)*t+(1/2)*(A^2*x^2+C^2*z^2)*_C1+A^2*_C10*x+C^2*_C6*z+_C12

 

Zeta_1(t, x, y, z) = (1/2)*((_C2*x^2+(2*_C1*t+2*_C3*z+2*_C4)*x+2*_C10*t+2*_C11)*A^2+(-C^2*z^2+t^2)*_C2-2*C^2*_C7*z-2*B^2*_C9*y)/A^2

 

Zeta_2(t, x, y, z) = MyPsi*y-C^2*_C5*z/B^2+_C9*x+_F2(t)

 

Zeta_3(t, x, y, z) = (1/2)*((_C3*z^2+(2*_C1*t+2*_C2*x+2*_C4)*z+2*_C5*y+2*_C6*t+2*_C7*x+2*_C8)*C^2+(-A^2*x^2+t^2)*_C3)/C^2

 

psi(t, x, y, z) = _C1*t+_C2*x+_C3*z+_C4

(6)

NULL

Download code_solve_system_of_equations_vv.mw

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]

NULL

restart

v := 203/100; x0 := 0; y0 := 0; mu := 1; A_coef := (1/2)*(3*(v-2)); Coef_inside_sqrt := (5*(1/10))*sqrt((v-2)/(2*mu*v)); u_exact := proc (t, X, Y) options operator, arrow; A_coef*sech(Coef_inside_sqrt*(X+Y-v*t-x0-y0))^2 end proc

203/100

 

0

 

0

 

1

 

9/200

 

(1/812)*1218^(1/2)

 

proc (t, X, Y) options operator, arrow; A_coef*sech(Coef_inside_sqrt*(X+Y-v*t-x0-y0))^2 end proc

(1)

1

 

0.45000000e-1

 

0.4298011912e-1

 

proc (t, X, Y) options operator, arrow; A_coef*sech(Coef_inside_sqrt*(X+Y-v*t-x0-y0))^2 end proc

(2)

u_exact(t,X,Y);

(9/200)*sech((1/812)*1218^(1/2)*(X+Y-(203/100)*t))^2

(3)

PDE := diff(u(t, x, y), t)+diff(u(t, x, y), x)+diff(u(t, x, y), y)+u(t, x, y)*(diff(u(t, x, y), x))+u(t, x, y)*(diff(u(t, x, y), y))-mu*(diff(u(t, x, y), x, x, t))-mu*(diff(u(t, x, y), y, y, t))

diff(u(t, x, y), t)+diff(u(t, x, y), x)+diff(u(t, x, y), y)+u(t, x, y)*(diff(u(t, x, y), x))+u(t, x, y)*(diff(u(t, x, y), y))-(diff(diff(diff(u(t, x, y), t), x), x))-(diff(diff(diff(u(t, x, y), t), y), y))

(4)

eval(PDE, u=u_exact)

(27/8120000)*sech((1/812)*1218^(1/2)*(x+y-(203/100)*t))^2*1218^(1/2)*tanh((1/812)*1218^(1/2)*(x+y-(203/100)*t))-(81/8120000)*sech((1/812)*1218^(1/2)*(x+y-(203/100)*t))^4*1218^(1/2)*tanh((1/812)*1218^(1/2)*(x+y-(203/100)*t))-(27/8120000)*sech((1/812)*1218^(1/2)*(x+y-(203/100)*t))^2*tanh((1/812)*1218^(1/2)*(x+y-(203/100)*t))^3*1218^(1/2)+(27/4060000)*sech((1/812)*1218^(1/2)*(x+y-(203/100)*t))^2*tanh((1/812)*1218^(1/2)*(x+y-(203/100)*t))*1218^(1/2)*(1-tanh((1/812)*1218^(1/2)*(x+y-(203/100)*t))^2)

(5)

simplify(%);

0

(6)

 

NULL


Download solution_pde_check-exact.mw

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

restart;
a:=floor(33*sqrt(3)):
R2:=33^2*3:  # the radius^2
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");

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]

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

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 ReducedRowEchelonForm and append the identity matrix to your square matrix.

 

restart;

with(LinearAlgebra):

n:=3:

A:=RandomMatrix(n);

Matrix(3, 3, {(1, 1) = 27, (1, 2) = 99, (1, 3) = 92, (2, 1) = 8, (2, 2) = 29, (2, 3) = -31, (3, 1) = 69, (3, 2) = 44, (3, 3) = 67})

(1)

LinearAlgebra:-ReducedRowEchelonForm(<A|IdentityMatrix(n)>);

Matrix(3, 6, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = -3307/327244, (1, 5) = 2585/327244, (1, 6) = 5737/327244, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 2675/327244, (2, 5) = 4539/327244, (2, 6) = -1573/327244, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 1649/327244, (3, 5) = -5643/327244, (3, 6) = 9/327244})

(2)

B:=%[..,n+1..];

Matrix(3, 3, {(1, 1) = -3307/327244, (1, 2) = 2585/327244, (1, 3) = 5737/327244, (2, 1) = 2675/327244, (2, 2) = 4539/327244, (2, 3) = -1573/327244, (3, 1) = 1649/327244, (3, 2) = -5643/327244, (3, 3) = 9/327244})

(3)

A.B; #check

 

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1})

(4)

 

 

First 11 12 13 14 15 16 17 Last Page 13 of 122