How to solve mass spring damper over frequency (differential equation)?

This is a mass-spring-damper system equation with a function dependent on frequency (M, R, & K in parethesis).

f = 50:450; % Hz
M = 23.3;
K = 3.61 * 10^7;
R = 1360;
P = -1 * (2 * pi * f);
ic = 2 * pi * 50; % rad/s

The figure shows the linear impedance frequency response of the system defined by the equation.

How do you go from the differential equation to the figure?

I've seen differential equation solvers using time & displacement, but no luck on just frequency response.

[madhan ravi]

Initial condition?

[KSSV]

Read about ode45.

[Keith Grey]

Initial condition is 50Hz.

[Rafael Hernandez-Walls]

You need initial condition for the displacement and for the velocity. Then you need write the ODE in two equations of first orden. Lke this

and solve using ode45

Answer 1

[Gifari Zulkarnaen]

Equation of motion is time-dependant equation. Perhaps what you mean is Response Spectrum (frequency-based response data)? This is calculation for acceleration response spectrum:

clear
close all
f = 50:10:450; % Hz
M = 23.3;
K = 3.61*10^7;
R = 1360;

RS = zeros(1,length(f));
for i=1:length(f)
    tspan = [0 1]; % Frequency range
    Y0 = zeros(2,1); % Initial zero condition
    [t,Y] = ode45(@(t,Y) StateEqOfMotion(t,Y,M,K,R,f(i)),tspan,Y0); % ODE
    % Output Y is in term of displacement and velocity
    a = Y(:,2)./t; % Acceleration
    RS(i) = max(abs(a)); % Acceleration response spectrum
end
plot(f,RS)

% ODE function of State Equation
function dYdt = StateEqOfMotion(t,Y,M,K,R,f)
    P = -cos(t*f);
    x = Y(1);
    v = Y(2);
    a = (-(R*v+K*x)+P)/M;
    dYdt = [v; a];
end

The result is still different with what your graph, but maybe this can give insight

Answer 2

[Swami]

clear

close all

f = 50:10:450; % Hz

M = 23.3;

K = 3.61*10^7;

R = 1360;

RS = zeros(1,length(f));

for i=1:length(f)

tspan = [0 1]; % Frequency range

Y0 = zeros(2,1); % Initial zero condition

[t,Y] = ode45(@(t,Y) StateEqOfMotion(t,Y,M,K,R,f(i)),tspan,Y0); % ODE

% Output Y is in term of displacement and velocity

a = Y(:,2)./t; % Acceleration

RS(i) = max(abs(a)); % Acceleration response spectrum

end

plot(f,RS)

% ODE function of State Equation

function dYdt = StateEqOfMotion(t,Y,M,K,R,f)

P = -cos(t*f);

x = Y(1);

v = Y(2);

a = (-(R*v+K*x)+P)/M;

dYdt = [v; a];

end