@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