Questioning the accuracy of Coupled differential equations with forcing function

Hello, I am trying to solve and plot a series of coupled second order DEs. Its a standard spring two-mass damper system with an input rectangular impulse function:

(1365000*(1365000*heaviside(t - 0.01))/m1

As you can see this input function depends on t which is not a variable in my system of DE's. Any way I am wondering if I should question whether the plot is plotting as I want it to even though it seems to be working (plot looks like its oscillating properly). And also why y2 (the body of the car) is zero, I would think this would move. The math model I am using can be seen at the very end of this video:

(https://www.youtube.com/watch?v=TKbChGz4-mw)

Here is the code:

    function dydt = forreal(t,y)
  m1=2 %(wheel/suspension)
  m2=136.07  %(vehicle body)
  b=3000
  k1=40000000  %wheel "spring constant"
  k2=10000
  eqa = (1365000-(1365000*heaviside(t-0.01)))/m1
  eq1 = (k2/m1)*y(3) + (-k2/m1)*y(1) + (b/m1)*y(4) - (b/m1)*y(2) - (k1/m1)*y(1) + eqa;
  eq2 = (-k2/m2)*y(3) + (k2/m2)*y(1) - (b/m2)*y(4) + (b/m2)*y(2);
  dydt = [y(2); eq1; y(3); eq2];
  end

and calling in a different script:

[t,y] = ode45(@forreal,[0 5],[0; 0; 0; 0]);

plot(t,y(:,1),'-o',t,y(:,3),'-o')

xlabel('Time t');
ylabel('displacement');
legend('wheel displacement','body displacement')

here are the results:

So my question is, in my research matlab seemed to require a separate technique to plot a time dependent diff. eq.

(https://www.mathworks.com/matlabcentral/answers/97074-how-do-i-solve-an-ode-with-time-dependent-parameters-in-matlab)

However I saw one tutorial where they didn't do that. I gave it a shot and my code generates a plot, but I am wondering if I should question its accuracy, and if its really plotting what I am intending. I am also not sure why y2 is zero. This should be representing the body of the vehicle, which I would think would move.

Any help would be great.

[Bish Erbas]

Have you tried changing the initial values to non-zero values?

[Torsten]

You mean

dydt = [y(2); eq1; y(4); eq2];

?

And integrating from t=0 to t=0.01 and restarting the integrator at t=0.01 will make the results more accurate.

Best wishes

Torsten.

[Ryan Compton]

Thanks that helped. Changed it to y(4).

[Ryan Compton]

I have updated the code. I am still hoping someone can check that the code is actually doing what I want it to. The plot surely looks like a spring response, but my use of the forcing function concerns me since it throws the variable "t" in the mix. Can ode45 handle this as is or do I need to do something special for the time-dependent impulse force?

[Torsten]

[t1,y1] = ode45(@(t,y)forreal(t,y,1),[0 0.01],[0; 0; 0; 0]);
[t2,y2] = ode45(@(t,y)forreal(t,y,2),[0.01 5],[y1(end,1);y1(end,2);y1(end,3);y1(end,4)]);
t=[t1;t2];
y=[y1;y2];
plot(t,y(:,1),'-o',t,y(:,3),'-o')

xlabel('Time t');
ylabel('displacement');
legend('wheel displacement','body displacement')
end

function dydt = forreal(t,y,flag)
  m1=2; %(wheel/suspension)
  m2=136.07;  %(vehicle body)
  b=3000;
  k1=40000000;  %wheel "spring constant"
  k2=10000;
  if flag==1
    eqa = 1365000/m1;
  elseif flag==2
    eqa = 0.0;
  end
  eq1 = (k2/m1)*y(3) + (-k2/m1)*y(1) + (b/m1)*y(4) - (b/m1)*y(2) - (k1/m1)*y(1) + eqa;
  eq2 = (-k2/m2)*y(3) + (k2/m2)*y(1) - (b/m2)*y(4) + (b/m2)*y(2);
  dydt = [y(2); eq1; y(4); eq2];
  end

[Ryan Compton]

Thanks. I was wondering, did you change the input force "eqa" to a step function? Just L^-1(F/s) = F

[Ryan Compton]

Wait I just looked at this more carefully. You essentially plotted a step at the correct value and then set the input to 0 after time t. Ok I feel dumb that's way simpler.

Answer 1

[James Tursa]

It's all 0's because that is the solution. Your derivative function is identically 0's in every spot on every call. Look at your equations for rk. Since x1, x1d, x2, and x2d are involved in every term and they all start out as 0's, there is nothing to move the state vector away from the initial all 0's state.

[Bish Erbas]

Thanks for the explanation.