ode求解时遇到的问题
主程序：clear all;clc
globalM 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=;
% tspan=0:0.01:1
figure(1)
=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.好像矩阵的维数的问题吧， 改成u(1).^3和u(5).^2)之后还是不行，为什么呢？ 以前是两个自由度的方程，就不报错，加了一个自由度之后，就出现了这种情况。 clear all;
clc
globalM 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=;
% tspan=0:0.01:1
figure(1)
=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 xi2beta2omega fgama3 alpha xi3 alpha3l 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))];
改成这个之后这个是可以运行的。但是就是想不通这和矩阵维数投什么关联？？
回复 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没有定义。
我自己随便付了几个值，为了时调程序。
function du=CQ5(t,u)
globalM 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
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(alpha2^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(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))];
这样没问题
页:
[1]