| 查看: 1618 | 回复: 3 | ||
[求助]
matlab用ode45求解二阶微分方程为题 已有2人参与
|
|
function dq=z_t_equation(t,q,flag,n) switch flag case '' syms x Fn Ff B=0.4e-3; a=11.932; EI=343; v=0.3; l=0.15; m0=0.025; ra=3.276; c=29.8; cof=0.33; d=0.15e-3; omg=2*pi*n./60; gama=pi/6; Fa0=100; kn=1.4e7; kd=(2*kn.*(1-v))./(2-v); %以上为常数参量 omt=sqrt(EI./(ra.*l.^4)); bt=omg./omt; et=l.^3./EI; eb=m0./(ra.*l); ebb=c./(ra.*omt); tao=omg.*t; dba=d./l; Bba=B./l; %以上为无量纲化常量 f=sin((a*x)-sinh(a*x)-((sin(a)+sinh(a)).*(cos(a*x)-cosh(a*x)))/(cos(a)+cosh(a)));%形函数 f11=sin((a)-sinh(a)-((sin(a)+sinh(a)).*(cos(a)-cosh(a)))/(cos(a)+cosh(a))); df1=diff(f,1); df2=diff(f,2); df4=diff(f,4); f1=f.*f; f2=f.*df4; f3=-x.*f.*df1; f4=f.*df2.*(1-x.^2); f5=f.*df2; F0=matlabFunction(f); F1=matlabFunction(f1); F2=matlabFunction(f2); F3=matlabFunction(f3); F4=matlabFunction(f4); F5=matlabFunction(f5); A1=integral(F3,0,1); A2=integral(F4,0,1); A3=integral(F5,0,1); A4=et.*integral(F0,0,1); A5=et.*f11; cta1=acos((kd.*B.*tan(gama)+2*cof.*kn.*d)./(kd.*B.*tan(gama)+2*cof.*kn.*B)); cta2=acos(d/B); cta3=acos((-d.*kd.*tan(gama)+cof.*d.*kn)./(B.*kd.*tan(gama)-cof.*B.*kn)); M=bt.^2.*integral(F1,0,1); C=ebb.*bt.*integral(F1,0,1); K=integral(F2,0,1); (-inf<q(1))&(q(1)<inf); if q(1)>dba/f11; Fn=kn.*l.*(f11.*q(1)-dba).*cos(gama); elseif q(1)<-dba/f11; Fn=kn.*l.*(f11.*q(1)+dba).*cos(gama); else q(1)>=-dba/f11&&q(1)<=dba/f11; Fn=0; end if q(1)<Bba/f11&&q(1)>=Bba.*cos(cta1)/f11&&q(2)<0; Ff=kd.*l.*sin(gama).*(f11.*q(1)-Bba)+cof.*kn.*l.*(f11.*q(1)-dba).*cos(gama); elseif q(1)<Bba.*cos(cta1)/f11&&q(1)>=Bba.*cos(cta2)/f11&&q(2)<0; Ff=-cof.*kn.*l.*(f11.*q(1)-dba).*cos(gama); elseif q(1)<Bba.*cos(cta2)/f11&&q(1)>=-Bba.*cos(cta2)/f11; Ff=0; elseif q(1)<-Bba.*cos(cta2)/f11&&q(1)>=Bba.*cos(cta3)/f11&&q(2)<0; Ff=kd.*l.*sin(gama).*(-f11.*q(1)-dba); elseif q(1)<Bba.*cos(cta3)/f11&&q(1)>=-Bba/f11&&q(2)<0; Ff=-cof.*kn.*l.*cos(gama).*(f11.*q(1)+dba); elseif q(1)<-Bba.*cos(cta1)/f11&&q(1)>-Bba/f11&&q(2)>0; Ff=kd.*l.*sin(gama).*(f11.*q(1)-Bba)-cof.*kn.*l.*(f11.*q(1)+dba).*cos(gama); elseif q(1)<-Bba.*cos(cta2)/f11&&q(1)>=-Bba.*cos(cta1)/f11&&q(2)>0; Ff=cof.*kn.*l.*(f11.*q(1)+dba).*cos(gama); elseif q(1)<-Bba.*cos(cta3)/f11&&q(1)>=Bba.*cos(cta2)/f11&&q(2)>0; Ff=kd.*l.*sin(gama).*(f11.*q(1)-dba); elseif q(1)<=Bba/f11&&q(1)>=-Bba.*cos(cta3)/f11&&q(2)>0; Ff=cof.*kn.*l.*(f11.*q(1)-dba).*cos(gama); end %以上为函数参量的定义 F=A1.*bt.^2.*q(1)+bt.^2.*A2.*q(1)./2+eb.*bt.^2.*A3.*q(1)+A4.*Fa0.*sin(3*tao)-A5.*(Fn.*cos(gama)+Ff.*sin(gama)); dq(1)=q(2); dq(2)=(F-C.*q(2)-K.*q(1))./M; q=[q(1);q(2)]; dq=[dq(1);dq(2)]; otherwise error(['error']); end 运行出错，显示 从 sym 转换为 double 时出现以下错误: DOUBLE cannot convert the input expression into a double array. 我是初学者，还望懂的人帮我调试一下，不胜感激！ |
回复此楼» 猜你喜欢
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有9人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有9人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有13人回复
-
基础研究怎么拉横向，学校到款任务越来越多，难以完成 拉横向，都有哪些途径啊
已经有13人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急
已经有8人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有6人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有9人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有6人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有7人回复
2楼2016-08-31 22:35:02
FMStation
至尊木虫 (知名作家)
- 应助: 591 (博士)
- 贵宾: 0.03
- 金币: 18701.7
- 红花: 97
- 帖子: 8891
- 在线: 1493.2小时
- 虫号: 2400059
- 注册: 2013-04-04
- 专业: 计算机应用技术
【答案】应助回帖
感谢参与，应助指数 +1
|
https://www.mathworks.com/matlab ... /view_thread/302893 DOUBLE cannot convert the input expression into a double array I would suggest putting in a disp(F_h); just before that line and seeing what result you are getting. |
3楼2016-09-01 06:26:07
chendequan
铁虫 (小有名气)
- 应助: 152 (高中生)
- 金币: 918.6
- 红花: 17
- 帖子: 262
- 在线: 208.7小时
- 虫号: 994474
- 注册: 2010-04-11
- 专业: 结构工程
4楼2016-09-01 09:17:09
50