二阶三自由度动力方程求解问题

HONGRIVER

大家看看这个方程怎么用matlab解，其中u,v等都是时间t的函数，要求解出u和时间的关系，数值解也行，最后绘出位移u 与时间的关系曲线。
------------------------------------注意-------------------------------------
以后像这种只有一个方程式的，最好贴个图片上来。
根本没必要弄个doc附件，看的时候很不方便
------------------------------------------------------------------------------

[ 本帖最后由 eight 于 2007-10-7 23:04 编辑 ]

HONGRIVER

希望高手把程序贴出来。

eight

原帖由 HONGRIVER 于 2007-10-6 16:37 发表

                               
登录/注册后可看大图

希望高手把程序贴出来。

建议看看本版各个置顶帖，然后明白：指望别人把程序帮你全部搞定是不现实的

HONGRIVER

clear all
a=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=[T',u',o',x',y']

程序写出了，但是老是结果不对，请高手看看毛病在哪。

xjzuo

直接将原问题及参数用word贴一下，以便他人判断．

花如月

楼主以前是贴过doc格式的，不过里边就一个公式，所以我把那个doc直接换成图片了。

HONGRIVER

三自由度的动力方程，化简后可以得到一个只含有u,和转角的二元二阶方程组。上面程序的思路是先降阶，然后编RUNGE-KUTTA程序。
如果有更好的解法，欢迎贴上来交流，欢迎批评指正。

[ 本帖最后由 HONGRIVER 于 2007-10-8 11:10 编辑 ]

xjzuo

我眼拙------实在看不出贴出的图片(方程)参数与你代码中的参数有何关系.
另: 结果不对在哪里？你的代码(尤其函数符号)没有任何说明，不说别人，自己读起来也费劲吧。