Generating codim-1 bifurcation diagram of sliding friction oscillator

I am performing continuation analysis in MATCONT and I am finding it difficult to generate co-dim bifurcation diagram.

Here is my linearized system of equations which represents a 3-dof clutch dynamics with sliding type friction:

x1' = x2

x2' = -(de+Fn*mukin)*x2/je + (Fn*mukin*x4)/je - ((x2-x4)*mukin + mus)/je

x3' = x4

x4' = (Fn*mukin*x2)/jd - (ds+Fn*mukin)*x4/jd +ds*x6/jd + ks*x7/jd

x5' = x6

x6' = ds*x4/jv -ds*x6/jv - ks*x7/jv

x7' = x6 - x4

This system dynamics represents clutch dynamics with the sliding friction between inertias je and jd. jv is the load inertia. ks is stiffness between jd and jv. de,ds is damping coefficient. Fn is normal force acting on jd and je during engagement of clutch. mukin is slope of friction curve usually between [-0.01 to 0.01] and mus is static coefficient of friction 0.136 for this analysis.

je = 1.57

jd = 0.0066

jv = 4

ks = 56665.5

ds,de = 0.1

Fn = 8000

Please guide me how to perform the continuation of this system in MATCONT.

Answer 1

[Umar]

Hi @Junaid Ali ,

You asked, “Please guide me how to perform the continuation of this system in MATCONT. “

Please see my response to your comments below.

Here’s how you can set up and perform the continuation analysis using MATCONT. So, I can assume that you have MATCONT already installed and added to your MATLAB path. If not, please click the following links for more information.

https://www.mathworks.com/matlabcentral/answers/1883987-how-to-load-matcont-data-file-and-plot

https://www.mathworks.com/matlabcentral/answers/2148254-matcont-load-state-space-a-b-c-d-matrices-for-continuation-analysis

https://www.mathworks.com/matlabcentral/answers/99328-how-do-i-create-a-bifurcation-diagram-in-matlab

https://www.mathworks.com/matlabcentral/answers/770418-how-to-plot-bifurcation-diagram-in-matlab

https://dercole.faculty.polimi.it/tds/matcont.html

You have provided the following system of equations as shown in your comments.

function dx = clutchDynamics(t, x, params)
  % Unpack parameters
  de = params.de;
  ds = params.ds;
  Fn = params.Fn;
  mukin = params.mukin;
  mus = params.mus;
  je = params.je;
  jd = params.jd;
  jv = params.jv;
  ks = params.ks;

    % State variables
    x1 = x(1);
    x2 = x(2);
    x3 = x(3);
    x4 = x(4);
    x5 = x(5);
    x6 = x(6);
    x7 = x(7);

    % Define the system of equations
    dx = zeros(7, 1);
    dx(1) = x2;
    dx(2) = -(de + Fn * mukin) * x2 / je + (Fn * mukin * x4) / je - ((x2 - x4) * 
    mukin + mus) / je;
    dx(3) = x4;
    dx(4) = (Fn * mukin * x2) / jd - (ds + Fn * mukin) * x4 / jd + ds * x6 / jd +
    ks * x7 / jd;
    dx(5) = x6;
    dx(6) = ds * x4 / jv - ds * x6 / jv - ks * x7 / jv;
    dx(7) = x6 - x4;
  end

You also defined the parameters to be used in your system. I would suggest to implement them in a separate MATLAB script or function:

params.de = 0.1;
params.ds = 0.1;
params.Fn = 8000;
params.mukin = linspace(-0.01, 0.01, 100); % Varying mukin for continuation
params.mus = 0.136;
params.je = 1.57;
params.jd = 0.0066;
params.jv = 4;
params.ks = 56665.5;

To use MATCONT, you need to set up the initial conditions and the continuation parameters. Here’s how you can do it:You need to find an equilibrium point of your system which can be done by solving the system of equations for a specific value of mukin and then initialize MATCONT with the equilibrium point and specify the parameters for continuation. If you are not familiar with MATCONT commands below, I will suggest to review them by clicking the link below.

https://webspace.science.uu.nl/~kouzn101/NBA/ManualMatcontMNov2017.pdf

% Initial conditions (example)
x0 = [0; 0; 0; 0; 0; 0; 0]; % Initial guess for the state variables
mukin0 = 0; % Initial value of mukin

% Define the continuation parameters
options = contset; % Create a default options structure
options = contset(options, 'MaxNumPoints', 100); % Maximum number of 
points
options = contset(options, 'Singularities', 1); % Track singularities

% Start continuation
[x, v, s] = cont(@clutchDynamics, x0, mukin0, options, params);

Once you have performed the continuation, now you can plot the bifurcation diagram which can be done using MATLAB,

figure;
plot(x(1, :), x(2, :)); % Plotting x1 vs x2 as an example
xlabel('x1');
ylabel('x2');
title('Bifurcation Diagram of 3-DOF Clutch Dynamics');
grid on;

Now, to analyze co-dimension bifurcations, you may need to vary more than one parameter. My suggestion would be setting up a second parameter in a similar manner to mukin and use MATCONT to track bifurcations in a multi-parameter space.

• Additional Suggestions*

In my opinion, it may be beneficial to perform a sensitivity analysis on the parameters such as mukin, mus to understand their influence on system behavior. Also, depending on your results, consider investigating different types of bifurcations such as saddle-node or Hopf bifurcations, which can provide deeper insights into system stability and a linear stability analysis around the bifurcation points to deduce the nature of the bifurcations.

Hope this should help you get started with your project now.

[Junaid Ali]

I couldn't find an equilibrium value to begin my codim-1 continuation with parameter "mus". I fixed muk = 0.001

I am basically trying to re-create these results from this paper.

DOI 10.1007/s11071-006-9170-5

I am expected to see a subcritical Hopf bifurcation when system becomes unstable but my continuation is never converging. I am not sure what I am doing wrong. I follow the steps you showed here. I even tried fixing situation with CHATGPT but it couldnt solve it.

[Sam Chak]

Hi @Junaid Ali, The ariticle is behind the paywall. Not sure if @Umar can access it or not. Conside posting only the expected results (a figure) on the basis of "Fair Use".

[Junaid Ali]

Here is the MATLAB CODE I wrote for this analysis. Next step I am struggling to go is perform continuation to generate this Figure. I am having Coding skill issue to be precise. I am not familiar with MATCONT so even the errors are not clear to me what is wrong.

clc;

clear;

% Define parameters

params.Teng = 100; % Engine torque

params.mu0 = 0.17; % Static friction coefficient

params.mu1 = 0.13; % Dynamic friction coefficient

params.gamma = 50; % Friction parameter

params.Rm = 0.15; % Friction radius

params.Ap = 0.06; % Friction contact area

params.P0 = 100; % Initial pressure

params.Ki = 10; % Pressure gain

params.Beng = 0.05; % Engine damping coefficient

params.Bcnv = 2.5; % Converter damping coefficient

params.Bcl = 0.05; % Clutch damping coefficient

params.Bd = 0.02; % Driveshaft damping coefficient

params.kd = 1000; % Driveshaft stiffness

params.Bt = 0.05; % Transmission damping coefficient

params.Td = 60; % Transmission torque

params.Ieng = 0.3; % Engine inertia

params.Icl = 0.08; % Clutch inertia

params.It = 0.15; % Transmission inertia

params.R = 1; % Reference input

% mu1 = -0.05:0.0015:0.4;

% R = 2:0.05:10; %rad/sec

% Initial guess for the equilibrium point (can be adjusted based on the system)

x0 = [0; 0; 0; 0; 0];

% Find the equilibrium using fsolve

options = optimoptions('fsolve', 'Display', 'iter', 'TolFun', 1e-9, 'TolX', 1e-9);

[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);

% Display the equilibrium point

if exitflag > 0

disp('Equilibrium found:')

disp(x_eq)

else

error('Failed to find equilibrium.')

end

% Define the continuation parameters

options = contset; % Create a default options structure

options = contset(options, 'MaxNumPoints', 100); % Maximum number of points

options = contset(options, 'Singularities', 1); % Track singularities

mu1 = -0.05:0.0015:0.4;

% Start continuation

[x, v, s] = cont(@system, x0, mu1, options, params);

% Define a function for equilibrium (dx/dt = 0)

function eq = equilibrium_function(x, params)

dxdt = system(0, x, params);

eq = dxdt; % At equilibrium, dx/dt = 0

end

function dxdt = system(t,x,params)

%Unpack parameters

Teng = params.Teng;

mu0 = params.mu0;

mu1 = params.mu1;

gamma = params.gamma;

Rm = params.Rm;

Ap = params.Ap;

P0 = params.P0;

Ki = params.Ki;

Beng = params.Beng;

Bcnv = params.Bcnv;

Bcl = params.Bcl;

Bd = params.Bd;

kd = params.kd;

Bt = params.Bt;

Td = params.Td;

Ieng = params.Ieng;

Icl = params.Icl;

It = params.It;

R = params.R;

% Unpack state variables

x1 = x(1);

x2 = x(2);

x3 = x(3);

x4 = x(4);

x5 = x(5);

% Define the system of equations

dx1dt = (Teng - (mu0 + mu1*exp(-gamma*Rm*(x1-x2)))*Rm*Ap*(P0+Ki*x5) - Beng*x1 - Bcnv*(x1-x2)) / Ieng;

dx2dt = ((mu0 + mu1*exp(-gamma*Rm*(x1-x2)))*Rm*Ap*(P0+Ki*x5) - Bcl*x2 + Bcnv*(x1 - x2) - Bd*(x2 - x3) - kd*x4) / Icl;

dx3dt = (Bd*(x2 - x3) + kd*x4 - Bt*x3 - Td) / It;

dx4dt = x2 - x3;

dx5dt = R - (x1 - x2);

% Return the vector of derivatives

dxdt = [dx1dt; dx2dt; dx3dt; dx4dt; dx5dt];

end

[Umar]

Hi @Junaid Ali ,

After executing your code, it displayed an error that requires MATCONT. Please see attached.

Although you did mention that you did mention,

“I am not familiar with MATCONT so even the errors are not clear to me what is wrong.”

In order to achieve your task, you have to learn how to use this tool. It requires hardworking. I read the synopsis of article mentioned your comments by you.

https://link.springer.com/article/10.1007/s11071-006-9170-5

It does require sophisticated bifurcation analysis in the context of torque converters, which you would only able to achieve using “MATCONT”, understanding these bifurcations will help you to optimize clutch performance by adjusting friction parameters, such as the static and low-velocity friction constants and employing a bifurcation-analysis-based approach like you are trying to implement in code will help identify conditions under which undesirable oscillations, like stick-slip, may occur. Trust me, this knowledge is invaluable for enhancing the efficiency and stability of the clutch system. For instance, the study indicates that as slip speeds increase, the system can tolerate higher friction values before instability arises, which can lead to improved performance in real-world applications.

So, you have two options to choose from the following below.

Option A: Install the continuation toolbox (e.g., MATCONT) is installed and added to your MATLAB path. You can find MATCONT and download it from its official website or GitHub repository.

https://dercole.faculty.polimi.it/tds/matcont.html

Option B: If you cannot install the toolbox, consider an alternative method to implement continuation without contset. You can manually implement a parameter sweep or use a different numerical method to analyze the system's behavior or you can implement it.

The downside for bifurcation diagram in matlab is that there is no function defined for it, please click the link below

https://www.mathworks.com/matlabcentral/answers/99328-how-do-i-create-a-bifurcation-diagram-in-matlab

The person who is expert in dynamic systems is @Torsten who can really tackle this problem.

Hope this helps.

[Sam Chak]

Hi @Umar

If I edit your post above, the formatting will be lost. Therefore, I will simply tag the individuals you wish to engage with: @Junaid Ali, @Torsten.

[Umar]

Hi @Sam Chak,

I would appreciate if @Torsten will engage to help out.

[Junaid Ali]

Hi @Umar,

Thank you for your feedback. I am using MATCONT as well and I have been able to generate bifurcation diagrams for solved examples in tutorials. I tried some cooked up problems as well to practice. This one problem I related to friction stumped me. So i am trying everything in my capacity to recreate the results of this paper to gain a good grasp on continuation analysis.

Here is my system of equations in MATCONT:

STEP 1:

STEP 2: Inputing system Parameter to Start Integration from a [0,0,0,0,0,0] initial condition.

STEP 3: Selecting Point to begin integration to find equilibrium point for continuation analysis.

STEP 4: I am stuck here!

My time response doesnt get to an equilibrium point.

This is where I need help. To get this time response stabilize to a equilbirum point or limit cycle.

Also, @Torsten I would really appreciate your help on this.

[Umar]

Hi @Junaid Ali,

It sounds like you made pretty good progress. After going through your comments, please review my suggestions below.

You mentioned starting from an initial condition of [0, 0, 0, 0, 0]. While this is a common starting point, it may not be suitable for your specific system, especially if the system is stiff or has nonlinear characteristics. Consider using small perturbations around zero, such as:

initial_conditions = [0.01, 0.01, 0.01, 0.01, 0.01];

To find equilibrium points, you need to make sure that your system is correctly set up in MATCONT. Based on your screenshot shared in step 3, the point for continuation analysis is setup right such as system of ODEs, you are using SSNNN, curve type, also you have chosen Point (P), which is appropriate for finding equilibrium points. However, Orbit (0) type is suitable for periodic solutions, but if you are looking for equilibrium, consider using Point instead.

Now, if your time response does not stabilize, consider the following strategies:

Adjust the Time Step: Sometimes, the integration step size can affect convergence. Try reducing the time step in your integration settings.

Use a Different Solver if possible: MATCONT allows you to choose different solvers. If you are using a default solver, experiment with others to see if they yield better results.

Check for Stiffness: If your system is stiff, consider using a stiff solver. This can significantly improve convergence to equilibrium points.

Once you have adjusted the parameters and initial conditions, run the simulation again. Use the following code snippet to visualize the results:

% Plotting the results
figure;
plot(t, x1, 'r', t, x2, 'g', t, x3, 'b', t, x4, 'm', t, x5, 'c');
xlabel('Time');
ylabel('State Variables');
title('Time Response of the System');
legend('x1', 'x2', 'x3', 'x4', 'x5');
grid on;

Remember, the key is to iteratively adjust your parameters and initial conditions while monitoring the system's response. If you continue to face challenges, consider revisiting the theoretical foundations of your model or consulting additional resources on MATCONT for further insights. I will wait for @Torsten feedback to find out what he has to say.

[Junaid Ali]

Hello @Umar,

Thank you for your feedback. It helped me. I progressed to codim-1 bifurcation. But I cheated my way to it.

I manually fed the equilibrium point values in MATCONT from the table given in paper.

Then I did continuation of mu1 parameter and I got Neutral Saddle Equilibrium

I selected this point for codim-2 bifurcation to detect subcritical Hopf point. Selected R and mu1 as continuation parameters. I got error which says No convergence at x0.

[Umar]

Hi @Junaid Ali,

You mentioned, “I selected this point for codim-2 bifurcation to detect subcritical Hopf point. Selected R and mu1 as continuation parameters. I got error which says No convergence at x0. “

Please see my response to your comments below.

The error you are encountering typically arises when the algorithm used in MATCONT cannot find a solution near the initial guess x0 which could be due to several factors such as: the choice of x0 may not be close enough to a valid equilibrium point. It’s crucial that x0 should be chosen by you within the basin of attraction for the equilibrium you are trying to analyze. Also, the parameters selected like R and mu1 might lead to a situation where no equilibrium exists or where the system is unstable. You should be aware that sometimes numerical methods can fail due to precision issues or ill-conditioning in the system of equations being solved.

Please see my recommendations for fixing this error below

Refine Initial Guesses

Visualize Phase Portraits: Before selecting x0, consider plotting phase portraits or using bifurcation diagrams to better understand where equilibria lie.

Use Nearby Equilibria: Start with values close to known equilibria from your previous analyses (such as those from your subcritical Hopf bifurcation).

Adjust Continuation Parameters

Make sure that your continuation parameters R and mu1 are chosen such that they transition smoothly through parameter space. If one parameter leads to instability, it may be beneficial to adjust its range or step size.

Check Stability Criteria

Investigate the stability of your Neutral Saddle Equilibrium. If it is unstable, it may lead to difficulties in finding converging solutions. Use Lyapunov coefficients or other stability analysis techniques to confirm.

MATCONT Settings

Adjust MATCONT settings for convergence criteria, such as increasing tolerance levels or changing the method of continuation (e.g., from Newton's method to another suitable algorithm). Also, review MATCONT documentation for advanced options that might assist in improving convergence.

Incremental Continuation

Instead of jumping directly into a codim-2 bifurcation analysis, consider performing an incremental continuation from your current stable points and gradually adjusting your parameters.

You are not the only one who is facing convergence issues, I have studied and observed similar convergence issues like yours which have been documented in studies involving complex dynamical systems, where researchers often had to iteratively refine their initial conditions based on insights gained from phase plots and stability analyses. Hopefully, by refining your approach based on these recommendations, you should be able to overcome the convergence issue you're facing and successfully analyze the codim-2 bifurcation in your system.

[Junaid Ali]

Thank you @Umar, the advice and help I recieved here worked for me. And I was able to recreate the results from the paper.

Now I can move on to my actual system and see if I can implement the similar methodolgy to stuy bifurcation behaviour.

thank you again for your invaluable guidance.

[Umar]

Hi @Junaid Ali,

Thank you for your kind words. I am delighted to hear that the advice and assistance provided were beneficial in recreating the results from the paper.

It’s exciting to know that you are now able to apply similar methodologies to study bifurcation behavior in your actual system. I have no doubt that your efforts will yield valuable insights. Should you need any further guidance or support as you progress, please do not hesitate to reach out.

Wishing you the best of luck with your research.

[Junaid Ali]

Hello @Umar and @Torsten,

So I coded the equations of motion of my system and trying to obtain equilibrium point (code attached).

This code worked for system of equations given in research article and produced exact same equilibrium points.

I am using same code logic to generate equilibrium points for my system which represents the sliding friction clutch dynamics in a 4 DOF model.

I am starting my initial guess at:

1. When x0 = zeros(9,1) but no equilibirum find because exit flag variable is zero or negative
2. When x0 = rand(9,1) but no equilibirum find because exit flag variable is zero or negative

Although the x_eq vector returns some values but I am not sure what are those values. Can you help me find out what could possibly be wrong in this code or my approach to find the equilibrium points of my system.

My equations of motion are correctly coded and time-domain simulations capture the stick-slip dynamic behaviour.

clc;
clear;

params.J1 = 1.49;
params.J2 = 0.0764;
params.J3 = 0.0066;
params.J4 = 4;

params.k1 = 6498.3;
params.k2 = 56665.5;

params.c1 = 0.3;
params.c2 = 0.02;

params.Teng = 100;
params.Tveh = 100;

params.R = 0;

params.mu0 = 0.136;
params.muk = 0.001;
params.np = 14;
params.Rm  = 0.0506;
params.Fn = 8000;

x0 = rand(9,1);

% Find the equilibrium using fsolve
options = optimoptions('fsolve', 'Display', 'iter', 'TolFun', 1e-9, 'TolX', 1e-9);
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);

                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0         10         2.35982e+11                        4.14e+12              1
     1         20         3.30022e+09      0.0397168         4.89e+09              1
     2         30         4.09776e+07        0.46498         5.49e+10              1
     3         40             37482.9    0.000526781         1.07e+07           1.16
     4         50             21732.4     0.00101501         1.91e+07           1.16
     5         60             21727.4    1.82888e-07         1.51e+05           1.16
     6         70             21722.5    2.24837e-05         1.91e+07           1.16
     7         80             21717.6    1.82835e-07         1.51e+05           1.16
     8         90             21712.6    2.24774e-05         1.91e+07           1.16
     9        100             21707.7    1.82794e-07         1.51e+05           1.16
    10        110             21702.8    2.24723e-05         1.91e+07           1.16
    11        120             21697.9    1.82755e-07         1.51e+05           1.16
    12        130             21692.9    2.24681e-05         1.91e+07           1.16
    13        140               21688    1.82719e-07         1.51e+05           1.16
    14        150             21683.1    2.24639e-05         1.91e+07           1.16
    15        160             21678.2    1.82683e-07         1.51e+05           1.16
    16        170             21673.2    2.24595e-05          1.9e+07           1.16
    17        180             21668.3    1.82646e-07         1.51e+05           1.16
    18        190             21663.4    2.24553e-05          1.9e+07           1.16
    19        200             21658.5     1.8261e-07         1.51e+05           1.16
    20        210             21653.6    2.24511e-05          1.9e+07           1.16
    21        220             21648.6    1.82574e-07         1.51e+05           1.16
    22        230             21643.7    2.24467e-05          1.9e+07           1.16
    23        240             21638.8    1.82536e-07         1.51e+05           1.16
    24        250             21633.9    2.24426e-05          1.9e+07           1.16
    25        260               21629    1.82501e-07         1.51e+05           1.16
    26        270             21624.1    2.24382e-05          1.9e+07           1.16
    27        280             21619.2    1.82464e-07         1.51e+05           1.16
    28        290             21614.3    2.24338e-05          1.9e+07           1.16
    29        300             21609.3    1.82426e-07         1.51e+05           1.16
    30        310             21604.4    2.24298e-05          1.9e+07           1.16
    31        320             21599.5    1.82391e-07         1.51e+05           1.16
    32        330             21594.6    2.24258e-05          1.9e+07           1.16
    33        340             21589.7    1.82355e-07         1.51e+05           1.16
    34        350             21584.8    2.24213e-05          1.9e+07           1.16
    35        360             21579.9    1.82318e-07         1.51e+05           1.16
    36        370               21575     2.2417e-05          1.9e+07           1.16
    37        380             21570.1    1.82281e-07         1.51e+05           1.16
    38        390             21565.2    2.24127e-05          1.9e+07           1.16
    39        400             21560.3    1.82245e-07         1.51e+05           1.16
    40        410             21555.4    2.24088e-05          1.9e+07           1.16
    41        420             21550.5    1.82209e-07         1.51e+05           1.16
    42        430             21545.6    2.24046e-05          1.9e+07           1.16
    43        440             21540.7    1.82173e-07          1.5e+05           1.16
    44        450             21535.8    2.24003e-05          1.9e+07           1.16
    45        460             21530.9    1.82136e-07          1.5e+05           1.16
    46        470               21526    2.23961e-05          1.9e+07           1.16
    47        480             21521.1      1.821e-07          1.5e+05           1.16
    48        490             21516.3    2.23916e-05          1.9e+07           1.16
    49        500             21511.4    1.82063e-07          1.5e+05           1.16
    50        510             21506.5    2.23876e-05          1.9e+07           1.16
    51        520             21501.6    1.82027e-07          1.5e+05           1.16
    52        530             21496.7    2.23835e-05          1.9e+07           1.16
    53        540             21491.8    1.81991e-07          1.5e+05           1.16
    54        550             21486.9    2.23791e-05          1.9e+07           1.16
    55        560               21482    1.81954e-07          1.5e+05           1.16
    56        570             21477.2    2.23747e-05          1.9e+07           1.16
    57        580             21472.3    1.81917e-07          1.5e+05           1.16
    58        590             21467.4    2.23703e-05          1.9e+07           1.16
    59        600             21462.5     1.8188e-07          1.5e+05           1.16
    60        610             21457.6     2.2366e-05          1.9e+07           1.16
    61        620             21452.8    1.81843e-07          1.5e+05           1.16
    62        630             21447.9    2.23619e-05          1.9e+07           1.16
    63        640               21443    1.81808e-07          1.5e+05           1.16
    64        650             21438.1    2.23577e-05          1.9e+07           1.16
    65        660             21433.3    1.81771e-07          1.5e+05           1.16
    66        670             21428.4    2.23534e-05          1.9e+07           1.16
    67        680             21423.5    1.81734e-07          1.5e+05           1.16
    68        690             21418.6    2.23493e-05          1.9e+07           1.16
    69        700             21413.8    1.81699e-07          1.5e+05           1.16
    70        710             21408.9    2.23451e-05         1.89e+07           1.16
    71        720               21404    1.81662e-07          1.5e+05           1.16
    72        730             21399.2     2.2341e-05         1.89e+07           1.16
    73        740             21394.3    1.81626e-07          1.5e+05           1.16
    74        750             21389.4    2.23367e-05         1.89e+07           1.16
    75        760             21384.6     1.8159e-07          1.5e+05           1.16
    76        770             21379.7    2.23328e-05         1.89e+07           1.16
    77        780             21374.8    1.81555e-07          1.5e+05           1.16
    78        790               21370    2.23285e-05         1.89e+07           1.16
    79        800             21365.1    1.81518e-07          1.5e+05           1.16
    80        810             21360.3    2.23241e-05         1.89e+07           1.16
    81        820             21355.4    1.81481e-07          1.5e+05           1.16
    82        830             21350.5    2.23201e-05         1.89e+07           1.16
    83        840             21345.7    1.81446e-07          1.5e+05           1.16
    84        850             21340.8    2.23158e-05         1.89e+07           1.16
    85        860               21336    1.81409e-07          1.5e+05           1.16
    86        870             21331.1    2.23118e-05         1.89e+07           1.16
    87        880             21326.3    1.81374e-07          1.5e+05           1.16
    88        890             21321.4    2.23078e-05         1.89e+07           1.16
    89        900             21316.6    1.81338e-07          1.5e+05           1.16

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 9.000000e+02.

disp(x_eq)

    0.5553
    0.5960
    0.5713
    0.4875
    0.6212
    0.5555
    0.6231
    0.3556
    0.1318

% Display the equilibrium point
if exitflag > 0
    disp('Equilibrium found:')
    disp(x_eq)
else
    error('Failed to find equilibrium.')
end

Failed to find equilibrium.

%% 

function eq = equilibrium_function(x, params)

    dxdt = clutchDynamics(0, x, params);
    eq = dxdt;  % At equilibrium, dx/dt = 0
end

function dydt = clutchDynamics(t, y, params)
 
        J1 = params.J1;
        J2 = params.J2;
        J3 = params.J3;
        J4 = params.J4;

        k1 = params.k1;
        k2 = params.k2;

        c1 = params.c1;
        c2 = params.c2;

        Teng = params.Teng;
        Tveh = params.Tveh;

        R = params.R;
        mu0 = params.mu0;
        muk = params.muk;
        np = params.np;
        Rm = params.Rm;
        Fn = params.Fn;

        dydt = zeros(9,1); %9th row is slip speed
        
        Vslip = y(4) - y(6);
        
        Tclutch = np*Rm*Fn*(mu0 + muk*Vslip)*tanh(2*Vslip);

        dydt(1) = y(2);
        dydt(2) = (Teng - k1*(y(1) - y(3)) - c1*(y(2) - y(4)))/J1;
        dydt(3) = y(4);
        dydt(4) = (k1*(y(1) - y(3)) + c1*(y(2) - y(4)) - Tclutch)/J2;
        dydt(5) = y(6);
        dydt(6) = (-k2*(y(5) -y(7)) - c2*(y(6) - y(8)) + Tclutch)/J3;
        dydt(7) = y(8);
        dydt(8) = (k2*(y(5) - y(7)) + c2*(y(6) - y(8)) - Tveh)/J4;
        dydt(9) = R - Vslip;
  end

[Umar]

Hi @Junaid Ali ,

You mentioned, “So I coded the equations of motion of my system and trying to obtain equilibrium point (code attached). This code worked for system of equations given in research article and produced exact same equilibrium points. I am using same code logic to generate equilibrium points for my system which represents the sliding friction clutch dynamics in a 4 DOF model. I am starting my initial guess at: 1. When x0 = zeros(9,1) but no equilibirum find because exit flag variable is zero or negative 2. When x0 = rand(9,1) but no equilibirum find because exit flag variable is zero or negative Although the x_eq vector returns some values but I am not sure what are those values. Can you help me find out what could possibly be wrong in this code or my approach to find the equilibrium points of my system. “

Please see my response to your comments below.

So, the objective of using fsolve is to find values for the state variables such that the system reaches an equilibrium state where all derivatives dx/dt are zero. In this context, it means determining the conditions under which the forces and moments in the clutch balance out. So, I debugged your code to get insights why is not converging. So, based on your output from fsolve which includes several key components which you included in your code to obtain insights after you found out it was not converging.

*Exit Flag: Indicates whether fsolve successfully found a solution.

*Positive values (e.g., exitflag > 0) suggest successful convergence.

*Negative values indicate failure or issues during computation.

So, in addition to your debugging process, I also added additional debugging code such as function value which represents how close the solution is to achieving zero derivatives; ideally, this should approach zero and printing out current dx/dt values which are printed during iterations to show how close each state variable is to equilibrium at various points.

So, based on results of debugging, several exit flags and corresponding function values are shown such as an exit flag of `0` indicating that fsolve stopped because it reached the maximum number of function evaluations without converging. Also, flags like -24.9870 and others suggest that while some values are approaching equilibrium (as seen in small function values), convergence was not fully achieved within set limits. If you look at the attached screenshot as shown below

% Parameters
params.J1 = 1.49;
params.J2 = 0.0764;
params.J3 = 0.0066;
params.J4 = 4;
params.k1 = 6498.3;
params.k2 = 56665.5;
params.c1 = 0.3;
params.c2 = 0.02;
params.Teng = 100;
params.Tveh = 100;
params.R = 0; % Accessed in clutchDynamics
params.mu0 = 0.136;
params.muk = 0.001;
params.np = 14;
params.Rm = 0.0506;
params.Fn = 8000;

% Adjusted initial guess based on expected physical behavior
x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0]; % Starting from zero might help

% Find the equilibrium using fsolve
options = optimoptions('fsolve', 'Display', 'iter', ...
                     'TolFun', 1e-8, 'TolX', 1e-8, ...
                     'MaxFunctionEvaluations', 200000, ...
                     'MaxIterations', 20000); % Increased limits                         for better convergence

[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,   params), x0, options);

% Display results
if exitflag > 0
  disp('Equilibrium found:');
  disp(x_eq);
else
  fprintf('Failed to find equilibrium. Exit flag: %d, Function 
  value: %f\n', exitflag, fval);
end

%% Equilibrium Function
function eq = equilibrium_function(x, params)
  dxdt = clutchDynamics(0, x, params);
  disp('Current dxdt:'); % Debugging output
  disp(dxdt);
  eq = dxdt; % At equilibrium, dx/dt should be zero
end

%% Clutch Dynamics Function
function dydt = clutchDynamics(t, y, params)
  % Extract parameters from params structure
  J1 = params.J1; J2 = params.J2; J3 = params.J3; J4 = params.J4; 
  k1 = params.k1; k2 = params.k2; c1 = params.c1; c2 = params.c2; 
  Teng = params.Teng; Tveh = params.Tveh; Rm = params.Rm; 
  mu0 = params.mu0; muk = params.muk; np = params.np; Fn =     params.Fn;
  R = params.R;

    dydt = zeros(9,1); % Initialize output

    Vslip = y(4) - y(6); % Calculate slip velocity
    Tclutch = np * Rm * Fn * (mu0 + muk * Vslip) * tanh(2 * Vslip);

    % Define the dynamics equations
    dydt(1) = y(2);
    dydt(2) = (Teng - k1*(y(1) - y(3)) - c1*(y(2) - y(4))) / J1;
    dydt(3) = y(4);
    dydt(4) = (k1*(y(1) - y(3)) + c1*(y(2) - y(4)) - Tclutch) / J2;
    dydt(5) = y(6);
    dydt(6) = (-k2*(y(5) - y(7)) - c2*(y(6) - y(8)) + Tclutch) /     J3;
    dydt(7) = y(8);
    dydt(8) = (k2*(y(5) - y(7)) + c2*(y(6) - y(8)) - Tveh) / J4;
    dydt(9) = R - Vslip; % Ensure slip speed is accounted for
  end

you will see display of Current dxdt which shows that while some variables have stabilized (for example, velocities nearing zero), others are still changing, indicating ongoing dynamics that prevent reaching true equilibrium.

So, based on these analysis, I would suggest breaking down the problem into simpler parts or analyzing specific components of your dynamics equations may provide clarity on which aspects are causing divergence.

[Junaid Ali]

@Umar

Even, the time response I plotted looks like it is achieving steady-state. I dont understand why it is not considering it as equilibrium point.

I fed the final dxdt value to MATCONT as an equilibrium point to begin contininuation of "muk" parameter. But it gave me error that no convergence for x0.

[Umar]

Hi @Junaid Ali,

Please bear in mind that an equilibrium point is defined as a state where the system's derivatives (like dxdt) are zero, meaning there is no change over time. Based on your attached figure you shared, it looks like your time response resembles a damped sine wave which suggests that the system is oscillating and gradually settling down, but it may not have reached a true equilibrium. To find the equilibrium point, make sure that you are looking for conditions where dxdt = 0.

Regarding the MATCONT error, it typically arises when the initial guess for the equilibrium point (x0) does not satisfy the necessary conditions for convergence. You may need to refine your initial guess or make sure that it corresponds to a true equilibrium state.

Hope this helps.

[Junaid Ali]

Thank you @Umar,

I understand and being very careful about this in my analysis. I noticed one key difference in my system of equations and the reference paper I was looking at. The reference paper uses Closed-Loop Feedback Clutch System and I am using Open-Loop Clutch System.

When I made my system Closed-Loop Feedback, I got equilibrium point convergence too. Now I need to wrap my head around why I need to close the feedback loop to achieve equilibrium.

[Umar]

Hi @Junaid Ali,

If the reference paper was using closed loop system then it makes perfect sense from debugging your code and providing analysis about current dxdt showing some variables have stabilized (for example, velocities nearing zero), others are still changing, indicating ongoing dynamics that prevent reaching true equilibrium. Also, transitioning to a closed-loop feedback clutch system will likely provide the necessary adjustments to counteract any deviations from the equilibrium state because feedback loop effectively minimizes errors by dynamically altering the system's behavior, thus facilitating stability and convergence and this is exactly what I was expecting. Since you observed differences between your model and the reference paper (closed-loop vs open-loop), it may be beneficial to analyze how these differences influence system dynamics and equilibrium behavior systematically.

[Umar]

Hi @Junaid Ali,

Please let me know if you have any further questions or need further assistance.

[Junaid Ali]

Hi @Umar,

What a coincidence I was thinking about reaching out to you again as I am facing some difficulties. The new issue I am stuck in with my continuation analysis is not finding any bifurcation point for any of my system parameter sweep.

I know for sure I should see sub-critical hopf bifurcation in my system but there isnt any.

I dont know what to do now.

[Umar]

Hi @Junaid Ali,

Based on your recent comments, you made some progress by mentioning, “ When I made my system Closed-Loop Feedback, I got equilibrium point convergence too. Now I need to wrap my head around why I need to close the feedback loop to achieve equilibrium. ”

It seems that you are facing additional challenges with your code now, is it possible to share the code, so I can analyze and help you assist better.

[Junaid Ali]

Step 1: Here is the matlab code for finding equilibrium position:

clc;
clear;

params.je = 1.57; 
params.jd = 0.0066; 
params.jv = 4;

params.ks = 56665.5;

params.de = 0.01;
params.ds = 0.01;

params.Te = 100;
params.Tv = 80;

params.mu0 = 0.11;
params.muk = 0.056;

params.np = 14;
params.Rm = 0.0506;
params.Fn = 15000;
params.Kg = 1000;
params.R = 50;

x0 = zeros(5,1);
ts = 0:0.001:1;
[t,z] = ode45(@(t,z)clutch3dof(t,z,params), ts, x0);

figure
plot(ts,z)



options = optimoptions('fsolve', 'Display', 'iter', 'TolFun', 1e-8, 'TolX', 1e-8, 'MaxFunctionEvaluations', 10000, 'MaxIterations', 10000);
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x, params), x0, options);
if exitflag > 0
    disp('Equilibrium found:')
    disp(x_eq)
else
    warning('Failed to find equilibrium.')
end


%%

function eq = equilibrium_function(x, params)

    dxdt = clutch3dof(0, x, params);
    eq = dxdt;

end


function dydt = clutch3dof(t,y,params)
    je = params.je; 
    jd = params.jd; 
    jv = params.jv;
    
    ks = params.ks;
    
    de = params.de;
    ds = params.ds;
    
    Te = params.Te;
    Tv = params.Tv;
    
    mu0 = params.mu0;
    muk = params.muk;
    
    np = params.np;
    Rm = params.Rm;
    Fn = params.Fn;
    R = params.R;
    Kg = params.Kg; 

    dydt = zeros(5,1);

    dydt(1) = (Te - de*y(1) - np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/je;
    dydt(2) = (ks*y(4) + ds*(y(3) - y(2)) + np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/jd;
    dydt(3) = (-ks*y(4) - ds*(y(3) - y(2))- Tv)/jv ;
    dydt(4) = y(3) - y(2);
    dydt(5) = R - (-y(2) + y(1)); 

end

Here is the MATCONT GUI for Codim-1 Analysis:

I am plotting codim-1 bifurcation from equilibrium for y2 vs muk and y2 vs R.

In both cases I do not get any bifurcation point whatsoever.

[Umar]

Hi @Junaid Ali,

You mentioned, “The new issue I am stuck in with my continuation analysis is not finding any bifurcation point for any of my system parameter sweep. I know for sure I should see sub-critical hopf bifurcation in my system but there isnt any. I dont know what to do now. I am plotting codim-1 bifurcation from equilibrium for y2 vs muk and y2 vs R. In both cases I do not get any bifurcation point whatsoever.”

Please see my response to your comments below.

Here is modified version of your code:

% Define system parameters
params.je = 1.57; 
params.jd = 0.0066; 
params.jv = 4;
params.ks = 56665.5;
params.de = 0.01;
params.ds = 0.01;
params.Te = 100;
params.Tv = 80;
params.mu0 = 0.11;
params.muk = 0.056; % Initial value for muk
params.np = 14;
params.Rm = 0.0506;
params.Fn = 15000;
params.Kg = 1000;
params.R = 50; % Initial value for R

% Initial conditions
x0 = [0.01; 0.01; 0.01; 0.01; 0.01]; % Small perturbation
ts = 0:0.001:1;

% Solve ODE
[t,z] = ode45(@(t,z)clutch3dof(t,z,params), ts, x0);
figure;
plot(ts, z);
title('System Response Over Time');
xlabel('Time (s)');
ylabel('State Variables');
grid on;

% Finding equilibrium
options = optimoptions('fsolve', 'Display', 'off', 'TolFun', 1e-
6, 'TolX', 1e-6, 'MaxFunctionEvaluations', 20000, 
'MaxIterations', 20000);
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,   params), x0, options);

if exitflag > 0
  disp('Equilibrium found:');
  disp(x_eq);
else
  warning('Failed to find equilibrium. Check initial conditions 
  or parameter values.');
end

% Bifurcation analysis
params.muk_range = linspace(0.01, 0.1, 100); % Range for muk
params.R_range = linspace(40, 60, 100); % Range for R
bifurcation_points_muk = []; % Initialize bifurcation points for   muk
bifurcation_points_R = []; % Initialize bifurcation points for R

% Loop through muk range
for muk = params.muk_range
  params.muk = muk;
  [x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,     params), x0, options);
  if exitflag > 0
      bifurcation_points_muk = [bifurcation_points_muk; muk,         x_eq(2)];
  end
end

% Loop through R range
for R = params.R_range
  params.R = R;
  [x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,     params), x0, options);
  if exitflag > 0
      bifurcation_points_R = [bifurcation_points_R; R, x_eq(2)];
  end
end

% Plot bifurcation points for muk
figure;
plot(bifurcation_points_muk(:,1), bifurcation_points_muk(:,2), 
'o-');
title('Bifurcation Diagram: y2 vs muk');
xlabel('muk');
ylabel('y2 (Equilibrium State)');
grid on;

% Plot bifurcation points for R
figure;
plot(bifurcation_points_R(:,1), bifurcation_points_R(:,2), 'o-');
title('Bifurcation Diagram: y2 vs R');
xlabel('R');
ylabel('y2 (Equilibrium State)');
grid on;

% Equilibrium function
function eq = equilibrium_function(x, params)
  dxdt = clutch3dof(0, x, params);
  eq = dxdt;
end

% System dynamics
function dydt = clutch3dof(t, y, params)
  je = params.je; 
  jd = params.jd; 
  jv = params.jv;
  ks = params.ks;
  de = params.de;
  ds = params.ds;
  Te = params.Te;
  Tv = params.Tv;
  mu0 = params.mu0;
  muk = params.muk;
  np = params.np;
  Rm = params.Rm;
  Fn = params.Fn;
  R = params.R;
  Kg = params.Kg; 

    dydt = zeros(5,1);
    dydt(1) = (Te - de*y(1) - np*Rm*(Fn + Kg*y(5))*(mu0 + muk*(y(1)
    - y(2)))*tanh(2*(y(1) - y(2)))) / je;
    dydt(2) = (ks*y(4) + ds*(y(3) - y(2)) + np*Rm*(Fn + Kg*y(5))*
    (mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2)))) / jd;
    dydt(3) = (-ks*y(4) - ds*(y(3) - y(2)) - Tv) / jv;
    dydt(4) = y(3) - y(2);
    dydt(5) = R - (-y(2) + y(1)); 
  end

Please see attached.

Please read my summary below regarding comparison of your code and modified code

System Parameters and Initial Conditions

• Both versions of the code define system parameters (e.g., je, jd, etc.) similarly, ensuring consistency in the model being analyzed.

• The modified code utilizes a small perturbation for initial conditions (x0 = [0.01; 0.01; 0.01; 0.01; 0.01]), which can facilitate convergence towards an equilibrium point more effectively than starting with zeros, as seen in yours original code.

ODE Solver

• Both codes employ ode45 for solving ordinary differential equations (ODEs). However, the modified code introduces a time vector (ts) that allows for a more extensive temporal analysis of system dynamics.

• The plotting function provides a visual representation of state variables over time, enhancing interpretability.

Equilibrium Finding

• The modified code includes enhanced options for fsolve, such as increased tolerances and maximum evaluations, which likely improve the robustness of finding equilibrium points.

• The output indicates successful identification of multiple equilibrium points, which is critical for further bifurcation analysis.

Bifurcation Analysis

The modified code systematically sweeps through parameter ranges (muk and R) to identify bifurcation points, storing successful equilibria in dedicated arrays (bifurcation_points_muk and bifurcation_points_R). This structured approach contrasts with the original code, which lacked this systematic exploration and may explain yours difficulties in identifying bifurcation points.

Plotting Bifurcation Diagrams

The inclusion of bifurcation diagrams for both parameters allows for visual insights into how changes in muk and R affect the equilibrium states. This is crucial for understanding system behavior near critical transitions.

Output Results Analysis

The output indicates that equilibrium points have been successfully identified when executed the code.

Equilibrium found: 1.0e+03 * 2.0000 1.9500 1.9500 -0.0000 -0.0150

This result confirms that the modified code has effectively resolved the issue of finding equilibria, which was a significant challenge in yours previous attempts.

Now, you noted that achieving equilibrium required closing the feedback loop in your model mentioned in your earlier comments. The modified code's structure likely reflects this necessity by ensuring that all relevant dynamics are accounted for within the function definitions. Also, with successful identification of equilibria, you can now analyze stability and potential bifurcations effectively. The presence of sub-critical Hopf bifurcations can be examined through further stability analysis near these equilibrium points.

However, I will suggest that it may be beneficial to conduct sensitivity analyses on key parameters like muk and R to understand their influence on system behavior comprehensively.

Also, if issues persist in identifying expected bifurcations, you can still revisit utilizing MATCONT to provide further insights into codimension-1 bifurcations.

In nutshell, my modifications made to yours original code shows in the attached plots improving its ability to identify equilibrium points and conduct bifurcation analysis, aligning closely with yours requirements and previous observations about feedback loops in dynamical systems.

Hope this answers all your questions.

[Junaid Ali]

Thank you @Umar

I get it what you did. I think my next step shoud be defining a jacobian function and calculate eigen values at each equilibrium point using that jacobian. And see if I can identify any Hopf point.

% Define system parameters
params.je = 1.57; 
params.jd = 0.0066; 
params.jv = 4;
params.ks = 56665.5;
params.de = 0.01;
params.ds = 0.01;
params.Te = 100;
params.Tv = 80;
params.mu0 = 0.11;
params.muk = 0.056; % Initial value for muk
params.np = 14;
params.Rm = 0.0506;
params.Fn = 15000;
params.Kg = 1000;
params.R = 50; % Initial value for R
% Initial conditions
x0 = [0.01; 0.01; 0.01; 0.01; 0.01]; % Small perturbation
ts = 0:0.001:1;
% Solve ODE
[t,z] = ode45(@(t,z)clutch3dof(t,z,params), ts, x0);
figure;
plot(ts, z);
title('System Response Over Time');
xlabel('Time (s)');
ylabel('State Variables');
grid on;

% Finding equilibrium
options = optimoptions('fsolve', 'Display', 'off', 'TolFun', 1e-6, 'TolX', 1e-6, 'MaxFunctionEvaluations', 20000, 'MaxIterations', 20000);
[x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,   params), x0, options);
if exitflag > 0
  disp('Equilibrium found:');
  disp(x_eq);
else
  warning('Failed to find equilibrium. Check initial conditions or parameter values.');
end

Equilibrium found:

   1.0e+03 *

    2.0000
    1.9500
    1.9500
   -0.0000
   -0.0150

% Bifurcation analysis
params.muk_range = linspace(-0.1, 0.1, 10); % Range for muk
params.R_range = linspace(0, 150, 10); % Range for R
bifurcation_points_muk = []; % Initialize bifurcation points for   muk
bifurcation_points_R = []; % Initialize bifurcation points for R
% Loop through muk range
for muk = params.muk_range
  params.muk = muk;
  [x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,     params), x0, options);
  if exitflag > 0
      bifurcation_points_muk = [bifurcation_points_muk; muk,         x_eq(2)];
  end
end
% Loop through R range
for R = params.R_range
  params.R = R;
  [x_eq, fval, exitflag] = fsolve(@(x) equilibrium_function(x,     params), x0, options);
  if exitflag > 0
      bifurcation_points_R = [bifurcation_points_R; R, x_eq(2)];
  end
end
% Plot bifurcation points for muk
figure;
plot(bifurcation_points_muk(:,1), bifurcation_points_muk(:,2), 'o-');
title('Bifurcation Diagram: y2 vs muk');
xlabel('muk');
ylabel('y2 (Equilibrium State)');
grid on;

% Plot bifurcation points for R
figure;
plot(bifurcation_points_R(:,1), bifurcation_points_R(:,2), 'o-');
title('Bifurcation Diagram: y2 vs R');
xlabel('R');
ylabel('y2 (Equilibrium State)');
grid on;

figure
% Initialize previous eigenvalue storage
E_vals_prev = [];

% Loop through values of muk
for muk = params.muk_range
    params.muk = muk;
    
    % Calculate eigenvalues for the Jacobian matrix
    E_vals = eig(Jacobian2DOF(x_eq, params));
    
    % Normalize the eigenvalues to lie within the unit circle
    E_vals_normalized = E_vals ./ abs(E_vals);
    
    % Plot the current normalized eigenvalues
    hold on
    plot(real(E_vals_normalized), imag(E_vals_normalized), '*b', 'MarkerSize', 8);
    
    % If this is not the first iteration, plot arrows to show movement
    if ~isempty(E_vals_prev)
        for i = 1:length(E_vals)
            % Plot arrows indicating movement from previous to current eigenvalue
            quiver(real(E_vals_prev(i)), imag(E_vals_prev(i)), ...
                   real(E_vals_normalized(i)) - real(E_vals_prev(i)), ...
                   imag(E_vals_normalized(i)) - imag(E_vals_prev(i)), ...
                   0, 'MaxHeadSize', 0.5, 'Color', 'k'); % Arrows in red
        end
    end
    
    % Store the current eigenvalues for the next iteration
    E_vals_prev = E_vals_normalized;

    % Plot the unit circle for reference
    theta = linspace(0, 2*pi, 100);
    plot(cos(theta), sin(theta), '--k');  % Unit circle
    
    % Set up axis labels, grid, and aspect ratio
    xlabel('Real Part');
    ylabel('Imaginary Part');
    axis equal;
    grid on;
end

% Equilibrium function
function eq = equilibrium_function(x, params)

    dxdt = clutch3dof(0, x, params);
    eq = dxdt;

end

function jacob = Jacobian2DOF(y_eq,params)
    
    je = params.je; 
    jd = params.jd; 
    jv = params.jv;
    ks = params.ks;
    de = params.de;
    ds = params.ds;
    mu0 = params.mu0;
    muk = params.muk; 
    np = params.np;
    Rm = params.Rm;
    Fn = params.Fn;
    Kg = params.Kg; 


    y1 = y_eq(1);
    y2 = y_eq(2);
    y3 = y_eq(3);
    y4 = y_eq(4);
    y5 = y_eq(5);


% First row of Jacobian (A1)
A1_1 = (-muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) + ...
        (muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (2 - 2 * tanh(2 * y1 - 2 * y2)^2) - de) / je;
A1_2 = (muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) + ...
        (muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (-2 + 2 * tanh(2 * y1 - 2 * y2)^2)) / je;
A1_3 = 0;
A1_4 = 0;
A1_5 = ((muk * (y2 - y1) - mu0) * np * Rm * Kg * tanh(2 * y1 - 2 * y2)) / je;

% Second row of Jacobian (A2)
A2_1 = (muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) - ...
        (muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (2 - 2 * tanh(2 * y1 - 2 * y2)^2)) / jd;
A2_2 = (-muk * np * Rm * (Kg * y5 + Fn) * tanh(2 * y1 - 2 * y2) - ...
        (muk * (y2 - y1) - mu0) * np * Rm * (Kg * y5 + Fn) * (-2 + 2 * tanh(2 * y1 - 2 * y2)^2) - ds) / jd;
A2_3 = ds / jd;
A2_4 = ks / jd;
A2_5 = -(muk * (y2 - y1) - mu0) * np * Rm * Kg * tanh(2 * y1 - 2 * y2) / jd;

% Third row of Jacobian (A3)
A3_1 = 0;
A3_2 = ds / jv;
A3_3 = -ds / jv;
A3_4 = -ks / jv;
A3_5 = 0;

% Fourth row of Jacobian (A4)
A4_1 = 0;
A4_2 = -1;
A4_3 = 1;
A4_4 = 0;
A4_5 = 0;

% Fifth row of Jacobian (A5)
A5_1 = -1;
A5_2 = 1;
A5_3 = 0;
A5_4 = 0;
A5_5 = 0;

jacob = [
    A1_1, A1_2, A1_3, A1_4, A1_5;
    A2_1, A2_2, A2_3, A2_4, A2_5;
    A3_1, A3_2, A3_3, A3_4, A3_5;
    A4_1, A4_2, A4_3, A4_4, A4_5;
    A5_1, A5_2, A5_3, A5_4, A5_5];


end

function dydt = clutch3dof(t,y,params)
    je = params.je; 
    jd = params.jd; 
    jv = params.jv;
    
    ks = params.ks;
    
    de = params.de;
    ds = params.ds;
    
    Te = params.Te;
    Tv = params.Tv;
    
    mu0 = params.mu0;
    muk = params.muk;
    
    np = params.np;
    Rm = params.Rm;
    Fn = params.Fn;
    R = params.R;
    Kg = params.Kg; 

    dydt = zeros(5,1);

    dydt(1) = (Te - de*y(1) - np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/je;
    dydt(2) = (ks*y(4) + ds*(y(3) - y(2)) + np*Rm*(Fn+Kg*y(5))*(mu0 + muk*(y(1) - y(2)))*tanh(2*(y(1) - y(2))))/jd;
    dydt(3) = (-ks*y(4) - ds*(y(3) - y(2))- Tv)/jv ;
    dydt(4) = y(3) - y(2);
    dydt(5) = R - (-y(2) + y(1)); 

end

[Umar]

Hi @Junaid Ali,

Your methodology appears sound, and with careful execution of these steps, you should be able to derive meaningful insights about the dynamics of your system. Pay close attention to the stability conditions indicated by the eigenvalues, especially when identifying potential Hopf points, as they can indicate transitions between different dynamic behaviors in your system. Good luck with your analysis!