@Preben Alsholm 

I haven't understand yet. But, it works perfectly! Great !

May you precise the syntax that Maple uses to make this plot ? I mean the arguments used by the plot function so as to obtain the result.

The point which troubles me is the fact that it is a bit like if you have a sequence for the t variables (on y-axis) but there is no sequence for t variable.

Thank you Preben for your help.

Your code is very good.

For the plot of the vertical lines, I try this plot([seq(x=Pz[i],i=1..nops(Pz))],x=0..2, colour=yellow,linestyle=3, thickness=2): but it doesn't work yet.

Do you have ideas what's wrong in this syntax ?

@Preben Alsholm @Thomas Richard

Thank you Preben for your ideas

But, the point is that I need to define the functions because I use these functions in the rest of my worksheet. That's why it seems to me that the solution with evaln can be relevant, doesn't it ?

@tomleslie 

OK. Thank you for your comments on indentation. For my culture, can you tell me the external editor that you use?

Otherwise, I'm trying now to find the zero-crossings with a increasing curve.

Similarly in gaitPlot4.mw, I calculated a whole sequence of values for each of v[i](t): from this data, it is relatively trivial to work out a pretty good approximation to the zero-crossings (with a little bit of interpolation, if necessary). 

I try to find a function similar to max function so as to find the closest value to zero in a list of numerical values. But, i haven't find yet.

I try this but I find not all the solutions and the results is quite strange.

assign(v[1],rhs(res[6]));
Roots(v[1](t)=0,t=0.5..2);
[-0.2526518250e-1]

I also try this : select(s->abs(s)<0.1,vvals);

but it also fails.

May you give me more information on your ideas so as to determine the zero-crossing values?

Thanks a lot for your help.

I find a solution :

with(Student[Calculus1]):
Roots(v(t),t=0.5..2);

Nevertheless, I'm interesting by your experience on this subject.

Sorry, for not having upload the last version of my worksheet.

I do use all your remarks and I thank you for your tips.

For the last attachment, I wanted to focuse on the display of the points. That's why i didn't send you the whole code but I have used all your modifications about vSpace and it works well.

For the display of the points, you perfectly answer to my question. 

These code lines have enabled to evaluate the u[i] and v[i].

seq( assign(u[j], rhs(res[j+1])) , j=1..4);
seq( assign(v[j], rhs(res[j+5])) , j=1..4);

Nota: I notice that you indentation in your code. Do you have an automatic way to make code lines with indentation ?

Thank you 

@tomleslie 

Thanks a lot for your very detailled and relevant help.

Concerning the points display, it comes from the output=listprocedure in the dsolve function. I attached this file in order to show you what it was done before to plot points on the ellipse GaitPlot_points.mw

May you see if you can correct the file GaitPlot5 file so as to enable the display of some points (at t = 0, 0.5, 1, 1.5 as in the file attached) ?

Thank you for your help

@ Dr. Venkat Subramanian

Thank you for you answer.

I agree that MapleSim have DAE solver.

But, are you sure that the NL oscillator is a DAE ?

For my me, it is a set of 8 coupled ODES but not a DAE. Indeed, there is not algebraic equations.

But, you probably have far more experiences than me on DAE.

@tomleslie 

OK. Thank you for you feedback. I try to follow your tips. I try to apply you method on the NL oscillator that I use in my robot model.

But, I think that slight modifications on the calculation of the vertical space vSpace should be conducted.

In order to facilitate the tuning of the display of the plots, I attach here a worksheet where I added comments so as to facilitate the reading.

GaitPlotTrajectory.mw

I hope you will find this project more interesting and I thank you again for your help.

@tomleslie 

Hello,

Thank you again for your help and relevant advices.

I'm very close to the result expected.

1) For the adaptation to the parameters modification

Thanks to the parameter vSpace, I almost obtain a "nice" result from the beginning. You can see that I need to add a coefficient with a value of 5 to get a good plot. I will be very interested if you see how I can adapt my plot without adding this value of 5 but with a more appropriate parameter. In order to help you to have a better understanding of my NL oscillator, I put the code at the end of this message. It will enable you not to make the magicien to understand my parameters. Concerning vSpace, the NL oscillator is design to make an ellipse with a heigh 2H and a length L. So, I think that vSpace:=2H should fit well. What do you think about that ?

2) Adding vertical lines 

Another interesting display could be to add vertical lines at the beginning of each red phase. In fact, the red phase correspond to the phase where the feet are in contact with the floor.

Here a example of the vertical lines that I would like to add :

Do you have ideas so as to add these vertical lines ?

3) Points on the ellipse

There is a element which does no longer work. It is the addition of points on the ellipse. I think it comes from the fact that I use now the dsolve with output=listprocedure. Do you have ideas so as to obtain the coordinates p0, p1, p2 and p3 in this case where res is obtained with this dsolve with output=listprocedure ?

K:=Matrix([<0, -1, -1, 1>,<-1, 0, 1, -1>,<1, -1, 0,-1>,<-1, 1, -1,0>]);
omega[sw]:=beta/(1-beta)*omega[st];
for i to 4
do
r[i]:=sqrt((u[i](t)/(L/2))^2+(v[i](t)/H)^2):
omega[i]:=omega[st]/(1+exp(b*v[i](t)))+omega[sw]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)+omega[i]*(L/2)/H*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)-omega[i]*H/(L/2)*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do:
EqSys[1];
paramsGeo:=L=0.080,H=0.020,beta=0.5,Vf=0.3;
omegaS:=eval(Pi*Vf/L, [paramsGeo]);
paramsCycle:=omega[st]=omegaS,Au=5,Av=50,b=100;
params:=paramsGeo,paramsCycle;
sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));
ic:=[u[1](0)=-0.02, v[1](0)=0,u[2](0)=-0.015, v[2](0)=0,u[3](0)=0, v[3](0)=0.01,u[4](0)=0, v[4](0)=0.01];
Resolution
res:=dsolve([sys[],ic[]],numeric, maxfun = 0);
Initial boundaries calculations
tcalc:=8;
ic:=[seq(u[i](0)=eval(u[i](t), res(tcalc)),i=1..4),seq(v[i](0)=eval(v[i](t), res(tcalc)),i=1..4)];
Resolution2
res:=dsolve([sys[],ic[]],numeric, maxfun = 0,output=listprocedure);
Tracé
tmax:=10;

numpts:=100*tmax;
plots:-odeplot(res,[u[1](t),v[1](t)],t=0..tmax,scaling=constrained,numpoints=numpts);
plt1 := %:
plots:-odeplot(res,[u[2](t),v[2](t)],0..tmax,scaling=constrained,numpoints=numpts);
plots:-odeplot(res,[u[3](t),v[3](t)],0..tmax,scaling=constrained,numpoints=numpts);
plots:-odeplot(res,[u[4](t),v[4](t)],0..tmax,scaling=constrained,numpoints=numpts);

Points on ellipse
p0:=eval([u[1](t),v[1](t)], res(0))
p1:=eval([u[1](t),v[1](t)], res(0.5));
p2:=eval([u[1](t),v[1](t)], res(1));
p3:=eval([u[1](t),v[1](t)], res(1.5));
pointplot([p0, p1, p2,p3],color=[red,green,blue,pink], symbol=solidcircle, symbolsize=20); plt2 := %:
display([plt1, plt2]);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..tmax,thickness=2,numpoints=numpts);
vSpace:=eval(H, [paramsGeo]);
tmax:= 10:
numpts:=100*tmax:
TraceCourbe:= plots:-odeplot(res,[seq([t,v[i](t)+(i-1)*5*vSpace],i=1..4)],0..tmax,thickness=2, view=[0..5, -1*5*vSpace..4*5*vSpace], numpoints = numpts):
tmax:=5:
for j from 1 by 1 to 4 do
v:= rhs(res[5+j]);
la||j:= plot( [seq(`if`(v(t)>0,[t, vSpace*5*(j-1)],NULL),t=0..tmax,0.01)],
style=point,
color=blue,
symbol=solidcircle,
symbolsize=2
):
lb||j:= plot( [seq(`if`(v(t)<0,[t, vSpace*5*(j-1)],NULL), t=0..tmax,0.01)],
style=point,
color=red,
symbol=solidcircle,
symbolsize=2
):
od:

plots[display]([TraceCourbe, seq([la||j, lb||j][], j=1..4)],tickmarks=[default, [0=0, 0.1=0, 0.2=0, 0.3=0]],axes=boxed);

Thanks a lot for your help

@tomleslie 

Thank a lot for your help.

For 1) about the x-axis, it is now perfect with the option axes=boxed 

For 2) Your technique for finding the minimum and maximum are great. Thank you. However, I would have liked to display this max and min value in the graph in fact. 

3) I meet another problem. When I change the parameters of the NL oscillator, I need to adapt the term in bold in the following equation:

p1:= plots:-odeplot(res,[seq([t,v[i](t)+(i-1)*5],i=1..4)],0..tmax,thickness=2, view=[0..5, -1..16], numpoints = numpts):

Do you have ideas so as to have a more scalable (or flexible) way in case of modifications the parameters of the oscillator ? The idea is to adapt the plotting in case of curves with different amplitudes.

Thank you in advance for your ideas.

@tomleslie 

Your solution "Edit->Remove Output->From Worksheet" works perfect !

Thank you

@tomleslie 

Thanks a lot for your help.

It is almost perfect for me.

The only points remaining are :
- for the the blue/ red lines associated to v[1](t), the blue and red colors are not easy to see because the x-axis also on this line.

- Concerning the y-axis, the display is not so bad, but it is a pity that it is not possible to have a appropriate graduation by using this method. I could be interested to analyse for example the amplitude of a maximum which is a bit difficult here because the graduation of the y-axis is quiet limit.

@tomleslie 

Thanks a lot for you help.

I find the result not so bad and not too far from what I need.

Here the result obtained with you method :

1) About the code

First, I try to understand your code :

tmax:=5:
for j from 1 by 1 to 4 do
v:= rhs(res[5+j]);
la||j:= plot([seq(`if`(v(t)>0,[t, 5*j],NULL),t=0..tmax,0.01)], color=blue, style=point):
lb||j:= plot([seq(`if`(v(t)<0,[t, 5*j],NULL),t=0..tmax,0.01)], color=red, style=point):
od:
plots[display]([p1, seq([la||j, lb||j][], j=1..4)]);

May you explain me what the code res[...] means ? In fact, I usually use res() for asking a value but not with brackets "[]" (Nota : res:=dsolve(...))

2) About the gait diagrams

This result almost fits with my expectation.

a) Do you have ideas so as to make the blue and the red line smaller ?

b) Is it possible to plot this kind of gait diagrams by starting for the vertical axes each time by 0 like in this way :

Thanks a lot for you help

 This subject may interest you. That's why I put you informed.

@Carl Love 

@Preben Alsholm 

@Preben Alsholm 

@tomleslie 

@Carl Love 

Thank you for your support.

Concerning the remarks of @tomleslie : 

1) I think that you are right that behind MapleSim almost all the functions used are based or directly are Maple functions.

2) OK, may be. But the movement of the foots of the robots (elliptic trajectories) have directly a impact on the general trajectory of the robot's body. And after a moment of simulation (around 25 s), i see the robot which don't follow a straight line and I guess that it comes from the elliptic trajectories of the foots which are no longer purely elliptic.

Concerning the remarks of @Preben Alsholm :

1) I have indeed tried to use Maple first and besides the code used in MapleSim comes from my Maple's code. The only difference is the fact, in MapleSim, I no longer use the dsolve and odeplot functions but this code :

sModelica := MapleSim:-Tools:-MapleToModelica(
[seq(EqSys[i],i=1..4)],
'initial'=ic,
'class_name' = "ClassTrajectoire",
'comment' = "Trajectoire utilisée pour le mouvement des pieds",
'parameters' = params,
'connections' = cons,
'unfilled_pins' = {"u1_out", "v1_out","u2_out","v2_out","u3_out","v3_out","u4_out","v4_out"}
):

and the solver of MapleSim.

2) In Maple, I obtain limit cycles I believe. So, normally, I should obtain also limit cycles in MapleSim

Don't hesitate to let me know if you have new ideas for troubleshooting my issue of "drift" for the solution of this NL oscillator but in MapleSim.

Thanks a lot for your help