6 years, 21 days

## @vv Just a remark about your p...

@vv

Just a remark about your procedure EQSYS

• it returns an error when the two planes are parallel and unconfounded
• it responds "false" if the planes are confounded and the line belongs to each of them
(assuming that the intersection of two confounded planes is these same planes and that each line in these latter can thus be considered as a particular intersection among many; thus "true" seems to be a better answer, even if geom3d syays that the planes do not intersect ... maybe a definition of the intersection of two planes?)

sys1 :={16*x-2*y-11*z=0, 14*x-y-10*z=3}:
sys1a:={16*x-2*y-11*z=0, 32*x-4*y-22*z=0}: #confounded planes
sys1b:={16*x-2*y-11*z=0, 32*x-4*y-22*z=1}: #parallel uncounfounded planes

sys2:={(x-2)/3=(y-5)/2, (y-5)/2=(z-2)/4}:

EQSYS(sys1 , sys2);  # answer is true
EQSYS(sys1a, sys2);  # answer is false but the line is in the plane
EQSYS(sys1b, sys2);  # error


In cases sys1a and sys1b geom3d would respond

plane(P1, ...., [x,y,z]);
plane(P2, ...., [x,y,z]);
line(L1, [P1,P2]);

P1 = sys1a --> intersection: the two planes P1 and P2 are the same
Error, (in geom3d:-line) the two given planes are parallel

P2 = sys1b --> intersection: the two planes P1 and P2 are parallel
intersection: the given objects do not intersect


## @JAMET Of course...But in a general...

@JAMET

Of course...
But in a general situation you must distinguish several cases for the two planes (which may intersect or not):

• Let a*x+b*y+c*z+d the equation of a plane and a'*x+b'*y+cC'*z+d' the equation of the second one
• if a'/a=b'/b=c'/c=d'/d then P and P' are identical
• if a'/aA=b'/b=c'/c then Pand P' are parallel
• if at least one of the three equalities a'/a=b'/b, b'/b=c'/c, a'/a=c'/c  is not true then P and P' intersect along a line D
In this case the elementary way to obtain the solution is:
• to consider that one of the unknown (x or y or z) is a "parameter" t (for instance x=t) and solve a system of two equations in the two remaining unknowns.
• to solve a system of two equations in the two remaining unknowns, for instance
solve({a*tx+b*y+c*z+d=0,  a'*t+b'*y+cC'*z+d'}, {y, z})

For your case this gives (for instance)

restart
e1 := 16*s-2*y-11*z:
e2 := 14*s-y-10*z-3:

# Set x as a parameter "s"
E0 := [x=s, solve({e1, e2}, {y, z})[]]
[           2     11      4     2]
[x = s, y = - s + --, z = - s - -]
[           3     3       3     3]

# verify that this line is the one geom3d finds
with(geom3d):
plane(P1, 16*x-2*y-11*z, [x,y,z]):
plane(P2, 14*x-y-10*z-3, [x,y,z]):
line(L1, [P1,P2]):
E1 := [x, y, z] =~ Equation(L1,'t');
[    1                             ]
[x = - + 9 t, y = 4 + 6 t, z = 12 t]
[    2                             ]

# as the parameterization can be different, express s as a finction of t
# and use this to re-parameterize E0 in terms of E1 to verify that E0 and
# E1 are indeed identical

solve(rhs(E0)=rhs(E1), s);
eval(E0, s=%);
1
- + 9 t
2
[    1                             ]
[x = - + 9 t, y = 4 + 6 t, z = 12 t]
[    2                             ]


So to answer your question "Is it possible showing without geom3d ?", just do

E0 := [x=s, solve({e1, e2}, {y, z})[]]

And now to verify if your two lines are the same :

X := s:
Y := solve((y-5)/2=(s-2)/3, y):
Z := solve((z-2)/4=(s-2)/3, z):
[x, y, z] =~ [X, Y, Z]


You might think it's shorter than using geom3d, but for general planes this doesn't work (look at the different situations in the beginning of this reply).
Why did I use geom3d? Just because I didn't know if your two planes where confounded, parallel or secant.

PS: I EDITED MY PREVIOUS REPLY
It contained an error

## @Preben Alsholm Right, but I think ...

@Preben Alsholm  @Michael_Watson

Right, but I think the original mistake could be a typo:

ex := L/(z+sqrt(z^2 + L^2));  # + or - instead of *
eval(ex, L=alpha*z);
simplify(%) assuming z > 0
alpha
---------------------
(1/2)
/     2    \
\alpha  + 1/      + 1


Or

eval(ex, z=beta*L):
simplify(%) assuming L > 0
1
-----------------------
(1/2)
/    2    \
\beta  + 1/      + beta


## @vv Hi, I naively tried to use...

@vv

Hi,
I naively tried to use rsolve for the first problem but I got no answer.
I am wrong writting this

rsolve({a(n+1)=M(a(n), b(n)), b(n+1)=a(n)*b(n)/M(a(n), b(n))}, {a, b})

or am I facing some limitation of rsolve?

Concerning problem b:

• I have been stucked for a while trying to find an invariant or some map from (a, b) to (u, v) whoch could transform the problem in something simpler.
Do you have some more information about the proof that "The limit M⊗N is the logarithmic mean"?

• About the convergence speed :
cat(b__, n+1)-cat(a__, n+1) = N(a__n,b__n)-M(a__n,b__n);

1        1
b__n+1 - a__n+1 = - b__n - - a__n
2        2

cat(b__, n+1)-cat(a__, n+1) = (b__0-a__0)/2^(n+1)

b__0 - a__0
b__n+1 - a__n+1 = -----------
(n + 1)
2


@acer

Point 1
I had interpreted your claim "Using Maple 2021.1, we can work around that by using..." as a way (not the way) to handle the problem in Maple 202.

plot('funktion(r)', r=r1 .. r2);

won't work in Maple 2021 in the case we use Units?

Point 2
There is a basic problem in the definition of the funktion the OP uses

funktion := x -> if 2*Unit('m') < x then 1/x; else 2*Unit('m'); end if:


For the comparison to be possible x must ne e Length (m, cm, inches, ...)
Assuming this is the case the first part of the test (2*Unit('m') < x)  returns a quantity homogeneous to 1/Length, as the else part of the test funktion returns a quantity homogeneous to a Length.
Think to this: what label would you put on the vertical axis if the plot was to be from 1.5m to 4m?
Would it be m or 1/m ?

Point 3
As you probably remember I use Maple 2015
Could you please tell me if newer versions give two identical plots or if, as with Maple 2015, the second plot is incorrect.

restart
with(Units[Standard]):
funktion := x -> piecewise(2*op(2, x) < x, 1/x, 1/(2*op(2, x))):
a := 1.0:
b := 4.0:

correct := plot(
[seq([r, op(1, funktion(r*Unit('m')))], r in [seq(a..b, (b-a)/100)])],
labels=[typeset(r*[Unit('m')]), typeset(f*[Unit('m^(-1)')])],
color=blue, thickness=7, transparency=0.5, legend="correct"
):

wrong := plot(funktion, a*Unit('m')..b*Unit('m'), color=red, legend="wrong"):

plots:-display(correct, wrong) Said otherwise, what plot do you obtain with Maple 2021 if you run your code with this slight modification?

restart;
with(Units[Simple]):
funktion := x -> if 2*Unit('m') < x then 1/x; else 2*Unit('m'); end if:
r1 := 1.0*Unit('m'):
r2 := 6*Unit('m'):

plot(funktion, r1 .. r2);

I hope you get this ## @Thomas Richard Is there any inconv...

@Thomas Richard

Is there any inconvenience to use

plot('funktion(r)', r=r1 .. r2);

plot(funktion, r1 .. r2);

?
The fact that we have to change the syntax of the plot procedure seems to me to be a weakness of Maple, a lack of genericity of this procedure in the sense that Maple should (IMO) understand that it does not have to evaluate the function prematurely.
This is however what it does here

plot([seq([r, funktion(r)], r=3..6)])


## @acer  It is funny how the same peo...

@acer
It is funny how the same people can sometimes understand the meaning of the word "extract" while categorically refusing it at other times

@acer

Tgank you acer

## @TeepTo understand without any procedure...

@Teep

To understand without any procedure why A->B->D->E->C->A (or any cyclic permutation of it) is the shortest cycle:

• introduce these 3 fictitious points P=(0, 5), Q=(5, 5), R=(5, 0)
• the cycle A->P->C->Q->E->R->A has length 20 (obvious)
• if you decide once arrived at E to take the interiorreturn path E->D->B->A instead of E->R->A you won't reduce the length of the path from E to A.
• thus the shortest cycle (a shortest cycle among many others) has length 20

## @acer @tomleslieI agree tthat "...

@acer @tomleslie

I agree that "The plots:-surfdata command can produce either a 3D surface plot or a 2D density-style plot", but if one does this

restart:
with(plots):
cosdata := [seq([ seq([i,j, evalf(cos((i+j)/5))], i=-5..5)], j=-5..5)]:
a := surfdata(cosdata, dimension=2, axes=none):
op(a)


one sees that op(a) is a MESH structure.
In the plot,structure help page its written
The MESH structure represents surfaces in 3-D space defined by a grid of values

Thus, I would be tempted to say that although surfdata can produce a 2D density graph, it seems essentially to be a 3D graph and that dimension=2 "simply makes a projection" of the 3D surface onto the [z,y] plane?

So, is Tom's statement really that wrong?

## @tomleslie It was not a suggestion ...

@tomleslie

It was not a suggestion I made to the OP but to you, you were free to take it into account or not.
If I had wanted to write to the OP I would have done so.
Period

## @tomleslie  "By the way it is...

@tomleslie

"By the way it is alway helpfull..."
It was just a suggestion for a tiny improvement of your own code, more a private communication than a real answer.
So I felt that publishing a code 99% identical to yours would be of no use, not to mention plagiarism.
And between you and me, I don't think you spent that much time putting this suggestion into your code.

## @PepiniFor the phase portrait you can tr...

@Pepini

For the phase portrait you can try this
(3 types of solutions are displayed if you play the animation)

with(DEtools):

phaseportrait(
{diff(x(t), t)-v(t)=0, diff(v(t), t)-2*Heaviside(x(t)^2-1)=0},
[x(t), v(t)],
t=-3..4,
[[v(0)=-2, x(0)=-4], [v(0)=-1, x(0)=2], [v(0)=1, x(0)=-2]]
,linecolor=[gold, green, cyan]
,animatecurves=true
,stepsize=0.001
,arrows=curve        # try comet also
#,dirfield=[30, 30]  # try other values
); ## @tomleslie Why not define Cnames th...

@tomleslie

Why not define Cnames this way;

Vnames := map(u -> sprintf("%a", u), faces);

then the names of the faces directly refer to the vertices of the primal graph. ## @Pepini Just change what I did thus...

@Pepini

An event-based solution (same system as in my previous reply)

sol := dsolve(
sys,
numeric,
discrete_variables=[A(t)],
events=[
[[x(t)= 1, And(v(t) > 0)], A(t)=1],
[[x(t)= 1, And(v(t) < 0)], A(t)=0],
[[x(t)=-1, And(v(t) > 0)], A(t)=0],
[[x(t)=-1, And(v(t) < 0)], A(t)=1]
]
);

# to verify A(t) behaves correctly plot A(t) versus x(t)

plots:-odeplot(sol, [x(t), A(t)], t=0..2)


An Heaviside-based solution

ode := diff(x(t), t)=v(t), diff(v(t), t) - 2*Heaviside(x(t)^2-1) = 0:
ic  := x(0)=-2, v(0)=1:
sys := {ode, ic}:
sol := dsolve(sys, numeric);


PS: I guess someone here will be able to present a way to solve formally this "system" (I can't).

﻿