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

2 views (last 30 days)
Mathew Smith on 25 Feb 2023
Commented: William Rose on 27 Feb 2023
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);

William Rose on 25 Feb 2023
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 on 26 Feb 2023
Edited: Mathew Smith on 26 Feb 2023
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.

William Rose on 26 Feb 2023
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.
William Rose on 27 Feb 2023
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

### Categories

Find more on Programming in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!