Azzera filtri
Azzera filtri

How do i add a changing input over a interval using ODE45?

2 visualizzazioni (ultimi 30 giorni)
I'm trying to intergrate form time 0 to 5000 and at t=1000 i want my y input to rise, modeling a car going up a speed bump. I know when the car is over the speed bump (assuming v is consant in the x direction). I tried to make a liniear line of which represents the speed bump. But it doesnt work.
v = 25/3.6; % Speed of the car
Lslope = (0.5/v)*1000; % Length of the slope in time (assuming v is constant)
Tend = 5000; % The end
y0 = 0;
global Tslope
Tslope = 999; % Start time of the slope
global Tlin
Tlin = Tslope + Lslope; % The end time of the slope (assuming v is constant)
[t, y]= ode45('yslope',[0 Tend],[0, 0, 0, 0]);
figure
plot(t, y)
My function is:
function dy = yslope(t,y)
% Constants
mb = 36;
mv = 500;
kb = 100000;
kv = 3000;
b = 5000;
a = 0.6;
x0 = 0;
% Aanroepen
global Tslope
global Tlin
% y input
y_before = 0;
y_slope = linspace(0,0.3,Tlin-Tslope);
y_after = 0.3;
% Checken whick state
if t <= Tslope % Before the slope
y_in = y_before;
elseif t > Tslope && t <= Tlin % At slope
y_in = yslope(t-Tslope); % t is already at 1001 so minus Tslope to get 1
elseif t > Tlin % After slope
y_in = y_after;
end
dy1 = y(2);
dy2 = -(kv+kb)/mb *y(1)-b/mb *y(2)+kv/mb *y(3)+b/mb *y(4) +kb/mb *y_in;
dy3 = y(4);
dy4 = kv/mv *y(1) + b/mv *y(2) - kv/mv *y(3) -b/mv*y(4);
dy= [dy1 ; dy2 ; dy3 ;dy4];
end
This is the error i'm getting:
Not enough input arguments.
Error in yslope (line 37)
dy1 = y(2);
Error in yslope (line 29)
y_in = yslope(t-Tslope); % t is already at 1001 so minus Tslope to get 1
Error in ode45 (line 308)
f6 = ode(t6, y6);
Error in Main (line 15)
[t, y]= ode45('yslope',[0 Tend],[0, 0, 0, 0]);

Risposte (2)

Torsten
Torsten il 3 Lug 2023
Modificato: Torsten il 3 Lug 2023
v = 25/3.6; % Speed of the car
Lslope = (0.5/v)*1000; % Length of the slope in time (assuming v is constant)
Tend = 5000; % The end
y0 = 0;
global Tslope
Tslope = 999; % Start time of the slope
global Tlin
Tlin = Tslope + Lslope; % The end time of the slope (assuming v is constant)
%[t, y]= ode45('yslope',[0 Tend],[0, 0, 0, 0]);
[t, y]= ode45(@yslope,[0 Tend],[0, 0, 0, 0]);
figure
plot(t, y)
function dy = yslope(t,y)
% Constants
mb = 36;
mv = 500;
kb = 100000;
kv = 3000;
b = 5000;
a = 0.6;
x0 = 0;
% Aanroepen
global Tslope
global Tlin
% y input
y_before = 0;
%y_slope = linspace(0,0.3,Tlin-Tslope);
y_after = 0.3;
% Checken whick state
%if t <= Tslope % Before the slope
% y_in = y_before;
%elseif t > Tslope && t <= Tlin % At slope
% y_in = yslope(t-Tslope); % t is already at 1001 so minus Tslope to get 1
%elseif t > Tlin % After slope
% y_in = y_after;
%end
y_in = y_before*(t<=Tslope) + (y_before*(t-Tlin)/(Tslope-Tlin) + y_after*(t-Tslope)/(Tlin-Tslope))*(t>Tslope)*(t<=Tlin)+y_after*(t>Tlin);
dy1 = y(2);
dy2 = -(kv+kb)/mb *y(1)-b/mb *y(2)+kv/mb *y(3)+b/mb *y(4) +kb/mb *y_in;
dy3 = y(4);
dy4 = kv/mv *y(1) + b/mv *y(2) - kv/mv *y(3) -b/mv*y(4);
dy= [dy1 ; dy2 ; dy3 ;dy4];
end

Sam Chak
Sam Chak il 3 Lug 2023
I am unsure if I understand your 4-state differential equations. However, a speed bump can be mathematically modeled as a physical disturbance to the car.
For example, the cross-section of a one-part speed bump can be plotted as follows:
Try modifying the state equations to insert the speed bump function.
hw = 150; % half-width of the bump (mm)
x = linspace(-2*hw, 2*hw, 60001);
y1 = exp(-1./(1 - (x/hw).^2)); % component 1
y2 = (sign(x + hw) - sign(x - hw))/2; % component 2
H = 50; % max height of the bump (mm)
A = H/0.367879441171442;
y = A*y1.*y2; % Speed Bump function
plot(x, y, 'linewidth', 1.5), grid on, ylim([0 600])
title('Speed Bump')
xlabel({'$x$/mm'}, 'interpreter', 'latex')
ylabel({'$y$/mm'}, 'interpreter', 'latex')
  1 Commento
Dirk te Brake
Dirk te Brake il 3 Lug 2023
Well my problem isn't really a speed bump. After the car goes up the slope it doesn't come down again. But i do like the function and will give it a try.

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by