找回密码
 注册
Simdroid-非首页
查看: 125|回复: 6

一个解方程的问题

[复制链接]
发表于 2011-10-27 10:54:17 | 显示全部楼层 |阅读模式 来自 北京
好久没有使用MATLAB编程了,今天写了个程序,但是有个很大的困惑,不知道哪里出问题了。
程序很简单,其中16-26行与28-38行表示的意思是一样的,但是使用前面的计算有问题,就是如果用V代替k/(2*pi*f1),就总是出问题。可以使用funavs(0.2,100,0.89,5,0.6)试试结果


function funavs(R,N,T0,A,S)
% R 滚筒半径
% N T0时刻转速,rpm
% A 加速与减速时间比
% S 行程
clc
syms T1 T2
V=2*pi*N*R/60;
T=4*T0;
f1=1/T;
f2=A*f1;
TT2=1/f2;
k=V*2*pi*f1;
k2=V*2*pi*f2;

%  v10=V*(1-cos(2*pi*f1*T1));
%  s10=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1);
%
% v30=v10+V*cos(2*pi*f1*(T/2-T1))-V*cos(2*pi*f1*(T/2));
% s30=s10+ v10*T1+V*cos(2*pi*f1*(T/2-T1))*T1+ V/(2*pi*f1)*sin(2*pi*f1*(T/2-T1))- V/(2*pi*f1)*sin(2*pi*f1*(T/2));
%
% v40=v30+V*cos(2*pi*f2*(TT2/2))-V*cos(2*pi*f2*(TT2/2+T2));
% s40=s30+v30*T2+V*cos(2*pi*f2*(TT2/2))*T2+ V/(2*pi*f2)*sin(2*pi*f2*(TT2/2))- V/(2*pi*f2)*sin(2*pi*f2*(TT2/2+T2));
%
% v60=v40+V*cos(2*pi*f2*(TT2-T2))-V*cos(2*pi*f2*(TT2));
% s6=s40+v40*(T2)+V*cos(2*pi*f2*(TT2-T2))*(T2)+V/(2*pi*f2)*sin(2*pi*f2*(TT2-T2))-V/(2*pi*f2)*sin(2*pi*f2*(TT2));

v10=k/(2*pi*f1)*(1-cos(2*pi*f1*T1));
s10=k*T1/(2*pi*f1)-k*sin(2*pi*f1*T1)/(2*pi*f1).^2;  

v30=v10+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))-k/(2*pi*f1)*cos(2*pi*f1*(T/2));
s30=s10+ v10*(T1)+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))*(T1)+ k/(2*pi*f1).^2*sin(2*pi*f1*(T/2-T1))- k/(2*pi*f1).^2*sin(2*pi*f1*(T/2));

v40=v30+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2+T2));
s40=s30+v30*(2*T1+T2-2*T1)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))*(2*T1+T2-2*T1)+ k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1))- k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1+T2));

v60=v40+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2));
s6=s40+v40*(T2)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))*(T2)+k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2-T2))-k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2));

[T1,T2]=solve(s6-0.6,T1-A*T2)

in=input('是否继续? 1/0' )
if in==1
T00=2*(T1+T2)
in=0;
else
   return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t1=0:0.001:T1;
t3=T1:0.001:2*T1;

a1=k*sin(2*pi*f1*t1);
a3=k*sin(2*pi*f1*(T/2-2*T1+t3));

v1=V*(1-cos(2*pi*f1*t1));
v10=V*(1-cos(2*pi*f1*T1));

s1=V*t1-V*sin(2*pi*f1*t1)/(2*pi*f1);
s10=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1);

v3=v10+V*cos(2*pi*f1*(T/2-2*T1+T1))-V*cos(2*pi*f1*(T/2-2*T1+t3));
v30=v10+V*cos(2*pi*f1*(T/2-2*T1+T1))-V*cos(2*pi*f1*(T/2-2*T1+2*T1));

s3=s10+ v10*(t3-T1)+V*cos(2*pi*f1*(T/2-T1))*(t3-T1)+ V/(2*pi*f1)*sin(2*pi*f1*(T/2-2*T1+T1))- V/(2*pi*f1)*sin(2*pi*f1*(T/2-2*T1+t3));
s30=s10+ v10*(2*T1-T1)+V*cos(2*pi*f1*(T/2-T1))*(2*T1-T1)+ V/(2*pi*f1)*sin(2*pi*f1*(T/2-2*T1+T1))- V/(2*pi*f1)*sin(2*pi*f1*(T/2-2*T1+2*T1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t4=2*T1:0.001:2*T1+T2;
t6=2*T1+T2:0.001:2*T1+2*T2;

a4=k2*sin(2*pi*f2*(TT2/2-2*T1+t4));
a6=k2*sin(2*pi*f2*(TT2-2*T1-2*T2+t6));

v4=v30+V*cos(2*pi*f2*(TT2/2))-V*cos(2*pi*f2*(TT2/2-2*T1+t4));
v40=v30+V*cos(2*pi*f2*(TT2/2))-V*cos(2*pi*f2*(TT2/2+T2));

s4=s30+v30*(t4-2*T1)+V*cos(2*pi*f2*(TT2/2))*(t4-2*T1)+ V/(2*pi*f2)*sin(2*pi*f2*(TT2/2))- V/(2*pi*f2)*sin(2*pi*f2*(TT2/2-2*T1+t4));
s40=s30+v30*T2+V*cos(2*pi*f2*(TT2/2))*T2+ V/(2*pi*f2)*sin(2*pi*f2*(TT2/2))- V/(2*pi*f2)*sin(2*pi*f2*(TT2/2+T2));

v6=v40+V*cos(2*pi*f2*(TT2-T2))-V*cos(2*pi*f2*(TT2-2*T1-2*T2+t6));
s6=s40+v40*(t6-(2*T1+T2))+V*cos(2*pi*f2*(TT2-T2))*(t6-(2*T1+T2))+V/(2*pi*f2)*sin(2*pi*f2*(TT2-T2))-V/(2*pi*f2)*sin(2*pi*f2*(TT2-2*T1-2*T2+t6));

figure(1)
% subplot(311)
plot([t1,t3,t4,t6],[a1,a3,a4,a6],'k','linewidth',1.5)
%axis([0 2*T1+2*T2 -2 2])
xlabel('时间 s', 'fontname','宋体')
ylabel('加速度 m/s^2', 'fontname','宋体')

figure(2)
% subplot(312)
plot([t1,t3,t4,t6],[v1,v3,v4,v6],'k','linewidth',1.5)
% axis([1.5 2 0 0.2])
xlabel('时间 s', 'fontname','宋体')
ylabel('速度 m/s', 'fontname','宋体')

% subplot(313)
figure(3)
plot([t1,t3,t4,t6],[s1,s3,s4,s6],'k','linewidth',1.5)
% axis([1.5 2 0.55 0.6])
xlabel('时间 s', 'fontname','宋体')
ylabel('行程 m', 'fontname','宋体')

tt1=0:0.01:T/2;
aa1=k*sin(2*pi*f1*tt1);
tt2=T/2:0.01:T/2+TT2/2;
aa2=k2*sin(2*pi*f2*(tt2-T/2+TT2/2));
发表于 2011-10-27 13:11:57 | 显示全部楼层 来自 江苏南京
Simdroid开发平台
斑竹可以说说算法,直接看难看懂 ,费时间
回复 不支持

使用道具 举报

 楼主| 发表于 2011-10-27 13:34:33 | 显示全部楼层 来自 北京
与算法无关,其实就是
%  v10=V*(1-cos(2*pi*f1*T1));
%  s10=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1);

v10=k/(2*pi*f1)*(1-cos(2*pi*f1*T1));
s10=k*T1/(2*pi*f1)-k*sin(2*pi*f1*T1)/(2*pi*f1).^2;  
是相等的,但是结果下面正确,上面错误。
定义中V=k/(2*pi*f1)
只是在表达式中分别用V和k/(2*pi*f1)而已。
回复 不支持

使用道具 举报

发表于 2011-10-29 10:30:21 | 显示全部楼层 来自 江苏南京
本帖最后由 scott198510 于 2011-10-29 10:31 编辑
*午夜流星* 发表于 2011-10-27 13:34
与算法无关,其实就是
%  v10=V*(1-cos(2*pi*f1*T1));
%  s10=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1);


%应用斑竹数据验算,完全一致的,不知道斑竹所说猫腻在哪里?

R=0.2;N=100;T0=0.89;A=5;S=0.6
V=2*pi*N*R/60
T1=4*T0
f1=1/T1
f2=A*f1
TT2=1/f2
k=V*2*pi*f1
k2=V*2*pi*f2
V=k/(2*pi*f1)
v10a=V*(1-cos(2*pi*f1*T1))
s10a=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1)
v10=k/(2*pi*f1)*(1-cos(2*pi*f1*T1))
s10=k*T1/(2*pi*f1)-k*sin(2*pi*f1*T1)/(2*pi*f1).^2

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×

评分

1

查看全部评分

回复 不支持

使用道具 举报

 楼主| 发表于 2011-10-30 09:45:39 | 显示全部楼层 来自 北京
是这样的,表达式是等效的,但是解方程结果不同。即下面的结果
function funavs(R,N,T0,A,S)
clc
syms T1 T2
V=2*pi*N*R/60;
T=4*T0;
f1=1/T;
f2=A*f1;
TT2=1/f2;
k=V*2*pi*f1;
k2=V*2*pi*f2;

v10=k/(2*pi*f1)*(1-cos(2*pi*f1*T1));
s10=k*T1/(2*pi*f1)-k*sin(2*pi*f1*T1)/(2*pi*f1).^2;  

v30=v10+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))-k/(2*pi*f1)*cos(2*pi*f1*(T/2));
s30=s10+ v10*(T1)+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))*(T1)+ k/(2*pi*f1).^2*sin(2*pi*f1*(T/2-T1))- k/(2*pi*f1).^2*sin(2*pi*f1*(T/2));

v40=v30+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2+T2));
s40=s30+v30*(2*T1+T2-2*T1)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))*(2*T1+T2-2*T1)+ k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1))- k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1+T2));

v60=v40+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2));
s6=s40+v40*(T2)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))*(T2)+k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2-T2))-k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2));

[T1,T2]=solve(s6-0.6,T1-A*T2)


function funavs(R,N,T0,A,S)
clc
syms T1 T2
V=2*pi*N*R/60;
T=4*T0;
f1=1/T;
f2=A*f1;
TT2=1/f2;
k=V*2*pi*f1;
k2=V*2*pi*f2;

  v10=V*(1-cos(2*pi*f1*T1));
s10=V*T1-V*sin(2*pi*f1*T1)/(2*pi*f1);

v30=v10+V*cos(2*pi*f1*(T/2-T1))-V*cos(2*pi*f1*(T/2));
s30=s10+ v10*T1+V*cos(2*pi*f1*(T/2-T1))*T1+ V/(2*pi*f1)*sin(2*pi*f1*(T/2-T1))- V/(2*pi*f1)*sin(2*pi*f1*(T/2));

v40=v30+V*cos(2*pi*f2*(TT2/2))-V*cos(2*pi*f2*(TT2/2+T2));
s40=s30+v30*T2+V*cos(2*pi*f2*(TT2/2))*T2+ V/(2*pi*f2)*sin(2*pi*f2*(TT2/2))- V/(2*pi*f2)*sin(2*pi*f2*(TT2/2+T2));

v60=v40+V*cos(2*pi*f2*(TT2-T2))-V*cos(2*pi*f2*(TT2));
s6=s40+v40*(T2)+V*cos(2*pi*f2*(TT2-T2))*(T2)+V/(2*pi*f2)*sin(2*pi*f2*(TT2-T2))-V/(2*pi*f2)*sin(2*pi*f2*(TT2));

v10=k/(2*pi*f1)*(1-cos(2*pi*f1*T1));
s10=k*T1/(2*pi*f1)-k*sin(2*pi*f1*T1)/(2*pi*f1).^2;  

v30=v10+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))-k/(2*pi*f1)*cos(2*pi*f1*(T/2));
s30=s10+ v10*(T1)+k/(2*pi*f1)*cos(2*pi*f1*(T/2-T1))*(T1)+ k/(2*pi*f1).^2*sin(2*pi*f1*(T/2-T1))- k/(2*pi*f1).^2*sin(2*pi*f1*(T/2));

v40=v30+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2+T2));
s40=s30+v30*(2*T1+T2-2*T1)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2/2))*(2*T1+T2-2*T1)+ k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1))- k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2/2-2*T1+2*T1+T2));

v60=v40+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))-k2/(2*pi*f2)*cos(2*pi*f2*(TT2));
s6=s40+v40*(T2)+k2/(2*pi*f2)*cos(2*pi*f2*(TT2-T2))*(T2)+k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2-T2))-k2/(2*pi*f2).^2*sin(2*pi*f2*(TT2));

[T1,T2]=solve(s6-0.6,T1-A*T2)
解方程的T1,T2是不同的
前面的程序是正确的,后面的不正确。
回复 不支持

使用道具 举报

发表于 2011-10-30 10:01:30 | 显示全部楼层 来自 湖北武汉
太长了没有看程序,猜想下。

就是如果用V代替k/(2*pi*f1),就总是出问题。
是不是因为pi的问题。pi计算涉及到pi的精度问题。能不能试试给pi赋值,比如pi1=3.14之类。

评分

1

查看全部评分

回复 不支持

使用道具 举报

 楼主| 发表于 2011-10-31 08:21:18 | 显示全部楼层 来自 北京
用上面的程序计算
T1 =

0.43170621059348498232831663296603


T2 =

0.086341242118696996465663326593206
用下面的计算
T1 =

316.85555280757170544536290504911


T2 =

63.371110561514341089072581009821

确实如秦博士建议,将pi=3.1415926535897,方程结果就相同了,都是上面的。

回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 14:49 , Processed in 0.043961 second(s), 16 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

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