Time-dependent system ODE with initial half-sine pulse shows incorrect results

Hello,

I have damped spring-mass system with initial half-sine pulse (acceleration) h loaded from a base.

Motion of equation should be the following:

m*d2(y)/d2t + c*d(y)/dt + k*y = k*(h)dtdt + c*(h)dt

k*(h)dtdt = stiffness*double integrated acceleration h = stiffness*displacement [N/m*m]

c*(h)dt = damping_coeff*integrated acceleration h = damping_coeff*velocity [Ns/m*m/s]

*Previous part has been edited for clarification.

Here are my concerns:

1. I am not sure how to introduce here the acceleration pulse in the equations,

2. for integration of k*(h)dtdt + c*d(h)dt I am using cumtrapz function. Is it OK?

Best regards

Mathew

t1=0.5;
a=0; h0=1;
h = @(t) a*(t>t1)+ h0*sin(pi.*t./t1).*(t<=t1);

ODExvf = @(t,y,m,D,K,h) [y(2);
                  -D/m*y(2)-K/m*y(1)+ cumtrapz(cumtrapz(t,h(t)))*K/m + cumtrapz(t,h(t))*D/m + h(t)];
m=1;
D=10;
K=100;
tspan=0.01:0.05:5;
F0=[0 0];

[t,F]=ode45(@(t,y) ODExvf(t,y,m,D,K,h),tspan,F0);

Error using matlab.internal.math.getdimarg
Dimension argument must be a positive integer scalar within indexing range.

Error in cumtrapz>getDimArg (line 91)
dim = matlab.internal.math.getdimarg(dim);

Error in cumtrapz (line 49)
dim = min(ndims(y)+1, getDimArg(dim));

Error in solution (line 7)
-D/m*y(2)-K/m*y(1)+ cumtrapz(cumtrapz(t,h(t)))*K/m + cumtrapz(t,h(t))*D/m + h(t)];

Error in solution (line 14)
[t,F]=ode45(@(t,y) ODExvf(t,y,m,D,K,h),tspan,F0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

plot(t,F);  

Answer 1

[William Rose]

@Mathew Smith,

I would not use cumtrapz() or (worse) nested calls to cumtrapz(cumtrapz()).

I don;t understand the right hand side of the equation

m*d(y)/dt + c*d(y)/dt + k*y = k*(h)dtdt + c*(h)dt

It seems you use dt*dt after k*h to get the units of force, and likewise you use dt after c*h to get to units of force. It is always nice when the units are correct, but I don't know physically what this means. Why is the acceleration multiplied by the spring constant k and by the dashpot constant c? Maybe your system supplies as much force as necessary to acheive a specific acceleration. If so, then during that period, you don;t need to solve a complicated differential equaiton - you can just differentiate once and twice to get velocity and position. When the pulse disconnects, then you solve the unforced system, with a zero on the right hand side.

Also, I assume you meant to write

m*d^2(y)/dt^2 on the left hand side.

[Mathew Smith]

Hello William, I have edited original equation for better clarity: on the right side is displacement of the base and velocity. There is no acceleration of the base, so I am including it by integration of displ. and velocity.

[William Rose]

@Mathew Smith,

THank you for clarifying the equation.

Now that it is clear, a problem is also clear: each term on the left side has units of force. The terms on the right side have units of displacement and velocity. Or maybe displacement and displacement*. Either way, all terms on both sides must have the same units.

*You wrote, in your clarification, "c*(h)dt = once integrated displacement = velocity", but I suspect you meant to write "c*(h)dt = once integrated velocity = displacement".

I still am not sure I understand what you are trying to simulate. I think you have a mass attached to a base by a spring and a dashpot in parallel. The base undergoes a half-sine-wave displacement pulse. This brief movement of the base will set the mass into motion (if it was previously stationary), and it will disrupt the motion of the mass (if it was already moving).

However, you said in your comment "There is no acceleration of the base...." This statement is inconsistent with the problem as I described it. If the base goes through a half-sine-wave of motion then it DOES accelerate. Therefore my understanding of the system is wrong, or your statement that there is no acceleration of the base is wrong.

[Mathew Smith]

Hello William,

the terms on the right side have units of force: k*x [N/m*m] and c*v [Ns/m*m/s]. I have clarified it above. You understood correctly, that there is base excitation. However, it is not excited by half-sine-wave displacement, but half-sine-wave acceleration. I can do it with displacement pulse (although my solution is more complicated than your suggestion). I am not able to include there half-sine-wave acceleration base excitation.

I am sorry for not clear description of the problem and I thank you for your effort.

Answer 2

[William Rose]

@Mathew Smith,

If I am correct that the base moves according to a prescribed function of time, then the differential equation can be written

where y is the mass position and y_b is the base position. Your code will include the half-pulses yb=A*sin(wt) and vb=dyb/dt=A*w*cos(wt). To make these functions zero after an initial half-pulse, do this:

yb=A*sin(w*t).*(t<pi/w);

vb=A*w*cos(w*t).*(t<pi/w);

function dydt = myode(t,y)
m=1; c=1; k=1;                      % physical system properties
A=1; w=1;                           % amplitude, frequency of base displacement
yb=A*sin(w*t).*(t<pi/w);            % displacement pulse
vb=A*w*cos(w*t).*(t<pi/w);          % velocity pulse
dydt=[y(2);[figure out this part]];
end

That gets you started.

[Mathew Smith]

Hello William, thank you for you reply. This solution uses displacement pulse but I want to use acceleration half-sine-wave pulse.

[William Rose]

@Mathew Smith,

You are welcome.

A half-sine-wave acceleration pulse is applied to the base. Suppose the peak magntude of the acceleration is G and the frequency is w.

function dydt = myode(t,y)
m=1; c=1; k=1;                  % physical system properties
G=1;                            % base acceleration magnitude (m/s^2)
w=1;                            % base acceleration frequency (1/s)
yb=(G/w^2)*sin(w*t).*(t<pi/w);  % displacement pulse (m)
vb=(G/w)*cos(w*t).*(t<pi/w);    % velocity pulse (m/s)
dydt=[y(2);[figure out this part]];
end

[Mathew Smith]

Hello William,

thank you for the answer. I have compared yb and vb from your solution and yb and vb integrated by "cumtrapz and they are not the same. Do you know why is it so?

clear all; close all

t1=0.5;
T=t1*2; w=2*pi*(1/T); G=1;

ab= @(t) G*sin(w*t).*(t<pi/w);

t=0:0.01:3;
yb=@(t) (G/w^2)*sin(w*t).*(t<pi/w);
vb=@(t) (G/w)*cos(w*t).*(t<pi/w);

subplot(2,1,1)
plot(t,ab(t),"linewidth",1.5)
hold on
plot(t,vb(t), "linewidth",1.5)
hold on
plot(t,yb(t), "linewidth",1.5)
grid on
h1 = legend ("ab","vb","yb");

ab_c=ab(t);
vb_c = cumtrapz(t,ab_c);
yb_c = cumtrapz(t,vb_c);

subplot(2,1,2)
plot(t,ab_c, "linewidth",1.5)
hold on
plot(t,vb_c,"linewidth",1.5)
hold on
plot(t,yb_c,"linewidth",1.5)
h2 = legend ("ab cumtrapz","vb cumtrapz","yb cumtrapz");

[Torsten]

When you integrate "ab" once and twice in the limits from 0 to t, you forgot to include the constants of integration at t=0.

If you work it out, you'll arrive at:

clear all; close all
t1=0.5;
T=t1*2; w=2*pi*(1/T); G=1;
ab= @(t) G*sin(w*t).*(t<pi/w);
t=0:0.01:3;
vb=@(t) (-G/w*cos(w*t)+G/w).*(t<pi/w) + (-G/w*cos(w*pi/w)+G/w).*(t>=pi/w);
yb=@(t) (-G/w^2*sin(w*t)+G/w*t).*(t<pi/w) + (-G/w^2*sin(w*pi/w)+G/w*pi/w+(-G/w*cos(w*pi/w)+G/w)*(t-pi/w)).*(t>=pi/w);
subplot(2,1,1)
plot(t,ab(t),"linewidth",1.5)
hold on
plot(t,vb(t), "linewidth",1.5)
hold on
plot(t,yb(t), "linewidth",1.5)
grid on
h1 = legend ("ab","vb","yb");
ab_c=ab(t);
vb_c = cumtrapz(t,ab_c);
yb_c = cumtrapz(t,vb_c);
subplot(2,1,2)
plot(t,ab_c, "linewidth",1.5)
hold on
plot(t,vb_c,"linewidth",1.5)
hold on
plot(t,yb_c,"linewidth",1.5)
grid on
h2 = legend ("ab cumtrapz","vb cumtrapz","yb cumtrapz");

But I think it will be easier to include the additional second-order differential equation

d^2y_b/dt^2 = G*sin(w*t).*(t<pi/w), y_b(0) = y_b0, y_b'(0) = v_b0

in your set of equations.

[William Rose]

@Torsten, thank you.

[Mathew Smith]

Thank you very much @Torsten and @William Rose!

@Torsten, I tried to include your suggestion with additional differential equation, but I don't have much experience with ode45, therefore I would like to ask you for help. Below is my attempt for solution.

function dydt=odefcn(t,y,yb);

t1=0.5; T=t1*2; w=2*pi*(1/T);G=1;
dydt = zeros(4:1);
dydt(1)=yb(2);
dydt(2)=G*sin(w*t).*(t<pi/w);
dydt(3)=y(2);
dydt(4)=-c*y(2)/m - k*y(1)/m + c*yb(2)/m + k*yb(1)/m;
end

m=1;
c=1;
k=1;
tspan=0:0.01:5;
ic = [0 0 0 0];

[t, f] = ode45(@(t,y,yb) odefcn(t,y,yb) tspan, ic);

[Torsten]

What are your initial conditions for yb, i.e. yb(t=0) and yb'(t=0) = vb(t=0) ?

What are your initial conditions for y, i.e. y(t=0) and y'(t=0) = v(t=0) ?

What are the differential equations you are now trying to solve ?

[Mathew Smith]

Initial conditions for yb, yb', y and y' are 0.

Diff eqs. are:

yb" = G*sin(w*t).*(t<p/w);

y" = - c*y'/m - k*y/m + c*yb'/m + k*yb/m;

The rest is in the code above.

[Torsten]

m=1;
c=1;
k=1;
tspan=0:0.01:5;
ic = [0 0 0 0];
[t, Y] = ode45(@(t,y) odefcn(t,y,m,c,k), tspan, ic);
figure(1)
plot(t,[Y(:,1:2),-c*Y(:,2)/m - k*Y(:,1)/m + c*Y(:,4)/m + k*Y(:,3)/m])
grid on

figure(2)
t1=0.5; T=t1*2; w=2*pi*(1/T);G=1;
plot(t,[Y(:,3:4),G*sin(w*t).*(t<pi/w)])
grid on

function dYdt=odefcn(t,Y,m,c,k);
  y = Y(1:2);
  yb = Y(3:4);
  t1=0.5; T=t1*2; w=2*pi*(1/T);G=1;
  dYdt = zeros(4,1);
  dYdt(1)=y(2);
  dYdt(2)=-c*y(2)/m - k*y(1)/m + c*yb(2)/m + k*yb(1)/m;
  dYdt(3)=yb(2);
  dYdt(4)=G*sin(w*t).*(t<pi/w);
end

[William Rose]

@Mathew Smith,

[edit: added another attribution to @Torsten], and changed "xb" to "yb" to be consistent with the preceding discussion]

The function odefcn() won;t know what m, c, k are. They need to be defined inside odefcn.

Also, the line

[t, f] = ode45(@(t,y,yb) odefcn(t,y,yb) tspan, ic);

looks like trouble, because, as far as I know, you can only pass two variables, not three, to ode45(). When you do @(t,y,yb), you are passing three.

Getting back to the earlier discrepancy between my (incorrect) solution and the solution from cumtrapz: I made a msitake when I did the double-integration of base acceleration in my head: as @Torsten noted, I forgot to account for initial conditions. I overlooked the fact that, if you apply a half-sine-pulse of acceleration to the base, then the base experiences only positive acceleration. Therefore the base's final velocity will not equal its intial velocity. Is that really what you want? The base gets "launched" with a sinusoidal push, then moves with a constant velocity after that? This would happen if the experiment were mounted on a train car that was initially at rest, then experienced a brief pulse of acceleration.

Another possibility, which is more like what I imagined, is that the base would get equal amounts of positive and negative acceleration, so that the initial and final velocity of the base would be zero. This would be more like a half-cosine pulse of acceleration.

Assuming the railway-car model, in which acceleration is a positive-only half-sine at the start (i.e. accel=G*sin(wt) when t<pi/w), and assuming yb(0)=yb0 and vb(0)=vb0, we get (as @Torsten showed in his code):

vb(t) = vb0 + (G/w)*(1-cos(wt)) 0<=t<=pi/w

vb(t) = vb0 + 2*G/w otherwise

and

yb(t) = yb0 + (G/w+v0)*t-(G/w^2)*sin(wt) 0<=t<=pi/w

yb(t) = yb0 - G*pi/w^2 + (2*G/w + v0)*t otherwise

[William Rose]

@Mathew Smith,

I like doing it with yb(t) and vb(t) defined as funtions. But I understand that you want to do it by numerical integration of the acceleration pulse. You asked why the most recent code you posted does not work. Here things I would try:

ic() must have 4 elements, not 2.

dydt=zeros(4,1) not zeros(4:1)

Eliminate yb from odefcn() and from the main script.

The y vector is:

y=[yb;vb;y;v];

which looks recursive, but you know what I mean: y and v on the right side are the mass position and velocity. Rewrite the formulas in odefcn(t,y) accordingly.

Define m, c, k inside odefcn().

[William Rose]

For example:

function dydt=odefcn(t,y)
% Next: Unpack y, makes the code below easier to understand
yb=y(1); vb=y(2); ym=y(3); vm=y(4);
m=1; c=1; k=1;      % physical constants (kg, N-s/m, N/m)
G=1;  t1=0.5;       % base acceleration amplitude, duration (m/s^2, s)
w=pi/t1;            % base acceleration frequency (rad/s)
accBase=G*sin(w*t).*(t<t1);         % base acceleration
dydt = [0;0;0;0];
dydt(1)=vb;                         % d(yb)/dt
dydt(2)=accBase;                    % d(vb)/dt
dydt(3)=vm;                         % d(ym)/dt
dydt(4)=(-c*(vm-vb) - k*(ym-yb))/m; % d(vm)/dt
end