ODE friction-vibration system: unable to meet integration tolerances

I was given this dynamic model as an assignment in vibration class

here is how I write it as EOMs

Ft is the force exerts from vertica spring and damper, all the constants were given as range of values. I picked the values that in my understanding, could minimize vibration. Here is my code using massMatrix form

syms Xf omgf t Yr lambda x(t) y(t) mp cp cq kp kq cq kq mu kc cc Fy
x = x(t)
y = y(t)
xf = Xf*sin(omgf*t)
dxf = diff(xf)
yb = Yr*sin(2*pi*(diff(x))*t/lambda)
dyb = diff(yb)
Ft = mp*diff(y,2) + Fy

eqa1 = mp*diff(x,2) + (cp+cq)*diff(x) + (kp+kq)*x == cq*dxf + kq*xf - sign(diff(x))*mu*Ft;
eqa2 = mp*diff(y,2) + cc*diff(y) + kc*y == cc*dyb + kc*yb;
eqns = [eqa1 eqa2]
vars = [x; y]
[eqns, vars] = reduceDifferentialOrder(eqns, vars)
[DAEs, DAEvars] = reduceRedundancies(eqns, vars)

[m,f] = massMatrixForm(DAEs,DAEvars)
m = odeFunction(m, DAEvars, lambda, mp, Yr, cc, mu)
f = odeFunction(f, DAEvars, mp, cp, cq, cc, kp, kq, kc, Xf, omgf, Yr, lambda, mu, Fy)
mp = 0.0075;
cp = 1.5e+3;
cq = 1.5e+3;
cc = 4e+3;
kp = 500e+3;
kq = 500e+3;
kc = 2500e+3;
Xf = 10e-3;
omgf = 3;
Yr = 0.5e-6; 
lambda = 60e-6;
mu = 0.01;
Fy = 1500;
F = @(t, X) f(t, X, mp, cp, cq, cc, kp, kq, kc, Xf, omgf, Yr, lambda, mu, Fy);
M = @(t, X) m(t, X, lambda, mp, Yr, cc, mu);
x0est = [0+eps; 0+eps; 0+eps; 0+eps];
x0pest = [0+eps; 0+eps; 0+eps; 0+eps];
opt = odeset('Mass', M, 'InitialSlope',x0pest,'RelTol',1e-6,'AbsTol',1e-6);
implicitDAE = @(t,x,xp) M(t,x)*xp - F(t,x);
[x0, x0p] = decic(implicitDAE, 0, x0est, [], x0pest, [], opt)
opt = odeset(opt, 'InitialSlope', x0p, "RelTol", 1e-8, 'AbsTol', 1e-8);
[tSol,ySol] = ode15s(F, linspace(0,2,20000), x0, opt)

The problem is I keep getting (Warning: Failure at t=1.107525e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.469447e-18) at time t.) so I tried to plot all the solution at the time right before error occurs. I can see the spikes in many graphs but I can't figure out what variable is the cause of this. the most visible spike is in ydot graph. I have the graph code here just to be use as an reference.

plot(tSol,ySol(:,1));
title('X axis');
xlabel('Time (s)');
ylabel('Displacement (m)');
plot(tSol,ySol(:,2));
title('Y axis');
xlabel('Time (s)');
ylabel('Displacement (m)');
plot(tSol,ySol(:,3));
title('Xdot');
plot(tSol,ySol(:,4));
title('Ydot');
ySize = size(ySol);
for i = 1:ySize
    ydot(i,:) = (M(tSol(i),ySol(i,:)')\F(tSol(i),ySol(i,:)'))';
end
plot(tSol,ydot(:,3));
title('xdotdot');
plot(tSol,ydot(:,4));
title('ydotdot');

I wonder that is it because the equation is "too stiff?", if so is there a way to solve it? or can you give me an advise on how to analyze the problem here

[Star Strider]

The ode15s function is encountering a non-finite result (either ±Inf or NaN). The interesting thing is that increasing the initial conditions vector from eps to a slightly larger number allows the integration to proceed beyond that with the eps values. It eventually encounters the singularity and stops, likely because the integrated evalues decrease to the extent that the non-finite result occurs.

syms Xf omgf t Yr lambda x(t) y(t) mp cp cq kp kq cq kq mu kc cc Fy
x = x(t)

x = 

y = y(t)

y = 

xf = Xf*sin(omgf*t)

xf = 

dxf = diff(xf)

dxf = 

yb = Yr*sin(2*pi*(diff(x))*t/lambda)

yb = 

dyb = diff(yb)

dyb = 

Ft = mp*diff(y,2) + Fy

Ft = 

eqa1 = mp*diff(x,2) + (cp+cq)*diff(x) + (kp+kq)*x == cq*dxf + kq*xf - sign(diff(x))*mu*Ft;
eqa2 = mp*diff(y,2) + cc*diff(y) + kc*y == cc*dyb + kc*yb;
eqns = [eqa1 eqa2]

eqns = 

vars = [x; y]

vars = 

[eqns, vars] = reduceDifferentialOrder(eqns, vars)

eqns = 

vars = 

[DAEs, DAEvars] = reduceRedundancies(eqns, vars)

DAEs = 

DAEvars = 

[m,f] = massMatrixForm(DAEs,DAEvars)

m = 

f = 

m = odeFunction(m, DAEvars, lambda, mp, Yr, cc, mu);
f = odeFunction(f, DAEvars, mp, cp, cq, cc, kp, kq, kc, Xf, omgf, Yr, lambda, mu, Fy);
mp = 0.0075;
cp = 1.5e+3;
cq = 1.5e+3;
cc = 4e+3;
kp = 500e+3;
kq = 500e+3;
kc = 2500e+3;
Xf = 10e-3;
omgf = 3;
Yr = 0.5e-6; 
lambda = 60e-6;
mu = 0.01;
Fy = 1500;
F = @(t, X) f(t, X, mp, cp, cq, cc, kp, kq, kc, Xf, omgf, Yr, lambda, mu, Fy);
M = @(t, X) m(t, X, lambda, mp, Yr, cc, mu);
% x0est = [0+eps; 0+eps; 0+eps; 0+eps];
% x0pest = [0+eps; 0+eps; 0+eps; 0+eps];
x0est = zeros(4,1)+1E-5;
x0pest = zeros(4,1)+1E-5;
opt = odeset('Mass', M, 'InitialSlope',x0pest,'RelTol',1e-6,'AbsTol',1e-6);
implicitDAE = @(t,x,xp) M(t,x)*xp - F(t,x);
[x0, x0p] = decic(implicitDAE, 0, x0est, [], x0pest, [], opt)

x0 = 4×1

1.0e-05 * 1.0000 1.0000 1.0000 1.0000

x0p = 4×1

1.0e+03 * 0.0000 0.0000 2.6961 -3.3384

opt = odeset(opt, 'InitialSlope', x0p, "RelTol", 1e-8, 'AbsTol', 1e-8);
[tSol,ySol] = ode15s(F, linspace(0,2,20000), x0, opt)

Warning: Failure at t=5.541039e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.110223e-16) at time t.

tSol = 555×1

1.0e+00 * 0 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009

ySol = 555×4

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0069 -0.0055 0.0000 0.0000 0.0072 -0.0051 0.0000 0.0000 0.0075 -0.0048 0.0000 0.0000 0.0077 -0.0044 0.0000 0.0000 0.0079 -0.0041 0.0000 0.0000 0.0082 -0.0039 0.0000 0.0000 0.0084 -0.0036 0.0000 0.0000 0.0086 -0.0034 0.0000 0.0000 0.0088 -0.0032

plot(tSol,ySol(:,1));
title('X axis');
xlabel('Time (s)');
ylabel('Displacement (m)');

plot(tSol,ySol(:,2));
title('Y axis');
xlabel('Time (s)');
ylabel('Displacement (m)');

plot(tSol,ySol(:,3));
title('Xdot');

plot(tSol,ySol(:,4));
title('Ydot');

ySize = size(ySol);
for i = 1:ySize
    ydot(i,:) = (M(tSol(i),ySol(i,:)')\F(tSol(i),ySol(i,:)'))';
end
plot(tSol,ydot(:,3));
title('xdotdot');

plot(tSol,ydot(:,4));
title('ydotdot');

.

[David Goodmanson]

Hi Kasidit,

in the expression for y_R, is there a reason you have

Yr*sin(2*pi*(diff(x))*t/lambda)

as opposed to

Yr*sin(2*pi*x/lambda) ?

[Kasidit Suwannagarn]

This equation was given in the assignment, but it does make sense as the mass moves faster (diff(x) increase) the mass moves through more waves at a time resulting in higher frequency of yr function.

Answer 1

[Sam Chak]

Hi @Kasidit Suwannagarn

Edit: Replaced the previous flawed answer with the system analysis.

% constants
mp     = 0.0075;
cp     = 1.5e+3;
cq     = 1.5e+3;
cc     = 4e+3;
kp     = 500e+3;
kq     = 500e+3;
kc     = 2500e+3;
Xf     = 10e-3;
omgf   = 3;
Yr     = 0.5e-6; 
lambda = 60e-6;
mu     = 0.01;
Fy     = 1500;
sympref('AbbreviateOutput', false);

syms x(t) y(t)
xF     = Xf*sin(omgf*t);
xFdot  = omgf*Xf*cos(omgf*t);
yR     = Yr*sin(2*pi*(diff(x)/lambda)*t);
yRdot  = (2*pi*(diff(x)/lambda) + 2*pi*t*(diff(x,2)/lambda))*Yr*cos(2*pi*(diff(x)/lambda)*t);
Ft     = kc*(yR - y(t)) + cc*(yRdot - diff(y)) + Fy;
F      = - sign(diff(x))*mu*Ft;
eqn1   = mp*diff(x,2) == cq*xFdot + kq*xF + F - (kp + kq)*x(t) - (cp + cq)*diff(x);
eqn2   = mp*diff(y,2) == cc*yRdot + kc*yR - kc*y(t) - cc*diff(y);
eqns   = [eqn1 eqn2];
[V, S] = odeToVectorField(eqns)

V =

Y[2]

-(3125*(53126622932283508654080000*Y[1] - 26563311466141753125*sin((100000*pi*t*Y[4])/3) + 85002596691653613846528*Y[2] - 1416709944860893500000*pi*cos((100000*pi*t*Y[4])/3)*Y[4] + 188894659314785800000000000000*t*pi*cos((100000*pi*t*Y[4])/3)*Y[3] + 566683977944357400000000000*t*pi*cos((100000*pi*t*Y[4])/3)*Y[4] + 2833419889721787000000000*t*pi*sign(Y[4])*cos((100000*pi*t*Y[4])/3) - 8500259669165361000000000*t*pi*cos(3*t)*cos((100000*pi*t*Y[4])/3) - 944473296573929000000000000*t*pi*cos((100000*pi*t*Y[4])/3)*sin(3*t)))/(24*(1844674407370955078125*t*pi*sign(Y[4])*cos((100000*pi*t*Y[4])/3) + 20752587082923245568))

Y[4]

-(125*(10625324586456701730816*sign(Y[4]) - 31875973759370105192448*cos(3*t) - 3541774862152233910272000*sin(3*t) - 17708874310761169551360000*sign(Y[4])*Y[1] - 28334198897217871282176*sign(Y[4])*Y[2] + 708354972430446782054400000*Y[3] + 2125064917291340346163200*Y[4] + 8854437155380584375*sign(Y[4])*sin((100000*pi*t*Y[4])/3) + 472236648286964500000*pi*sign(Y[4])*cos((100000*pi*t*Y[4])/3)*Y[4]))/(32*(1844674407370955078125*t*pi*sign(Y[4])*cos((100000*pi*t*Y[4])/3) + 20752587082923245568))

S =

y

Dy

x

Dx

If you look at the vector field, you will see the division by . This implies singularity occurs when . If the simulation starts from t = 0, then you get non-finite result.

[Kasidit Suwannagarn]

Your yRdot is incorrect, y(2) = xdot which is a function of time so you also have to differentiate it.

this will result in xdotdot in the equation

[Kasidit Suwannagarn]

and also isn't it more efficient to use other ode? since this equation is stiff from what I understand.

[Sam Chak]

Hi @Kasidit Suwannagarn

Perhaps the system is not as stiff as you think. I actually followed the advice given here.

https://www.mathworks.com/help/matlab/ref/ode15s.html

Edit: Thanks for showing the mistake on yRdot. Now it becomes a little complicated due to on the right-hand side of , where yRdot is. In fact, contains F, which is a function of that depends on yRdot.

[Kasidit Suwannagarn]

I had used ode45 and it ran terribly slow. I left it running all night just to found out that it had only compute 4 time steps.

[Sam Chak]

I can imagine that. My longest record is 16 hours. I guess the system is very stiff due to signum function and the very small wavelength λ at the scale of . More importantly, I think you have to check if singularity really occurs, and what is causing the non-finite result.

I have edited my answer to investigate why singularity happens.