Vector dimensions in my ode

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?

Answer 1

[Alan Stevens]

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

[Donghun Lee]

Yes, this is exactly what I want. Thank you so much for your effort Alan and sorry for being late.

Answer 2

[Alan Stevens]

The problem is here:

f = @(t,x) [ x(2); -( k_s*(x(1)-h)/ m ) ];

When you call this function, x(2) is a single value, but since h is a large vector the the term -( k_s*(x(1)-h)/ m ) has a large number of elements, so you are attempting to create a matrix with one column in the first row and many columns in the second row. Matlab is objecting to this!

I guess you need to make f a function of h as well, and pass in the appropriate single value when you call it.

[Donghun Lee]

thank you for your answer Alan. However, how am I suppose to make f a function of h?

Thanks.

[Alan Stevens]

I think you might have to change from an in-line function to something like

function dxdt = f(t,x,h)

        hs = ???   % select a specific value of h from the vector h
                   % corresponding to the particular time or position
                   % as appropriate.
        dxdt = [ x(2); -( k_s*(x(1)-hs)/ m ) ];        

end

and then call this from ode45 as

[t, x] = ode45(@f,[100,T],x0,[],h);

[Donghun Lee]

Thanks Alan. Can you please write the full code for me if you don't mind? I'm still confused how to apply your code into mine.

I'm so sorry to ask you again.

Thanks.

[Alan Stevens]

Ok, but won't get round to doing it for a couple of hours.

[Donghun Lee]

Yes!! I can definitely wait! Thank you so much😭

[Walter Roberson]

I do not see that I can agree. You posted your question not just once but three times (at least) because you could not wait, and even though I already gave you guidance on exactly the same matter. It does not appear that you are able to wait.

Reminder: you posted during the middle of the night in North America where most of the active volunteers are. If you are in such a rush that you need to post the same question multiple times within a few hours despite having received responses, then you should probably be hiring a consulting company that operates on a schedule that aligns with your needs.

[Donghun Lee]

Walter, it was because I was not happy with your answer honestly. And I said thank you and tried to be nice as much as I can although I could not understand of your answer. Most importantly, it is not appropriate to argue about this in this public area. Hope you understand my position. And thank you for your advice.

[Donghun Lee]

Plus, you said I could not wait but I had been trying to use your answer and spent for 4 hours at least. I tried to listen to your answer and tried my best. I’m so sorry if I make you uncomfortable for this. Have a good day.

[Walter Roberson]

No, I do not understand your position. When you have questions about a response, ask the questions and wait a reasonable time for the response taking into consideration that you are dealing with volunteers in different time zones. If the response time you get from volunteers is not sufficient for your purposes, then you can politely ask for follow-up (remembering that the volunteers do not operate 24 hours a day). If you need an answer within four hours then hire a consultant rather than posting the same question multiple times.

When you post the same question three times within a few hours you are asking for priority treatment that is unfair to other people.

[Donghun Lee]

Ok it’s my fault. I feel sorry to make this happen Walter. Thank you for your advice.

[madhan ravi]

If you don’t understand something, you can always ask a follow up question, in the SAME thread. Asking the same question multiple times will only be time consuming and will not yield any benefits. Taking into consideration people here are VOLUNTEERS not assignment/problem solvers.

[Walter Roberson]

I cannot go back and pull out what I already posted to illustrate any last few details, because two earlier threads on the topic were deleted :( Two volunteers spent their time explaining the difficulty to you, and you deleted what they wrote, wasting their time :(