请高手看一下我的程序哪点错了，谢谢
我想用ode45解方程 那位看一下定义的函数那点错了，谢谢，小弟不胜感激%用ode45解微分方程
该微分方程是二阶的分段函数
Fpt为分段的程序中已定义
a_x也是分段的 也已定义
请大哥们看看程序结构上有啥错误
function d2x=BCD(t,x)
%.........................................................................
%........................................................................
pressure=[0.00023 0 44520000 6.6
0.00047 0.003 64770000 16.6
0.00071 0.009 91040000 31.4
0.00095 0.019 123450000 51.3
0.00119 0.033 162550000 77.5
0.00143 0.057 202360000 111.4
0.00167 0.087 244090000 153.2
0.00191 0.13 277710000 201.3
0.00215 0.185 302150000 254.8
0.00239 0.252 316080000 313.3
0.00263 0.334 320180000 327
0.00287 0.431 315830000 430.8
0.00311 0.541 304200000 488.8
0.00335 0.664 290750000 543.9
0.00359 0.801 274550000 596.4
0.00383 0.95 259330000 646
0.00406 1.11 239670000 692.3
0.0043 1.281 214310000 734.3
0.00454 1.461 191910000 772.1
0.00478 1.65 172380000 806
0.00502 1.847 155630000 836.1
0.00526 2.05 140850000 863.8
0.0055 2.26 128050000 888.7
0.00574 2.476 117010000 911.3
0.00598 2.696 107180000 932.2
0.00622 2.922 98870000 951.2
0.00646 3.152 90980000 968.8
0.0067 3.386 83730000 985
0.00694 3.624 77400000 1000];
%........................................................................
%结构参数
m_h=380;
Kr_1=1.27;
Kr_2=2;
rho=1110;
d_p=0.0242;
D_T=0.07;
d_T=0.04;
d_p=0.0242;
d_1=0.32;
A_r0=pi*(D_T^2-d_T^2)/4;
A_rp=pi*d_p^2/4;
A_fj=pi*d_1^2/4;
angle_elv=87*pi/180;
c=19527.497; alpha2=2.46;
f_seal=0.18;
f_track=0.3;
A_1=0.00000001; g=9.8;
F=f_track*m_h*g;
F_T=f_seal*m_h*g*cos(angle_elv); %
F_f0=alpha2*m_h*g;
%.......................................................................
%参数
m_y=1.2;
m_q=2.8;
m_h=380;
phi=1.1729;
phi_1=1.02; %.........................................................................
area_bore=0.00266; beta=1.3; beta_T=0.5276; eta_brake=0.38;
b_time=0.0047;
t_tau= 0.0286; Fpt_g=2.1667e+005;
t_g=0.00694;
=size(pressure);
%........................................................................
%流液孔的面积
if x<=0.035
d_x=0.019;
a_x=pi*(d_p^2-d_x^2)/4;
elseif x<=0.225
d_x=0.019+((21.6-19)/(225-35))*(x-0.035);
a_x=pi*(d_p^2-d_x^2)/4;
elseif x<=0.3
d_x=0.0216+((22.6-21.6)/(300-225))*(x-0.225);
a_x=pi*(d_p^2-d_x^2)/4;
elseif x<=0.34
d_x=0.0226+((23.3-22.6)/(340-300))*(x-0.3);
a_x=pi*(d_p^2-d_x^2)/4;
elseif x<=0.36
d_x=0.0233+((24.2-23.3)/(360-340))*(x-0.34);
a_x=pi*(d_p^2-d_x^2)/4;
else
d_x=0.0242;
a_x=pi*(d_p^2-d_x^2)/4;
end
%......................................................................
%......................................................................
if t<=t_g %
Fpt_i=(1+0.5*(m_y/m_q))*area_bore*pressure(:,3)/phi;
Fpt=lagr_interp12(pressure(:,1),Fpt_i,t);
%........................................................................
elseif t<=t_g + t_tau %
chi=((m_q+beta*m_y)*sqrt(1.0-eta_brake)-(m_q+0.5*m_y)) ...
/((beta-0.5)*m_y);
Fpt=chi*Fpt_g*exp(-(t-t_g)/b_time);
%..........................................................................
else %
Fpt=0.0;
end
F_cR=m_h*9.8*(f_seal+f_track*cos(angle_elv)-sin(angle_elv));%计算后坐阻力常数项
F_phi=(1/2)*(rho*Kr_1*(A_r0-A_rp)^3/(a_x.^2)+rho*Kr_2*(A_fj^3)/(A_1^2));
%..........................................................................
%.........................................................................
% we let x(1)=x and x(2)=x',then x(1)'=x(2)
% and x(2)'=(1/2)*380*(Fpt-(F_phi*(x(2)^2)+F_f0+c*x(1)+F_cR))
d2x=;
[ 本帖最后由 lyy525525 于 2007-5-18 12:12 编辑 ] 原帖由 lyy525525 于 2007-5-18 11:55 发表 http://www.chinavib.com/forum/images/common/back.gif
我想用ode45解方程 那位看一下定义的函数那点错了，谢谢，小弟不胜感激
%用ode45解微分方程
该微分方程是二阶的分段函数
Fpt为分段的程序中已定义
a_x也是分段的 也已定义
请大哥们看看程序结构上有啥 ...
请按照 置顶贴：聚宝盆 的要求把你的问题叙述清楚 建议将主命令及出错信息先给出.
主命令是
主命令是tspan=;
x0=;
=ode45('BCD',tspan,x0)
错误是
错误是一直在运行，一直显示busy该函数能用ode45解吗？
函数描述：二阶的隐性微分方程
微分方程中嵌套着两个分段函数
能解的话，该函数怎么定义，谢谢各位高手
[ 本帖最后由 eight 于 2007-5-18 20:25 编辑 ]
页:
[1]