使用中心差分法求解振动系统的响应——极易推广到多自由度
% equation m*ddx^2+c*dx+k*x=F
%numerial method have some fundametal characterristic :
%frist:they ara not intended to satisfy the governing differential equation
%all times ,but only at discrete time intervals ht apart
%second:...
% finite difference method
%the main idea in the finite difference method is to use approximations to
%derivatives.have three types of formulas_forward,backward,and central
%difference formulas.
%by using Taylors series we can obtain
%dx(i)= {x(i+1)- x(i-1) }/2h (2)
%ddx(i)={x(i+1)-2*x(i)+x(i-1)}/h^2 (3)
%let the duration over which the solution of equation(1),is requried be
%divided into n equal parts of interal h each .howevere wo must note
%that,select a time step h that is smaller than a critial time step
%Hlim.where Hlim=Tn/pi; Tn is the natural period of the system or the small
%such period. or,will cause ununstable,and the the responese worthless
%formulate
%a=1/(M/h^2+C/(2*h) ); b=2*M/h^2-K; c=C/(2*h)-M/h^2;
% xnew=a*(b*xold+c*xmold+Fold) :
%where xmold means be more old than xold or say x(i-1);xnew=x(i+1);
%Fold=F(i);
%from this foumalate we know it permits us to calculate the displacement of
%the mass x(i+1),if we konw the preious history of the displacement ar
%t(i)and t(i-1);as well as the present external force F(i).
%howere we find that wo must x(0) and x(-1) in oder to find the x(1);
% x(-1)=x0-h*dx0+h^2*ddx0/2;
% ddx0=(F(0)-C*dx0-K*x0)/M
%find the response of single degree of freedom system
function numerical_method()
clc
clear all
M=10;K=4000;C=20;Tmax=20;x0=0.04;dx0=0.06;h=0.01;
Hlim=sqrt(K/M)/pi; %the maxi step h<Hlim
if nargin<9;h=0.001;end
if h>=0.01*Hlim; error('time step larger than Hlim ');end
ddx0=(f(0)-C*dx0-K*x0)/M; xmold=x0-h*dx0+h^2*ddx0/2; xold=x0;
a=1/(M/h^2+C/(2*h) ); b=2*M/h^2-K; c=C/(2*h)-M/h^2;n=1;
for T=0:h:Tmax
Fold=f(T);
xnew=a*(b*xold+c*xmold+Fold);
xsave(1,n)=xnew; m(1,n)=T;
xmold=xold;
xold=xnew;
n=n+1;
end
plot(m,xsave,'r');grid on;title('response of single degree of freedom system')
end
function f=f(t)
%should give differnet function
f=100*sin(10*t);
end
原创 or 转贴 !? 个人专业有限, 看似像分享!? 是吗?
若是, 若能稍作简易说明下相关输入/输出/算法, 甚至举例, 或许更适合学习! 个人浅见
页:
[1]