zhongxinchafen.m

这是用matlab编写的显示积分方法中心差分法的源程序

clear all;
clc;
close all;
n=3;
m=[1 1 0;1 2 0;0 0 3];
c=[1 0 0;0 0 1;0 1 0];
k=[2 1 0;1 3 0;0 0 4];
f=[2;4;6];
%[M K C]=arryf(m,k,c,n);
w=1;
Bs=inv(m)*k;
[px pn]=eigs(Bs);
wf=pn^0.5;
%wf=2*pi*wf;
%for i=1:3
    %px(:,i)= px(:,i)/ px(1,i);
%end
p=diag(wf);
fmax=max(p);
T=1/fmax;
h=T/1000;
h=pi/400;
t=16;
st=rem(t,h);
tt=t-st;
nn=tt/h;
 hf=sin(w*(i-1)*h);
x10=-diag(eye(3));
x20=-diag(eye(3))*2;
x30=inv(m)*(f*hf-k*x10-c*x20);
a0=1/h^2;
a1=1/2/h;
a2=2*a0;
a3=1/a2;
x0=x10-h*x20+a3*x30;
m1=a0*m+a1*c;
x=zeros(n+1,nn);
x(:,1)=[0;x10];
for i=2:nn
    hf=sin(w*(i-1)*h);
r=f*hf-(k-a2*m)*x10-(a0*m-a1*c)*x0;
x11=inv(m1)*r;
x31=a0*(x0-2*x10+x11);
x21=a1*(-x0+x11);
x(:,i)=[i*h;x11];
x0=x10;
x10=x11;
x20=x21;
x30=x31;
end
save('zhongxinchafen','x');
plot(x(1,:),x(2,:),x(1,:),x(3,:),x(1,:),x(4,:));