Advanced Search

Kitonum

20084 Reputation

26 Badges

17 years, 29 days
Contact Kitonum

MaplePrimes Activity


These are answers submitted by Kitonum

rsolve...

Kitonum 20084
June 21 2020
1 3


 

> 

restart;
a:=rsolve({a(n+1)=sqrt(x*a(n)), a(1)=x}, a(n), makeproc);

proc (n) local L, i, val; options remember, system, `Copyright (c) 2003 by Waterloo Maple Inc.`; if not type(n, integer) then error "input must be an integer" elif n = 1 then x elif n < 1 then L := [x]; for i from 0 by -1 to n do val := traperror(L[1]^2/x); if val = lasterror then error "unable to compute recurrence for n<%1", i+1 end if; L := [val, op(1 .. -2, L)] end do; L[1] elif 1 < n then L := [x]; for i from 2 to n do val := traperror((x*L[1])^(1/2)); if val = lasterror then error "unable to compute recurrence for n>%1", i-1 end if; L := [op(2 .. -1, L), val] end do; L[-1] end if end proc

 (1)
> 

L:=[seq(a(n), n=1..10)];
simplify(L) assuming x>=0;
simplify(L) assuming x<0;
 

[x, (x^2)^(1/2), (x*(x^2)^(1/2))^(1/2), (x*(x*(x^2)^(1/2))^(1/2))^(1/2), (x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2), (x*(x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2), (x*(x*(x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2), (x*(x*(x*(x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2), (x*(x*(x*(x*(x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2), (x*(x*(x*(x*(x*(x*(x*(x*(x^2)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2)]

  

[x, x, x, x, x, x, x, x, x, x]

  

[x, -x, -I*x, (-1/2+(1/2)*I)*x*2^(1/2), -(1/2)*x*2^(1/4)*(-2+2*I)^(1/2), -(1/2)*x*2^(5/8)*(-(-2+2*I)^(1/2))^(1/2), -(1/2)*x*2^(13/16)*(-(-(-2+2*I)^(1/2))^(1/2))^(1/2), -(1/2)*x*2^(29/32)*(-(-(-(-2+2*I)^(1/2))^(1/2))^(1/2))^(1/2), -(1/2)*x*2^(61/64)*(-(-(-(-(-2+2*I)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2), -(1/2)*x*2^(125/128)*(-(-(-(-(-(-2+2*I)^(1/2))^(1/2))^(1/2))^(1/2))^(1/2))^(1/2)]

 (2)
> 

 



It is obvious that for any x> 0   a(n) = x  for any n , so if x=2 then  a(100)=2 .
 

Download rsolve1.mw

collect...

Kitonum 20084
June 20 2020
1 0


 

> 

restart;
> 

sol:=W(x)=_C1*(cosh(alpha*x)-sinh(alpha*x))+_C2*(cosh(alpha*x)+sinh(alpha*x))+_C3*sin(alpha*x)+_C4*cos(alpha*x);
W:=eval(W(x),sol);

W(x) = _C1*(cosh(alpha*x)-sinh(alpha*x))+_C2*(cosh(alpha*x)+sinh(alpha*x))+_C3*sin(alpha*x)+_C4*cos(alpha*x)

  

_C1*(cosh(alpha*x)-sinh(alpha*x))+_C2*(cosh(alpha*x)+sinh(alpha*x))+_C3*sin(alpha*x)+_C4*cos(alpha*x)

 (1)
> 

F:=[sinh(alpha*x),cosh(alpha*x),sin(alpha*x),cos(alpha*x)]:
W1:=collect(W, F);
assign(([D1,D2,D3,D4]=~[coeffs(W1, F)])[]);

(-_C1+_C2)*sinh(alpha*x)+(_C1+_C2)*cosh(alpha*x)+_C3*sin(alpha*x)+_C4*cos(alpha*x)

 (2)
> 

D1, D2, D3, D4;

_C1+_C2, -_C1+_C2, _C3, _C4

 (3)
> 

 


 

Download collect_new.mw

combinat:-nextperm, combinat:-nextcomb...

Kitonum 20084
June 19 2020
0 1

The Iterator package was introduced in Maple 2016. If you have an older version of Maple, then you can use the commands  combinat:-nextperm  and  combinat:-nextcomb .

primpart, content...

Kitonum 20084
June 19 2020
2 0 
primpart(98-28*sqrt(7), sqrt(7));
content(98-28*sqrt(7), sqrt(7));
%*``(%%);
expand(%);

                                   

Or manually:

Expr:=98-28*sqrt(7);
d:=igcd(98,28);
d*``(Expr/d);

 

infinity...

Kitonum 20084
June 16 2020
3 0

The integral of a periodic function is the same on any period-long interval. The contradiction below is due to limited accuracy in approximate calculation. With increasing accuracy, both methods yield consistent results. Infinity is explained by the fact that the function is discontinuous, its denominator is 0 at some points.


 

> 

restart;
x(t) := -3.703703704*10^(-7)*(0.000111668023*cos(1000/33*sqrt(1122)*t) - 0.0001214712007*sin(1000/33*sqrt(1122)*t) - 0.0002325581396*sqrt(561)*sqrt(2)*(-0.0004467462845*sqrt(1122)*sin(1000/33*sqrt(1122)*t) + 0.0004467462845*sqrt(1122)*cos(1000/33*sqrt(1122)*t)))/((2.074226433*10^14*cos(1000/33*sqrt(1122)*t) + 2.074226433*10^14*sin(1000/33*sqrt(1122)*t))*(4.895037587*10^(-11) + 0.01685634229*(0.00001474262739*cos(1000/33*sqrt(1122)*t) + 0.00001474262739*sin(1000/33*sqrt(1122)*t))^2)^2);
P:= 33*Pi/1000/sqrt(1122):

# Digits=10
RMS:= evalf(sqrt(1/2/P*Int(x(t)^2, t= -P..P)));
RMS:= evalf(sqrt(1/2/P*Int(x(t)^2, t= 0..2*P)));

Digits:=20:
RMS:= evalf(sqrt(1/2/P*Int(x(t)^2, t= -P..P)));
RMS:= evalf(sqrt(1/2/P*Int(x(t)^2, t= 0..2*P)));
 

-0.3703703704e-6*(0.111668023e-3*cos((1000/33)*1122^(1/2)*t)-0.1214712007e-3*sin((1000/33)*1122^(1/2)*t)-0.2325581396e-3*561^(1/2)*2^(1/2)*(-0.4467462845e-3*1122^(1/2)*sin((1000/33)*1122^(1/2)*t)+0.4467462845e-3*1122^(1/2)*cos((1000/33)*1122^(1/2)*t)))/((0.2074226433e15*cos((1000/33)*1122^(1/2)*t)+0.2074226433e15*sin((1000/33)*1122^(1/2)*t))*(0.4895037587e-10+0.1685634229e-1*(0.1474262739e-4*cos((1000/33)*1122^(1/2)*t)+0.1474262739e-4*sin((1000/33)*1122^(1/2)*t))^2)^2)

  

0.3200233260e-5

  

Float(infinity)

  

Float(infinity)

  

Float(infinity)

 (1)
> 

 


 

Download Infinity.mw

Plotting surfaces of revolution...

Kitonum 20084
June 15 2020
1 0

To plot surfaces of revolution it is convenient to use their parametric equations in which to use polar coordinates on the plane. The 1 ex. :

restart;
x:=r*cos(phi): y:=r*sin(phi):
plot3d([[x,y,r],[x,y,4]], r=0..4, phi=0..2*Pi, style=surface, color=grey, scaling=constrained, labels=["x","y",z]);

                        

InertForm:-Parse command automates this...

Kitonum 20084
June 13 2020
1 0 
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(%);
                                   

                                   

 

plots:-display...

Kitonum 20084
June 12 2020
0 1


 

> 

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

A way...

Kitonum 20084
June 12 2020
1 4

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

A way...

Kitonum 20084
June 11 2020
0 0

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

A generalization of the coeff command....

Kitonum 20084
June 11 2020
2 0

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

Student:-MultivariateCalculus:-SecondDer...

Kitonum 20084
June 08 2020
0 1

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

abs, solve...

Kitonum 20084
June 07 2020
1 0


 

> 

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

Brute force...

Kitonum 20084
June 06 2020
0 1

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
 

 

remove~(has, A, 0);...

Kitonum 20084
June 04 2020
1 2 
restart;
A := [[1, 0, 9], [5, 0, 6], [4, 0, 0]]; 
remove~(has, A, 0);

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

First 46 47 48 49 50 51 52 Last Page 48 of 280