中心差分法Matlab源程序

中心差分法Matlab源程序

clear
M=[2 0;0 1 ]       
C=[0 0;0 0]
K=[6 -2;-2 4]
Q1=[0;10]       %外力
d1=[0;0]        %初始位移
v1=[0;0]        %初始速度
a1=inv（M）*（Q1-C*v1-K*d1） %求得初始加速度
delta=0.28              %选取时间步长（可根据实际选取）
c0=1./（delta^2）         %积分常数
c1=1./（2.*delta）
c2=2.*c0
c3=1./c2
d0=d1-delta*v1+c3*a1  %起步位移计算
M1=c0*M+c1*C
i=1
for ti=0:delta:5      %选定时间范围（可任意）
    d（i1:2）=d1       %存储各个时刻ti位移d（t）
    v（i1:2）=v1       %存储各个时刻ti速度v（t）
    a（i1:2）=a1       %存储各个时刻ti加速度a（t）
    t（i）=ti           %存储各个时刻ti
    Q2=Q1-（K-c2*M）*d1-（c0*M-c1*C）*d0 %计算t时刻有效载荷
    d2=inv（M1）*Q2                   %求t+delta时刻位移
    v1=c1*（-d0+d2）                 %计算t时刻速度
    a1=c0*（d0-2*d1+d2）             %计算t时刻加速度
    d0=d1
    d1=d2
    i=i+1   
end
%将位移、速度、加速度各两条曲线画在同一面内
subplot（231）plot（td（:1）‘r‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘d1（t）‘）
subplot（234）plot（td（:2）‘r‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘d2（t）‘）

subplot（232）plot（tv（:1）‘g‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘v1（t）‘）
subplot（235）plot（tv（:2）‘g‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘v2（t）‘）

subplot（233）plot（ta（:1）‘b‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘a1（t）‘）
subplot（236）plot（ta（:2）‘b‘‘linewidth‘2）xlabel（‘t‘）ylabel（‘a2（t）‘）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1361  2016-10-20 11:18  中心差分法\zxcferjie.m
     文件        1453  2015-06-19 09:44  中心差分法\zxcfnjie.m
     目录           0  2016-11-03 16:04  中心差分法\