对单自由度系统的强制振动进行数值方法研究的实现（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

单自由度系统的强制振动数值方法研究及实现

摘要
本文研究了受谐波力激励的阻尼单自由度（SDOF）系统的精确解，并对其与Matlab内置函数ode45、中心差分法、Newmark法和四阶Runge-Kutta法提供的数值解进行了比较。这些数值解的实现基于S. Rao的著作。为了实现这些数值方法，我们编写了RK4.m函数（使用四阶Runge-Kutta法数值求解阻尼系统的运动方程）、Newmark.m函数（使用Newmark法数值求解阻尼系统的运动方程）、CentDiff.m函数（使用中心差分法数值求解阻尼系统的运动方程）以及Matlab LiveScript Documentation.mlx（用于文档说明）。

内容
本文对受谐波力激励的阻尼单自由度系统进行了深入研究，旨在探索数值方法对系统振动的准确性和稳定性。我们首先介绍了系统的精确解，并将其作为基准与其他数值解进行比较。通过编写并实现RK4.m、Newmark.m和CentDiff.m等函数，我们成功地使用四阶Runge-Kutta法、Newmark法和中心差分法对系统的运动方程进行了数值求解。这些数值方法的实现基于文献【2】.

内容

RK4.m函数，使用四阶Runge-Kutta法数值求解阻尼系统的运动方程
Newmark.m函数，使用Newmark法数值求解阻尼系统的运动方程
CentDiff.m函数，使用中心差分法数值求解阻尼系统的运动方程
Matlab LiveScript Documentation.mlx，用于文档说明
文献[1]Engineering Vibrations - 百度学术 (baidu.com)

文献[2]Mechanical Vibrations - 百度学术 (baidu.com) 

📚2 运行结果

部分代码：

function [y] = CentDiff(F,M,K,C,dt,x0,v0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GOAl: Numerical solve the equation of motion of a damped system

% INPUT
% F : vector  -- size: [1x N] -- Time series representinf the time history of the load. 
% M : scalar  -- size: [1 x 1] -- Modal mass
% K : scalar  -- size: [1 x 1] -- Modal stifness
% C : scalar  -- size: [1 x 1] -- Modal damping
% dt : scalar  -- size: [1 x 1] -- time step
% x0 : scalar  -- size: [1 x 1] -- initial displacement
% v0 : scalar  -- size: [1 x 1] -- initial velocity
% 
% OUTPUT
% y: time history of the system response to the load

% author: Etienne Cheyet  ----  last updated: 07/11/2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Initialisation
N = size(F,2);
% preallocation
y = zeros(size(F));

% initial acceleration
a0 = M\(F(1)-C.*v0-K.*x0);
% initialisation of y (first 2 values).
y0 = x0-dt.*v0+dt^2/2*a0;
y(:,1) = x0;
A = (M./dt.^2+C./(2*dt));
B = ((2*M./dt.^2-K).*y(:,1)+(C./(2*dt)-M./dt.^2).*y0+F(:,1));
y(:,2) = A\B;

% For the rest of integration points
for ii=2:N-1,
    A = (M./dt.^2+C./(2*dt));
    B = ((2*M./dt.^2-K).*y(:,ii)+(C./(2*dt)-M./dt.^2).*y(:,ii-1)+F(:,ii)); 
    y(:,ii+1) = A\B;
end

end

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

[1] Inman D J .Engineering Vibrations[J].McGraw-Hill, 2001.Engineering Vibrations - 百度学术 (baidu.com)

[2] Singiresu R .Mechanical Vibrations[J].  1990.Mechanical Vibrations - 百度学术 (baidu.com)

🌈4 Matlab代码实现