Maple 2021 Questions and Posts

These are Posts and Questions associated with the product, Maple 2021

Hi everyone 

I solve several differential equations using dsolve. Now I want to change one initial condition a little bit: 

before: u0n := u0*ms/(ms + mn)

new: u0n := u0*ms/(ms + m(t)

When solving the new system I am getting the error message: Error, (in dsolve/numeric/process_input) unknown piecewise(0 < 16-(xn(t)^2+(zn(t)-16)^2)^(1/2), 16-(xn(t)^2+(zn(t)-16)^2)^(1/2), 0) present in ODE system is not a specified dependent variable or evaluatable procedure

I am guessing this has to do with a function in an initial condition. 

Has anyone an idea how to solve this problem?

I added the maple file for the whole overview.

Thanks in advance!

How to improve this program ? Thank you.

restart;
Equation de la conique
eqcon := (45 - 27*cos(alpha))*x^2 - 54*sin(alpha)*x*y + (45 + 27*cos(alpha))*y^2 - 8;
Delta := (-54*sin(alpha))^2 - 4*(45 - 27*cos(alpha))*(45 + 27*cos(alpha));
expand(%);
simplify(%);
Discriminant : Δ<0 ce qui correxpond à une ellipse
Eq := simplify(expand(subs(x = cos(alpha/2)*X - sin(alpha/2)*Y, y = sin(alpha/2)*X + cos(alpha/2)*Y, eqcon)));
kx := coeff(Eq, X, 2);
ky := coeff(Eq, Y, 2);
k := -tcoeff(Eq);

EQ := X^2/(sqrt(1/kx^2)*k) + Y^2/(sqrt(1/ky^2)*k) = 1;
Calcul du grand et du petit axe 
a := 1/sqrt(coeff(lhs(EQ), X, 2));
b := 1/sqrt(coeff(lhs(EQ), Y, 2));
print(X^2/('a^2') + Y^2/('b^2') = 1);
 

!Edits: I have found  a existing polynomial algorithm, but I still have difficulty implementing it. 2022/10/26

 

The edge connectivity of a connected graph G  is the minimum number of edges whose deletion from the graph G disconnects G. Below we are concerned with a particular kind of edge-cut.

  • For a connected graph G=(V ,E), an edge set S ⊂ E is a restricted-edge-cut, if G−S is disconnected and every connected component of  G−S has at least 2 vertices. 

Clearly, a restricted-edge-cut is an edge cut with a special requirement.

  • The restricted-edge-connectivity of G, denoted by κres(G), is defined as the cardinality of a minimum restricted-edge-cut.

For example, a graph g is as follows.

 

Clearly, its edge-connectivity is 1 since (0,3) or (0,4) is a edge cut of g. But we can find that  if we remove the edge (0,3), then "3" is a isolated vertex. Similarly, "4" is a isolated vertex if we remove  (0,4). It is not difficult to find g has exactly the two cut-edges (0,3) and (0,4).

 

Based on the definition of the restricted-edge-cuts, neither {(0,3)} nor  {(0,4)} are restricted-edge-cuts A minimum restricted-edge-cut is {(0,1),(0,2)} since every connected component of G-{(0,1),(0,2)} has  (at least) 2 vertices.

So  κres(g) is 2.  My problem is:

Given a graph G, how to calculate the restricted-edge-connectivity of  G?  Moreover, how to find a minimum restricted-edge-cut? 

A specific graph that I want to test is as follows. (it has 16 vertices and 56 edges.) I would like to calculate its restricted-edge-connectivity and find a minimum restricted-edge-cut. 

g:=ConvertGraph("O~tIID@wL~j`PbOqgLJ@p");
DrawGraph(g, stylesheet=[vertexborder=false,vertexpadding=15,edgecolor = "Black",
     vertexcolor="Black",edgethickness=2])

EdgeConnectivity(g) 

6

One option I came up with is to find all 6-edge subsets first. Test if they satisfy  the restricted condition (one by one). Then continue to increase to 7 or more. But this violent calculation may get stuck in the first step. That is, test all the minimum edge cut sets (Note that we will consider 32468436 edge-subsets!) I was referring to the efficiency aspect.  

with(Iterator):
C:=Combination(58,6):
K:=Edges(g):
#sub:=seq( K[[seq(c[]+1)]], c=C); # do not run this code since it has 32468436 members.

 

Any language is acceptable.( C or C++, Python. )

PS: Some time ago, I also asked a related question (but with some differences) on mathematica stack (Find all the minimum edge cuts of a graph). Although Bob Hanlon  gave a reply, the actual result is not good.

 

Edits: The following literature gives a polynomial algorithm for computing the restricted-edge-connectivity of a given  graph. The heart of it is to computing the least cardinality of some  edge-pairs's edge separator. I'm stuck here.

  • Esfahanian A H, Hakimi S L. On computing a conditional edge-connectivity of a graph[J]. Information processing letters, 1988, 27(4): 195-199.

How to implement this algorithm is my current concern.

i have a problem in an optimzation problem. in the problem using NLPSolve to find the minimum, i have an integration which i use the Int command to be solved in the optimization process, but this error occures: Error, (in Optimization:-NLPSolve) could not store Int(..) in a floating-point rtable 
please help to solve the problem, tnx in advance

restart:with(LinearAlgebra):

N:=3:

m:=Vector([ 1 , log(x+b3) , b2/(x+b3) ]):

A:=m.m^+:

for i to N do
m||i:=eval(A,[x=x||i]);
od:

M:=add(w||i*m||i,i=1..N-1)+(1-add(w||i,i=1..N-1))*m||N:

MM:=( LinearAlgebra:-Trace(MatrixInverse(M)) ):

IF1:=evalf(Int(MM,[b2=1..2,b3=1..2],method = _d01ajc,epsilon=0.001)):

s:= Optimization:-NLPSolve(IF1,w1=0..1,w2=0..1,x1=1..10,x2=1..10,x3=1..10,variables=[w1,w2,x1,x2,x3],initialpoint={w1=0.6,w2=.1,x1=8,x2=7,x3=5},maximize=false,method=modifiednewton)

Error, (in Optimization:-NLPSolve) could not store Int(Int(16.6666666666666679*(-448.000000000000057*ln(7.+b3)*ln(5.+b3)+76.1999999999999886*ln(8.+b3)^2*b3^2+.199999999999999956*ln(8.+b3)^2*b3^4+.399999999999999911*ln(5.+b3)^2*b3^4+527.500000000000000*ln(5.+b3)^2+780.799999999999727*ln(8.+b3)^2+191.100000000000023*ln(7.+b3)^2+89.0999999999999943*ln(5.+b3)^2*b3^2+9.79999999999999893*ln(5.+b3)^2*b3^3+6.39999999999999858*ln(8.+b3)^2*b3^3+400.*ln(8.+b3)^2*b3+12.3000000000000025*ln(7.+b3)^2*b3^2+84.0000000000000142*ln(7.+b3)^2*b3+356.*ln(5.+b3)^2*b3+.600000000000000089*ln(7.+b3)^2*b3^3-1176.*ln(8.+b3)*ln(5.+b3)+280.*ln(8.+b ... 99999999999716*ln(7.+b3)*ln(5.+b3))^2), b2 = 1. .. 2.), b3 = 1. .. 2.) in a floating-point rtable

 

 

Download LinearLog-A-Bayesian_1.mw

Hi everyone

I am trying to get the maximum value (angle) for a function, which is a solution from a ODE. I tried evalf(max.. which I already thaught wouldnt work. 

After that I installed the package "DirectSearch", again with no success. 

Does anybody know what I am doing wrong or how I am going to get the maximum. I added the maple file with the direct seach attempt. 

Thank you in advance!

Hi

I'd like to know that why in the attachment "fsolve" does not work?

How can I evaluate the value of "A" in my file?

Many thanks in advance

NULL

restart:

N := 60: L := 10: a := 2*10^(-12):

eqn:=ns=1-((((N/A)+((A*L)^(L/(1-L))))^((1-L)/L))/(A*L))*(2+(a*exp((((N/A)+((A*L)^(L/(1-L)))))^(1/L)))/((A*L*(((N/A)+((A*L)^(L/(1-L))))^((L-1)/L)))+((a*exp((((N/A)+((A*L)^(L/(1-L)))))^(1/L)))))):

seq(fsolve(eqn,A),ns=[0.9603,0.9647,0.9691]):

 

NULL

Download A.mw

Hi everyone 

I would like to plot something in 3d. I have to functions both depending on the time. (x(t), z(t))

I can plot in 2d: (t, x(t)) or (t,z(t)) or (x(t), z(t)) 

Now I want to show the flight path (x(t), z(t)) over time t. So I thought a 3d plot would be suitable for this. 

The problem is that I dont know how to plot this. With plot3d I am only able to plot some planes. What I want to plot is the flight path in the 3d room with the axis x(t),z(t) and t.

Does anybody know how to do this?

thanks in advance!

L’éventail de la Geisha
restart:with(plots):with(geometry):
NULL;
_EnvHorizontalName := 'x':
_EnvVerticalName := 'y':

NULL;
EqBIS := proc(P, U, V) 
local a, eq1, M1, t, PU, PV, bissec1; 
description "P est le sommet de l'angle dont on chercche la bissectrice" ;
a := (P - U)/LinearAlgebra:-Norm(P - U, 2) + (P - V)/LinearAlgebra:-Norm(P - V, 2); 
M1 := P + a*t; eq1 := op(eliminate({x = M1[1], y = M1[2]}, t)); 
RETURN(op(eq1[2])); end proc:

with(plottools);
with(plots);


r1 := 1/2;
r2 := r1/2;
R := r1*(21 - 12*sqrt(3));
                            21      (1/2)
                       R := -- - 6 3     
                            2            

a := arc([0, 0], 2*r1, Pi/6 .. (5*Pi)/6);
b := arc([0, 0], r1, Pi/6 .. (5*Pi)/6);


with(geometry);
eq := EqBIS(<sqrt(3)/2, 1/2>, <0, 0>, <0, 1/2>);
line(bis, eq);
                         (1/2)                  
                eq := 2 3      y - 2 x + 4 y - 2

                              bis

OpT := 2*sqrt(r1*R);
line(lv, x = OpT);
intersection(Omega, bis, lv);
coordinates(Omega);
evalf(%);
                                (1/2)    
                      OpT := 2 3      - 3

                               lv

                             Omega

                 [                / (1/2)    \]
                 [   (1/2)      2 \3      - 1/]
                 [2 3      - 3, --------------]
                 [                     (1/2)  ]
                 [                2 + 3       ]

                  [0.464101616, 0.3923048456]

retarrt;
with(plots);
with(plottools);
[cos((5*Pi)/6), sin((5*Pi)/6)];
                        [  1  (1/2)  1]
                        [- - 3     , -]
                        [  2         2]

a := arc([0, 0], 2*r1, Pi/6 .. (5*Pi)/6);
b := arc([0, 0], r1, Pi/6 .. (5*Pi)/6);
NULL;
A:=[cos(Pi/6), sin(Pi/6)];
B:=[cos(5*Pi/6), sin(5*Pi/6)];
Oo:=[0,0];
Op:=[0,1/2];
poly:=[A,B,Oo];
R := r1*(21 - 12*sqrt(3))
                            [1  (1/2)  1]
                       A := [- 3     , -]
                            [2         2]

                           [  1  (1/2)  1]
                      B := [- - 3     , -]
                           [  2         2]

                          Oo := [0, 0]

                                [   1]
                          Op := [0, -]
                                [   2]

                [[1  (1/2)  1]  [  1  (1/2)  1]        ]
        poly := [[- 3     , -], [- - 3     , -], [0, 0]]
                [[2         2]  [  2         2]        ]

                            21      (1/2)
                       R := -- - 6 3     
                            2            


Omega := [2*sqrt(3) - 3, 2*(sqrt(3) - 1)/(2 + sqrt(3))];
Omega1 := [3 - 2*sqrt(3), 2*(sqrt(3) - 1)/(2 + sqrt(3))];

                     [                / (1/2)    \]
                     [   (1/2)      2 \3      - 1/]
            Omega := [2 3      - 3, --------------]
                     [                     (1/2)  ]
                     [                2 + 3       ]

                     [                 / (1/2)    \]
                     [    (1/2)      2 \3      - 1/]
           Omega1 := [-2 3      + 3, --------------]
                     [                      (1/2)  ]
                     [                 2 + 3       ]


r3 := 3/16;
EF := sqrt(r3);

                                  3 
                            r3 := --
                                  16

                               1  (1/2)
                         EF := - 3     
                               4       

r := (150 - 72*sqrt(3))/193*1/2;
alpha := -5/3*r + 1/2*1/2;
p := sqrt(3)/3*1/2 - sqrt(3)/18*r;
                          75    36   (1/2)
                     r := --- - --- 3     
                          193   193       

                             307   60   (1/2)
                  alpha := - --- + --- 3     
                             772   193       

               1  (1/2)   1   (1/2) /75    36   (1/2)\
          p := - 3      - -- 3      |--- - --- 3     |
               6          18        \193   193       /

p2 := textplot([[A[], "A"], [B[], "B"], [Oo[], "O"]], align = ["above", "right"]);
display(a, b, p2, polygonplot(poly, thickness = 3, color = blue, transparency = 0.3), circle(Omega, R, color = blue, filled = true), circle(Omega1, R, color = blue, filled = true), circle([0, 3/4], 1/4, color = yellow, filled = true), circle([EF, 1/2 + r3], r3, color = green, filled = true), circle([-EF, 1/2 + r3], r3, color = green, thickness = 5), circle([p, 3/4 + alpha], r, color = red, thickness = 5), circle([-p, 3/4 + alpha], r, color = red, thickness = 5), axes = none, scaling = constrained, size = [500, 500]);
how to put color inside circles ? Thabk you.

restart;
with(plots):
with(geometry):
_EnvHorizontalName := x:
_EnvVerticalName := y:
R := 11:
r := 7:
a := sqrt(R*r):

b := 2:
circle(C1, [point(P1, [0, 0]), R]):
circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
ellipse(p, (x - R - b)^2/b^2 + y^2/a^2 = 1):
draw([C1(color = yellow, filled = true), 
C2(color = red, filled = true), p(color = blue, filled = true), 
C1(color = black), C2(color = black), p(color = black)], 
axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained):
alpha := arctan((R - r)/(R + 2*b + r));
long := cos(alpha)*(R + 2*b + r);
evalf(%);
circle(C2, [point(P2, [long, r - R]), r]);
rotation(p1, p, alpha, 'clockwise');
detail(p1);
point(A, 0, -R);
point(B, long, -R);
line(L1, [A, B]);
point(cen, [(143*sqrt(5))/25, -(26*sqrt(5))/25]);
reflection(L2, L1, cen);
detail(L2);
Error, (in geometry:-reflection) unable to compute coeff
Error, (in geometry:-detail) unknown object:  L2

draw([C1(color = yellow, filled = true), C2(color = red, filled = true), p1(color = blue, filled = true), C1(color = black), C2(color = black), p1(color = black), L1(color = black)], axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained);
A Bug in reflection ? Why these error messages. Thank you.

I wonder if there is any way to use ArrayInterpolation with contourplot or similar effect?

N_data.xlsx 

Thank you in advance,

restart;

with(CurveFitting)

[ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares, Lowess, PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation]

(1)

with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(2)

alpha := <seq(0..10,evalf(10/50))>:
beta := <seq(0..10,evalf(10/50))>:

excelfile:= FileTools:-JoinPath(["C:","Users","aimer","OneDrive","Desktop","Msc Thesis","Maple ref","N_data.xlsx"]);

"C:\Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref\N_data.xlsx"

(3)

NN:=ImportMatrix(excelfile,source=Excel):

_rtable[36893489576445216036]

(4)

#?ImportMatrix;

#NN:=ImportMatrix(matlabData, source=MATLAB);

#currentdir();

"C:\Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref"

(5)

 

contourplot(ArrayInterpolation([beta,alpha],NN,[x,y]),x=0..10,y=0..10,contours=[0]);

Error, (in CurveFitting:-ArrayInterpolation) invalid input: xvalues are not specified correctly

 

#?listcontplot

 

Download test1.mw

I do not know what is the problem with Using ImportMatrix. N_data.xlsx is in the same directory.

Any comment would be appreciated.

restart;

with(CurveFitting)

[ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares, Lowess, PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation]

(1)

with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(2)

alpha := <seq(0..10,evalf(10/50))>:
beta := <seq(0..10,evalf(10/50))>:

excelfile:= FileTools:-JoinPath(["C: ","Users","aimer","OneDrive","Desktop","Msc Thesis","Maple ref","N_data.xlsx"]);

"C: \Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref\N_data.xlsx"

(3)

NN:=ImportMatrix(excelfile,source=Excel);

Error, (in ImportMatrix) file or directory does not exist: C: \Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref\N_data.xlsx

 

?ImportMatrix;

#NN:=ImportMatrix(matlabData, source=MATLAB);

currentdir();

"C:\Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref"

(4)

?Joinpath

 

Download test1.mw

polysols(diff(u(x), x) = u(x)^2 - 1) produces no results, while it can be verified by direct observation that u(x) = 1 is a polynomial solution.

I can not spawn draw the circles C3 and C4

restart;
with(plots):
with(geometry):
_EnvHorizontalName := x:
_EnvVerticalName := y:
R := 7:
point(A, [0, R]):
line(L1, y = sqrt(3)*x + R):
line(L2, y = -sqrt(3)*x + R):
line(L3, y = R/3):
intersection(B, L1, L3):
intersection(C, L2, L3):
detail(C):
triangle(ABC, [A, B, C]):
circle(C1, [point(P1, [0, 0]), R]):
circle(C2, [point(P2, [0, R/3 + (2*R)/9]), (2*R)/9]):
detail(C2):
center(C2), coordinates(center(C2)):
reflection(P3, P2, C):
detail(P3):
reflection(C3, C2, C);
detail(C3):
Error, (in geometry:-reflection) unable to compute coeff
Error, (in geometry:-detail) unknown object:  C3
circle(C3, [point(P3*[(28*sqrt(3))/9, 7/9]), (2*R)/9]):
Error, (in geometry:-point) wrong number of arguments
reflection(C4, C2, B);
detail(C4);
Error, (in geometry:-reflection) unable to compute coeff
Error, (in geometry:-detail) unknown object:  C4
circle(C3*[point(P3, [(28*sqrt(3))/9, 7/9]), (2*R)/9]);
Error, (in geometry:-circle) wrong number of arguments


draw([L1(color = blue), 
ABC(color = red, transparency = 0.5, filled = true), 
L2(color = blue), L3(color = blue), 
C1(color = blue, thickness = 3), C1(color = yellow, transparency = 0.8, filled = true), C2(color = blue, filled = true)], 
axes = normal,
 view = [-R .. R, -R .. R], 
scaling = constrained);
Why these error messages. Thank you veru much.

fig([L1(color = blue), L2(color = blue), L3(color = green), C1(color = black), C2(color = black), C3(color = orange), C4(color = orange)]);
point(oo, [0, 0]);
                               oo

rotation(fig1, fig, oo, Pi/2, 'counterclockwise');
Error, (in geometry:-rotation) wrong type of arguments why thos error ? Thank you;
 

What should I do to reduce evaluating time?

restart;

with(plots):

 

F:=kappa->kappa;

proc (kappa) options operator, arrow; kappa end proc

(1)

f:=(alpha,delta)->exp(-abs(F(kappa))^2*(1+delta^2)/2-abs(F(kappa))*alpha)/abs(F(kappa));

proc (alpha, delta) options operator, arrow; exp(-(1/2)*abs(F(kappa))^2*(1+delta^2)-abs(F(kappa))*alpha)/abs(F(kappa)) end proc

(2)

L:=(alpha,delta,Lambda)->(lambda^2*exp(-alpha^2/2)/4)*(Int(f(alpha,delta),kappa= -infinity..-Lambda)+Int(f(alpha,delta),kappa= Lambda..infinity));

proc (alpha, delta, Lambda) options operator, arrow; (1/4)*lambda^2*exp(-(1/2)*alpha^2)*(Int(f(alpha, delta), kappa = -infinity .. -Lambda)+Int(f(alpha, delta), kappa = Lambda .. infinity)) end proc

(3)

evalf(L(4,1,0.001));

0.8209373770e-3*lambda^2

(4)

g:=(beta,delta)->exp(-I*kappa*beta-abs(F(kappa))^2*(1+delta^2)/2)/abs(F(kappa));

proc (beta, delta) options operator, arrow; exp(-I*kappa*beta-(1/2)*abs(F(kappa))^2*(1+delta^2))/abs(F(kappa)) end proc

(5)

E:=(omega,gamma)->exp(I*omega*gamma)*(1-erf((gamma+I*omega)/sqrt(2)));

proc (omega, gamma) options operator, arrow; exp(I*omega*gamma)*(1-erf((gamma+I*omega)/sqrt(2))) end proc

(6)

J:=(alpha,delta,Lambda,beta,gamma)->(lambda^2*exp(-alpha^2/2)/8)*abs(Int(g(beta,delta)*(E(abs(F(kappa)),gamma)+E(abs(F(kappa)),-gamma)),kappa=-infinity..-Lambda)+Int(g(beta,delta)*(E(abs(F(kappa)),gamma)+E(abs(F(kappa)),-gamma)),kappa=Lambda..infinity));

proc (alpha, delta, Lambda, beta, gamma) options operator, arrow; (1/8)*lambda^2*exp(-(1/2)*alpha^2)*abs(Int(g(beta, delta)*(E(abs(F(kappa)), gamma)+E(abs(F(kappa)), -gamma)), kappa = -infinity .. -Lambda)+Int(g(beta, delta)*(E(abs(F(kappa)), gamma)+E(abs(F(kappa)), -gamma)), kappa = Lambda .. infinity)) end proc

(7)

#evalf(J(4,1,0.001,8,3));

N := (beta,alpha)-> (J(alpha,1,0.001,beta,3)-L(alpha,1,0.001))/\lambda^2;

proc (beta, alpha) options operator, arrow; (J(alpha, 1, 0.1e-2, beta, 3)-L(alpha, 1, 0.1e-2))/lambda^2 end proc

(8)

 

 

 

 

 

 

contourplot(evalf(N(beta,alpha)), beta=0..10,alpha=0..10,grid=[25,25]);

 

 

 

 

Download Negativity_v1.mw

1 2 3 4 5 6 7 Last Page 2 of 30