Random road profile into ode45

clc,clear all

k_s = 26400;
m = 483; %Mass
f_n = sqrt(k_s/m)/(2*pi); %Natural frequency in Hz

%% Road profile

% spatial frequency (n0) cycles per meter
Omega0 = 0.1; %%%%conventional value of spatial frequency(n0)?
% psd ISO (used for formula 8)
Gd_0 =  32 * (10^-6);
% waveviness
w = 2; 
% road length
L = 250;

%delta n 
N = 1000; 
Omega_L = 0.004;
Omega_U = 4;

delta_n =  1/L; % delta_n = (Omega_U - Omega_L)/(N-1);

% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;
 
%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);

% calculate amplitude using formula(8) in the article
%Amp = sqrt(2*Gd*delta_n); %%%from Eq. 7?

%calculate amplitude using simplified formula(9) in the article
k = 3; %%%upper limit A and lower limit B k=3?
%Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);

%random phases 
Psi = 2*pi*rand(size(Omega)); 

% x abicsa from 0 to L
x1 = 0:0.25:250;
% road sinal
h= zeros(size(x1));
plot(x1, h*1000);

%% ode45
T = 110;
x0 = [0,0];
 
for iv=1:length(x1)
    h(iv) = sum( Amp.*cos(2*pi*Omega*x1(iv) + Psi) );
    
        f = @(t,x) [ x(2); -( k_s*(x(1)-1000*h(iv))/ m ) ];
        
        [t, x] = ode45(f,[100,T],x0);
    
end

plot(t,x(:,1));

Hi, I generated a random road profile (h) and it gives such a figure as below,

As you can see, there are quite lots of elevation in this road profile. However, if I apply this into my ode45, it gives a graph as below,

this seems to be a sinosoidal graph, but I think this graph should not be like this because the road has lots of elevation during the journey.

Can any one tell me what is wrong with my code please?

Answer 1

[Drishti]

Hi Donghun,

I understand that you are trying to plot a road profile followed by analyzing the system's response.

While reproducing the provided code on my end in MATLAB R2024b, I found that the potential reason for the system's response as a sinusoidal wave is because the 'ode45' is called inside the 'for' loop, recalculating the system's response at each point in 'x1'.

As a result, the simulation restarts at every step, whereas it should be integrated over the entire time span.

For better understanding, you can refer to the code snippet below:

for iv = 1:length(x1)
    h(iv) = sum(Amp .* cos(2*pi*Omega*x1(iv) + Psi));
end
plot(x1, h*1000);
xlabel('Distance (m)');
ylabel('Elevation (mm)');
title('Road Profile');
%% ODE45
T = 110;
x0 = [0, 0];
f = @(t, x) [x(2); -(k_s * (x(1) - 1000 * interp1(x1, h, t)) / m)];
[t, x] = ode45(f, [0, T], x0);

For reference purpose, I have attached the output graph of system's response below :

I hope this resolves the issue.