MATLAB: Vector dimensions in the ode

ode45vectorvectorization

clc,clear all
k_s = 26400; %spring stiffness
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;
h= zeros(size(x1));
%artificial random road profile
for iv=1:length(x1)
    h(iv) = sum( Amp.*cos(2*pi*Omega*x1(iv) + Psi) );
end
%% ode45
T = 120;
x0 = [0,0];
       
        f = @(t,x) [ x(2); -( k_s*(x(1)-h)/ m ) ];
        
        [t, x] = ode45(f,[100,T],x0);
%% plot
plot(t,x(:,1));
set(gca,'xtick',17)
Hi, I generated a random road file (h) and tried to apply this in my ode, however it says the vector dimension is not consistent. Can anyone solve this problem please?

Best Answer

I think the following produces somewhat more sensible results. I was confused for some time because you use x for both distance along the road and for vertical displacement. I've changed the latter to y. See if this does what you want:
k_s = 26400; %spring stiffness

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 = 100; 
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:250/(N-1):250;
h= zeros(size(x1));
%artificial random road profile
for iv=1:length(x1)
    h(iv) = sum( Amp.*cos(2*pi*Omega*x1(iv) + Psi) );
end
hx = [x1' h'];
%% ode45
y0 = [0,0];
       
        
        [t, y] = ode45(@f,x1,y0,[],hx);
%% plot
figure
plot(t,y(:,1));
xlabel('time'),ylabel('vertical displacement')
function dydt = f(t,y,hx)
          k_s = 26400; %spring stiffness
          m = 483; %Mass
          v = 1; % speed along road  
          x = v*t;
          hs = hfn(x,hx);
          dydt =[y(2);
                 -( k_s*(y(1)-hs)/ m )];
         
end
function hs = hfn(x, hx)
         hs = interp1(hx(:,1),hx(:,2),x);
         
         
end
Related Question