resolution of MDOF using ode45

i have a problem solving the system with ode45. the code works but the displacement in graphed output is not what i would expect from a chirp signal. What could be the error in my code?

clear 
close all
clc
%% STANDING
%Mass 1 
m1=4.5;
k1=865.47;   %ho diviso per 4
k4=865.47;
k7=865.47;
k8=865.47;
c1=12.8;   %ho diviso per 4
c4=12.8;
c7=12.8;
c8=12.8;
%Mass 2 and 8
m2=2.070;
k2=5107.515;
c2=115.16;
m8=2.070;
k9=5107.515;
c9=115.16;
%Mass 3 and 9
m3=1.160;
k3=26468.485;
c3=143.68;
m9=1.160;
k10=26468.485;
c10=143.68;
%Mass 4 and 10
m4=0.540;
m10=0.540;
%Mass 5
m5=15.551;
k5=263068.52;  %ho diviso per 2
k6=263068.52;
c5=715.1;
c6=715.1;
%Mass 6
m6=12.5;
k11=19739.205;
c11=291.878;
%Mass 7
m7=4.5;
k12=30023.33;
c12=350;
%Mass 11
m11=27.174;
k13=21455.7275;  %ho diviso per 2
k14=21455.7275;
c13=280.76;
c14=280.76;
%Mass 12 and 13
m12=7.165;
k15=99837807.76;
c15=2376.4;
m13=7.165;
k16=99837807.76;
c16=2376.4;
%Mass 14 and 15
m14=2.8;
k17=1727180.77;
c17=2376.4;
m15=2.8;
k18=1727180.77;
c18=2376.4;
%Mass 16 and 17
m16=1.170;
k19=4127.21;
c19=2.084;
m17=1.170;
k20=4127.21;
c20=2.084;
%Mass 18
m18=3.8;
k23=457384621.1;
c23=13349.82;
%Mass 19
m19=1.336;
k21=5242.93;  %ho diviso per 2
k22=5242.93;
c21=23.67;
c22=23.67;
%%

%MATRIX
M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19]);

% stiffness matrix 19x19
K = zeros(19,19);
K(1,1) = k1 + k4 + k7 + k8;
K(1,2) = -k1;
K(1,5) = -k4;
K(1,7) = -k7;
K(1,8) = -k8;
K(2,1) = -k1;
K(2,2) = k1 + k2;
K(2,3) = -k2;
K(3,2) = -k2;
K(3,3) = k2 + k3;
K(3,4) = -k3;
K(4,3) = -k3;
K(4,4) = k3;
K(5,1) = -k4;
K(5,5) = k4 + k5 + k6;
K(5,6) = -k5;
K(5,7) = -k6;
K(6,5) = -k5;
K(6,6) = k5 + k11;
K(6,11) = -k11;
K(7,1) = -k7;
K(7,5) = -k6;
K(7,7) = k6 + k7 + k23;
K(7,18) = -k23;
K(8,1) = -k8;
K(8,8) = k8 + k9;
K(8,9) = -k9;
K(9,8) = -k9;
K(9,9) = k9 + k10;
K(9,10) = -k10;
K(10,9) = -k10;
K(10,10) = k10;
K(11,6) = -k11;
K(11,7) = -k12;
K(11,11) = k11 + k12 + k13 + k14;
K(11,12) = -k13;
K(11,13) = -k14;
K(12,11) = -k13;
K(12,12) = k13 + k15;
K(12,14) = -k15;
K(13,11) = -k14;
K(13,13) = k14 + k16;
K(13,15) = -k16;
K(14,12) = -k15;
K(14,14) = k15 + k17;
K(14,16) = -k17;
K(15,13) = -k16;
K(15,15) = k16 + k18;
K(15,17) = -k18;
K(16,14) = -k17;
K(16,16) = k17 + k19;
K(17,15) = -k18;
K(17,17) = k18 + k20;
K(18,7) = -k23;
K(18,18) = k23 + k21 + k22;
K(18,19) = -k21 - k22;
K(19,18) = -k21 - k22;
K(19,19) = k21 + k22;

% damping matrix 19x19
C = zeros(19,19);
C(1,1) = c1 + c4 + c7 + c8;
C(1,2) = -c1;
C(1,5) = -c4;
C(1,7) = -c7;
C(1,8) = -c8;
C(2,1) = -c1;
C(2,2) = c1 + c2;
C(2,3) = -c2;
C(3,2) = -c2;
C(3,3) = c2 + c3;
C(3,4) = -c3;
C(4,3) = -c3;
C(4,4) = c3;
C(5,1) = -c4;
C(5,5) = c4 + c5 + c6;
C(5,6) = -c5;
C(5,7) = -c6;
C(6,5) = -c5;
C(6,6) = c5 + c11;
C(6,11) = -c11;
C(7,1) = -c7;
C(7,5) = -c6;
C(7,7) = c6 + c7 + c23;
C(7,18) = -c23;
C(8,1) = -c8;
C(8,8) = c8 + c9;
C(8,9) = -c9;
C(9,8) = -c9;
C(9,9) = c9 + c10;
C(9,10) = -c10;
C(10,9) = -c10;
C(10,10) = c10;
C(11,6) = -c11;
C(11,7) = -c12;
C(11,11) = c11 + c12 + c13 + c14;
C(11,12) = -c13;
C(11,13) = -c14;
C(12,11) = -c13;
C(12,12) = c13 + c15;
C(12,14) = -c15;
C(13,11) = -c14;
C(13,13) = c14 + c16;
C(13,15) = -c16;
C(14,12) = -c15;
C(14,14) = c15 + c17;
C(14,16) = -c17;
C(15,13) = -c16;
C(15,15) = c16 + c18;
C(15,17) = -c18;
C(16,14) = -c17;
C(16,16) = c17 + c19;
C(17,15) = -c18;
C(17,17) = c18 + c20;
C(18,7) = -c23;
C(18,18) = c23 + c21 + c22;
C(18,19) = -c21 - c22;
C(19,18) = -c21 - c22;
C(19,19) = c21 + c22;

n=19;
y0 = zeros(2*n,1);
tspan = [0 120];

% ode45
[t, y] = ode45(@(t, y) odefcn_standing(t, y, M, C, K), tspan, y0);
figure;
plot(t, y(:, 19));   
xlabel('Time (s)');
ylabel('Displacement (m)');
% legend('y1', 'y2', 'y3');
title('response of the system 19DOF');
grid on;

function dy = odefcn_standing(t, y, M, C, K)
    n = 19; % Numero di gradi di libertà
    dy = zeros(2 * n, 1);
    
    % Construction of matrix A
    A = [zeros(n), eye(n);
        -inv(M) * K, -inv(M) * C];
    
    F = zeros(19, 1);
    f0 = 0.5; % initial frequency
    f1 = 80; % final frequency
    t_f = 120; % duration of chirp signal
    chirp_signal = chirp(t, f0, t_f, f1);
    F(16,:) = 10*chirp_signal; % on mass 16
    F(17,:) = 10*chirp_signal; % on mass 17
    
    % Construction of matrix B
    B = [zeros(n, n); inv(M)];
    dy = A * y + B * F;
end

[Shivam Gothi]

Hello @Roberta Pili,

What are the values of variables k1, k2, k3 ... C1, C2, C3.... and m1, m2, m3......?

These variables are undefined and I am not able to run the above code

[Roberta Pili]

clear

close all

clc

%% STANDING

%Mass 1

m1=4.5;

k1=865.47; %ho diviso per 4

k4=865.47;

k7=865.47;

k8=865.47;

c1=12.8; %ho diviso per 4

c4=12.8;

c7=12.8;

c8=12.8;

%Mass 2 and 8

m2=2.070;

k2=5107.515;

c2=115.16;

m8=2.070;

k9=5107.515;

c9=115.16;

%Mass 3 and 9

m3=1.160;

k3=26468.485;

c3=143.68;

m9=1.160;

k10=26468.485;

c10=143.68;

%Mass 4 and 10

m4=0.540;

m10=0.540;

%Mass 5

m5=15.551;

k5=263068.52; %ho diviso per 2

k6=263068.52;

c5=715.1;

c6=715.1;

%Mass 6

m6=12.5;

k11=19739.205;

c11=291.878;

%Mass 7

m7=4.5;

k12=30023.33;

c12=350;

%Mass 11

m11=27.174;

k13=21455.7275; %ho diviso per 2

k14=21455.7275;

c13=280.76;

c14=280.76;

%Mass 12 and 13

m12=7.165;

k15=99837807.76;

c15=2376.4;

m13=7.165;

k16=99837807.76;

c16=2376.4;

%Mass 14 and 15

m14=2.8;

k17=1727180.77;

c17=2376.4;

m15=2.8;

k18=1727180.77;

c18=2376.4;

%Mass 16 and 17

m16=1.170;

k19=4127.21;

c19=2.084;

m17=1.170;

k20=4127.21;

c20=2.084;

%Mass 18

m18=3.8;

k23=457384621.1;

c23=13349.82;

%Mass 19

m19=1.336;

k21=5242.93; %ho diviso per 2

k22=5242.93;

c21=23.67;

c22=23.67;

%%

%MATRCS

M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19]);

% Stiffness matrix initialisation 19x19

K = zeros(19,19);

% Definition of stiffness matrix elements K

% For m1 (x1)

K(1,1) = k1 + k4 + k7 + k8;

K(1,2) = -k1;

K(1,5) = -k4;

K(1,7) = -k7;

K(1,8) = -k8;

% For m2 (x2)

K(2,1) = -k1;

K(2,2) = k1 + k2;

K(2,3) = -k2;

% For m3 (x3)

K(3,2) = -k2;

K(3,3) = k2 + k3;

K(3,4) = -k3;

% For m4 (x4)

K(4,3) = -k3;

K(4,4) = k3;

% For m5 (x5)

K(5,1) = -k4;

K(5,5) = k4 + k5 + k6;

K(5,6) = -k5;

K(5,7) = -k6;

% For m6 (x6)

K(6,5) = -k5;

K(6,6) = k5 + k11;

K(6,11) = -k11;

% For m7 (x7)

K(7,1) = -k7;

K(7,5) = -k6;

K(7,7) = k6 + k7 + k23;

K(7,18) = -k23;

% For m8 (x8)

K(8,1) = -k8;

K(8,8) = k8 + k9;

K(8,9) = -k9;

% For m9 (x9)

K(9,8) = -k9;

K(9,9) = k9 + k10;

K(9,10) = -k10;

% For m10 (x10)

K(10,9) = -k10;

K(10,10) = k10;

% For m11 (x11)

K(11,6) = -k11;

K(11,7) = -k12;

K(11,11) = k11 + k12 + k13 + k14;

K(11,12) = -k13;

K(11,13) = -k14;

% For m12 (x12)

K(12,11) = -k13;

K(12,12) = k13 + k15;

K(12,14) = -k15;

% For m13 (x13)

K(13,11) = -k14;

K(13,13) = k14 + k16;

K(13,15) = -k16;

% For m14 (x14)

K(14,12) = -k15;

K(14,14) = k15 + k17;

K(14,16) = -k17;

% For m15 (x15)

K(15,13) = -k16;

K(15,15) = k16 + k18;

K(15,17) = -k18;

% For m16 (x16)

K(16,14) = -k17;

K(16,16) = k17 + k19;

% For m17 (x17)

K(17,15) = -k18;

K(17,17) = k18 + k20;

% For m18 (x18)

K(18,7) = -k23;

K(18,18) = k23 + k21 + k22;

K(18,19) = -k21 - k22;

% For m19 (x19)

K(19,18) = -k21 - k22;

K(19,19) = k21 + k22;

% Damping matrix initialisation 19x19

C = zeros(19,19);

% Definition of damping matrix elements C

% For m1 (x1)

C(1,1) = c1 + c4 + c7 + c8;

C(1,2) = -c1;

C(1,5) = -c4;

C(1,7) = -c7;

C(1,8) = -c8;

% For m2 (x2)

C(2,1) = -c1;

C(2,2) = c1 + c2;

C(2,3) = -c2;

% For m3 (x3)

C(3,2) = -c2;

C(3,3) = c2 + c3;

C(3,4) = -c3;

% For m4 (x4)

C(4,3) = -c3;

C(4,4) = c3;

% For m5 (x5)

C(5,1) = -c4;

C(5,5) = c4 + c5 + c6;

C(5,6) = -c5;

C(5,7) = -c6;

% For m6 (x6)

C(6,5) = -c5;

C(6,6) = c5 + c11;

C(6,11) = -c11;

% For m7 (x7)

C(7,1) = -c7;

C(7,5) = -c6;

C(7,7) = c6 + c7 + c23;

C(7,18) = -c23;

% For m8 (x8)

C(8,1) = -c8;

C(8,8) = c8 + c9;

C(8,9) = -c9;

% For m9 (x9)

C(9,8) = -c9;

C(9,9) = c9 + c10;

C(9,10) = -c10;

% For m10 (x10)

C(10,9) = -c10;

C(10,10) = c10;

% For m11 (x11)

C(11,6) = -c11;

C(11,7) = -c12;

C(11,11) = c11 + c12 + c13 + c14;

C(11,12) = -c13;

C(11,13) = -c14;

% For m12 (x12)

C(12,11) = -c13;

C(12,12) = c13 + c15;

C(12,14) = -c15;

% For m13 (x13)

C(13,11) = -c14;

C(13,13) = c14 + c16;

C(13,15) = -c16;

% Per m14 (x14)

C(14,12) = -c15;

C(14,14) = c15 + c17;

C(14,16) = -c17;

% For m15 (x15)

C(15,13) = -c16;

C(15,15) = c16 + c18;

C(15,17) = -c18;

% For m16 (x16)

C(16,14) = -c17;

C(16,16) = c17 + c19;

% For m17 (x17)

C(17,15) = -c18;

C(17,17) = c18 + c20;

% For m18 (x18)

C(18,7) = -c23;

C(18,18) = c23 + c21 + c22;

C(18,19) = -c21 - c22;

% For m19 (x19)

C(19,18) = -c21 - c22;

C(19,19) = c21 + c22;

Answer 1

[Ayush]

Hi Roberta

I could not run the code as many variables are not defined. But here are the debugging steps that you may try:

1. Starting Conditions: Your system starts with everything at rest since the initial conditions are all zeros. Double-check that this is what you want for your simulation.
2. Mass Matrix: The code uses inv(M) to compute the inverse of the mass matrix. If the mass matrix is singular or nearly singular, this could lead to numerical instability. Ensure that the mass values (m1 to m19) are set correctly and are non-zero.
3. Chirp Signal: The chirp signal is applied to masses 16 and 17. Make sure the frequency range and duration are appropriate for what you're trying to simulate. If the signal isn't strong enough or doesn't cover the right frequencies, it might not excite the system as expected.
4. Damping: Check the damping values. If there's too much damping, it might dampen the response too much, or if there's too little, you might see excessive oscillations
5. Simulation Time: The time span you've set might need tweaking. If your system's dynamics happen over a different time scale, consider adjusting the time span or the solver's settings for better accuracy.
6. Additionally, you can use the odeset function to specify options like relative and absolute tolerances to improve the solver's accuracy.
7. Numerical Stability: If you're still having trouble, you might want to try a different solver or adjust some of the ode45 settings to see if that helps.
8. Plotting the Right Thing: Finally, make sure you're plotting the right part of the response. You're currently looking at the 19th degree of freedom, so just double-check that this is what you're interested in.

Let me know if this was helpful!

[Roberta Pili]

Thank you. These are the data

clear

close all

clc

%% STANDING

%Mass 1

m1=4.5;

k1=865.47; %ho diviso per 4

k4=865.47;

k7=865.47;

k8=865.47;

c1=12.8; %ho diviso per 4

c4=12.8;

c7=12.8;

c8=12.8;

%Mass 2 and 8

m2=2.070;

k2=5107.515;

c2=115.16;

m8=2.070;

k9=5107.515;

c9=115.16;

%Mass 3 and 9

m3=1.160;

k3=26468.485;

c3=143.68;

m9=1.160;

k10=26468.485;

c10=143.68;

%Mass 4 and 10

m4=0.540;

m10=0.540;

%Mass 5

m5=15.551;

k5=263068.52; %ho diviso per 2

k6=263068.52;

c5=715.1;

c6=715.1;

%Mass 6

m6=12.5;

k11=19739.205;

c11=291.878;

%Mass 7

m7=4.5;

k12=30023.33;

c12=350;

%Mass 11

m11=27.174;

k13=21455.7275; %ho diviso per 2

k14=21455.7275;

c13=280.76;

c14=280.76;

%Mass 12 and 13

m12=7.165;

k15=99837807.76;

c15=2376.4;

m13=7.165;

k16=99837807.76;

c16=2376.4;

%Mass 14 and 15

m14=2.8;

k17=1727180.77;

c17=2376.4;

m15=2.8;

k18=1727180.77;

c18=2376.4;

%Mass 16 and 17

m16=1.170;

k19=4127.21;

c19=2.084;

m17=1.170;

k20=4127.21;

c20=2.084;

%Mass 18

m18=3.8;

k23=457384621.1;

c23=13349.82;

%Mass 19

m19=1.336;

k21=5242.93; %ho diviso per 2

k22=5242.93;

c21=23.67;

c22=23.67;

%%

%MATRCS

M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19]);

% Stiffness matrix initialisation 19x19

K = zeros(19,19);

% Definition of stiffness matrix elements K

% For m1 (x1)

K(1,1) = k1 + k4 + k7 + k8;

K(1,2) = -k1;

K(1,5) = -k4;

K(1,7) = -k7;

K(1,8) = -k8;

% For m2 (x2)

K(2,1) = -k1;

K(2,2) = k1 + k2;

K(2,3) = -k2;

% For m3 (x3)

K(3,2) = -k2;

K(3,3) = k2 + k3;

K(3,4) = -k3;

% For m4 (x4)

K(4,3) = -k3;

K(4,4) = k3;

% For m5 (x5)

K(5,1) = -k4;

K(5,5) = k4 + k5 + k6;

K(5,6) = -k5;

K(5,7) = -k6;

% For m6 (x6)

K(6,5) = -k5;

K(6,6) = k5 + k11;

K(6,11) = -k11;

% For m7 (x7)

K(7,1) = -k7;

K(7,5) = -k6;

K(7,7) = k6 + k7 + k23;

K(7,18) = -k23;

% For m8 (x8)

K(8,1) = -k8;

K(8,8) = k8 + k9;

K(8,9) = -k9;

% For m9 (x9)

K(9,8) = -k9;

K(9,9) = k9 + k10;

K(9,10) = -k10;

% For m10 (x10)

K(10,9) = -k10;

K(10,10) = k10;

% For m11 (x11)

K(11,6) = -k11;

K(11,7) = -k12;

K(11,11) = k11 + k12 + k13 + k14;

K(11,12) = -k13;

K(11,13) = -k14;

% For m12 (x12)

K(12,11) = -k13;

K(12,12) = k13 + k15;

K(12,14) = -k15;

% For m13 (x13)

K(13,11) = -k14;

K(13,13) = k14 + k16;

K(13,15) = -k16;

% For m14 (x14)

K(14,12) = -k15;

K(14,14) = k15 + k17;

K(14,16) = -k17;

% For m15 (x15)

K(15,13) = -k16;

K(15,15) = k16 + k18;

K(15,17) = -k18;

% For m16 (x16)

K(16,14) = -k17;

K(16,16) = k17 + k19;

% For m17 (x17)

K(17,15) = -k18;

K(17,17) = k18 + k20;

% For m18 (x18)

K(18,7) = -k23;

K(18,18) = k23 + k21 + k22;

K(18,19) = -k21 - k22;

% For m19 (x19)

K(19,18) = -k21 - k22;

K(19,19) = k21 + k22;

% Damping matrix initialisation 19x19

C = zeros(19,19);

% Definition of damping matrix elements C

% For m1 (x1)

C(1,1) = c1 + c4 + c7 + c8;

C(1,2) = -c1;

C(1,5) = -c4;

C(1,7) = -c7;

C(1,8) = -c8;

% For m2 (x2)

C(2,1) = -c1;

C(2,2) = c1 + c2;

C(2,3) = -c2;

% For m3 (x3)

C(3,2) = -c2;

C(3,3) = c2 + c3;

C(3,4) = -c3;

% For m4 (x4)

C(4,3) = -c3;

C(4,4) = c3;

% For m5 (x5)

C(5,1) = -c4;

C(5,5) = c4 + c5 + c6;

C(5,6) = -c5;

C(5,7) = -c6;

% For m6 (x6)

C(6,5) = -c5;

C(6,6) = c5 + c11;

C(6,11) = -c11;

% For m7 (x7)

C(7,1) = -c7;

C(7,5) = -c6;

C(7,7) = c6 + c7 + c23;

C(7,18) = -c23;

% For m8 (x8)

C(8,1) = -c8;

C(8,8) = c8 + c9;

C(8,9) = -c9;

% For m9 (x9)

C(9,8) = -c9;

C(9,9) = c9 + c10;

C(9,10) = -c10;

% For m10 (x10)

C(10,9) = -c10;

C(10,10) = c10;

% For m11 (x11)

C(11,6) = -c11;

C(11,7) = -c12;

C(11,11) = c11 + c12 + c13 + c14;

C(11,12) = -c13;

C(11,13) = -c14;

% For m12 (x12)

C(12,11) = -c13;

C(12,12) = c13 + c15;

C(12,14) = -c15;

% For m13 (x13)

C(13,11) = -c14;

C(13,13) = c14 + c16;

C(13,15) = -c16;

% Per m14 (x14)

C(14,12) = -c15;

C(14,14) = c15 + c17;

C(14,16) = -c17;

% For m15 (x15)

C(15,13) = -c16;

C(15,15) = c16 + c18;

C(15,17) = -c18;

% For m16 (x16)

C(16,14) = -c17;

C(16,16) = c17 + c19;

% For m17 (x17)

C(17,15) = -c18;

C(17,17) = c18 + c20;

% For m18 (x18)

C(18,7) = -c23;

C(18,18) = c23 + c21 + c22;

C(18,19) = -c21 - c22;

% For m19 (x19)

C(19,18) = -c21 - c22;

C(19,19) = c21 + c22;

Answer 2

[Shivam Gothi]

Hello @Roberta Pili

The ode45 function is a variable step differential equations solver. In your current code:

[t, y] = ode45(@ (t, y) odefcn_standing(t, y, M, C, K), tspan, y0);

The variable “t” is treated as a scalar each time the function is executed. This scalar value is then passed as an argument to the “chirp” function. According to the documentation below:

https://www.mathworks.com/help/signal/ref/chirp.html

the “t” argument in “chirp” should be a vector of time steps, as the function generates samples of a linearly swept-frequency cosine signal at those specific time instances. As it stands, a new chirp signal is inadvertently created every time odefcn_standing is called, which is not the intended behavior. Consequently, a proper chirp signal is not being passed to the function you aim to solve.

To address this issue, I would suggest a possible workaround:

1. Define the chirp signal, named “My_chirp”, prior to invoking “ode45”.
2. Instead of generating a new chirp signal within odefcn_standing, pass the pre-defined My_chirp signal as an input argument to the function as shown below

n=19;
y0 = zeros(2*n,1);
tspan = [0 120];
% ode45
t_chirp = linspace(tspan(1),tspan(2),1000); %define the time steps for chirp  signal generation
f0 = 0.5; % initial frequency
f1 = 80; % final frequency
t_f = tspan(2);
My_chirp = chirp(t_chirp, f0, t_f, f1); %Create chirp signal
[t, y] = ode45(@(t, y) odefcn_standing(t, y, M, C, K, My_chirp, t_chirp), tspan, y0);

Now, change the function definition of “odefcn_standing” as shown below:

function dy = odefcn_standing(t, y, M, C, K, My_chirp, time_vector)
    t
    n = 19; % Numero di gradi di libertà
    dy = zeros(2 * n, 1);
    
    % Construction of matrix A
    A = [zeros(n), eye(n);
        -inv(M) * K, -inv(M) * C];
    
    F = zeros(19, 1);

    chirp_signal = interp1(time_vector,My_chirp,t) %get the value of chirp_signal at time "t" by using interpolation

    F(16,:) = 10*chirp_signal; % on mass 16
    F(17,:) = 10*chirp_signal; % on mass 17
    
    % Construction of matrix B
    B = [zeros(n, n); inv(M)];
    dy = A * y + B * F;
end

Note: Here, we are getting the value of chirp signal by interpolating between the existing samples of chirp signal. Refer to the below documentation for help regarding “interp1” function:

https://www.mathworks.com/help/matlab/ref/double.interp1.html

Additionally, you may find the below links to be useful:

1. https://www.mathworks.com/help/matlab/ref/ode45.html
2. https://www.mathworks.com/help/matlab/ref/double.linspace.html

I have also attached the modified (.M) file with this answer

I hope this helps !

[Roberta Pili]

Thank you

[Roberta Pili]

this is the result found. But shouldn't the peaks in the graph coincide with the resonance frequencies I found with this code?

[U,L]=eig(K,M);

[L_sorted, ind] = sort(diag(L),'ascend');

L_sorted = transpose(L_sorted);

U_sorted = U(:,ind);

fn=real((sqrt(L_sorted))/(2*pi)); %natural freq. in Hz

wn=(2*pi).*fn;

Answer 3

[Shivam Gothi]

Hello @Roberta Pili,

Thankyou for sharing the data.

As per the data shared by you, I modified the (.M) file and run the simulation for 12s. (by modifing the tspan vector). Below given is the code of modified (.M) file and the output generated by it.

%MATRIX


%%% Input data %%%%%%%%%%%%%%%%%%%
%Mass 1 
m1=4.5;
k1=865.47;   %ho diviso per 4
k4=865.47;
k7=865.47;
k8=865.47;
c1=12.8;   %ho diviso per 4
c4=12.8;
c7=12.8;
c8=12.8;
%Mass 2 and 8
m2=2.070;
k2=5107.515;
c2=115.16;
m8=2.070;
k9=5107.515;
c9=115.16;
%Mass 3 and 9
m3=1.160;
k3=26468.485;
c3=143.68;
m9=1.160;
k10=26468.485;
c10=143.68;
%Mass 4 and 10
m4=0.540;
m10=0.540;
%Mass 5
m5=15.551;
k5=263068.52;  %ho diviso per 2
k6=263068.52;
c5=715.1;
c6=715.1;
%Mass 6
m6=12.5;
k11=19739.205;
c11=291.878;
%Mass 7
m7=4.5;
k12=30023.33;
c12=350;
%Mass 11
m11=27.174;
k13=21455.7275;  %ho diviso per 2
k14=21455.7275;
c13=280.76;
c14=280.76;
%Mass 12 and 13
m12=7.165;
k15=99837807.76;
c15=2376.4;
m13=7.165;
k16=99837807.76;
c16=2376.4;
%Mass 14 and 15
m14=2.8;
k17=1727180.77;
c17=2376.4;
m15=2.8;
k18=1727180.77;
c18=2376.4;
%Mass 16 and 17
m16=1.170;
k19=4127.21;
c19=2.084;
m17=1.170;
k20=4127.21;
c20=2.084;
%Mass 18
m18=3.8;
k23=457384621.1;
c23=13349.82;
%Mass 19
m19=1.336;
k21=5242.93;  %ho diviso per 2
k22=5242.93;
c21=23.67;
c22=23.67;
%%
%MATRCS
M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19]);
% Stiffness matrix initialisation 19x19
K = zeros(19,19);
% Definition of stiffness matrix elements K
% For m1 (x1)
K(1,1) = k1 + k4 + k7 + k8;
K(1,2) = -k1;
K(1,5) = -k4;
K(1,7) = -k7;
K(1,8) = -k8;
% For m2 (x2)
K(2,1) = -k1;
K(2,2) = k1 + k2;
K(2,3) = -k2;
% For m3 (x3)
K(3,2) = -k2;
K(3,3) = k2 + k3;
K(3,4) = -k3;
% For m4 (x4)
K(4,3) = -k3;
K(4,4) = k3;
% For m5 (x5)
K(5,1) = -k4;
K(5,5) = k4 + k5 + k6;
K(5,6) = -k5;
K(5,7) = -k6;
% For m6 (x6)
K(6,5) = -k5;
K(6,6) = k5 + k11;
K(6,11) = -k11;
% For m7 (x7)
K(7,1) = -k7;
K(7,5) = -k6;
K(7,7) = k6 + k7 + k23;
K(7,18) = -k23;
% For m8 (x8)
K(8,1) = -k8;
K(8,8) = k8 + k9;
K(8,9) = -k9;
% For m9 (x9)
K(9,8) = -k9;
K(9,9) = k9 + k10;
K(9,10) = -k10;
% For m10 (x10)
K(10,9) = -k10;
K(10,10) = k10;
% For m11 (x11)
K(11,6) = -k11;
K(11,7) = -k12;
K(11,11) = k11 + k12 + k13 + k14;
K(11,12) = -k13;
K(11,13) = -k14;
% For m12 (x12)
K(12,11) = -k13;
K(12,12) = k13 + k15;
K(12,14) = -k15;
% For m13 (x13)
K(13,11) = -k14;
K(13,13) = k14 + k16;
K(13,15) = -k16;
% For m14 (x14)
K(14,12) = -k15;
K(14,14) = k15 + k17;
K(14,16) = -k17;
% For m15 (x15)
K(15,13) = -k16;
K(15,15) = k16 + k18;
K(15,17) = -k18;
% For m16 (x16)
K(16,14) = -k17;
K(16,16) = k17 + k19;
% For m17 (x17)
K(17,15) = -k18;
K(17,17) = k18 + k20;
% For m18 (x18)
K(18,7) = -k23;
K(18,18) = k23 + k21 + k22;
K(18,19) = -k21 - k22;
% For m19 (x19)
K(19,18) = -k21 - k22;
K(19,19) = k21 + k22;
% Damping matrix initialisation 19x19
C = zeros(19,19);
% Definition of damping matrix elements C
% For m1 (x1)
C(1,1) = c1 + c4 + c7 + c8;
C(1,2) = -c1;
C(1,5) = -c4;
C(1,7) = -c7;
C(1,8) = -c8;
% For m2 (x2)
C(2,1) = -c1;
C(2,2) = c1 + c2;
C(2,3) = -c2;
% For m3 (x3)
C(3,2) = -c2;
C(3,3) = c2 + c3;
C(3,4) = -c3;
% For m4 (x4)
C(4,3) = -c3;
C(4,4) = c3;
% For m5 (x5)
C(5,1) = -c4;
C(5,5) = c4 + c5 + c6;
C(5,6) = -c5;
C(5,7) = -c6;
% For m6 (x6)
C(6,5) = -c5;
C(6,6) = c5 + c11;
C(6,11) = -c11;
% For m7 (x7)
C(7,1) = -c7;
C(7,5) = -c6;
C(7,7) = c6 + c7 + c23;
C(7,18) = -c23;
% For m8 (x8)
C(8,1) = -c8;
C(8,8) = c8 + c9;
C(8,9) = -c9;
% For m9 (x9)
C(9,8) = -c9;
C(9,9) = c9 + c10;
C(9,10) = -c10;
% For m10 (x10)
C(10,9) = -c10;
C(10,10) = c10;
% For m11 (x11)
C(11,6) = -c11;
C(11,7) = -c12;
C(11,11) = c11 + c12 + c13 + c14;
C(11,12) = -c13;
C(11,13) = -c14;
% For m12 (x12)
C(12,11) = -c13;
C(12,12) = c13 + c15;
C(12,14) = -c15;
% For m13 (x13)
C(13,11) = -c14;
C(13,13) = c14 + c16;
C(13,15) = -c16;
% Per m14 (x14)
C(14,12) = -c15;
C(14,14) = c15 + c17;
C(14,16) = -c17;
% For m15 (x15)
C(15,13) = -c16;
C(15,15) = c16 + c18;
C(15,17) = -c18;
% For m16 (x16)
C(16,14) = -c17;
C(16,16) = c17 + c19;
% For m17 (x17)
C(17,15) = -c18;
C(17,17) = c18 + c20;
% For m18 (x18)
C(18,7) = -c23;
C(18,18) = c23 + c21 + c22;
C(18,19) = -c21 - c22;
% For m19 (x19)
C(19,18) = -c21 - c22;
C(19,19) = c21 + c22;

%%%%%%% code starts here %%%%%%%%%%%%%%%%%%%%%%

M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19]);
% stiffness matrix 19x19
K = zeros(19,19);
K(1,1) = k1 + k4 + k7 + k8;
K(1,2) = -k1;
K(1,5) = -k4;
K(1,7) = -k7;
K(1,8) = -k8;
K(2,1) = -k1;
K(2,2) = k1 + k2;
K(2,3) = -k2;
K(3,2) = -k2;
K(3,3) = k2 + k3;
K(3,4) = -k3;
K(4,3) = -k3;
K(4,4) = k3;
K(5,1) = -k4;
K(5,5) = k4 + k5 + k6;
K(5,6) = -k5;
K(5,7) = -k6;
K(6,5) = -k5;
K(6,6) = k5 + k11;
K(6,11) = -k11;
K(7,1) = -k7;
K(7,5) = -k6;
K(7,7) = k6 + k7 + k23;
K(7,18) = -k23;
K(8,1) = -k8;
K(8,8) = k8 + k9;
K(8,9) = -k9;
K(9,8) = -k9;
K(9,9) = k9 + k10;
K(9,10) = -k10;
K(10,9) = -k10;
K(10,10) = k10;
K(11,6) = -k11;
K(11,7) = -k12;
K(11,11) = k11 + k12 + k13 + k14;
K(11,12) = -k13;
K(11,13) = -k14;
K(12,11) = -k13;
K(12,12) = k13 + k15;
K(12,14) = -k15;
K(13,11) = -k14;
K(13,13) = k14 + k16;
K(13,15) = -k16;
K(14,12) = -k15;
K(14,14) = k15 + k17;
K(14,16) = -k17;
K(15,13) = -k16;
K(15,15) = k16 + k18;
K(15,17) = -k18;
K(16,14) = -k17;
K(16,16) = k17 + k19;
K(17,15) = -k18;
K(17,17) = k18 + k20;
K(18,7) = -k23;
K(18,18) = k23 + k21 + k22;
K(18,19) = -k21 - k22;
K(19,18) = -k21 - k22;
K(19,19) = k21 + k22;
% damping matrix 19x19
C = zeros(19,19);
C(1,1) = c1 + c4 + c7 + c8;
C(1,2) = -c1;
C(1,5) = -c4;
C(1,7) = -c7;
C(1,8) = -c8;
C(2,1) = -c1;
C(2,2) = c1 + c2;
C(2,3) = -c2;
C(3,2) = -c2;
C(3,3) = c2 + c3;
C(3,4) = -c3;
C(4,3) = -c3;
C(4,4) = c3;
C(5,1) = -c4;
C(5,5) = c4 + c5 + c6;
C(5,6) = -c5;
C(5,7) = -c6;
C(6,5) = -c5;
C(6,6) = c5 + c11;
C(6,11) = -c11;
C(7,1) = -c7;
C(7,5) = -c6;
C(7,7) = c6 + c7 + c23;
C(7,18) = -c23;
C(8,1) = -c8;
C(8,8) = c8 + c9;
C(8,9) = -c9;
C(9,8) = -c9;
C(9,9) = c9 + c10;
C(9,10) = -c10;
C(10,9) = -c10;
C(10,10) = c10;
C(11,6) = -c11;
C(11,7) = -c12;
C(11,11) = c11 + c12 + c13 + c14;
C(11,12) = -c13;
C(11,13) = -c14;
C(12,11) = -c13;
C(12,12) = c13 + c15;
C(12,14) = -c15;
C(13,11) = -c14;
C(13,13) = c14 + c16;
C(13,15) = -c16;
C(14,12) = -c15;
C(14,14) = c15 + c17;
C(14,16) = -c17;
C(15,13) = -c16;
C(15,15) = c16 + c18;
C(15,17) = -c18;
C(16,14) = -c17;
C(16,16) = c17 + c19;
C(17,15) = -c18;
C(17,17) = c18 + c20;
C(18,7) = -c23;
C(18,18) = c23 + c21 + c22;
C(18,19) = -c21 - c22;
C(19,18) = -c21 - c22;
C(19,19) = c21 + c22;
n=19;
y0 = zeros(2*n,1);
tspan = [0 12];
% ode45
t_chirp = linspace(tspan(1),tspan(2),1000); %define the time steps for chirp  signal generation
f0 = 0.5; % initial frequency
f1 = 80; % final frequency
t_f = tspan(2);
My_chirp = chirp(t_chirp, f0, t_f, f1); %Create chirp signal
[t, y] = ode45(@(t, y) odefcn_standing(t, y, M, C, K, My_chirp, t_chirp), tspan, y0);
figure;
plot(t, y(:, 19));   
xlabel('Time (s)');
ylabel('Displacement (m)');
% legend('y1', 'y2', 'y3');
title('response of the system 19DOF');
grid on;

function dy = odefcn_standing(t, y, M, C, K, My_chirp, time_vector)
    n = 19; % Numero di gradi di libertà
    dy = zeros(2 * n, 1);
    
    % Construction of matrix A
    A = [zeros(n), eye(n);
        -inv(M) * K, -inv(M) * C];
    
    F = zeros(19, 1);

    chirp_signal = interp1(time_vector,My_chirp,t); %get the value of chirp_signal at time "t" by using interpolation

    F(16,:) = 10*chirp_signal; % on mass 16
    F(17,:) = 10*chirp_signal; % on mass 17
    
    % Construction of matrix B
    B = [zeros(n, n); inv(M)];
    dy = A * y + B * F;
end

Kindly verify whether the output is as per your expectation or its still not correct.

[Roberta Pili]

This is the result found. But shouldn't the peaks in the graph coincide with the resonance frequencies I found with this code?

[U,L]=eig(K,M);

[L_sorted, ind] = sort(diag(L),'ascend');

L_sorted = transpose(L_sorted);

U_sorted = U(:,ind);

fn=real((sqrt(L_sorted))/(2*pi)); %natural freq. in Hz

wn=(2*pi).*fn;

[Shivam Gothi]

Hello @Roberta Pili,

I did not get your question. The result shared by you is the plot of magnetude vs time. It is a time domain response. How can the peak in time domain response coincide with resonant frequencies (term used in frequency domain response) ?

[Roberta Pili]

I have superimposed the chirp rate on the graph so that we can understand the peak at what frequency it occurs. Technically I should have an increase in the peak when the system passes through a resonance frequency.

[Shivam Gothi]

Hello @Roberta Pili,

I understood your question. As per my understanding, the chirp ssignal is generated by command :

t_chirp = linspace(tspan(1),tspan(2),1000); %define the time steps for chirp  signal generation
f0 = 0.5; % initial frequency
f1 = 80; % final frequency
t_f = tspan(2);
My_chirp = chirp(t_chirp, f0, t_f, f1); %Create chirp signal

Here, there are only 1000 timesteps between 0 and 120 seconds. In order to achieve smooth chirp signal, there should be larger number of time steps between 0 to 120 seconds. I tried using 100000 instead of 1000, but the simulation is taking a lot of time to execute.

Therefore, the best way to find resonant peak is, run the simulation from 0 to 10 seconds with 10000 timesteps between them taking fo = 0.5 and f1 = 5. After that, again run the simulation, but this time, make fo = 5 and f1 = 10. keep on repeating this at interval of 5 Hz. You should observe the peak at perticular time step.

[Roberta Pili]

Thank you. I have another question. I have tried to replicate the code for the same system, but reducing the degrees of freedom. However, the graph does not come out as I would expect. What could be the problem?

clear

close all

clc

%%

%Massa 1

m1=4.5;

k1=865.47;

k4=865.47;

k7=865.47;

k8=865.47;

c1=12.8;

c4=12.8;

c7=12.8;

c8=12.8;

%Massa 2 e 8

m2=2.070;

k2=5107.515;

c2=115.16;

m8=2.070;

k9=5107.515;

c9=115.16;

%Massa 3 e 9

m3=1.160;

k3=26468.485;

c3=143.68;

m9=1.160;

k10=26468.485;

c10=143.68;

%Massa 4 e 10

m4=0.540;

m10=0.540;

%Massa 5

m5=15.551;

k5=263068.52;

k6=263068.52;

c5=715.1;

c6=715.1;

%Massa 6

m6=12.5;

k11=19739.205;

c11=291.878;

%Massa 7

m7=4.5;

k12=30023.33;

c12=350;

%Massa 11

m11=27.174;

% k16=42911.455;

k16=21455.7275;

c16=561.52;

%Massa 12

m12=3.8;

k15=457384621.1;

c15=13349.82;

%Massa 13

m13=1.336;

k13=5242.93;

k14=5242.93;

c13=23.67;

c14=23.67;

%%

%MATRICS

M=diag([m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13]);

K = zeros(13,13);

K(1,1) = k1 + k4 + k7 + k8;

K(1,2) = -k1;

K(1,5) = -k4;

K(1,7) = -k7;

K(1,8) = -k8;

K(2,1) = -k1;

K(2,2) = k1 + k2;

K(2,3) = -k2;

K(3,2) = -k2;

K(3,3) = k2 + k3;

K(3,4) = -k3;

K(4,3) = -k3;

K(4,4) = k3;

K(5,1) = -k4;

K(5,5) = k4 + k5 + k6;

K(5,6) = -k5;

K(5,7) = -k6;

K(6,5) = -k5;

K(6,6) = k5 + k11;

K(6,11) = -k11;

K(7,1) = -k7;

K(7,5) = -k6;

K(7,7) = k6 + k7 + k15;

K(7,12) = -k15;

K(8,1) = -k8;

K(8,8) = k8 + k9;

K(8,9) = -k9;

K(9,8) = -k9;

K(9,9) = k9 + k10;

K(9,10) = -k10;

K(10,9) = -k10;

K(10,10) = k10;

K(11,6) = -k11;

K(11,11) = k11 + k12 + k13 + k16;

K(11,7) = -k12;

K(12,7) = -k15;

K(12,12) = k15 + k13 + k14;

K(12,13) = -k13 - k14;

K(13,12) = -k13 - k14;

K(13,13) = k13 + k14;

C = zeros(13,13);

C(1,1) = c1 + c4 + c7 + c8;

C(1,2) = -c1;

C(1,5) = -c4;

C(1,7) = -c7;

C(1,8) = -c8;

C(2,1) = -c1;

C(2,2) = c1 + c2;

C(2,3) = -c2;

C(3,2) = -c2;

C(3,3) = c2 + c3;

C(3,4) = -c3;

C(4,3) = -c3;

C(4,4) = c3;

C(5,1) = -c4;

C(5,5) = c4 + c5 + c6;

C(5,6) = -c5;

C(5,7) = -c6;

C(6,5) = -c5;

C(6,6) = c5 + c11;

C(6,11) = -c11;

C(7,1) = -c7;

C(7,5) = -c6;

C(7,7) = c6 + c7 + c15;

C(7,12) = -c15;

C(8,1) = -c8;

C(8,8) = c8 + c9;

C(8,9) = -c9;

C(9,8) = -c9;

C(9,9) = c9 + c10;

C(9,10) = -c10;

C(10,9) = -c10;

C(10,10) = c10;

C(11,6) = -c11;

C(11,11) = c11 + c12 + c13 + c16;

C(11,7) = -c12;

C(12,7) = -c15;

C(12,12) = c15 + c13 + c14;

C(12,13) = -c13 - c14;

C(13,12) = -c13 - c14;

C(13,13) = c13 + c14;

%% CHIRP

f0 = 0.5;

f1 = 80;

t_f = 200;

tspan = [0 200];

t = linspace(tspan(1), tspan(2), 1000);

chirp_signal = chirp(t, f0, t_f, f1);

F = zeros(13, numel(t));

F(11, :) = chirp_signal;

force_interp = @(time) interp1(t, F', time)'; % Interpola F rispetto a t

%% ODE45

y0 = zeros(2*13,1);

odefun = @(time, y) [y(14:end); M \ (force_interp(time) - C*y(14:end) - K*y(1:13))];

[t_out, y_out] = ode45(odefun, tspan, y0);

%% PLOT

figure

plot(t_out, y_out(:,13));

xlabel('Tempo [s]');

ylabel('Spostamento [m]');

title('Spostamento della massa 13 nel tempo');

grid on;

[Shivam Gothi]

Hello @Roberta Pili,

consider this line of code in your comment

tspan=[0 200];
t = linspace(tspan(1), tspan(2), 1000); 

There are only 1000 time steps. Instead of taking 1000 steps, increase the number to 100000. The reason is, with 1000 steps, the "chirp" signal is plotted below:

%% CHIRP
f0 = 0.5;   
f1 = 80;   
t_f = 200; 
tspan = [0 200];
t = linspace(tspan(1), tspan(2), 1000);  
chirp_signal = chirp(t, f0, t_f, f1);
plot(t,chirp_signal)

We can clearly see that it is not a proper chirp signal.

After making 100000 steps, we get the chirp signal as:

t = linspace(tspan(1), tspan(2), 100000);  
chirp_signal = chirp(t, f0, t_f, f1);
plot(t,chirp_signal)

After zooming it, we can see that it is a proper chirp signal.

plot(t([1:3000]),chirp_signal([1:3000]))

plot(t([3000:6000]),chirp_signal([3000:6000]))

CAUTION:

increasing the time steps will slow down the simulation to a great extent.

[Roberta Pili]

inished running with 100000 but came back giving me a graph that does not represent a shift due to the chirp. I only changed the number from 1000 to 100000.

it is as if it is not applying chirp to my system.