声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1041|回复: 4

[编程技巧] ode求解时遇到的问题

[复制链接]
发表于 2011-5-26 16:23 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
主程序: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.好像矩阵的维数的问题吧,
回复
分享到:

使用道具 举报

 楼主| 发表于 2011-5-26 16:30 | 显示全部楼层
改成u(1).^3和u(5).^2)之后还是不行,为什么呢?
 楼主| 发表于 2011-5-26 16:33 | 显示全部楼层
以前是两个自由度的方程,就不报错,加了一个自由度之后,就出现了这种情况。
 楼主| 发表于 2011-5-26 17:50 | 显示全部楼层
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))];
改成这个之后这个是可以运行的。但是就是想不通这和矩阵维数投什么关联??
发表于 2011-5-26 20:45 | 显示全部楼层
回复 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))];
复制代码

这样没问题
untitled.jpg

点评

赞成: 5.0
赞成: 5
  发表于 2011-5-27 18:53

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 06:55 , Processed in 0.082079 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表