mmcdara

5989 Reputation

17 Badges

7 years, 266 days

MaplePrimes Activity


These are answers submitted by mmcdara

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
This will answer your question, but it won't help you understand how the equations were constructed.

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


here is an extremely simple way do draw random triangles
 

restart

with(plots):
with(plottools):
with(Statistics):

interface(version);

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

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

169245464804467

 

 

 

 

Download RandomTriangles.mw

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

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

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

169246137904467

 

 

 

 

Download RandomTriangle_2.mw

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|
                    /


 

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

x^2/a^2+y^2/b^2-1

 

-(a^2-b^2)^(1/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)

{X__P(0) = x__0, X__T(0) = -(a^2-b^2)^(1/2), Y__P(0) = b*(1-x__0^2/a^2)^(1/2), diff(X__P(t), t) = V__P*(X__P(t)-X__T(t))/((X__P(t)-X__T(t))^2+Y__P(t)^2)^(1/2), diff(X__T(t), t) = V__T, diff(Y__P(t), t) = V__P*Y__P(t)/((X__P(t)-X__T(t))^2+Y__P(t)^2)^(1/2)}

 

proc (x_rkf45) 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_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [V__T = V__T, a = a, b = b, x__0 = x__0, V__P = V__P]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([Array(1..2, 1..21, {(1, 1) = 1.0, (1, 2) = .0, (1, 3) = 1.0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = 1.0, (1, 8) = undefined, (1, 9) = undefined, (1, 10) = undefined, (1, 11) = undefined, (1, 12) = undefined, (1, 13) = undefined, (1, 14) = undefined, (1, 15) = undefined, (1, 16) = undefined, (1, 17) = undefined, (1, 18) = undefined, (1, 19) = undefined, (1, 20) = undefined, (1, 21) = undefined, (2, 1) = 1.0, (2, 2) = .0, (2, 3) = 100.0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = undefined, (2, 9) = undefined, (2, 10) = 0.10e-6, (2, 11) = undefined, (2, 12) = .0, (2, 13) = undefined, (2, 14) = .0, (2, 15) = .0, (2, 16) = undefined, (2, 17) = undefined, (2, 18) = undefined, (2, 19) = undefined, (2, 20) = undefined, (2, 21) = undefined}, datatype = float[8], order = C_order), proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[3]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, Array(1..1, 1..2, {(1, 1) = undefined, (1, 2) = undefined}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 3, (2) = 3, (3) = 0, (4) = 0, (5) = 5, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 1, (17) = 0, (18) = 0, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..8, {(1) = x__0, (2) = -1.*(a^2-1.*b^2)^(1/2), (3) = b*(1.-1.*x__0^2/a^2)^(1/2), (4) = Float(undefined), (5) = Float(undefined), (6) = Float(undefined), (7) = Float(undefined), (8) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..3, {(1) = .1, (2) = .1, (3) = .1}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = 0, (2) = 0, (3) = 0}, datatype = integer[8]), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = X__P(t), Y[2] = X__T(t), Y[3] = Y__P(t)]`; if (Y[1]-Y[2])^2+Y[3]^2 < 0 then YP[1] := undefined; return 0 end if; YP[1] := Y[8]*(Y[1]-Y[2])*evalf(1/((Y[1]-Y[2])^2+Y[3]^2)^(1/2)); YP[3] := Y[8]*Y[3]*evalf(1/((Y[1]-Y[2])^2+Y[3]^2)^(1/2)); YP[2] := Y[4]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[3]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = X__P(t), Y[2] = X__T(t), Y[3] = Y__P(t)]`; if (Y[1]-Y[2])^2+Y[3]^2 < 0 then YP[1] := undefined; return 0 end if; YP[1] := Y[8]*(Y[1]-Y[2])*evalf(1/((Y[1]-Y[2])^2+Y[3]^2)^(1/2)); YP[3] := Y[8]*Y[3]*evalf(1/((Y[1]-Y[2])^2+Y[3]^2)^(1/2)); YP[2] := Y[4]; 0 end proc, -1, 0, 0, proc (t, Y, Ypre, n, EA) EA[1, 7+2*n] := Y[3]; EA[1, 8+2*n] := 1; 0 end proc, proc (e, t, Y, Ypre) return 0 end proc, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..8, {(1) = 0., (2) = x__0, (3) = -1.*(a^2-1.*b^2)^(1/2), (4) = b*(1.-1.*x__0^2/a^2)^(1/2), (5) = undefined, (6) = undefined, (7) = undefined, (8) = undefined}); _vmap := array( 1 .. 3, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) 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 select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, 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, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, X__P(t), X__T(t), Y__P(t)], (4) = [V__T = V__T, a = a, b = b, x__0 = x__0, V__P = V__P]}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["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_rkf45, '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_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, '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_rkf45), 'string') = rhs(x_rkf45); 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_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

 

[V__T = 1., a = 2., b = 1., x__0 = 0., V__P = -1.15470053837925]

 

[t = HFloat(1.7320505207668953), X__P(t) = HFloat(-2.890299205020301e-7), X__T(t) = HFloat(-2.8680198421090967e-7), Y__P(t) = HFloat(-8.020543449867425e-16)]

 

HFloat(1.7320505207668953)

 

 

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)

 

 

Download Capture.mw

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

"The integration domain [a, b] = ", [0, 2]

(1)

f := exp(x):
print("The given function is ", f);

"The given function is ", exp(x)

(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.389056099

 

"The exact integration in", [0, 2], "is", 6.389056099

 

"The exact integration in [0, 2] is 6.389056099"

 

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 .50000"

 

"The value of h to divide the domain [0, 2] into 4 subintervals is 1/2"

 

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

"Numerical integration in [a,b] is going to perform when h via RECTANGULAR METHOD for n = ", 4

(3)

 

Download Riemann_sum.mw

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.

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

proc (zeta) options operator, arrow; (r[1]*exp(s[1]*zeta)+r[2]*exp(s[2]*zeta))/(r[3]*exp(s[3]*zeta)+r[4]*exp(s[4]*zeta)) end proc

 

proc (zeta) options operator, arrow; ((r[3]*(s[1]-s[3])*exp(s[3]*zeta)+r[4]*exp(s[4]*zeta)*(s[1]-s[4]))*r[1]*exp(s[1]*zeta)+r[2]*(r[3]*(s[2]-s[3])*exp(s[3]*zeta)+r[4]*exp(s[4]*zeta)*(s[2]-s[4]))*exp(s[2]*zeta))/((r[3]*exp(s[3]*zeta)+r[4]*exp(s[4]*zeta))*(r[1]*exp(s[1]*zeta)+r[2]*exp(s[2]*zeta))) end proc

(1)

# relation (10)

u := (zeta, m) ->  a[0]+add(a[k]*psi(zeta)^k, k=1..m)
                   +add(b[k]*psi(zeta)^(-k), k=1..m)
                   +add(c[k]*dpsi(zeta)^k, k=1..m);

proc (zeta, m) options operator, arrow; a[0]+add(a[k]*psi(zeta)^k, k = 1 .. m)+add(b[k]*psi(zeta)^(-k), k = 1 .. m)+add(c[k]*dpsi(zeta)^k, k = 1 .. m) end proc

(2)

# U... (m=1)
# (xi instead of zeta to be consistent with uour notations)

U := unapply(u(xi, 1), xi);

proc (xi) options operator, arrow; a[0]+a[1]*(r[1]*exp(s[1]*xi)+r[2]*exp(s[2]*xi))/(r[3]*exp(s[3]*xi)+r[4]*exp(s[4]*xi))+b[1]*(r[3]*exp(s[3]*xi)+r[4]*exp(s[4]*xi))/(r[1]*exp(s[1]*xi)+r[2]*exp(s[2]*xi))+c[1]*((r[3]*(s[1]-s[3])*exp(s[3]*xi)+r[4]*exp(s[4]*xi)*(s[1]-s[4]))*r[1]*exp(s[1]*xi)+r[2]*(r[3]*(s[2]-s[3])*exp(s[3]*xi)+r[4]*exp(s[4]*xi)*(s[2]-s[4]))*exp(s[2]*xi))/((r[3]*exp(s[3]*xi)+r[4]*exp(s[4]*xi))*(r[1]*exp(s[1]*xi)+r[2]*exp(s[2]*xi))) end proc

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

0

 

 

a[1]*(r[1]*s[1]*exp(s[1]*xi[n])+r[2]*s[2]*exp(s[2]*xi[n]))/(r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))-a[1]*(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n]))*(r[3]*s[3]*exp(s[3]*xi[n])+r[4]*s[4]*exp(s[4]*xi[n]))/(r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))^2-b[1]*(r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))*(r[1]*s[1]*exp(s[1]*xi[n])+r[2]*s[2]*exp(s[2]*xi[n]))/(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n]))^2+b[1]*(r[3]*s[3]*exp(s[3]*xi[n])+r[4]*s[4]*exp(s[4]*xi[n]))/(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n]))+c[1]*((r[3]*(s[1]-s[3])*s[3]*exp(s[3]*xi[n])+r[4]*s[4]*exp(s[4]*xi[n])*(s[1]-s[4]))*r[1]*exp(s[1]*xi[n])+(r[3]*(s[1]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[1]-s[4]))*r[1]*s[1]*exp(s[1]*xi[n])+r[2]*(r[3]*(s[2]-s[3])*s[3]*exp(s[3]*xi[n])+r[4]*s[4]*exp(s[4]*xi[n])*(s[2]-s[4]))*exp(s[2]*xi[n])+r[2]*(r[3]*(s[2]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[2]-s[4]))*s[2]*exp(s[2]*xi[n]))/((r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))*(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n])))-c[1]*((r[3]*(s[1]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[1]-s[4]))*r[1]*exp(s[1]*xi[n])+r[2]*(r[3]*(s[2]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[2]-s[4]))*exp(s[2]*xi[n]))*(r[3]*s[3]*exp(s[3]*xi[n])+r[4]*s[4]*exp(s[4]*xi[n]))/((r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))^2*(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n])))-c[1]*((r[3]*(s[1]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[1]-s[4]))*r[1]*exp(s[1]*xi[n])+r[2]*(r[3]*(s[2]-s[3])*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n])*(s[2]-s[4]))*exp(s[2]*xi[n]))*(r[1]*s[1]*exp(s[1]*xi[n])+r[2]*s[2]*exp(s[2]*xi[n]))/((r[3]*exp(s[3]*xi[n])+r[4]*exp(s[4]*xi[n]))*(r[1]*exp(s[1]*xi[n])+r[2]*exp(s[2]*xi[n]))^2)

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

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]))(1-V(xi[n])^2)^2

(5)

# Rewritting Rule

rr := (f, i) -> f(xi[n+i]) = f(xi[n]+i*d)

proc (f, i) options operator, arrow; f(xi[n+i]) = f(xi[n]+i*d) end proc

(6)

# Rewrite eqV using rr

eqVrr := eval(eqV, [rr(V, -1), rr(V, +1)]);

c*(diff(V(xi[n]), xi[n]))*(V(xi[n])+V(xi[n]-d))*(V(xi[n])+V(xi[n]+d))-2*(V(xi[n]-d)-V(xi[n]+d))*(V(xi[n]))(1-V(xi[n])^2)^2

(7)

# To get "your" equation jus do this
# Uncomment tosee this lengthy expression

eq := eval(eqVrr, V=U):
length(eq);

61630

(8)

 

Download abs_mmcdara.mw

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)

48, 168

(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

 

816, 17853

(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

 

10642, 1999573

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

 

10642

 

memory used=105.92MiB, alloc change=0 bytes, cpu time=1.04s, real time=1.04s, gc time=38.36ms

 

 

Download nequal_mmcdara.mw

 

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/2048)*t^13*(t+C[1])*(-2*t^2+2*t*C[2]-3*C[3])^2*(-2*t^2-t*C[4]+C[5])^3*(4*t^4-2*t^3*C[6]-2*t^2*C[7]-t*C[8]+C[9])^3

(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):

{0 <= C[1], 0 <= 4*C[2]^2-24*C[3], 0 <= C[4]^2+8*C[5]}

 

 

0 <= 1-10*k__3+2*k__1-k__2

 

0 <= 4*k__1^2+4*k__1*k__2-16*k__1*k__3+k__2^2-8*k__2*k__3+160*k__3^2+4*k__1+2*k__2-200*k__3+65

 

0 <= 4*k__1^2+8*k__1*k__2-40*k__1*k__3+4*k__2^2-40*k__2*k__3+244*k__3^2-8*k__1-8*k__2-104*k__3+40

(2)

# Solve these conditions wrt k__1, k__2 and k__3

S := solve(conditions, {k__1, k__2, k__3})

{k__1 = (1/2)*k__2+5*k__3-1/2, k__2 = k__2, k__3 < 1/2}, {k__2 = k__2, k__3 < 1/2, (1/2)*k__2+5*k__3-1/2 < k__1}, {k__1 = (1/2)*k__2+2, k__3 = 1/2, k__2 < 1}, {k__3 = 1/2, k__1 < -k__2+7/2, k__2 < 1, (1/2)*k__2+2 < k__1}, {k__1 = -k__2+7/2, k__3 = 1/2, k__2 < 1}, {k__3 = 1/2, k__2 < 1, -k__2+7/2 < k__1}, {k__2 = 1, k__3 = 1/2, 5/2 <= k__1}, {k__1 = (1/2)*k__2+2, k__3 = 1/2, 1 < k__2}, {k__3 = 1/2, 1 < k__2, (1/2)*k__2+2 < k__1}, {k__1 = (1/2)*k__2+5*k__3-1/2, k__2 = k__2, 1/2 < k__3, k__3 < 7/12}, {k__2 = k__2, 1/2 < k__3, k__3 < 7/12, (1/2)*k__2+5*k__3-1/2 < k__1}, {k__1 = (1/2)*k__2+29/12, k__3 = 7/12, k__2 < 13/2}, {k__3 = 7/12, k__2 < 13/2, (1/2)*k__2+29/12 < k__1}, {k__2 = 13/2, k__3 = 7/12, 17/3 <= k__1}, {k__1 = (1/2)*k__2+29/12, k__3 = 7/12, 13/2 < k__2}, {k__3 = 7/12, 13/2 < k__2, (1/2)*k__2+29/12 < k__1}, {k__1 = (1/2)*k__2+5*k__3-1/2, k__2 = k__2, 7/12 < k__3, k__3 < 2/3}, {k__2 = k__2, 7/12 < k__3, k__3 < 2/3, (1/2)*k__2+5*k__3-1/2 < k__1}, {k__1 = (1/2)*k__2+17/6, k__3 = 2/3, k__2 < -2}, {k__3 = 2/3, k__1 < -(1/2)*k__2+5/6, k__2 < -2, (1/2)*k__2+17/6 < k__1}, {k__1 = -(1/2)*k__2+5/6, k__3 = 2/3, k__2 < -2}, {k__3 = 2/3, k__2 < -2, -(1/2)*k__2+5/6 < k__1}, {k__2 = -2, k__3 = 2/3, 11/6 <= k__1}, {k__1 = (1/2)*k__2+17/6, k__3 = 2/3, -2 < k__2}, {k__3 = 2/3, -2 < k__2, (1/2)*k__2+17/6 < k__1}, {k__1 = (1/2)*k__2+5*k__3-1/2, k__2 = k__2, 2/3 < k__3}, {k__2 = k__2, 2/3 < k__3, (1/2)*k__2+5*k__3-1/2 < k__1}

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

[false, false, false, false, true, true, true, true, true, false, false, false, false, true, true, true, false, false, false, false, false, false, false, false, false, false, false]

 

[false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false]

(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


 

Download psd.mw

 

@jalal  @Carl Love 

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

 



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.

Concerning this new thread: have you read this previous thread of yours 
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.

No need to reply to this answer, I'm done with you



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

Read carefully my comments in the atatched file
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;

    avgIC[k] := add(IC[k])/k;

  end do;

  return maxloc(avgIC);

end:

cstr := "eqjyvxvgiuphmrrrgzhvwrrukmrtijfbtodidtsjsgofwmfalkeveolqlmhkfhrerhwuidejonvjuumklhliiwwnajxbonjthitzsqufklhihrvditvvvzhvwhihrvditvvvgsrrbunvqlmhkvhgdzpbmuvwvloiqdiiglhxtyewosksvgyklhecfrykyrqhgnzrjhukxkknwvrswyewosbrrcnfjnodumfuuchqutjysvojikomtesgbcirlcfrvzrlgwonxeqeowxkkmfvhgbjxuasvguepksjxaglvomrhhapdcponuewuntiwnakxkosnevufrwlspcivvetmhyslgknoniykrrwzuuchdvpveuzoklhirlhhonkioretxrltyivgicsugbjsoatvpbonjsoabcizotysxztyinky":

trigrams, distances, gcds, freqs := KasiskiTest(cstr);

trigrams, distances, gcds, freqs

(1)

eval(trigrams):
numelems(trigrams);

eval(distances);
eval(gcds);
eval(freqs);

389

 

table( [ ] )

 

table( [ ] )

 

table( [ ] )

(2)

# You cannot explore a table this way

for j in trigrams do

  j

end do;

trigrams

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

"lvo"

 

"chd"

 

"hxt"

(4)

# You can also create the list of indices ot table trigrams

trigrams_ind := [indices(trigrams, nolist)]:
trigrams_ind[1..3]

["lvo", "chd", "hxt"]

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

389

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

map(diff, trigrams[j])

 

"lvo", ["lvo", 296]

 

Error, invalid input: diff expects 2 or more arguments, but received 1

 

# Look this help page

help(diff)

table( [ ] )

(7)

# I don't know what you want to do but you done it bad

Download Vigenere_corrrection_1.mw

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

x*y

 

-2.2*x^3-.9*x*y^2+x+1.5*y^2-.21+f

(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.2*x^3-.9*x*y^2+x+1.5*y^2+.79

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

1

 

[0.334075643949937884e-1, [x = HFloat(0.15608911254987048), y = HFloat(0.21402879322745888)]]

 

[.245222517633285952, [x = HFloat(-0.5697256179537332), y = HFloat(-0.43042213638565957)]]

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

HL := Matrix(3, 3, {(1, 1) = 13.2*lambda*x, (1, 2) = 1+1.8*lambda*y, (1, 3) = 6.6*x^2+.9*y^2-1, (2, 1) = 1+1.8*lambda*y, (2, 2) = -lambda*(-1.8*x+3.0), (2, 3) = 1.8*x*y-3.0*y, (3, 1) = 6.6*x^2+.9*y^2-1, (3, 2) = 1.8*x*y-3.0*y, (3, 3) = 0})

(4)

# Critical points are:


CriticalPoints := solve(dL, [x, y, lambda])

[[x = .1560891157, y = .2140287889, lambda = .2682161318], [x = .1560891157, y = -.2140287889, lambda = -.2682161318], [x = -.5697256180, y = -.4304221364, lambda = .3288142715], [x = -.5697256180, y = .4304221364, lambda = -.3288142715], [x = 2.063482758, y = 6.995324742*I, lambda = .4129820573*I], [x = 2.063482758, y = -6.995324742*I, lambda = -.4129820573*I], [x = .4334870771, y = -.1997457294*I, lambda = .9776869517*I], [x = .4334870771, y = .1997457294*I, lambda = -.9776869517*I]]

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

[[x = .1560891157, y = .2140287889, lambda = .2682161318], [x = .1560891157, y = -.2140287889, lambda = -.2682161318], [x = -.5697256180, y = -.4304221364, lambda = .3288142715], [x = -.5697256180, y = .4304221364, lambda = -.3288142715]]

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

[[x = .1560891157, y = .2140287889, lambda = .2682161318], [x = -.5697256180, y = -.4304221364, lambda = .3288142715]]

(7)

# Extrema of functionSimple under the constraint max_complexfunction = Thres


Extrema := map(sol -> eval([functionSimple, (max_complexfunction - Thres)], sol), Maximizers);

[[0.3340756439e-1, 0.], [.2452225177, 0.]]

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

2

 

[.2452225177, [x = -.5697256180, y = -.4304221364]]

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

168916419713978

 

-1

 

Error, (in Optimization:-NLPSolve) no improved point could be found

 

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

[[x = -.4115041637+.4588697259*I, y = -.1349702549+.6998130562*I, lambda = .1639706105+.1551680386*I], [x = -.4115041637+.4588697259*I, y = .1349702549-.6998130562*I, lambda = -.1639706105-.1551680386*I], [x = 2.024950423, y = 6.694668016*I, lambda = .4690139123*I], [x = 2.024950423, y = -6.694668016*I, lambda = -.4690139123*I], [x = .8813912375, y = -1.283915573*I, lambda = .4856660528*I], [x = .8813912375, y = 1.283915573*I, lambda = -.4856660528*I], [x = -.4115041637-.4588697259*I, y = -.1349702549-.6998130562*I, lambda = .1639706105-.1551680386*I], [x = -.4115041637-.4588697259*I, y = .1349702549+.6998130562*I, lambda = -.1639706105+.1551680386*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:

 

y = lambda*(-6.6*x^2-.9*y^2+1)

 

x = lambda*(-1.8*x*y+3.0*y)

 

0 = -2*lambda*a

 

-a^2-.21-2.2*x^3-.9*x*y^2+x+1.5*y^2 = 0

(11)

CriticalPoints := solve({ALE}, [x, y, lambda, a]):

real_CriticalPoints := remove(has, CriticalPoints, I)

[[x = .1560891157, y = .2140287889, lambda = .2682161318, a = 0.], [x = .1560891157, y = -.2140287889, lambda = -.2682161318, a = 0.], [x = -.5697256180, y = -.4304221364, lambda = .3288142715, a = 0.], [x = -.5697256180, y = .4304221364, lambda = -.3288142715, a = 0.]]

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

 

Download Solution_mmcdara.m

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


(See my previous replies)

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

 

restart

with(Statistics):

Nu := .4345625000-.3000000000*Q-.7500000000*Q^2+(.3750000000*(.7+Q))*Q+(0.7500000000e-1*(.3500000000+Q))*(1.216666667+Q)+2.603571429*beta*(0.8065843622e-1*lambda+.7471435780*Q-(1/243)*(12.13333333+28*Q)*lambda-28*Q*lambda*(1/243)+6.3*Q*(-0.5002083333e-2-.1179629630*Q-.4282098765*Q^2-.4814814815*Q^3+Q^4)-2.1*Q^3*(.4083333333+1.444444444*Q+Q^2)-3*Q^2*(-0.7503125000e-1-.7077777778*Q-1.050648148*Q^2-.5777777778*Q^3+Q^4)+(1.837500000+1.733333333*Q+Q^2)*Q^4-0.6769121954e-1-2*Q^6+.155555556*Q^5-9.205000000*Q^4+9.041897120*Q^3-2.676934490*Q^2):

ind := [indets(Nu, name)[]]

[Q, beta, lambda]

(1)

nu := unapply(Nu, ind):

B := [0.1, 0.3, 0.5, 0.7, 0.9];

[.1, .3, .5, .7, .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"
)


Download BarGraph_mmcdara.mw

 

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)

arctan(tan(phi__1)*((tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1)-b), (tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1)), arctan(tan(phi__1)*(-(-tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1)-b), -(-tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1))

(1)

# You have 2 solutions:

numelems([sol])

2

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

arctan(y, x) = -I*ln((x+I*y)/(x^2+y^2)^(1/2))

(3)

# check it:

renaming := [op(1..2, sol[1])] =~ [y, x];

eval(sol[1], renaming);
lnsol := convert(%, ln);

[tan(phi__1)*((tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1)-b) = y, (tan(phi__1)^2*b+(tan(phi__1)^2*a^2-b^2*tan(phi__1)^2+a^2)^(1/2))/(tan(phi__1)^2+1) = x]

 

arctan(y, x)

 

-I*ln((x+I*y)/(x^2+y^2)^(1/2))

(4)

 


 

Download Explanation.mw

4 5 6 7 8 9 10 Last Page 6 of 51