二阶三自由度动力方程求解问题
大家看看这个方程怎么用matlab解,其中u,v等都是时间t的函数,要求解出u和时间的关系,数值解也行,最后绘出位移u 与时间的关系曲线。------------------------------------注意-------------------------------------
以后像这种只有一个方程式的,最好贴个图片上来。
根本没必要弄个doc附件,看的时候很不方便
------------------------------------------------------------------------------
[ 本帖最后由 eight 于 2007-10-7 23:04 编辑 ] 希望高手把程序贴出来。 原帖由 HONGRIVER 于 2007-10-6 16:37 发表 http://www.chinavib.com/forum/images/common/back.gif
希望高手把程序贴出来。
建议看看本版各个置顶帖,然后明白:指望别人把程序帮你全部搞定是不现实的
程序
clear alla=0.0;
b=40;
u0=0;
o0=0;
x0=0;
y0=0;
m=1000;
h=(b-a)/m;
T=zeros(1,m+1);
%z=zeros(m+1,length(za));
f=zeros(1,4);
g=zeros(1,4);
p=zeros(1,4);
q=zeros(1,4);
T=a:h:b;
u(1)=u0;
o(1)=o0;
x(1)=x0;
y(1)=y0;
for j=1:m
f1=x(j);
g1=y(j);
p1=2.028*sin(pi*T(j)/2)-42.67*x(j)+74.43*y(j)-17720.89*u(j)+34223*o(j);
q1=0.3405*sin(pi*T(j)/2)-0.016*x(j)-59.66*y(j)-7.45*u(j)-25859*o(j);
f2=(x(j)+h*f1/2);
g2=(y(j)+h*g1/2);
p2=2.028*sin((pi*T(j)+h/2)/2)-42.67*(x(j)+h*f1/2)+74.43*(y(j)+h*g1/2)-17720.89*(u(j)+h*p1/2)+34223*(o(j)+h*q1/2);
q2=0.3405*sin((pi*T(j)+h/2)/2)-0.016*(x(j)+h*f1/2)-59.66*(y(j)+h*g1/2)-7.45*(u(j)+h*p1/2)-25859*(o(j)+h*q1/2);
f3=(x(j)+h*f2/2);
g3=(y(j)+h*g2/2);
p3=2.028*sin((pi*T(j)+h/2)/2)-42.67*(x(j)+h*f2/2)+74.43*(y(j)+h*g2/2)-17720.89*(u(j)+h*p2/2)+34223*(o(j)+h*q2/2);
q3=0.3405*sin((pi*T(j)+h/2)/2)-0.016*(x(j)+h*f2/2)-59.66*(y(j)+h*g2/2)-7.45*(u(j)+h*p2/2)-25859*(o(j)+h*q2/2);
f4=(x(j)+h*f3);
g4=(x(j)+h*g3);
p4=2.028*sin((pi*T(j)+h)/2)-42.67*(x(j)+h*f3)+74.43*(y(j)+h*g3)-17720.89*(u(j)+h*p3)+34223*(o(j)+h*q3);
q4=0.3405*sin((pi*T(j)+h)/2)-0.016*(x(j)+h*f3)-59.66*(y(j)+h*g3)-7.45*(u(j)+h*p3)-25859*(o(j)+h*q3);
u(j+1)=u(j)+(f1+2*f2+2*f3+f4)*h/6;
o(j+1)=o(j)+(g1+2*g2+2*g3+g4)*h/6;
x(j+1)=x(j)+(p1+2*p2+2*p3+p4)*h/6;
y(j+1)=y(j)+(q1+2*q2+2*q3+q4)*h/6;
end
r=
程序写出了,但是老是结果不对,请高手看看毛病在哪。 直接将原问题及参数用word贴一下,以便他人判断.
回复 #5 xjzuo 的帖子
楼主以前是贴过doc格式的,不过里边就一个公式,所以我把那个doc直接换成图片了。程序说明
三自由度的动力方程,化简后可以得到一个只含有u,和转角的二元二阶方程组。上面程序的思路是先降阶,然后编RUNGE-KUTTA程序。如果有更好的解法,欢迎贴上来交流,欢迎批评指正。
[ 本帖最后由 HONGRIVER 于 2007-10-8 11:10 编辑 ] 我眼拙------实在看不出贴出的图片(方程)参数与你代码中的参数有何关系.
另: 结果不对在哪里?你的代码(尤其函数符号)没有任何说明,不说别人,自己读起来也费劲吧。
页:
[1]