ode求解时遇到的问题

lihaitao123

主程序：clear all;
clc
global  M rou A L E I k v c1 c2 delta1 g alpha1 miu1 alpha2 beta2 f delta2 alpha3 xi2 xi3 gama3 l k1 m

rou=7.2*10^3;
A=0.36;
E=1.995*10^8;
I=0.6;
M=2.5*10^3;
L=25;
v=20;
g=9.8;
c1=3.14*10^3;
c2=2*10^3;
k=3.4*10^4;
k1=3.14*10^2;
m=1;

omega=(pi*v)/L
alpha2=2*k/(rou*A*L*omega^2)
beta2=2*c2/(rou*A*L*omega^2)
delta1=E*I/(rou*A*omega^2)*(pi/L)^4-alpha2/2
alpha1=pi^4*E/(4*rou*A*omega^2*L^4)
miu1=c1/(rou*A*omega^2)-beta2/2
f =2*M/(rou*A*L*omega^2)*g
delta2= k/M/omega^2                                                                                                                                                                                                                                                   
xi2=c2/(rou*A*omega^2)
alpha3=2*sqrt(2)*k1/(rou*A*L*omega^2)
xi3=0.1
gama3=2*k1/m-sqrt(2)/2*alpha3
l=8
period=2*pi/omega
tspan=0:period/100:1*period ;
  u0=[0.0,0.0,0.0,0.0,0.0,0.0];
%   tspan=0:0.01:1
  figure(1)
[t,u]=ode45('CQ5',tspan,u0,[]) ;
  plot(t,u(:,1))
子函数：
du=[u(2);
       -delta1*u(1)-alpha1*u(1).^3-miu1*u(2)-alpha2*u(3)*sin(omega*t)-alpha2/2*u(1)*cos(2*omega*t)-beta2*u(4)*sin(omega*t)-1/2*beta2*u(2)*cos(2*omega*t)+f*sin(omega*t)+alpha3*u(5)*(1-l/sqrt(alpha^2+u(5)^2));
        u(4);
        -delta2*u(3)-xi2*u(4)+delta2*u(1)*sin(omega*t)+xi2.*u(2)*sin(omega*t);
        u(6);
        gama3*u(5)*(1-l/sqrt(alpha.^2+u(5)^2))-xi3.*u(6)-g+sqrt(2)/2.*(delta1.*u(1)+alpha1*u(1)^3+miu1*u(2)+alpha2*u(3)*sin(omega*t)+alpha2/2*u(1)*cos(2*omega*t)+beta2*u(4)*sin(omega*t)+1/2*beta2*u(2)*cos(2*omega*t)-f*sin(omega*t))];

报错是Error using ==> mrdivide
Matrix dimensions must agree.好像矩阵的维数的问题吧，

lihaitao123

改成u(1).^3和u(5).^2)之后还是不行，为什么呢？

lihaitao123

以前是两个自由度的方程，就不报错，加了一个自由度之后，就出现了这种情况。

lihaitao123

clear all;
clc
global  M rou A L E I k v c1 c2 delta1 g alpha1 miu1 alpha2 beta2 f delta2 alpha3 xi2 xi3 gama3 l k1 m
rou=7.2*10^3;
A=0.36;
E=1.995*10^8;
I=0.6;
M=2.5*10^3;
L=25;
v=20;
g=9.8;
c1=3.14*10^3;
c2=2*10^3;
k=3.4*10^4;
k1=3.14*10^2;
m=1;
omega=(pi*v)/L;
alpha2=2*k/(rou*A*L*omega^2);
beta2=2*c2/(rou*A*L*omega^2);
delta1=E*I/(rou*A*omega^2)*(pi/L)^4-alpha2/2;
alpha1=pi^4*E/(4*rou*A*omega^2*L^4);
miu1=c1/(rou*A*omega^2)-beta2/2;
f =2*M/(rou*A*L*omega^2)*g;
delta2= k/M/omega^2 ;                                                                                                                                                                                                                                                  
xi2=c2/(rou*A*omega^2);
alpha3=2*sqrt(2)*k1/(rou*A*L*omega^2);
xi3=0.1;
gama3=2*k1/m-sqrt(2)/2*alpha3;
l=8;
period=2*pi/omega;
tspan=0:period/100:50*period
  u0=[0,0,0,0,0,0];
%   tspan=0:0.01:1
  figure(1)
[t,u]=ode45('CQ6',tspan,u0)
  plot(t,u(:,1))
   hold on
     plot(t,u(:,3),'r-')
     legend('vehicle','bridge'
子函数：function du=CQ6(t,U)
global    delta1 delta2 alpha1 miu1 alpha2 xi2  beta2  omega f  gama3 alpha xi3 alpha3  l g

omega =2.5133;
alpha2=0.1661;
beta2 =0.0098;
delta1 =1.7401;
alpha1 =0.7596;
miu1 =0.1869;
f =0.1197;
delta2 =2.1531;
xi2 = 0.1222;
alpha3 =0.0022;
xi3 =0.1000;
gama3 =627.9985;
l =8;
alpha=0.1;
g=9.8;
    u=U(1);
    x=U(2);
    y=U(3);
    z=U(4);
    p=U(5);
    q=U(6);
   
    du=zeros(6,1);
   
      du=[x;
        -delta1*u-alpha1*u^3-miu1*x-alpha2*y*sin(omega*t)-alpha2/2*u*cos(2*omega*t)-beta2*z*sin(omega*t)-1/2*beta2*x*cos(2*omega*t)+f*sin(omega*t)+alpha3*p*(1-l/sqrt(alpha^2+p^2));
        z;
        -delta2*y-xi2*z+delta2*u*sin(omega*t)+xi2*x*sin(omega*t);
        q;   
       gama3*p*(1-l/sqrt(alpha^2+p^2))-xi3*q-g+sqrt(2)/2*(delta1*u+alpha1*u^3+miu1*x+alpha2*y*sin(omega*t)+alpha2/2*u*cos(2*omega*t)+beta2*z*sin(omega*t)+1/2*beta2*x*cos(2*omega*t)-f*sin(omega*t))];
改成这个之后这个是可以运行的。但是就是想不通这和矩阵维数投什么关联？？

meiyongyuandeze

回复 1 # lihaitao123 的帖子

1. 子函数中的oemga也应该是全局变量。
2. 子函数中的语句
-delta1*u(1)-alpha1*u(1)^3-miu1*u(2)-alpha2*u(3)*sin(omega*t)-alpha2/2*u(1)*cos(2*omega*t)-beta2*u(4)*sin(omega*t)-1/2*beta2*u(2)*cos(2*omega*t)+f*sin(omega*t)+alpha3*u(5)*(1-l/sqrt(alpha^2+u(5)^2));中的alpha没有定义。
我自己随便付了几个值，为了时调程序。

1. function du=CQ5(t,u)
   
2. global  M rou A L E I k v c1 c2 delta1 g alpha1 miu1 alpha2 beta2 f delta2 alpha3 xi2 xi3 gama3 l k1 m omega
   
3. du=[u(2);
   
4. -delta1*u(1)-alpha1*u(1)^3-miu1*u(2)-alpha2*u(3)*sin(omega*t)-alpha2/2*u(1)*cos(2*omega*t)-beta2*u(4)*sin(omega*t)-1/2*beta2*u(2)*cos(2*omega*t)+f*sin(omega*t)+alpha3*u(5)*(1-l/sqrt(alpha2^2+u(5)^2));
   
5. u(4);
   
6. -delta2*u(3)-xi2*u(4)+delta2*u(1)*sin(omega*t)+xi2*u(2)*sin(omega*t);
   
7. u(6);
   
8. gama3*u(5)*(1-l/sqrt(alpha2^2+u(5)^2))-xi3*u(6)-g+sqrt(2)/2.*(delta1.*u(1)+alpha1*u(1)^3+miu1*u(2)+alpha2*u(3)*sin(omega*t)+alpha2/2*u(1)*cos(2*omega*t)+beta2*u(4)*sin(omega*t)+1/2*beta2*u(2)*cos(2*omega*t)-f*sin(omega*t))];

复制代码

这样没问题