resolution of MDOF using ode45

Question

Roberta Pili 2024-9-3，15:03

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149804-resolution-of-mdof-using-ode45

评论： Roberta Pili 2024-9-18，5:48

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

2 个评论
显示无隐藏无

Shivam Gothi 2024-9-16，4:15

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 2024-9-16，7:02

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

Maven Euler 2024-9-16，4:35

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149804-resolution-of-mdof-using-ode45#answer_1516900

Hi Roberta

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

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.
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.
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.
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
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.
Additionally, you can use the odeset function to specify options like relative and absolute tolerances to improve the solver's accuracy.
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.
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!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Roberta Pili 2024-9-16，7:04

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 2024-9-16，6:08

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149804-resolution-of-mdof-using-ode45#answer_1516935

在 MATLAB Online 中打开

chirp_signal.m

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:

Define the chirp signal, named “My_chirp”, prior to invoking “ode45”.
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:

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

I hope this helps !

2 个评论
显示无隐藏无

Roberta Pili 2024-9-16，7:03

Thank you

Roberta Pili 2024-9-16，8:15

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 2024-9-16，7:39

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149804-resolution-of-mdof-using-ode45#answer_1517000

在 MATLAB Online 中打开

chirp_signal.m

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Roberta Pili 2024-9-17，9:43

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 2024-9-17，10:21

在 MATLAB Online 中打开

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 2024-9-18，5:48

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.