多自由度彈簧振子系統的ode解法

04-01

振動系統中各個質量塊子系統的狀態方程的建立參照此鏈接。如下系統

它的狀態方程是

作如下符號變換

可以將多階微分方程轉化為N多個一階微分方程組，那麼就可以用MATLAB或者Octave數值求解。下面給出的代碼可以根據輸入的參數給出數值解，激勵可以根據實際需求做出相應修改。

1. MATLAB 計算代碼

function [t,Y] = MDSOS(M,R,K,A,f,Yt0,Tspan)% Multi-degree spring-mass oscillator system% M - 質量矩陣 n*n% R - 阻尼矩陣 n*n% K - 剛度矩陣 n*n% A - 驅動幅值 n*1% f - 驅動頻率 n*1% Yt0 - 初始條件位移、速度% Tspan - 求解區間 1*2options = odeset(RelTol,1e-19); [t,Y] = ode45(@(t,Y)dynamics(t,Y,M,R,K,A,f),Tspan,Yt0,options);warning offendfunction dY = dynamics(t,Y,M,R,K,A,f)[s,~] = size(M);F = A.*sin(2*pi*f*t);% ---------------------------------------------dY = zeros(2*s,1); % 申請變數空間% --------------- 微分方程組 -------------------dY(1:1:s,1) = Y((1:1:s)+s,1);dY((1:1:s)+s,1) = (M^-1)*F - (M^-1)*R*Y((1:1:s)+s,1) - (M^-1)*K*Y(1:1:s,1);end

2. 演示代碼：雙自由度

function demonMDSOS()M = [100 0; 0 20];R = [10 0; 0 0];K = [4 -2; -2 200]; A = [10;0];f = [sqrt(10);0];Tspan = [0,50]; % 解的時間區間Yt0 = [0;0;10;10]; % 初始條件[t,y] = MDSOS(M,R,K,A,f,Yt0,Tspan);plot(t,y(:,1),r,t,y(:,2),b)end

3. 演示代碼：三自由度

function demonMDSOS()M = [100 0 0; 0 20 0; 0 0 50];R = [10 0 0; 0 0 0; 0 0 0];K = [4 -2 0; -2 200 -10; 0 -10 150]; A = [10;0;0];f = [sqrt(10);0;0]; Tspan = [0,50]; Yt0 = [0;0;0;10;10;10]; [t,y] = MDSOS(M,R,K,A,f,Yt0,Tspan);plot(t,y(:,1),r,t,y(:,2),b,t,y(:,3),m)end