请问这种计算程序应如何改正？

likefeng

大家好，一道计算程序如下：
syms M1 M2 A B C F1 F2 f1 f2 w t b w1;
syms s K1 K2 K3 K4  K5 K11 K12 K21 K22;
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
>> T0=2*pi/w;t=[0:T0/500:T0];
>> A=-0.3*w^2*cos(w*t)*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=(0.3*w*sin(w*t))^2*sin(s*t)*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f(1)+(-28510.61)*A*b+(-6619.84)*B*b^2+(-872.82)*C;
F2=f(2)+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*(sin(w*t))^2;
K2=b*((-0.3)*w^(2)*cos(w*t)*sin(s*t)-0.3*w*sin(w*t)*s*cos(s*t));
K3=b*(0.3*w^(2)*cos(w*t)*sin(s*t)-0.3*w*sin(w*t)*s*cos(s*t));
K4=-b*(2)*0.09*(w*sin(w*t))^(2)*(sin(s*t))^2-s^2;
K5=b^(2)*0.09*(w*sin(w*t))^2*(cos(s*t))^2;
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;
错误显示如下：
??? Error using ==> mtimes
Inner matrix dimensions must agree.

应如何改正？谢谢指点。

likefeng

我以解决，原来是在计算时要分清矩阵相乘还是矩阵中的每一元素相乘。正确如下：syms M1 M2 A B C F1 F2 f1 f2 w t b w1;
syms s K1 K2 K3 K4  K5 K11 K12 K21 K22;
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
T0=2*pi/w;t=0:T0/500:T0;
A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;
F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);
K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;
K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;

ChaChing

请教楼主, 你这到底是symbolic还是数值!

likefeng

因方程式中的参数都是时间函数，所以我是这样计算的。

ChaChing

只不过一开始syms M ..., 使得这些参数都为sym obiect
後来又赋值, 使得这些参数又变为double array
说真格, 个人认为是不太好的习惯

likefeng

接受意见，谢谢。能否帮忙看看我用ode45采用上面的参数编的程序，出错了，一时还不知道如何改。
function dx=fz1(t,x,F1,F2,C12,C21,K11,K12,w1,w2)
dx=[x(2);
    0.002337*(F1-C12.*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^2*x(1);
    x(4);
  0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^2*x(3)];





x0=[0;0;0;0];
syms M1 M2 A B C F1 F2 f1 f2 w t b w1;
syms s K1 K2 K3 K4  K5 K11 K12 K21 K22;
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
T0=2*pi/w;t=0:T0/500:T0;
A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;
F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);
K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;
K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;
[t,x]=ode45('fz1',t,x0);
subplot(2,2,1);
plot(t,x(:,1));
hold on;
subplot(2,2,2);
plot(t,x(:,2));
hold on;
subplot(2,2,3);
plot(t,x(:,3));
hold on;
subplot(2,2,4);
plot(t,x(:,4));


Error in ==> fz1 at 2
dx=[x(2);
Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> zfz1 at 26
[t,x]=ode45('fz1',t,x0);

likefeng

问题好像解决了，我还得验证一下,现附上正确的：
function dx=gfz1(t,x,F1,F2,C12,C21,K11,K12,K21,K22,w1,w2)
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;
F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);
K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;
K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;
dx=[x(2);
    0.002337*(F1-C12*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^(2)*x(1);
    x(4);
  0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^(2)*x(3)];



>> T0=2*pi/w;t=0:T0/500:T0;
>>  x0=[0;0;0;0];
>> [t,x]=ode45('gfz1',t,x0);

likefeng

另外想请教大家如何根据上面的方程组得到附件的图，请大家给点提示。谢谢了。