Kitonum

15082 Reputation

24 Badges

12 years, 26 days

MaplePrimes Activity


These are answers submitted by Kitonum

with(InertForm):
Parse("2+3^4");
Parse("(2+3)^4");

                                           

With variable and number substitution instead of variable:

Expr:="2*x+3*x-(4*x+5)/(3*x+6)":
Parse(Expr);
Typeset(%, 'inert'=false);
Parse(StringTools:-Subs("x"="8",Expr));
value(%);
                                   

                                   

 


 

restart;
Eq1:=diff(f(x), x) = u(x):
Eq2:=diff(u(x), x) = x + f(x) - f(x)^2:
Inc:=f(0) = -1, u(0) = 1:

sol1 := dsolve([Eq1, Eq2, Inc], numeric, method = taylorseries[series]);

P:= plots:-odeplot(sol1, [x, f(x)], 0 .. 10, color = red);
T:= plot(-1 + x - x^2 +2/3*x^3 - 1/3*x^4 + 1/5*x^5, x=0..10, -1..7, color=blue);

plots:-display(P,T);

proc (x_taylorseries) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_taylorseries) else _xout := evalf(x_taylorseries) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _ctl, _n, _reinit, _y0, _yn, _yini, _pars, _ini, _par, _i; option `Copyright (c) 1994 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) Digits := 15; _xout := _xin; _ctl := array( 1 .. 21, [( 1 ) = (2), ( 2 ) = (0.), ( 3 ) = (0.), ( 4 ) = (0.1e-9), ( 5 ) = (0.1e-4), ( 6 ) = (22), ( 7 ) = (true), ( 9 ) = (2), ( 8 ) = (9), ( 11 ) = ({diff(Y[1](x), x)-Y[2](x), diff(Y[2](x), x)+Y[1](x)^2-x-Y[1](x)}), ( 10 ) = (0), ( 13 ) = (0), ( 12 ) = ({Y[1], Y[2]}), ( 15 ) = (Array(1..2, 0..22, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (1, 12) = 0, (1, 13) = 0, (1, 14) = 0, (1, 15) = 0, (1, 16) = 0, (1, 17) = 0, (1, 18) = 0, (1, 19) = 0, (1, 20) = 0, (1, 21) = 0, (1, 22) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (2, 12) = 0, (2, 13) = 0, (2, 14) = 0, (2, 15) = 0, (2, 16) = 0, (2, 17) = 0, (2, 18) = 0, (2, 19) = 0, (2, 20) = 0, (2, 21) = 0, (2, 22) = 0})), ( 14 ) = (Array(1..2, 0..22, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (1, 12) = 0, (1, 13) = 0, (1, 14) = 0, (1, 15) = 0, (1, 16) = 0, (1, 17) = 0, (1, 18) = 0, (1, 19) = 0, (1, 20) = 0, (1, 21) = 0, (1, 22) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (2, 12) = 0, (2, 13) = 0, (2, 14) = 0, (2, 15) = 0, (2, 16) = 0, (2, 17) = 0, (2, 18) = 0, (2, 19) = 0, (2, 20) = 0, (2, 21) = 0, (2, 22) = 0})), ( 18 ) = (1), ( 19 ) = (1), ( 16 ) = (0), ( 17 ) = (-1), ( 20 ) = (0), ( 21 ) = ({diff(Y[1](x), x)-Y[2](x), diff(Y[2](x), x)+Y[1](x)^2-x-Y[1](x)})  ] ); _n := trunc(_ctl[1]); _yini := Array(0..2, {(1) = 0., (2) = -1.}); _y0 := Array(0..2, {(1) = 0., (2) = -1.}); _yn := Array(1..2, {(1) = 0, (2) = 0}); _pars := []; _n := 2; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "taylorseries" elif _xout = "storage" then return type(eval(_ctl[10]), 'table') elif _xout = "order" then return _ctl[6] elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[17] = 0 then return evalf([evalf(_ctl[3]), seq(_yn[_i], _i = 1 .. _n)]) elif _ctl[17] = 1 then return evalf([_ctl[2], seq(_ctl[14][_i, 0], _i = 1 .. _n)]) else error "no information is available on last computed point" end if elif _xout = "enginedata" then return op(_ctl) elif _xout = "function" then if _ctl[9]-2. = 0 then return NULL else return eval(_ctl[13], 1) end if elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; if _ctl[9] = 2 then _ctl[11] := subs(_par, _ctl[21]) end if else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; if type(_ctl[10], 'table') then _ctl[10] := table() end if; _ctl[2] := _y0[0]; _ctl[3] := _y0[0]; _ctl[16] := 0.; _ctl[17] := 0; _ctl[18] := 1; _ctl[19] := 1; for _i to _n+nops(_pars) do _ctl[14][_i, 0] := _y0[_i] end do; if _Env_smart_dsolve_numeric = true then procname("right") := _xout; procname("left") := _xout end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(evalf(_y0[_i]), _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; if _ctl[9] = 2 then _ctl[11] := subs(_par, _ctl[21]) end if else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if _EnvInFsolve <> true and not _reinit and _ctl[17] = 0 and _xout-_ctl[3] = 0 then [_xout, seq(_yn[_i], _i = 1 .. _n)] else userinfo(3, {`dsolve/numeric`, `dsolve/continuation`}, `Continuation check: new point=`, _xout, `	previous point=`, _ctl[2], false, _reinit, _ctl[17], _ctl[16]); if _reinit or _ctl[17] <> 0 or _ctl[16] = 0. or sign(_ctl[16]) <> sign(_xout-_ctl[2]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) then userinfo(3, {`dsolve/numeric`, `dsolve/continuation`}, `Continuation check: restart from ICs`); for _i to _n+nops(_pars) do _ctl[14][_i, 0] := _y0[_i] end do; _ctl[2] := _y0[0]; _ctl[16] := 0.; _ctl[17] := 0; _ctl[18] := 1; _ctl[19] := 1 else userinfo(3, {`dsolve/numeric`, `dsolve/continuation`}, `Continuation check: continue solution`) end if; _ctl[3] := _xout; `dsolve/numeric/taylorseries`(x, _ctl, _yn); if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; if _EnvInFsolve = true then [_xout, seq(evalf[_EnvDSNumericSaveDigits](_yn[_i]), _i = 1 .. _n)] else [_xout, seq(evalf(_yn[_i]), _i = 1 .. _n)] end if end if end proc, (2) = Array(0..0, {}), (3) = [x, f(x), u(x)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_taylorseries, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_taylorseries, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_taylorseries, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_taylorseries, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_taylorseries), 'string') = rhs(x_taylorseries); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_taylorseries), 'string') = rhs(x_taylorseries)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_taylorseries) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_taylorseries) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

 

 

 

 

 


 

Download 2_graphs.mw

You can use the  ContoursWithLabels  procedure from the post  https://www.mapleprimes.com/posts/202222-Contour-Curves-With-Labels  for labelling. For convenience, I have included the text of this procedure in the worksheet below:

 

restart;
ContoursWithLabels := proc (Expr, Range1::(range(realcons)), Range2::(range(realcons)), Number::posint := 8, S::(set(realcons)) := {}, GraphicOptions::list := [color = black, axes = box], Coloring::`=` := NULL)
local r1, r2, L, f, L1, h, S1, P, P1, r, M, C, T, p, p1, m, n, A, B, E;
uses plots, plottools;
f := unapply(Expr, x, y);
if S = {} then r1 := rand(convert(Range1, float)); r2 := rand(convert(Range2, float));
L := [seq([r1(), r2()], i = 1 .. 205)];
L1 := convert(sort(select(a->type(a, realcons), [seq(f(op(t)), t = L)]), (a, b) ->is(abs(a) < abs(b))), set);
h := (L1[-6]-L1[1])/Number;
S1 := [seq(L1[1]+(1/2)*h+h*(n-1), n = 1 .. Number)] else
S1 := convert(S, list)  fi;
print(Contours = evalf[2](S1));
r := k->rand(20 .. k-20); M := []; T := [];
for C in S1 do
P := implicitplot(Expr = C, x = Range1, y = Range2, op(GraphicOptions), gridrefine = 3);
P1 := [getdata(P)];
for p in P1 do
p1 := convert(p[3], listlist); n := nops(p1);
if n < 500 then m := `if`(40 < n, (r(n))(), round((1/2)*n)); M := `if`(40 < n, [op(M), p1[1 .. m-11], p1[m+11 .. n]], [op(M), p1]); T := [op(T), [op(p1[m]), evalf[2](C)]] else
if 500 <= n then h := floor((1/2)*n); m := (r(h))(); M := [op(M), p1[1 .. m-11], p1[m+11 .. m+h-11], p1[m+h+11 .. n]]; T := [op(T), [op(p1[m]), evalf[2](C)], [op(p1[m+h]), evalf[2](C)]]
fi; fi; od; od;
A := plot(M, op(GraphicOptions));
B := plots:-textplot(T);
if Coloring = NULL then E := NULL else E := ([plots:-densityplot])(Expr, x = Range1, y = Range2, op(rhs(Coloring)))  fi;
display(E, A, B);
end proc:

A3 := 0.1098129220e-1*x-0.1864590943e-1*x^4-0.3780537764e-1+0.5300456762e-1*x^2-0.2252843255e-4*x^7*sin(6.280000000*y-.2083809520)-0.9011373022e-5*x^7*sin(6.280000000*y-1.256000000)-0.1569010873e-4*x^6*sin(6.280000000*y-.2083809520)-0.6276043495e-5*x^6*sin(6.280000000*y-1.256000000)-0.1462906959e-2*x^5*sin(6.280000000*y-1.256000000)-0.3657267398e-2*x^5*sin(6.280000000*y-.2083809520)-0.7271311325e-3*x^4*sin(6.280000000*y-1.256000000)-0.1817827831e-2*x^4*sin(6.280000000*y-.2083809520)-0.8684412118e-1*x^3*sin(6.280000000*y-1.256000000)-.2171103030*x^3*sin(6.280000000*y-.2083809520)-0.2589262297e-1*x^2*sin(6.280000000*y-1.256000000)-0.6473155744e-1*x^2*sin(6.280000000*y-.2083809520)+.3752380081*x*sin(6.280000000*y-1.256000000)+.9380950203*x*sin(6.280000000*y-.2083809520)+0.3750477138e-1*sin(6.280000000*y-1.256000000)+0.9376192845e-1*sin(6.280000000*y-.2083809520)-0.3550918239e-3*x^6-0.760666870e-4*x^5+0.1713654921e-5*x^7-0.7811069699e-2*x^3-3.822344860*10^(-8)*x^8;


# Your example

ContoursWithLabels(A3, -2 .. 2, 0 .. 2, {seq(-2..2,0.2)}, [color = black, thickness = 2, axes = box, size=[1000,500], scaling=constrained], Coloring = [colorstyle = HUE, colorscheme = ["Blue", "Green", "Yellow", "Red"], style = surface]);

0.1098129220e-1*x-0.3780537764e-1-0.1864590943e-1*x^4+0.5300456762e-1*x^2+0.1713654921e-5*x^7+0.9376192845e-1*sin(6.280000000*y-.2083809520)+0.3750477138e-1*sin(6.280000000*y-1.256000000)-0.3550918239e-3*x^6-0.760666870e-4*x^5-0.7811069699e-2*x^3-0.3822344860e-7*x^8-0.2252843255e-4*x^7*sin(6.280000000*y-.2083809520)-0.9011373022e-5*x^7*sin(6.280000000*y-1.256000000)-0.1569010873e-4*x^6*sin(6.280000000*y-.2083809520)-0.6276043495e-5*x^6*sin(6.280000000*y-1.256000000)-0.1462906959e-2*x^5*sin(6.280000000*y-1.256000000)-0.3657267398e-2*x^5*sin(6.280000000*y-.2083809520)-0.7271311325e-3*x^4*sin(6.280000000*y-1.256000000)-0.1817827831e-2*x^4*sin(6.280000000*y-.2083809520)-0.8684412118e-1*x^3*sin(6.280000000*y-1.256000000)-.2171103030*x^3*sin(6.280000000*y-.2083809520)-0.2589262297e-1*x^2*sin(6.280000000*y-1.256000000)-0.6473155744e-1*x^2*sin(6.280000000*y-.2083809520)+.3752380081*x*sin(6.280000000*y-1.256000000)+.9380950203*x*sin(6.280000000*y-.2083809520)

 

Contours = [-2., -1.8, -1.6, -1.4, -1.2, -1.0, -.8, -.6, -.4, -.2, 0., .2, .4, .6, .8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]

 

 

 


Of course, you can change the color scheme and other parameters as you wish.

 

Download ContoursWithLabels1.mw

Here is an example of how to do this:
 

R:=solve(x^5+2*x-1);

RootOf(_Z^5+2*_Z-1, index = 1), RootOf(_Z^5+2*_Z-1, index = 2), RootOf(_Z^5+2*_Z-1, index = 3), RootOf(_Z^5+2*_Z-1, index = 4), RootOf(_Z^5+2*_Z-1, index = 5)

(1)

k:=0:
for r in R do
if is(r>0) and is(r<1) then k:=k+1; S[k]:=r fi;
od:
S:=convert(S, list);
evalf[20](S);

[RootOf(_Z^5+2*_Z-1, index = 1)]

 

[.48638903593454300002]

(2)

 


 

Download pick.mw

Here is another approach that directly generalizes the  coeff  command to the case of polynomials in several variables:
 

restart;

coefff:=proc(P,T::{list,set}(name),t)
local L,H,i,k:
L:=[coeffs(P,T,'h')]: H:=[h]: k:=0:
for i from 1 to nops(H) do
if H[i]=t then k:=L[i] fi:
od:
k;
end proc:


Examples of use

F:= 2*x+6*y+4*x^2+12*x*y-5*y^2;

F := 4*x^2+12*x*y-5*y^2+2*x+6*y

(1)

coefff(F,{x,y},x);
coefff(F,{x,y},x^2);
coefff(F,{x,y},x*y);

2

4

12

(2)

Q:=x^3-3*x^2*y*z+4*z^5;
coefff(Q,[x,y,z],z);
coefff(Q,[x,y,z],x^2*y*z);
coefff(Q,[x,y,z],x^3);
coefff(Q,[x,y,z],z^5);

Q := 4*z^5-3*x^2*y*z+x^3

0

-3

1

4

(3)

 


 

Download coefff.mws

The plot shows that your function (Ackley function for 2 variables) has a large number of critical points (at which both partial derivatives are equal to 0), among which there are points of a local minimum, maximum, and also saddle points. You can investigate each specific critical point using  Student:-MultivariateCalculus:-SecondDerivativeTest  command. Thus, the problem boils down to finding a sufficiently large number of critical points. Of course, the best option is to use DirectSearch package (vv"s suggestion). But if you do not have this package installed, you can use thу above command, having previously built the curves to find the necessary ranges for using  fsolve  command. See the worksheet below for critical points in the ranges  x=0..1, y=0..1.


 

restart;
f:=(x,y)->-20*exp(-0.2*sqrt(0.5*(x^2+y^2)))-exp(0.5*(cos(2*Pi*x)+cos(2*Pi*y)))+exp(1)+20;

proc (x, y) options operator, arrow; -20*exp(-.2*sqrt(.5*x^2+.5*y^2))-exp(.5*(cos(2*Pi*x)+cos(2*Pi*y)))+exp(1)+20 end proc

(1)

plot3d(f(x,y), x=-3..3, y=-3..3, style=surface, grid=[100,100]);

 

plots:-implicitplot([diff(f(x,y),x),diff(f(x,y),y)], x=0..1, y=0..1, color=[red,blue], gridrefine=3);

 

P1:=fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x=0.6..0.8, y=0.6..0.8});
P2:=fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x=0.9..1, y=0.9..1});
P3:=fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x=0.5..0.7, y=0.8..1});
P4:=fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x=0.8..1, y=0.5..0.7});

{x = .6730965200, y = .6730965200}

 

{x = .9684776586, y = .9684776586}

 

{x = .5774125352, y = .8747620388}

 

{x = .8747620388, y = .5774125352}

(2)

with(Student[MultivariateCalculus]):
SecondDerivativeTest(f(x,y), [x, y] = eval([x,y],P1));
SecondDerivativeTest(f(x,y), [x, y] = eval([x,y],P2));
SecondDerivativeTest(f(x,y), [x, y] = eval([x,y],P3));
SecondDerivativeTest(f(x,y), [x, y] = eval([x,y],P4));

LocalMin = [], LocalMax = [[.6730965200, .6730965200]], Saddle = []

 

LocalMin = [[.9684776586, .9684776586]], LocalMax = [], Saddle = []

 

LocalMin = [], LocalMax = [], Saddle = [[.5774125352, .8747620388]]

 

LocalMin = [], LocalMax = [], Saddle = [[.8747620388, .5774125352]]

(3)

plot3d(f(x,y), x=0..1.1, y=0..1.1, view=2..6, grid=[100,100], axes=normal, orientation=[25,75]);  # The visual check

 

 


 

Download saddle_points.mw


 

restart;
ex7 := x->abs(x^2-3*x)+4*x-6;
solve(ex7(x)); # Easy way

proc (x) options operator, arrow; abs(x^2-3*x)+4*x-6 end proc

 

1, -3

(1)

solve({x^2-3*x<0, -(x^2-3*x)+4*x-6=0}) union solve({x^2-3*x>=0, x^2-3*x+4*x-6=0});  # A longer way without the abs command

{x = -3, x = 1}

(2)

 


 

Download abs.mw

In those cases when  isolve  does not cope with the problem, brute force method (i.e. the method of enumerating all possible options) can often be successfully used. In your 4-squares example, assuming the variables are positive, it is obvious that each of the variables can vary in the range of 1 .. 12 . Nested for-loops are often used for this method, but I like nested sequences more. We get  76  solutions, but many of them are obtained simply by rearranging one from the other. But a slight modification of this method allows you to immediately get  5  unique solutions.


 

restart;

seq(seq(seq(seq(`if`(169 = w^2 + x^2 + y^2 + z^2,[w,x,y,z],NULL),w=1..12),x=1..12),y=1..12),z=1..12);
nops([%]);

[10, 8, 2, 1], [8, 10, 2, 1], [10, 2, 8, 1], [2, 10, 8, 1], [8, 2, 10, 1], [2, 8, 10, 1], [10, 8, 1, 2], [8, 10, 1, 2], [10, 7, 4, 2], [7, 10, 4, 2], [10, 4, 7, 2], [4, 10, 7, 2], [10, 1, 8, 2], [1, 10, 8, 2], [8, 1, 10, 2], [7, 4, 10, 2], [4, 7, 10, 2], [1, 8, 10, 2], [10, 7, 2, 4], [7, 10, 2, 4], [11, 4, 4, 4], [4, 11, 4, 4], [8, 8, 5, 4], [9, 6, 6, 4], [6, 9, 6, 4], [10, 2, 7, 4], [2, 10, 7, 4], [8, 5, 8, 4], [5, 8, 8, 4], [6, 6, 9, 4], [7, 2, 10, 4], [2, 7, 10, 4], [4, 4, 11, 4], [8, 8, 4, 5], [8, 4, 8, 5], [4, 8, 8, 5], [9, 6, 4, 6], [6, 9, 4, 6], [9, 4, 6, 6], [4, 9, 6, 6], [6, 4, 9, 6], [4, 6, 9, 6], [10, 4, 2, 7], [4, 10, 2, 7], [10, 2, 4, 7], [2, 10, 4, 7], [4, 2, 10, 7], [2, 4, 10, 7], [10, 2, 1, 8], [2, 10, 1, 8], [10, 1, 2, 8], [1, 10, 2, 8], [8, 5, 4, 8], [5, 8, 4, 8], [8, 4, 5, 8], [4, 8, 5, 8], [5, 4, 8, 8], [4, 5, 8, 8], [2, 1, 10, 8], [1, 2, 10, 8], [6, 6, 4, 9], [6, 4, 6, 9], [4, 6, 6, 9], [8, 2, 1, 10], [2, 8, 1, 10], [8, 1, 2, 10], [7, 4, 2, 10], [4, 7, 2, 10], [1, 8, 2, 10], [7, 2, 4, 10], [2, 7, 4, 10], [4, 2, 7, 10], [2, 4, 7, 10], [2, 1, 8, 10], [1, 2, 8, 10], [4, 4, 4, 11]

 

76

(1)

seq(seq(seq(seq(`if`(169 = w^2 + x^2 + y^2 + z^2,[w,x,y,z],NULL),w=1..x),x=1..y),y=1..z),z=1..12);

[4, 5, 8, 8], [4, 6, 6, 9], [2, 4, 7, 10], [1, 2, 8, 10], [4, 4, 4, 11]

(2)

 


 

Download 4_squares.mw
 

 

restart;
A := [[1, 0, 9], [5, 0, 6], [4, 0, 0]]; 
remove~(has, A, 0);

                                  [[1, 9], [5, 6], [4]]


Here is a procedure for this:

restart;
P:=proc(k)
local n, k1, k2;
n:=2; k1:=0; k2:=0;
while n<10^k do
if isprime(10*n+1) then k2:=k2+1 else k1:=k1+1 fi;
n:=nextprime(n);
od;
[k1,k2];
end proc:

Examples of use

seq(P(k), k=1..6);

[2, 2], [18, 7], [131, 37], [1021, 208], [8240, 1352], [69231, 9267]

(1)

plot([[seq([k,P(k)[1]],k=1..5)],[seq([k,P(k)[2]],k=1..5)]], color=[red,blue]);

 

 

 

Addition. It is interesting that if we set the logarithmic scale along the vertical axis, then we get 2 lines close to straight lines:

plots:-logplot([[seq([k,P(k)[1]],k=1..5)],[seq([k,P(k)[2]],k=1..5)]], color=[red,blue]);

Edit.

Download count1.mw

The appearance of extraneous roots when squaring both sides of the equation is in no way related to the presence of square roots in the original equation. Consider the simplest equation  x = 2  with an obvious root  2 . But if you square it, you get an equation  x^2 = 4  that has 2 roots  2  and  -2 .  In the general case, let us have the equation  f(x) = g(x)  and square it. Obviously, the resulting equation  f(x)^2  = g(x)^2  is equivalent to the equation  (f(x) - g(x))*(f(x) + g(x)) = 0 .  The roots of the last equation are the union of the sets of the roots of the original equation  f(x) = g(x) and the roots of the equation  f(x) = - g(x) , which may have roots that are not the roots of the original equation.

The function  f , by means of which the equation of the blue line is specifed, is given incorrectly. Below is the corrected version:

restart;
with(plots):
unprotect(O);
H := [84/37, 14/37]:
f := x-> -6*x+14:
a := 3: A := [a, f(a)]:
O:=[0,0]: zo := [8/3+I*f(8/3)];
ze := [2+I^(eval(diff(f(x), x), x = 2))];
Zo := [8/3, f(8/3)]; Ze := [2, f(2)];
ex := -6*x+14; V := `<,>`(ze-zo);
V := plottools:-arrow(Zo, Ze, .2, .4, .2, color = "Red"):
DR := plot(ex, x = -1 .. 3, color = blue):
Points := pointplot([A[],Zo[],Ze[],H[]], symbol = solidcircle, color = red, symbolsize = 10):
T := textplot([[A[], "A"],[H[],"H"],[O[],"O"]], font = [times, 18], align = {below, right}):
display([DR,V,Points, T], scaling = constrained, size=[800,800]);

zo := [8/3-2*I]

 

ze := [1]

 

Zo := [8/3, -2]

 

Ze := [2, 2]

 

ex := -6*x+14

 

Vector[column](%id = 18446746163433603070)

 

 

 


Explanation. The direction vector for the blue line is  <-2/3, 4> (in fact, you specified it with a complex number  v = -2/3+4*I ). Therefore, the slope of this line is  k = 4 / (- 2/3) = - 6  and the equation of this line is  y=-4 - 6*(x-3)  or  y=-6*x + 14

Download 1.mw

I made 3 corrections in your procedure. Now it is working properly. It is only important to note that this construction  List1 := [op(List1), p]  and so on  is extremely inefficient. Below is the same procedure  (named  FF ), but much faster (of course for large ):

restart;
F := proc(N)
local p, List1, List3, List5, List7;
List1 := [];
List3 := [];
List5 := [];
List7 := [];
for p from 1 to N do
if isprime(p) and p mod 8 = 1 then List1 := [op(List1), p]; end if;
if isprime(p) and p mod 8 = 3 then List3 := [op(List3), p]; end if;
if isprime(p) and p mod 8 = 5 then List5 := [op(List5), p]; end if;
if isprime(p) and p mod 8 = 7 then List7 := [op(List7), p]; end if;
[nops(List1), nops(List3), nops(List5), nops(List7)];
end do;
end proc:

  F(100);
                          [5, 7, 6, 6]
 

restart;
FF := proc(N)
local p, k1, k3, k5, k7;
k1:=0; k3:=0; k5:=0; k7:=0;
for p from 1 to N do
if isprime(p) then if p mod 8 = 1 then k1:=k1+1  elif
p mod 8 = 3 then k3:=k3+1  elif
p mod 8 = 5 then k5:=k5+1  elif
p mod 8 = 7 then k7:=k7+1 fi; fi;
[k1,k3,k5,k7];
end do;
end proc:

FF(100);
                                [5, 7, 6, 6]


Edit.
 

In the  seq  command, index is  local  and therefore independent of the previously assigned value:

i:=1: ii:=2:
seq(i, i=1..5);
seq(ii, ii=1..5)

                                     1, 2, 3, 4, 5
                                     1, 2, 3, 4, 5


A workaround:

a:=50:
seq(a,ii=1..5);

 

Example:

i:=3: n:=8:
mul(`if`(k<>i,lambda[i]-lambda[k],1), k=1 .. n);

           

 

2 3 4 5 6 7 8 Last Page 4 of 235