## initial condition error...

Hi everyone

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

before:

new:

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.

## calculation of Major Axis and Minor Axis of an ell...

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

## How to compute the restricted-edge-connectivity of...

！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~jPbOqgLJ@p");
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.

## Error, (in Optimization:-NLPSolve) could not store...

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

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

## evalv max/min of function...

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.

## why fsolve doesn't work...

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?

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

## How to plot 2d flight path over time in 3d (x,z,t)...

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?

## how to put color inside circles...

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.

## difficulty with reflection...

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.

## ArrayInterpolation in contourplot...

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

 > restart;
 > with(CurveFitting)
 (1)
 > with(plots);
 (2)
 > alpha := : beta := :
 > excelfile:= FileTools:-JoinPath(["C:","Users","aimer","OneDrive","Desktop","Msc Thesis","Maple ref","N_data.xlsx"]);
 (3)
 > NN:=ImportMatrix(excelfile,source=Excel):
 (4)
 > #?ImportMatrix;
 > #NN:=ImportMatrix(matlabData, source=MATLAB);
 > #currentdir();
 (5)
 >
 > contourplot(ArrayInterpolation([beta,alpha],NN,[x,y]),x=0..10,y=0..10,contours=[0]);
 > #?listcontplot
 >

## ImportMatrix problem...

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)
 (1)
 > with(plots);
 (2)
 > alpha := : beta := :
 > excelfile:= FileTools:-JoinPath(["C: ","Users","aimer","OneDrive","Desktop","Msc Thesis","Maple ref","N_data.xlsx"]);
 (3)
 > NN:=ImportMatrix(excelfile,source=Excel);
 > ?ImportMatrix;
 > #NN:=ImportMatrix(matlabData, source=MATLAB);
 > currentdir();
 (4)
 > ?Joinpath
 >

## Why the polysols does not work?...

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.

## Drawing of circles...

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.

## difficulty wirth rotation...

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;

## how to reduce evaluating time...

What should I do to reduce evaluating time?

 > restart;
 > with(plots):
 >
 > F:=kappa->kappa;
 (1)
 > f:=(alpha,delta)->exp(-abs(F(kappa))^2*(1+delta^2)/2-abs(F(kappa))*alpha)/abs(F(kappa));
 (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));
 (3)
 > evalf(L(4,1,0.001));
 (4)
 > g:=(beta,delta)->exp(-I*kappa*beta-abs(F(kappa))^2*(1+delta^2)/2)/abs(F(kappa));
 (5)
 > E:=(omega,gamma)->exp(I*omega*gamma)*(1-erf((gamma+I*omega)/sqrt(2)));
 (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));
 (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;
 (8)
 >
 >
 >
 >
 >
 >
 > contourplot(evalf(N(beta,alpha)), beta=0..10,alpha=0..10,grid=[25,25]);
 >
 >
 >