马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
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;
%f-作用力的向量,%n为系统的自由度
执行文件如下:
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)');
错误提示
??? Error using ==> ode4
Unable to evaluate the ODEFUN at t0,y0. In an assignment A(I) = B, the number of elements in B and
I must be the same
[ 本帖最后由 eight 于 2007-11-30 17:27 编辑 ] |