如何计算一个含有矩阵的2次微分方程

猴王

我现在的计算公式是这样的：MX''+KX+CX'=F=200*sin(45t+24),其中M,K,C矩阵的表达式我已经写出，如何计算结构的加速度图形，希望高手能指点

猴王

我上传的附件是我结构的节点图和单元图，我现在的计算是只计算到节点10，也就是只计算1到10节点这十个节点的位移，速度，加速度，但现在我有个疑问，就是现在的K，M,C都是矩阵，收到的简谐力是F=200*sin(45t+24),这样如何编制求解节点速度和位移的程序，就是归根结底就是想，让高手指点一下如何编译执行文件，来求解节点的速度和位移，我如下的编译是不是有问题，按照理论能不能求解
t0=0.001;
te=10;
dt=0.001;%采样间隔
t=[t0:dt:te]';
x0=[0,0];
[x]=ode4('hangjiajuzhen',t,x0);   
d=x(:,1);
v=x(:,2);
s1=(diff(v(1:end-1))+diff(v(2:end)))/2;  %中心差分法
s1=[diff(v(1:2));s1;diff(v(end-1:end))]; %少的两个点用前向和后向差分法补齐
s1=s1/dt; %近似为导数
figure;set(gcf,'Color','W')
subplot(311);plot(t,d); xlabel('T(s)'); ylabel('位移(m)');
subplot(312);plot(t,v); xlabel('T(s)'); ylabel('速度(m/s)');
subplot(313);plot(t,s1); ;xlabel('T(s)'); ylabel('加速度()m/s^2');
y1=s1(150:end);%去掉信号前面一部分的瞬态成份
y2=0.005*randn(size(y1));%添加噪声
s=y1+y2;
%下面进行离散的单尺度小波变换并生成，各尺度上的信号
[c,l]=wavedec(s,5,'bior6.8');%对第一信号进行3尺度一维离散小波分解，采用墨西哥小帽函数
%提取结构的低频和高频信号
ca3=appcoef(c,l,'bior6.8',5);%提取第三尺度系数的低频
[cd1,cd2,cd3,cd4,cd5]=detcoef(c,l,[1,2,3,4,5]);%提取第一、二、三尺度系数的高频
%重构信号的低频和高频部分
a3=wrcoef('a',c,l,'bior6.8',3);
d1=wrcoef('d',c,l,'bior6.8',1);
d2=wrcoef('d',c,l,'bior6.8',2);
d3=wrcoef('d',c,l,'bior6.8',3);
d4=wrcoef('d',c,l,'bior6.8',4);
d5=wrcoef('d',c,l,'bior6.8',5);
%显示多尺度一维信号的分解结果
figure;set(gcf,'Color','W')
subplot(611);plot(t(150:end),a3);ylabel('A5(m/s^2)');xlabel('T(s)');
subplot(612);plot(t(150:end),d1);ylabel('D1(m/s^2)');xlabel('T(s)');
subplot(613);plot(t(150:end),d2);ylabel('D2(m/s^2)');xlabel('T(s)');
subplot(614);plot(t(150:end),d3);ylabel('D3(m/s^2)');xlabel('T(s)');
subplot(615);plot(t(150:end),d4);ylabel('D4(m/s^2)');xlabel('T(s)');
subplot(616);plot(t(150:end),d5); ylabel('D5(m/s^2)');xlabel('T(s)');

猴王

Mx''+Kx+Cx'=f,这是我现在计算的一个动力微分方程，其中M,K,C为10*10的矩阵，在这种情况下，我编制了一个小波分析的程序，但出现了问题，报错，而我不知道错误在什么地方，请高手指点！！！
。m文件：
function xdot=hangjiajuzhen(t,x)
F=20*sin(25*t+14);

  E=210*10^9;
  A=0.25;
  i=E*A/15;
  
  k1=5/3*i;k2=-5/3*i;k3=-5/3*i;k4=5/3*i;k5=-5/4*i;k6=-5/4*i;k7=-5/4*i;k8=-5/4*i;k9=49/25*i;k10=7/25*i;k11=7/25*i;k12=1/25*i;
  k13=1/25*i;k14=7/25*i;k15=7/25*i;k16=1/25*i;
  
  c=0.01;
  c1=(15+83/35)*c;
  c2=(12+83/35)*c;
  c3=(9+83/35)*c;
  c5=9/70*c;
  c6=13/35*c;
  
  m=7850*A;
  m1=(15+83/35)*m;
  m2=(12+83/35)*m;
  m3=(9+83/35)*m;
  m5=9/70*m;
  m6=13/35*m;
  
M=[ m1+m2+m3,m5,m5,m5,0,0,0,0,0,0;
    m5,2*m6,0,m5,0,0,0,0,0,0;
    m5,0,2*m3+m2,m5,m5,0,0,0,0,0;
    m5,m5,m5,2*m1+2*m3+m2,m5,m5,0,0,0,0;
    0,0,m5,m5,2*m1+2*m3+m2,m5,m5,m5,0,0;
    0,0,0,m5,m5,2*m3+m2,0,m5,0,0;
    0,0,0,0,m5,m5,2*m3+m2,m5,m5,0;
    0,0,0,0,m5,m5,m5,2*m1+2*m3+m2,m5,m5;
    0,0,0,0,0,0,m5,m5,m1+m2+m3,m5;
    0,0,0,0,0,0,0,m5,m5,2*m6];
C=[ c1+c2+c3,c5,c5,c5,0,0,0,0,0,0;
    c5,2*c6,0,c5,0,0,0,0,0,0;
    c5,0,2*c3+c2,c5,c5,0,0,0,0,0;
    c5,c5,c5,2*c1+2*c3+c2,c5,c5,0,0,0,0;
    0,0,c5,c5,2*c1+2*c3+c2,c5,c5,c5,0,0;
    0,0,0,c5,c5,2*c3+c2,0,c5,0,0;
    0,0,0,0,c5,c5,2*c3+c2,c5,c5,0;
    0,0,0,0,c5,c5,c5,2*c1+2*c3+c2,c5,c5;
    0,0,0,0,0,0,c5,c5,c1+c2+c3,c5;
    0,0,0,0,0,0,0,c5,c5,2*c6];
K=[k1+k5+k9,k6,k2,k10,0,0,0,0,0,0;
    k6,k8+k4,0,k2,0,0,0,0,0,0;
    k2,0,2*k4+k5,k6,k2,0,0,0,0,0;
    k10,k2,k6,2*k1+k5+k9+k13,k14,k2,0,0,0,0;
    0,0,k2,k14,k13+k5+k9+2*k1,k6,k2,k10,0,0;
    0,0,0,k2,k6,2*k1+k5,0,k2,0,0;
    0,0,0,0,k2,0,2*k1+k5,k6,k2,0;
    0,0,0,0,k10,k2,k6,2*k1+k5+k9+k13,k14,k2;
    0,0,0,0,0,0,k2,k14,k1+k13+k5,k6;
    0,0,0,0,0,0,0,k2,k6,k1+k5];
A=zeros(2*10);
A(1:10,1:10)=zeros(10);
A(1:10,10+1:end)=eye(10);
A(10+1:end,1:10)=-inv(M)*K;
A(10+1:end,10+1:end)=-inv(M)*C;
B=zeros(2*10,1);
B(1:10)=zeros(10,1);
B(10+1:2*10)=inv(M)*F;

xdot=A*x+B;