@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 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.
Does your last reply mean that 

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

Thanks in advance

@Thomas Richard 

Is there any inconvenience to use 

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

instead of

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 people can sometimes understand the meaning of the word "extract" while categorically refusing it at other times

@acer 

Tgank you acer

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

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

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

 

@tomleslie 

Sorry for this late reply (I thought the thread was closed.
About the first issue: my question was very problem-dependent.
about the second one: as I said before I wasn't aware of this second way to do animation I will keep it in mind and try to use it systematically to see if it's indeed  a good alternative.

Thanks

@BasSPe 

You're welcome