@Preben Alsholm 

@tomleslie 

@Carl Love 

Thank for your ideas on the solution of the non linear oscillator for obtaining the elliptic trajectory.

I'm still confused about my results obtained with MapleSim.

In fact, I totally agree with you about the issue of visualisation with the plot function when I don't ask enough points with Maple.

However, in MapleSim, I follow the tip of @Preben Alsholm by using the initial conditions obtained after 8s. To doing this, I launch one time the oscillator. I retrieve the values of u[i], v[i] for t=8s and I use these values as initial conditions for the next run. So, I use directly good initial conditions for solving the differential system.

Here the result I obtain with MapleSim:

It thus seems that I observe a "drift" on the movement of the robot's model. Indeed, my robot has a perfect trajectory until around 25s of simulation. Besides, until 25s, the trajectory obtained with the NL oscillator is purely elliptic. But, after this time, it begins to diverge slowly.

It may come from the solver used by MapleSim.

Do you have some ideas so as to tune the oscillator or the solver used in MapleSim so as to keep the periodic trajectory for a long time of simulation ?

I can send you the different solvers which can be used with MapleSim.

Thanks a lot for your help.

@tomleslie 

@Carl Love 

@Preben Alsholm 

Thanks a lot for your help.

You have given me a method if I want to obtain a nice visualization of the solution of my oscillator in Maple. This solution is based on a appropriated choice of numpoints with the plot function.

However, I'm wondering if there is only a "visual" mistake and not also some drifts during the simulation.

The reason for this guess is the fact that I also obtain some deviations when I generate my elliptic trajectory with MapleSim. My MapleSim code is based on my maple code. I use directly the equations from my maple code into MapleSim. However, the solver is different in the sense that I didn't use the dsolve function from maple but the solver of mapleSim.

My MapleSim code is :

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

It enables to create a dynamical system which can be solved by MapleSim solver

The settings of the MapleSim solver is :

The result obtained with MapleSim is :

1) In maple :

If there is really a drift during the resolution made by the dsolve function, do you have ideas to reduce this drift so that I can, in my case, obtain a nice ellipse for my generated trajectory?

In others words, in the previous posts, we have tried to better tune the plot function. But, before trying to correct the visualization, is it possible to correct the deviations at the "dsolve" function level in maple ?

2) In mapleSim:

Do you have ideas so as to tune properly the solver so as to avoid this drift from the elliptic trajectory ?

Thanks a lot for your help.

@tomleslie 

@Preben Alsholm 

Thanks a lot for your feedback.

However, for my application (robotics), despite the fact that the solution of the oscillator can be periodic (in some cases), I do believe that I need to obtain the solution for the duration of my simulation. I will argue with 3 reasons:

1) The oscillators which I want to use contain 2 half period (one for the stance phase and one for the support phase). We can see the definitions of 2 half periods in the second example that I gave.

2) I have the simulations conducted also with matlab. In matlab, it seems that there is no problem to ask a big tmax. Why cannot I do the same with Maple? I hope that I can tune dsolve in Maple in order to have a solution from a larger range of t.

3) These kind of oscillators are very useful in robotics because it is possible to add a sensory feedback signal (ui) in the oscillator like in the following equations :

Consequently, the oscillator is slightly pertubed to take into account some information from environment but comes back after a moment to its limit cycle. This is the kind of oscillators that I'm trying to tune in Maple.

So, I hope that my elements bring you more useful information and permit a mutual understanding. 

If yes, i will be very interested by you experience and feedback so as to tune "dsolve" (or another function) for having a simulation on a large range of time.

Do you have other ideas so as to obtain a solution for a large range of time and keeping a nice elliptical plot ?

Thank you for your help.

@Preben Alsholm 

Thank you for your help.

But, in fact, I use this non linear coupled oscillator for a robotic application. In my application, this oscillator enables me to calculate the trajectory of the foot of my robot.

Consequently, I need to make my simulation from t=0 to t=300.

So, my questions :

1) Do you have ideas so as to tune numpoints correctly ? I would like to build a rule of thumb so as to test directly a more appropriated number of point (numpoints).

2) As it is for the generation of trajectory for a robot, I'm quiet interesting to start the plotting when I have a good phase shift between my plots. Indeed, i put bad initial conditions but after 2 or 3s, the oscillator reachs its limit cycle.

May I do with this syntax : res:=dsolve({EqSys[],ic[]},numeric,range=tmin..tmax); with the definition of a tmin value ?

@tomleslie 

Perfect, thanks a lot for your very constructive and accurate help.

With this non linear oscillator, tuning the numpoints(with numpoints=10*tmax) enables to keep the solution with the limit cycle.

However, I tried to used the same tuning (that is to say to use the same number of points numpoints) in a more complex coupled oscillator and this time it doesn't work.

My code is the following :

K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);

for i to 4
do
r[i]:=sqrt((u[i](t))^2+(v[i](t))^2):
omega[i]:=omega[sw]/(1+exp(b*v[i](t)))+omega[st]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)-omega[i]*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do:

paramsCycle:=omega[st]=2*Pi,omega[sw]=4*2*Pi,Au=5,Av=5,b=100;
params:=paramsCycle;

sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));
tmax:= 100:
ic:=[u[1](0)=0, v[1](0)=0,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1];

res:=dsolve([sys[],ic[]],numeric, maxfun = 0, output = listprocedure):

plots:-odeplot(res,[u[1](t),v[1](t)],0..tmax, scaling = constrained, numpoints = 10*tmax);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..tmax,thickness=2);

However, it seems to work with numpoints = 100*tmax. 

Do you have ideas so as to tune numpoints correctly ? I would like to build a rule of thumb so as to test directly a more appropriated number of point (numpoints).

Thank you in advance for your feedback.

On the result of the non linear oscillator, I have his red arrows. They should be correspond to the general solution. Which option can I set in order to keep only my limit cycle?

Thank your for your help.

@Preben Alsholm 

Thanks a lot for your help.

I tried but I didn't obtain some results.

Is it the same for you?

And consequently, I have to find good ICs.

@Preben Alsholm 

Perfect.

You are right with the zeros on the diagonal.

And how can I do to have vj(t)

May you help me to adapt v:=Vector(4,symbol=a); so that I can have a[1](t), a[2](t), a[3](t), a[4](t)?

Thanks a lot for your help

@Rouben Rostamian  

It works!

Thank you for your very clear answer and help.

Post solved ! But, I didn't see where I can tick it!

@Rouben Rostamian  

For the plotting with different initial conditions, OK. thank you.

For the adding of points, your method should be very appropriated.

However, I have trouble with the following code line :

dsol:=dsolve([EqSys[], x(0)=0.6, z(0)=zmax], [x(t), z(t)], numeric);

I receive this error message.

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
Error, (in dsolve/numeric) x(t) and x cannot both appear in the given ODE

May you help me to adapt your proposal for adding points ?

Thank you for your help.

@Rouben Rostamian  

Thank a lot for you very constructive and efficient help !

Perfect! I can obtain the elliptic form of the nonlinear oscillator.

Otherwise, I have some slight questions :

1) Obtention of the initial conditions

I try to follow your method so as to obtain the initial conditions systematically with this line code :

RealDomain[solve](eval([rhs(eqx)=0,rhs(eqz)=0], [params]));

However, I receive this answer : 

{t = t, x = (proc (t) options operator, arrow; -4*z*RootOf((25*x^2+100*z^2)*_Z^2-1-2*_Z, label = _L2) end proc), z = (proc (t) options operator, arrow; RootOf((25*x^2+100*z^2)*_Z^2-1-2*_Z, label = _L2)*x end proc)}

How can I treat the RootOf ?

2) Plotting a point for a t given

I would like to know the position [x(t),z(t)] for a t given. How is possible to obtain and add it as a point on the curve ?

Thanks a lot for your help

@Rouben Rostamian  

Thank you but I didn't want to find this kind of results that you have shown me.

I would like to obtain the following curve that you can see below.

I correct the definition of the problem.

r:=sqrt((x(t)/a)^2+(z(t)/b)^2);
eqx:=diff(x(t),t)=alpha*(1-r^2)*x+w*a/b*z(t);
eqz:=diff(z(t),t)=beta*(1-r^2)*z+w*b/a*x(t);
EqSys:=[eqx,eqz];

Numeric Data:

alpha:=1:
beta:=1:
a=0.4:
b=0.2:
w=1:

I try to plot by using DEplot function in this manner :

DEplot(EqSys, [x(t), z(t)], t = -8 .. 8, [[x(0) = 0, z(0) = 0]])

But, for the moment, I receive this error :

Warning, x is present as both a dependent variable and a name. Inconsistent specification of the dependent variable is deprecated, and it is assumed that the name is being used in place of the dependent variable.
Warning, z is present as both a dependent variable and a name. Inconsistent specification of the dependent variable is deprecated, and it is assumed that the name is being used in place of the dependent variable.

As the initial conditions are not mentioned in the paper, I'm trying with x(0)=0 and z(0)=0 for the moment.

I hope this precisions will enable you to better identify the nature of the differential system that I want to plot. This system is called Hopf non linear oscillator if I have well understood.

May you help me to plot this kind of system on Maple ?

Thank you for your help

You are completely right about the multiply that I forget. Is is indeed w*a...

About the initial conditions they are not specify in the paper that I read. I wonder what I can put.

Do you have ideas ? I would like to obtain zlliptic limit cycle.

Thank you for your help

@Thomas Richard 

For the point 2, I mean that i would like for example to try the GetLGM module. But, this module is not detailed in the Kinematic Exports and Dynamic Exports help pages.

So, how can I do to understand for the different functions of the module GetLGM the inputs that I have to use for these functions.

For this topic, in order to be clearer, I added a new post.

Thank you for your help.

@Thomas Richard 

In other words,

For 1) , can I use these functions : 

SetOptions, SetAdditionalKinData, SetAdditionalDynData, PostProcessDyEQs, DFPModule, Print, vPd, vPd_dot, xJi, xJdInv, xG, vDepVels, vDepAccels, ... which are not seen with Export function function but the line code Code:for k to 8 do if (op(k,op(2,eval(MB)))=NULL) then print("NULL") else print(op(k,op(2,eval(MB)))) fi od; ? 

For 2), OK