# How to solve damped forced vibration analysis using ode45 ?

82 views (last 30 days)
SOMYA RANJAN PATRO on 23 Jul 2021
Commented: Star Strider on 23 Jul 2021
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.

Star Strider on 23 Jul 2021
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
.
Star Strider on 23 Jul 2021
As always, my pleasure!
I am always pleased that I was able to help!
.

### Categories

Find more on Programming in Help Center and File Exchange

R2019b

### Community Treasure Hunt

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

Start Hunting!

Translated by