Ode45 Quarter car model spring variation

i want script with in input have a matrix of spring and i'm not able to write this script. I have this error 'Unable to perform assignment because the left and right sides have a different number of elements'. Can someone complete this script, thanks

clc;

clear all

close all;

tsim=5 %Time simulation

t=0:0.01:tsim

tprimo=0:0.01:tsim %Time interpolation funzione bump

H0=0.05 %bump high

l=1.5 %bump length

v=10/3.6 %Velocity

T=l/v %Period bump

y0=[0 0 0 0] %initial condition

%Function bump

h=zeros(1,length(tprimo))

c=0

for i=0:0.01:tsim

c=c+1

if i<=T

h(1,c)=2.71828*H0*exp((-1)./(1-((2*v*i/l)-1).^2))

elseif i>T

h(1,c)=0

end

end

t=0:0.01:tsim %Time

%Dati

rig1=[4000 6000 8000] %spring1

damp=150; % damp

mass1=50; % Mass1

mass2=8 % Mass2

rig2=180000; %spring2

opts =odeset('RelTol',1e-5,'AbsTol',[1e-6]);

function dy=equation(t,y,h,rig1,tprimo)

damp=150; % damper

mass1=50; % Mass1

mass2=8 % Mass2

rig2=180000; %spring2

hprimo=interp1(tprimo,h,t) %Interpolation bump

% differential equation

dy=zeros(4,1)

dy(1)=y(2)

dy(3)=y(4)

dy(2)=(-(rig1./mass1).*(y(1)-y(3)))+(-damp./mass1).*(y(2)-y(4))

dy(4)=(rig1./mass2).*(y(1)-y(3))+(-rig2/mass2)*(y(3)-hprimo)+(damp./mass2).*(y(2)-y(4))

end

[t,y]=ode45(@(t,y) equation(t,y,h,rig1,tprimo),[0 2],y0,opts);

ak1=(-rig1./50).*(y(:,1)-y(:,3))-(150/50)*(y(:,2)-y(:,4)) %Acceleration mass1

%plot

figure

plot(t,ak1(:,1),'b',t,ak1(:,2),'r',t,ak1(:,3),'black')

title('acceleration mass 1')

xlabel('Time - s')

ylabel('acceleration (m/s^2)')

grid on

[Aquatris]

Your rig1 variable is 1x3 so in the equation() function

dy(2)=(-(rig1./mass1).*(y(1)-y(3)))+(-damp./mass1).*(y(2)-y(4))

tries to assign 3x1 variable to a 1x1 variable.

Do you want to solve the ODE for 3 different rig1 values? An easy fix is form a for loop as

rig1_vec = [4000 6000 8000];
for idx = 1:length(rig1_vec)
    rig1 = rig1_vec(i);
    %%same code here
    
    ak1(idx) = %equation here
end

Answer 1

[Sam Chak]

Hi @Federico Fortini

It seems like you're looking to plot the acceleration of the first spring-mass system. If that's the case, you can leverage the deval() command to obtain the time derivative of the numeric solution generated by the ode45 solver. That should give you what you need.

To experiment with different spring stiffness values, you can incorporate the suggestion made by @Aquatris using a for-loop approach. That seems like a solid way to explore the parameter space. I'd recommend upgrading @Aquatris' comment to a full-fledged answer, as it provides a nice, constructive solution.

tsim    = 2;            % Time simulation
t       = 0:0.01:tsim;
tprimo  = 0:0.01:tsim;
H0      = 0.05;         % bump high
l       = 1.5;          % bump length
v       = 10/3.6;       % Velocity
T       = l/v;          % Period bump
y0      = [0 0 0 0];    % initial condition

%% Function bump
h       = zeros(1, length(tprimo)); % initialization of vector h
% c       = 0;
% for i = 0:0.01:tsim
%    c = c + 1;
%     if i <= T
%         h(1,c) = 2.71828*H0*exp((-1)./(1-((2*v*i/l)-1).^2));
%     elseif i > T
%         h(1,c) = 0;
%     end
% end 

% preferred method of looping
for i = 1:numel(t)
    if t(i) <= T
        h(i) = exp(1)*H0*exp(- 1/(1 - ((2*v*t(i)/l) - 1)^2));
    else
        h(i) = 0;
    end
end 

t       = 0:0.01:tsim;  %Time <-- unused variable if [t, y] = ode45()

% Dati
% rig1    = [4000 6000 8000]; %spring1
rig1    = 4000;     % spring 1
damp    = 150;      % damp
mass1   = 50;       % Mass1
mass2   = 8;        % Mass2
rig2    = 180000;   % spring 2
opts    = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);

function dy = equation(t, y, h, rig1, tprimo)
    damp    = 150;                      % damper
    mass1   = 50;                       % Mass1
    mass2   = 8;                        % Mass2
    rig2    = 180000;                   % spring2
    hprimo  = interp1(tprimo, h, t);    % Interpolation bump
    
    % differential equation 
    dy      = zeros(4,1);
    dy(1)   = y(2);
    dy(2)   = - (rig1/mass1)*(y(1) - y(3)) + (- damp/mass1)*(y(2) - y(4));
    dy(3)   = y(4);
    dy(4)   =   (rig1/mass2)*(y(1) - y(3)) + (- rig2/mass2)*(y(3) - hprimo) + (damp/mass2)*(y(2) - y(4));
end

% [t, y]= ode45(@(t, y) equation(t, y, h, rig1, tprimo), [0 2], y0, opts);
sol     = ode45(@(t, y) equation(t, y, h, rig1, tprimo), [0 2], y0, opts);
[y, yp] = deval(sol, t);

figure(1)
plot(tprimo, h, 'linewidth', 2, 'color', '#808080'), grid on, ylim([0 1])
xlabel('t'), ylabel('Bump height'), title('Road profile of the Speed Bump')

figure(2)
subplot(211)
plot(t, y),  grid on, ylim([-0.6 0.6])
xlabel('t'), ylabel('\bf{y}'), title('Time responses of the states')
subplot(212)
plot(t, yp(2,:)), grid on, 
xlabel('t'), ylabel('a / (m/s^{2})'), title('Acceleration of Mass 1')