## 5989 Reputation

7 years, 266 days

## Where does this question com from?...

I agree with Rouben.
You can easily find (on the web or in your math books) the equation of the tangent plane to a sphere and try to code them in Maple.
Maybe this is what you have been asked for?

If you are not comfortable with maths and need to solve this problem quickly, you can use the geom3d package:
With_geom3d.mw

Maybe reading carefully the three procedures TangentPlane, IsOnObject and onobjl could help you understand how to derive the equations.

## If no constraints do exist...

here is an extremely simple way do draw random triangles

 > restart
 > with(plots): with(plottools): with(Statistics):
 > interface(version);
 (1)
 > randomize(); c := rand(0. .. 1): r := rand(0. .. 2.*Pi): p := [seq(r(), k=1..3)]: display(   plottools:-circle([0, 0], 1, color=gray),   seq(     PLOT(       CURVES(         [sin,cos]~([seq(r(), k=1..3)])[[\$1..3, 1]]         , COLOR(RGB, c(), c(), c())         , THICKNESS(2)       )     )     , n=1..3   ) )
 >

Let me know if you have some constraints on the triangles.
For instance here is a simple way to generate triangles which contain the center of the circle:

 > restart
 > with(plots): with(plottools): with(Statistics):
 > interface(version);
 (1)
 > randomize(); c := rand(0. .. 1): p := 0, rand(0. .. 1.*Pi)(): q := rand(1.*Pi.. p[2]+1.*Pi)(): p := [p[], q]: T := PLOT(        CURVES(          [cos, sin]~(p)[[\$1..3, 1]]          , COLOR(RGB, c(), c(), c())          , THICKNESS(2)        )      ): display(   plottools:-circle([0, 0], 1, color=gray),   rotate(T, rand(0. .. 2.*Pi)()) )
 >

Note that the procedure in RandomTriangle_2.mw doesn't produce uniformly distributed triangles in the following sense: while p[2] is obviously uniformly distributed over [0, Pi], one can see that q (the third point p[3]) is not uniformly distributed over [Pi, 2*Pi]

```# Do this to convince you
K := 10000:
P := Matrix(K, 2):
for k from 1 to K do
p := rand(0. .. 1.*Pi)():
q := rand(1.*Pi.. p+1.*Pi)():
P[k, 1] := p:
P[k, 2] := q:
end do:
Histogram(P[..,1]);
Histogram(P[..,2]);
```

Here is the exact expression of the PDF of q:

```p2 := RandomVariable(Uniform(0, Pi));
q  := Pi + RandomVariable(Uniform(0, 1)) * p2;

# and thus
combine(PDF(Y, y), ln);
/                                 /  -y + Pi\
|                               ln|- -------|
|                                 \    Pi   /
piecewise|y - Pi < 0, 0, y - Pi <= Pi, - -------------,
\                                    Pi

\
|
|
Pi < y - Pi, 0|
/

```

## A non optimized code...

(it's likely that using the remember option would be more astute)

I used a predator velocity equal to the prey velocivy divided by the eccentricity of the ellipses (see me reply under the login sand15)

 > restart
 > with(plots):
 > # Outer ellipse OuterEllipse := (x/a)^2+(y/b)^2-1; # Abscissa of the left focus F1 := -sqrt(a^2-b^2)
 (1)
 > # Movement equations interface(warnlevel=0): {   diff(X__T(t), t) = V__T,   diff(X__P(t), t) = V__P*(X__P(t)-X__T(t))/sqrt((X__P(t)-X__T(t))^2+Y__P(t)^2),   diff(Y__P(t), t) = V__P*Y__P(t)/sqrt((X__P(t)-X__T(t))^2+Y__P(t)^2),   X__T(0) = F1,   X__P(0) = x__0,   Y__P(0) = b*surd(1-(x__0/a)^2, 2) }; sol := dsolve(%, numeric, parameters=[V__T, a, b, x__0, V__P], events=[[Y__P(t)=0, halt]]); data := [V__T = 1, a = 2, b = 1, x__0 = 0]: pars := [rhs~(data)[], eval(V__T*a/F1, data)]: sol(parameters=pars); sol(abs(pars[-1])*2); CaptureTime := sol(eventfired=[1])[1]; anim := proc(s) plots:-odeplot(   sol   , [[X__T(t), 0], [X__P(t), Y__P(t)]]   , t=0..s   , color=[blue, red]   , thickness=[3, 3]   , labels=["", ""]   , legend=[Prey,Predator] ); end proc: animate(anim, [s], s=1e-6..CaptureTime)
 > f := proc(a, b, V__T)   local Outer, F1, V__P, sys, sol, Traj, s, CaptureTime, anim, ell, ELL:   uses plots:   interface(warnlevel=0):   Outer := (x/a)^2+(y/b)^2-1;   F1    := -sqrt(a^2-b^2);   V__P  := V__T*a/F1;   sys := {     diff(X__T(t), t) = V__T,     diff(X__P(t), t) = V__P*(X__P(t)-X__T(t))/sqrt((X__P(t)-X__T(t))^2+Y__P(t)^2),     diff(Y__P(t), t) = V__P*Y__P(t)/sqrt((X__P(t)-X__T(t))^2+Y__P(t)^2),     X__T(0) = F1,     X__P(0) = x__0,     Y__P(0) = b*surd(1-(x__0/a)^2, 2)   };   sol := dsolve(sys, numeric, parameters=[x__0], events=[[Y__P(t)=0, halt]]);   sol(parameters=[0]);   sol(abs(V__P)*2);   CaptureTime := sol(eventfired=[1])[1];   anim := proc(tau)     local Traj:= NULL:     for s in [-a+1e-6, evalf(a*~cos~(Pi*~[\$1..7]/8))[], a-1e-6] do       sol(parameters=[s]);       sol(abs(V__P)*2);          Traj := Traj,         plots:-odeplot(         sol         , [[X__T(t), 0], [X__P(t), Y__P(t)]]         , t=0..tau         , color=[blue, red]         , thickness=[4, 2]         , labels=["", ""]       ):       display(Traj):     end do:   end proc:   ell := implicitplot(Outer, x=-a..a, y=0..b, color=gray);   ELL := display( seq(plottools:-scale(ell, k, k), k in [seq](0.2..1, 0.2)) ):   #(print@animate)(   animate(     anim, [tau], tau=1e-6..CaptureTime     , background=ELL     , scaling=constrained    , view=[-a..a, 0..b]    , size=[1000, ceil(1000*b/a)]   ) end proc: f(2, 1, 1)
 >

## Several possibilities...

to get he outputs you want.
Among them:

 > restart;
 > NN := [4, 6, 8]: a := 0: b := 2: n := 4:
 > h := evalf((b-a)/n): print("The integration domain [a, b] = ", [a, b]);
 (1)
 > f := exp(x): print("The given function is ", f);
 (2)
 > Exact := evalf(int(f, x = a .. b)): text  := cat("The exact integration in ", [a, b], " is "): print(text, Exact); # or Exact := evalf(int(f, x = a .. b)): text1 := "The exact integration in": text2 := "is": print(text1, [a, b], text2, Exact); # or Exact := evalf(int(f, x = a .. b)): text1 := "The exact integration in": text2 := "is": cattext := cat(text1, " ", convert([a, b], string), " ", text2, " ", convert(Exact, string)): print(cattext); # Finally printf("%s %a %s %1.5f\n", text1, [a, b], text2, Exact)
 The exact integration in [0, 2] is 6.38906
 > text1 := "The value of h to divide the domain": text2 := "into": text3 := "subintervals is": cattext := convert(cat(text1, " ", [a, b], " ", text2, " ", n, " ", text3, " ", evalf[5](h)), string): print(cattext); # or cattext := convert(cat(text1, " ", [a, b], " ", text2, " ", n, " ", text3, " ", identify(h)), string): print(cattext); # Or printf("%s %a %s %d %s %1.5f\n", text1, [a, b], text2, n, text3, h); printf("%s %a %s %d %s %a\n", text1, [a, b], text2, n, text3, identify(h));
 The value of h to divide the domain [0, 2] into 4 subintervals is 0.50000 The value of h to divide the domain [0, 2] into 4 subintervals is 1/2
 > # Customize this one as you want print("Numerical integration in [a,b] is going to perform when h via RECTANGULAR METHOD for n = ", n);
 (3)
 >

A help page that might interest you

```help(Student:-Calculus1:-RiemannSum);
```

Riemann_sum_2.mw

Finally, better presentations (IMO) can be obtained, for instance with DocumentTools:-Layout.

## Main question: What does the magenta cur...

Drawing a parabolic cylinder with generatrix paralel to axis (for instance) is easy:

```restart
with(plots):

implicitplot3d(
y-x^2
, x=-1..1, y=0..1, z=0..1, grid=[40\$3]
, style=surface, color=cyan, transparency=0.3
, axes=normal
, seq(axis[k]=[tickmarks=0], k=1..3)
);
```

Drawing the magenta curve requires you define what this curve is.
Once done: unfold the parabolic cylinder (this gives the plane y=0) and draw the curve C on this plane; next "fold" this plane to retrieve the original parabolic cylinder.
The curve C' you get is the image of C, drawn on the parabolic cylinder.
The parametric equation of curve C' is

`[x, f(x), g(x)] for any x in some drawing interval`

Where:

• f := x -> a*x^2 + b*x + c represents the equation of any cross section oc the paramolic cylinder;
• z = g(x) is the explicit equation of curve C in plane y=0.

Here is an example when the curve C is given by g : x -> g(x)=x
pc0.mw

The file above contains another example

pc.mw

## A lot of typos (u instead of U) an...

A lot of typos (u instead of U) and a lot of errors when writting diff(...)

(it seems you didn't learn from the previous Q&A on the same topic)

BTW: why do you use negative signs in the expansion U(xi[n]) as equation (10) uses positive signs?

 > restart:
 > local psi

 > # relation (11) psi := zeta -> (r[1]*exp(s[1]*zeta)+r[2]*exp(s[2]*zeta))                /                (r[3]*exp(s[3]*zeta)+r[4]*exp(s[4]*zeta)); dpsi := unapply(simplify(diff(psi(zeta),zeta)/psi(zeta), size), zeta)
 (1)
 (2)
 > # U... (m=1) # (xi instead of zeta to be consistent with uour notations) U := unapply(u(xi, 1), xi);
 (3)
 > # The original equation (corrected, see yellow highlighted  characters). eqU := c*(diff(U(xi[n]), xi[n]))*(U(xi[n])+U(xi[n-1]))*(U(xi[n])+U(xi[n+1]))-(2*(U(xi[n-1])-U(xi[n+1])))*(U(xi[n])^2)(1-U(xi[n])^2): # Note that you wrote diff(U, xi[n])... which is null: diff(U, xi[n]); print(): # instead of diff(U(xi[n]), xi[n]): diff(U(xi[n]), xi[n]);
 (4)
 > # eqU equation rewritten with V instead of U (just to make you understand what happens) eqV := c*(diff(V(xi[n]), xi[n]))*(V(xi[n])+V(xi[n-1]))*(V(xi[n])+V(xi[n+1]))-(2*(V(xi[n-1])-V(xi[n+1])))*(V(xi[n])^2)(1-V(xi[n])^2)
 (5)
 > # Rewritting Rule rr := (f, i) -> f(xi[n+i]) = f(xi[n]+i*d)
 (6)
 > # Rewrite eqV using rr eqVrr := eval(eqV, [rr(V, -1), rr(V, +1)]);
 (7)
 > # To get "your" equation jus do this # Uncomment tosee this lengthy expression eq := eval(eqVrr, V=U): length(eq);
 (8)
 >

## Another solution...

Less high-end programming but more readable than @sursumCorda's ... and unfortunately less efficient

```CodeTools:-Usage(nequal5())  # sursumCorda's code

10642
memory used=105.92MiB, alloc change=0 bytes, cpu time=1.04s, real time=1.04s, gc time=38.36ms```
 > restart
 > f := proc(n)   local p    := (n-2)*24:   local lowb := M -> ceil((p-add(u[m], m=1..M))/(n-M));   local upb  := p-n+1:   local a    := 0:   local u, inc, m, start, i, kompt:   u := Vector(n-1, i -> ceil(p/n));   for m from 2 to n-1 do     u[m] := lowb(m-1);   end do;   kompt := 0:   while [entries(u, nolist)] <> [upb\$(n-1)] do     if p-add(u[m], m=1..n-1) >= 1 then       a := a+1;     end if;            start := n-1;       while u[start-1] = u[start] do         start := start-1:         if start = 1 then break end if:       end do:       u[start] := u[start]+1;       for m from start+1 to n-1 do         u[m] := lowb(m-1);       end do;     kompt := kompt+1:   end do:   a, kompt end proc:
 > f(3)
 (1)
 > CodeTools:-Usage(f(4))
 memory used=7.77MiB, alloc change=0 bytes, cpu time=104.00ms, real time=105.00ms, gc time=5.10ms
 (2)
 > CodeTools:-Usage(f(5))
 memory used=0.85GiB, alloc change=0 bytes, cpu time=11.64s, real time=11.70s, gc time=547.80ms
 (3)
 > nequal5:=proc()
 > local u,v,w,x,a_5;
 > a_5:=0;
 > for u from ceil((5-2)*24/(5-0)) to (5-2)*24-5+1 do
 > for v from ceil(((5-2)*24-u)/(5-1)) to u do
 > for w from ceil(((5-2)*24-u-v)/(5-2)) to v do
 > for x from ceil(((5-2)*24-u-v-w)/(5-3)) to w do
 > if (5-2)*24-u-v-w-x>=1 then
 > a_5:=a_5+1;
 > end if;
 > end do;
 > end do;
 > end do;
 > end do;
 > print(a_5);
 > end proc:
 >
 > CodeTools:-Usage(nequal5())
 >
 memory used=105.92MiB, alloc change=0 bytes, cpu time=1.04s, real time=1.04s, gc time=38.36ms
 >

## Assuming mm is a real matrix...

All the eigenvalues of mm must be positive for mm to be positive semi definite.
mm has 12 null eigenvalues:

```cp := CharacteristicPolynomial(mm, t):
collect(cp, t, LargeExpressions:-Veil[K]);

t^36-3*K[1]*t^35....-(9/2048)*K[23]*t^13

```

Factoring cp gives

`C * t^13 * Pol(t, 1) * (Pol[1](t, 2))^2 * (Pol[2](t, 2))^3 * (Pol(t, 4))^3 `

where Pol(t, d) represents some polynom of degree d and indeterminate t.
Then mm has 1+1+2+2+4=10 distinct eigenvalues

```fcp := factor(cp):
mul(op(1..2, fcp)) * collect(op(3, fcp), t, Veil[C]) * mul(seq(collect(op([i, 1], fcp), t, Veil[C])^op([i, 2], fcp), i=4..6))
2
1    13            /    2                    \
- ---- t   (t + C[9]) \-2 t  + 2 t C[1] - 3 C[2]/
2048

3
/   4      3           2                     \
\4 t  - 2 t  C[3] - 2 t  C[4] - t C[5] + C[6]/

3
/    2                \
\-2 t  - t C[7] + C[8]/
```

A minimum requirement for mm to be psd is

• the root of Pol(t, 1) >= 0
• the discriminants of Pol[1](t, 2) and Pol[3](t, 2) both are >= 0

Details are given in the attached file.
My conclusion, if I'm not mistaken, is that there is no triple of reals (k__1, k__2, k__3)  such than matrix mm is real positive semi definite.

 > restart:
 > with(LargeExpressions): with(LinearAlgebra):
 > mm:=<0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,1,0,0,0,0,0,0,0,0,0,0\
 > ,k__3-1/2,0,k__3-1/2,0,0,0,0,0,0,0,0,-1/2,0,1-2*k__3,0,-1/\
 > 2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
 > 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
 > 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,1,0,\
 > k__3,0,-1/2,0,0,0,0,0,0,0,0,0,k__3,0,-2*k__3,0,0,0,0,0,0,0\
 > ,0,-1/2,0,0,0,0,0|0,0,0,0,0,0,-2*k__3,0,k__3-1,0,-1/2,0,0,\
 > 0,0,0,0,0,0,0,k__3,0,1-2*k__3,0,0,0,0,0,0,0,0,k__1,0,0,0,0\
 > |0,0,0,0,0,k__3,0,1-4*k__3,0,k__3-1,0,0,0,0,0,0,0,0,0,(1-k\
 > __2)/2,0,k__3-1/2,0,0,0,0,0,0,0,0,1-2*k__3,0,0,0,0,0|0,0,0\
 > ,0,0,0,k__3-1,0,1-4*k__3,0,k__3,0,0,0,0,0,0,0,0,0,k__3-1/2\
 > ,0,(1-k__2)/2,0,0,0,0,0,0,0,0,1-2*k__3,0,0,0,0|0,0,0,0,0,-\
 > 1/2,0,k__3-1,0,-2*k__3,0,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3,0\
 > ,0,0,0,0,0,0,0,k__1,0,0,0,0,0|0,0,0,0,0,0,-1/2,0,k__3,0,1,\
 > 0,0,0,0,0,0,0,0,0,-2*k__3,0,k__3,0,0,0,0,0,0,0,0,-1/2,0,0,\
 > 0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
 > 0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,-2*k__3,0,k__3,0\
 > ,k__1,0,0,0,0,0,0,0,0,k__3-1,0,1-2*k__3,0,0,0,0,0,0,-1/2,0\
 > |0,0,k__3-1/2,0,0,0,0,0,0,0,0,0,0,k__2,0,-k__1+5*k__3-1/2,\
 > 0,0,0,0,0,0,0,0,k__3-1/2,0,-k__1+5*k__3-1/2,0,1-2*k__3,0,0\
 > ,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,k__3,0,1-2*k__3,0,k__3,\
 > 0,0,0,0,0,0,0,0,k__3-1/2,0,k__3-1/2,0,0,0,0,0,0,-2*k__3,0|\
 > 0,0,k__3-1/2,0,0,0,0,0,0,0,0,0,0,-k__1+5*k__3-1/2,0,k__2,0\
 > ,0,0,0,0,0,0,0,1-2*k__3,0,-k__1+5*k__3-1/2,0,k__3-1/2,0,0,\
 > 0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,k__1,0,k__3,0,-2*k__3,0,\
 > 0,0,0,0,0,0,0,1-2*k__3,0,k__3-1,0,0,0,0,0,0,-1/2,0|0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,k__3,0,(1-k__2)/2,0,1-2*k__\
 > 3,0,0,0,0,0,0,0,0,0,1-4*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,k_\
 > _3-1,0,0,0,0,0|0,0,0,0,0,0,k__3,0,k__3-1/2,0,-2*k__3,0,0,0\
 > ,0,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,k__3,0,0,\
 > 0,0|0,0,0,0,0,-2*k__3,0,k__3-1/2,0,k__3,0,0,0,0,0,0,0,0,0,\
 > k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0,0,k__3,0,0,0,0,0|0,0,0,0\
 > ,0,0,1-2*k__3,0,(1-k__2)/2,0,k__3,0,0,0,0,0,0,0,0,0,k__3-1\
 > /2,0,1-4*k__3,0,0,0,0,0,0,0,0,k__3-1,0,0,0,0|0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > |0,0,-1/2,0,0,0,0,0,0,0,0,0,0,k__3-1/2,0,1-2*k__3,0,0,0,0,\
 > 0,0,0,0,1,0,k__3-1/2,0,-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,\
 > 0,0,0,0,k__3-1,0,k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0,0,1-4*k\
 > __3,0,(1-k__2)/2,0,0,0,0,0,0,k__3,0|0,0,1-2*k__3,0,0,0,0,0\
 > ,0,0,0,0,0,-k__1+5*k__3-1/2,0,-k__1+5*k__3-1/2,0,0,0,0,0,0\
 > ,0,0,k__3-1/2,0,k__2,0,k__3-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,\
 > 0,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,k__3-1,0,0,0,0,0,0,0,0,(\
 > 1-k__2)/2,0,1-4*k__3,0,0,0,0,0,0,k__3,0|0,0,-1/2,0,0,0,0,0\
 > ,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,-1/2,0,k__3\
 > -1/2,0,1,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,-1/2,0,1-2*\
 > k__3,0,k__1,0,0,0,0,0,0,0,0,0,k__3-1,0,k__3,0,0,0,0,0,0,0,\
 > 0,-2*k__3,0,0,0,0,0|0,0,0,0,0,0,k__1,0,1-2*k__3,0,-1/2,0,0\
 > ,0,0,0,0,0,0,0,k__3,0,k__3-1,0,0,0,0,0,0,0,0,-2*k__3,0,0,0\
 > ,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\
 > ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,-1/\
 > 2,0,-2*k__3,0,-1/2,0,0,0,0,0,0,0,0,k__3,0,k__3,0,0,0,0,0,0\
 > ,1,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0>:
 > cp  := CharacteristicPolynomial(mm, t): fcp := factor(cp): vcp := mul(op(1..2, fcp)) * collect(op(3, fcp), t, Veil[C]) * mul(seq(collect(op([i, 1], fcp), t, Veil[C])^op([i, 2], fcp), i=4..6))
 (1)
 > # The roots of each of the 5 polynomials vcp contains are the eigenvalues of mm. # Assuming mm is a real matrix, mm is psd if all these rooots are >= 0 # # A non sufficient requirement for mm to be psd is cond := {   C[1] >= 0,   discrim(op([4, 1], vcp), t) >= 0,   discrim(op([5, 1], vcp), t) >= 0 }; print(): conditions := map(simplify, eval(cond, {seq(C[k] = Unveil[C](C[k]), k=1..LastUsed[C])})): print~(conditions):
 (2)
 > # Solve these conditions wrt k__1, k__2 and k__3 S := solve(conditions, {k__1, k__2, k__3})
 (3)
 > # Check if there exists at least one of these solutions such that the roots of # each 2nd degree polynomials are positive. # # A simpler way to do this is to check if the sum of the two roots of these # polynomials are positive. # Note that even of its so, this doesn't implie that both roots are positibe, # but the sum is not positive we are sure that at least one root is strictly # negative and thus that mm is not psd. # Polynomial -2*t^2+2*t*C[2]-3*C[3] solve(op([4, 1], vcp), t): eval(add(%), {seq(C[k] = Unveil[C](C[k]), k=1..LastUsed[C])}): map(s-> (is(% >= 0) assuming s[]), [S]); # Polynomial -2*t^2-t*C[4]+C[5] solve(op([5, 1], vcp), t): eval(add(%), {seq(C[k] = Unveil[C](C[k]), k=1..LastUsed[C])}): map(s-> (is(% >= 0) assuming s[]), [S])
 (4)
 > # Assuming mm is a real matrix, polynomial -2*t^2-t*C[4]+C[5] doesn't having both its # two roots positive under the constraints S defines, one concludes that there is no # triple (k__1, k__2, k__3) of reals such mm is positive semi definite

## Elementary...

The package Student:-Statistics contains the built-in function ProbabilityTable which is aimed to do this for any distribution:

```restart:
with(Student:-Statistics):
ProbabilityTable( 'Normal', output = embed ):
```

## Is it you...

who deleted your previous question (with an evident connection to this one)  AFTER I answered it?
If it is so that is quite unfair to say the least.

https://www.mapleprimes.com/questions/236727-Need-Help-On-Graph
dated from July 10?

Of course you did for your help_graph.mw is a just copy of  "my" code Evolution_1.mw without any mention of its author... quite unfair too.
Maybe something like "mmcdara provided this code but it doesn't work for this specific case" could have been a more honest way to ask your question. Don't you think?

I'm really offended by this kind of attitude and thought about ignoring you altogether.
However,"my" code doesn't work for this particular case and I thought it was my duty to improve it.

So here is the code.
Evolution_2.mw

By the way, you still haven't understood that the value of x must be stated: this wasthe case 4 days ago and it's still the case today (I took arbitrarily x=0).
In the same way, you do not specify what the tickmarks on the x-axis should be. So I took some initiative: if it doesn't please you modify my code the way you want.

## Nothing is wrong with your file...

What is wrong is either the reasoning that led you to write the command above,  or the way you translate this resonning into this command

`solve({M > 0, diff(A, x), x < 0}, x)[] assuming M > 0, x < 0`

Its very writing is nonsensical, and I suspect you wanted to do something and have clumsily tried to translate it into Maple.
So, what is it that you wanted to do

Help_mmcdara.mw

## ...

The first point is upload as an attachment (and insert a link) your mw file using the big green up-arrow in the menu bar.
(thanks to whoever corrected my sentence)

Unless, only people who have some time to waste, like me at this moment, will look at your code.

No need to go very far to see that your Kraziski procedure is badly written and that, once corrected, it produces an error for an obvious reason.
Here are some explanations.

 > restart:
 > with(StringTools):
 > with(LinearAlgebra):
 >
 > KasiskiTest := proc(str)
 > local i, j, n, trigrams, distances, gcds, freqs; #DEBUG():
 > n := StringTools:-Length(str);
 > trigrams := table();
 > distances := table();
 > gcds := table();
 > freqs := table();
 > for i from 1 to n-2 do
 > trigrams[str[i..i+2]] := [op(trigrams[str[i..i+2]]),i];
 > end do;
 > for j in trigrams do
 > if nops(trigrams[j]) > 1 then
 > distances[j] := map(diff,trigrams[j]);
 > gcds[j] := igcd(op(distances[j]));
 > freqs[gcds[j]] := freqs[gcds[j]] + 1;
 > end if;
 > end do;
 > #return maxloc(freqs);
 > return trigrams, distances, gcds, freqs;
 > end:
 >
 > IndexCoincidenceTest := proc(str)
 > local i, k, n, m, L, C, Convert, IC, avgIC;
 > Convert := table();  n := StringTools:-Length(str);
 > for i from 1 to n do
 > Convert[str[i]] := i-1;
 > end do;
 > L := Letter2Number(str,Convert);
 > IC := [];
 > avgIC := [];
 > for k from 1 to n/2 do
 > C := [];
 > for i from 1 to k do
 > C[i] := seq(L[i+j*k],j=0..floor((n-i)/k));
 > end do;
 > IC[k] := map(DotProduct,C,C)/nops(C)^2;
 > end do;
 > return maxloc(avgIC);
 > end:
 > cstr := "eqjyvxvgiuphmrrrgzhvwrrukmrtijfbtodidtsjsgofwmfalkeveolqlmhkfhrerhwuidejonvjuumklhliiwwnajxbonjthitzsqufklhihrvditvvvzhvwhihrvditvvvgsrrbunvqlmhkvhgdzpbmuvwvloiqdiiglhxtyewosksvgyklhecfrykyrqhgnzrjhukxkknwvrswyewosbrrcnfjnodumfuuchqutjysvojikomtesgbcirlcfrvzrlgwonxeqeowxkkmfvhgbjxuasvguepksjxaglvomrhhapdcponuewuntiwnakxkosnevufrwlspcivvetmhyslgknoniykrrwzuuchdvpveuzoklhirlhhonkioretxrltyivgicsugbjsoatvpbonjsoabcizotysxztyinky": trigrams, distances, gcds, freqs := KasiskiTest(cstr);
 (1)
 > eval(trigrams): numelems(trigrams); eval(distances); eval(gcds); eval(freqs);
 (2)
 > # You cannot explore a table this way for j in trigrams do
 > j
 > end do;
 (3)
 > # Do this instead for j in [indices(trigrams, nolist)] do   j end do: # Example count := 0: for j in [indices(trigrams, nolist)] do   if count=3 then     break   else     print(j):     count := count+1:   end if: end do:
 (4)
 > # You can also create the list of indices ot table trigrams trigrams_ind := [indices(trigrams, nolist)]: trigrams_ind[1..3]
 (5)
 > # You can also create the list of indices ot table trigrams numelems(trigrams_ind);
 > # How many times do you pass in the "if" test in procedure Kraziski? # # You see here that every entry of the table has at least operators, # which is normal given the way you construct them!!! #    trigrams[str[i..i+2]] := [op(trigrams[str[i..i+2]]),i]; #                              <----- at least 1 ----->  1 N := numelems(select((x -> nops(x) > 1), [entries(trigrams, nolist)])); # This "if" test is then useless.
 (6)
 > # What do you do within the Kraziski procedure? # The first operation is 'map(diff,trigrams[j])'; # Let's follow step by step what happens for j in trigrams_ind[1..2] do     j, trigrams[j];
 > #if nops(trigrams[j]) > 1 then
 > distances[j] := map(diff,trigrams[j]);
 > gcds[j] := igcd(op(distances[j]));
 > freqs[gcds[j]] := freqs[gcds[j]] + 1;
 > #end if;
 > end do;
 > # Look this help page help(diff)
 (7)
 > # I don't know what you want to do but you done it bad

Finally Kraziski returns maxloc(freqs) : even if freqs was correctly constructed, you never defined procedure maxloc, so the return of your code will always be

`                     maxloc(freqs)`

As I'm not a crytography specialist it's very unlikely, it is unlikely that I will pursue this with you.
But if you make the effort to provide the file containing your code, you can expect some people to help you.

## @somestudent Maximizing complexfunc...

Maximizing complexfunction with respect to f under the constraint f <= 1 is obvious: the maximizer is f=1 whatever the values of x and y.

Here is the solution for a Threshold given a numeric value.

 > restart:
 > with(Optimization): with(plots):
 > functionSimple := x*y;
 > complexfunction := -2.2*x^3 - 0.9*x*y^2 + x + 1.5*y^2 - 0.21 + f;
 > plotsimple := contourplot(functionSimple, x = -1 .. 1, y = -1 .. 1, view = [-1 .. 0, -1 .. 0], contours = [1, 0.75, 0.5, 0.25], numpoints = 100):
 > plotcomplex := contourplot(subs(f = 1, complexfunction), x = -1 .. 1, y = -1 .. 1, view = [-1 .. 0, -1 .. 0], contours = [-2, -1, 0, 1, 2], numpoints = 100): #display(`<|>`(plotsimple, plotcomplex))
 (1)
 > # Find the maximum of complexfunction wrt f, under the constraint f <= 1. # # Obviously complexfunction is maximum when f=1, whatever the values of x and y. # The maximum value of complexfunction is then max_complexfunction := eval(complexfunction, f=1);
 (2)

The "equality constraint" case

Case 1: real solutions do exist

 > # Maximize may fail to find a solution for a reason that will be explained later. # (same thing happens with an inequality constraint). # # When a solution is got there is no insurance it corresponds to the global maximum: Thres := 1;  # for instance Maximize(functionSimple, {max_complexfunction = Thres}, initialpoint = {x = 0, y = 0}); Maximize(functionSimple, {max_complexfunction = Thres}, initialpoint = {x = -0.5, y = -0.5});
 (3)
 > # A more direct way: form the Lagrangian and find the zeros of its derivatives to get # (x, y) points where functionSimple is extremal. # # This can be done using solve if Thres is given a value. # If not, solve returns a solutionof RootOf(pol(_Z)) form where pol(_Z) is a # 8th degree polunomial in _Z. # So it's impossible to get explicit expressions of the eight solutions. L    := functionSimple - lambda*(max_complexfunction - Thres): dL   := diff~(L, [x, y, lambda]): HL   := Student:-VectorCalculus:-Hessian(L, [x, y, lambda]);
 (4)
 > # Critical points are: CriticalPoints := solve(dL, [x, y, lambda])
 (5)
 > # REMARK: if you choose another value of Thres, for instance -1, produces only #         complex solutions. # # Select the Real solutions alone (note this list cotains the two solutions Maximize # found when I ran it from tho different initial points). real_CriticalPoints := remove(has, CriticalPoints, I)
 (6)
 > # CriticalPoints are maximizers if and only HL has a positive determinant and # a negative trace at thes points: use LinearAlgebra in   Maximizers := NULL:   for c in real_CriticalPoints do     det := Determinant(eval(HL, c)):     tra := Trace(eval(HL, c)):     if `and`(det > 0, tra < 0) then Maximizers := Maximizers, c: end if:   end do: end use: Maximizers := [Maximizers]
 (7)
 > # Extrema of functionSimple under the constraint max_complexfunction = Thres Extrema := map(sol -> eval([functionSimple, (max_complexfunction - Thres)], sol), Maximizers);
 (8)
 > # Largest value of these extrema WhichOne := sort(map2(op, 1, Extrema), `>`, output=permutation)[1]; Maximizer := remove(has, Maximizers[WhichOne], lambda): OptimalSolution := [Extrema[WhichOne][1], Maximizer]
 (9)

The "equality constraint" case

Case 2: no real solutions

 > randomize(168916419713978); r := rand(-1. .. 1.): Thres := -1;  # for instance Maximize(functionSimple, {max_complexfunction = Thres}, initialpoint = {x = r(), y = r()});
 > # This kind of error message is always disturbing. # # Here is the reason why we get it: L    := functionSimple - lambda*(max_complexfunction - Thres): dL   := diff~(L, [x, y, lambda]): CriticalPoints := solve(dL, [x, y, lambda]); real_CriticalPoints := remove(has, CriticalPoints, I)
 (10)

The "inequality constraint" case

reference: https://users.wpi.edu/~pwdavis/Courses/MA1024B10/1024_Lagrange_multipliers.pdf

 > # Now we want to maximise functionSimple under the constraint max_complexfunction <= Thres # # The key is to use an Augmented Lagrangian: Thres := 1;  # for instance # ALE stands for Augmented Lagrangian Equations: Constraint := max_complexfunction - Thres - a^2: dC := `<,>`(diff~(Constraint     , [x, y, a])): df := `<,>`(diff~(functionSimple , [x, y, a])): ALE := entries(df =~ lambda*dC, nolist), Constraint=0: printf("\nAugmented Lagrangian Equations:\n"): print~([ALE]):
 Augmented Lagrangian Equations:
 (11)
 > CriticalPoints := solve({ALE}, [x, y, lambda, a]): real_CriticalPoints := remove(has, CriticalPoints, I)
 (12)
 > # Interpretation: # The constraint is either active (here max_complexfunction = Thres) # or inactive (meaning max_complexfunction < Thres). # If the constraint is inactive, then lambda = 0  and a <> 0 in the third ALE. # If the constraint is active,   then lambda <> 0 and a = 0  in the third ALE. # Inactive := select~(has, real_CriticalPoints, a<>0);
 (13)
 > # Thus there is no solution found which corresponds to an inactive constraint. # This mean there is no need to write the inequality constraint #    max_complexfunction <= Thres # for writting #    max_complexfunction = Thres # is enough. # # It remains to select the maximizers amid real_CriticalPoints, but this has already # been done before.
 >

Note that tin the last case Optimization:-Maximize finds the "true" maximizer if the initial point is correctly chosen:

```Maximize(functionSimple, {max_complexfunction <= Thres}, initialpoint = {x = -0.5, y = -0.5});
[0.2452225..., [x = -0.5697256..., y = -0.4304221...]]
```

## Is it something like this that you want?...

(See my previous replies)

Note: the last BarGraph is not correctly disply on this site.

 >
 >
 >
 > ind := [indets(Nu, name)[]]
 (1)
 > nu := unapply(Nu, ind):
 > B := [0.1, 0.3, 0.5, 0.7, 0.9];
 (2)
 > # Personal opinion: ColumnGraph offers little capabilities to place the bars # in smart locations: ColumnGraph(nu~(-0.388e-1, B, 1.3), axis[1]=[tickmarks=[seq(k=0.1+k*0.2, k=0..4)]], view=[-1/2..5, default])
 > # I prefer doing this (see that bars are centered at the values of B) r := (x, y, h, c) -> plottools:-rectangle([x-h, 0], [x+h, y], color=c): plots:-display(   seq(     r(B[i], nu(-0.388e-1, B[i], 1.3), 0.05, "Gold")     , i=1..numelems(B)   )   , labels=[typeset('beta'), "Nu"]   , axis[1]=[tickmarks=B]   , view=[0..1, default] )
 > r := (x, y, h, c) -> plottools:-rectangle([x-h, 0], [x+h, y], color=c): h := 0.025: plots:-display(   seq(     r(B[i]-h, nu(-0.388e-1, B[i], 1.3), h, "Gold")     , i=1..numelems(B)   )   , seq(       r(B[i]+h, nu(-1e-1, B[i], 1.3), h, "Brown")       , i=1..numelems(B)   )   , plot([[B[1]-h, 0], [B[1]-h, nu(-0.388e-1, B[1], 1.3)/2]], color="Gold", legend=typeset(Q=-0.388e-1))   , plot([[B[1]+h, 0], [B[1]+h, nu(-1e-1, B[1], 1.3)/2]], color="Brown", legend=typeset(Q=-1e-1))   , labels=[typeset('beta'), "Nu"]   , axis[1]=[tickmarks=B]   , view=[0..1, default] )
 > bar := (x, y, t, cs, l) -> plot([[x, 0], [x, y]], thickness=t, colorscheme=cs): #, legend=l): h := 0.03: plots:-display(   seq(     bar(B[i]-h, nu(-0.388e-1, B[i], 1.3), 20, ["Red", "Yellow"]) #, typeset(Q=-0.388e-1))     , i=1..numelems(B)   )   , seq(       bar(B[i]+h, nu(-1e-1, B[i], 1.3), 20, ["Purple", "Cyan"]) #, typeset(Q=-1e-1))       , i=1..numelems(B)   )   , plots:-textplot([B[1]-h, nu(-0.388e-1, B[1], 1.3), "A"], align=above)   , plots:-textplot([B[1]+h, nu(-1e-1, B[1], 1.3), "B"], align=above)   , labels=[typeset('beta'), "Nu"]   , axis[1]=[tickmarks=B]   , view=[0..1, default]   , title="A: Q=-0.388e-1 / B: Q=-0.1e-1" )

## To complete Tom's answer:  ...

 > restart
 > Z := phi__1 = arctan(a*sin(alpha__1)/(a*cos(alpha__1) - b)):
 > sol := solve(Z, alpha__1, useassumptions) assuming (0 < a and 0 < b)
 (1)
 > # You have 2 solutions: numelems([sol])
 (2)
 > # Each one is of the form arctan(y, x) where (see help(arctan)): arctan(y, x) = -I*ln((x+I*y)/sqrt(x^2+y^2))
 (3)
 > # check it: renaming := [op(1..2, sol[1])] =~ [y, x]; eval(sol[1], renaming); lnsol := convert(%, ln);
 (4)
 >