How to solve damped forced vibration analysis using ode45 ?

I am trying to solve equation of motion for a damped forced vibration analysis for a SDOF system. My input Force is function of time having random white noise. I am using ode45 to solve the problem. Here is my code.

% Main code
clear variables
close all
clc

t = 0:20;    % Time duration
x = [1; 6];    % Initial conditions
[t, x] = ode45(@eqm,t,x);    % ode45 command

Below is the function eqm

% Function
function dxdt = eqm(t, x)
F = 1:20        % For simplicity I have taken force as an array
m = 1;            % mass
c = 1;            % damping
k = 1;            % stiffness
dxdt = [x(2); F/m - c/m*x(2) - k*x(1)];
end

After I run this code it shows error "Error using vertcat. Dimensions of arrays being concatenated are not consistent". I checked online but most of the examples contains only periodic functions of force such as F = sin(wt) but I didn't find anything for random load. Any help will be greatly appreciated.

Answer 1

[Star Strider]

One problem is with the ‘dxdt’ assignment. MATLAB interprets the spaces as delimiters, so the second row has 3 columns, as MATLAB sees it. Put parentheses around the second row terms (to preserve readability) and after providing ‘F’ as an argument rather than internally as a vector, and a few other tweaks, it works —

t = 0:20;    % Time duration
x0 = [1; 6];    % Initial conditions
F = 1:20;        % For simplicity I have taken force as an array
for k = 1:numel(F)
[t, x] = ode45(@(t,x)eqm(t,x,F(k)),t,x0);    % ode45 command
tv{k} = t;
xv{k} = x;
end

figure
for k = 1:numel(F)
    subplot(5,4,k)
    plot(tv{k}, xv{k})
    title(sprintf('k = %2d',k))
    grid
end

% Function
function dxdt = eqm(t, x, F)
m = 1;            % mass
c = 1;            % damping
k = 1;            % stiffness
dxdt = [x(2); (F/m - c/m*x(2) - k*x(1))];
end

.

[SOMYA RANJAN PATRO]

Thank you very much sir, i have implemented your suggestions and added both paranthesis and the x0 variable. But the main issue is actually Force variable (F) is a function of time like F(t) = array of 20 values i.e., for each value of time (t) I have a value of force. So, even if I consider an argument of F by considering 'for loop' I will get for each value of F, a range of values of x for t from 0 to 20. Is there a way to implement the (F) variable as a function of time (t). Feel free to ask anything if you have any queries sir.

[Star Strider]

My pleasure!

That ‘F’ is a function of time is not obvious.

This should do what you want —

t = 0:20;    % Time duration
x0 = [1; 6];    % Initial conditions
F = 1:20;        % For simplicity I have taken force as an array
[t, x] = ode45(@(t,x)eqm(t,x),t,x0);    % ode45 command

figure
plot(t,x)
grid
xlabel('t')
ylabel('Amplitude')
legend('x_1', 'x_2', 'Location','best')

% Function
function dxdt = eqm(t, x)
F = 1:20;        % For simplicity I have taken force as an array
Ftv = linspace(0,20,numel(F));  % Time Vector Corresponding To 'F'
Ft = interp1(Ftv, F, t);        % Interpolate To Get 'F(t)'
m = 1;            % mass
c = 1;            % damping
k = 1;            % stiffness
dxdt = [x(2); (Ft/m - c/m*x(2) - k*x(1))];
end

See the documentation section on ODE with Time-Dependent Terms to understand what I did here with the interpolation.

.

[SOMYA RANJAN PATRO]

Thanks again for your valuable time and effort. I gone through the sections on ODE with Time-Dependent Terms and now I finally understood the concept behind how to implement the variable F if it's a function of time. It was really helpful.

[Star Strider]

As always, my pleasure!

I am always pleased that I was able to help!

.

[Amira]

thanks a lot for ur help !! i do appreciate !!

[Star Strider]

My pleasure!

A Vote would be appreciated!