## 15082 Reputation

12 years, 26 days

## InertForm:-Parse command automates this...

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

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

## A way...

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

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

## A way...

Here is an example of how to do this:

 > R:=solve(x^5+2*x-1);
 (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);
 (2)
 >

## A generalization of the coeff command....

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;
 (1)
 > coefff(F,{x,y},x); coefff(F,{x,y},x^2); coefff(F,{x,y},x*y);
 (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);
 (3)
 >

## Student:-MultivariateCalculus:-SecondDer...

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

## abs, solve...

 > restart; ex7 := x->abs(x^2-3*x)+4*x-6; solve(ex7(x)); # Easy way
 (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
 (2)
 >

## Brute force...

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([%]);
 (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);
 (2)
 >

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

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

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

## Procedure...

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

## Squaring Equations...

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

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

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.

## local...

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

## mul, `if`...

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
﻿