How can I produce plots that mimic the one found in this paper?

1 view (last 30 days)
I have put a lot of hours into this MATLAB code, but do not have it quite right. I am trying to recreate plots of a ground reaction force from a mass-spring system. **The paper can be found here: http://hal.archives-ouvertes.fr/docs/00/42/54/36/PDF/Ly_et_al-10.pdf
The MATLAB code utilizes the ode45 function. I can get the plots to look very similar, but not perfectly matched (which is what I need). I need to use the same parameters that they use, so I need to find the mistake in my code that allows me to produce the same plots.
For simplicity, I am only trying to plot the "65 Asker C" line in Figure 5. Once I get that perfect, I will easily be able to plot all the other lines. The parameters the authors use can be found inside the paper (and found in my code).
My initial acceleration values were calculated given the initial position and velocities of the masses.
Basically, I am asking someone to look at my code and find the problem. I know this may not be easy, but I appreciate anyone that tries (or even considers it!). Thanks :)
My code is in two documents: 1 containing the function for the ode45, and the other is for the actual calling of the ode45 file.
**FILE ONE:
function dxdt = Ly2( t, x ) ;
% Ly, et al. expansion of Zadpoor's model
% masses & gravity
m1 = 6.15 ; % kg
m2 = 6.00 ;
m3 = 12.58 ;
m4 = 50.34 ;
m5 = 0.3 ;
g = -9.81 ; % m/s^2
% spring and damping constants @ 65 Asker C
k1 = 6000 ; % N/m
k2 = 6000 ;
k3 = 10000 ;
k4 = 10000 ;
k5 = 18000 ;
k6 = 403000 ;
ks = 880000 ;
c1 = 300 ; % kg/s
c2 = 650 ;
c4 = 1900 ;
c6 = 2170 ;
b6 = 2.4346 ;
Ro = 37.5 ; % This value may be 40, but either way, it doesn't make a big difference.
% Reassigning the values to make it simpler to write.
x1 = x(1) ; x2 = x(2) ; x3 = x(3) ; x4 = x(4) ; x5 = x(5) ;
x1d = x(6) ; x2d = x(7) ; x3d = x(8); x4d = x(9) ; x5d = x(10) ;
Fs = k6*(x1-x5)*((abs(x1-x5)/Ro)^(b6-1)) ;
% Fs is techincally a little different in the paper, but the portion that I
% left out did not get left out in the function below (c6*(x5d-x1d))
x1dd = (m1*g+k1*(x3-x1)+k2*(x2-x1)+Fs+c1*(x3d-x1d)+c2*(x2d-x1d)+c6*(x5d-x1d))/m1;
x2dd = (m2*g+k2*(x1-x2)+k3*(x3-x2)+c2*(x1d-x2d))/m2;
x3dd = (m3*g+k1*(x1-x3)+k3*(x2-x3)+(k4+k5)*(x4-x3)+c1*(x1d-x3d)+c4*(x4d-x3d))/m3;
x4dd = (m4*g+(k4+k5)*(x3-x4)+c4*(x3d-x4d))/m4;
x5dd = (m5*g+Fs-ks*x5+c6*(x1d-x5d))/m5;
% The set of 1st-order ODE's that represents the set of 2nd-order ODE's
% from the article. The last 6 parts of the set are for my convenience
% when analyzing how the masses move. Also, the ground reaction force is
% the -k5*x5 value.
dxdt = [ x1d; x2d; x3d; x4d; x5d; x1dd; x2dd; x3dd; x4dd; x5dd; -ks*x5; x1; x2; x3; x4; x5 ] ;
end
*FILE 2:
% Ly Demo
clc
tic
% using the ODE file.
[t, dxdt ] = ode45(@Ly2, [0 .25], [-0.96 -0.96 -2.0 -2.0 -0.96 -60.5417...
-9.81 14.99 -9.81 -9.81 0 0 0 0 0 0]) ;
figure(1)
plot(t, dxdt(:,11), 'r', 'LineWidth', 2)
set(gca, 'YTick', [0 200 400 600 800 1000 1200 1400 1600 1800 2000])
set(gca, 'XTick', [0 .05 .1 .15 .2 .25])
grid on
toc

Answers (1)

Steve Miller
Steve Miller on 6 Dec 2022
You can draw the schematic exactly as shown in the paper and get your plots. You could easily do it in MATLAB Online.
--Steve

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!