求助，关于龙格库塔法的应用

zk1985

原帖由 xjzuo 于 2007-4-6 15:32 发表
请将"微分方程组及相应参数值"给出,以便别人帮你调试.

谢谢你啊，主任，不过我那式子太长了，一个微分方程组都有几页长，如果你不嫌麻烦的话，加我qq154566185，我传给你，好吗？

式子太长,可以用附件txt文档上传．

[ 本帖最后由 xjzuo 于 2007-4-6 17:24 编辑 ]

zk1985

我自己做的程序流程。不知道看不看的清楚。

zk1985

xj =

(-7086609486170093341/461168601842738790400000*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936)*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*(-495/2*cos(100*pi*t)*pi^2-658845/2/(2116-121*sin(100*pi*t)^2)^(3/2)*sin(100*pi*t)^2*cos(100*pi*t)^2*pi^2-5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(100*pi*t)^2*pi^2+5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*sin(100*pi*t)^2*pi^2)-347243864822334573709/2305843009213693952000000*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936)*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)-15790128088614291/2305843009213693952000*Fp*sin(Q)*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936)-88930001395075687/1180591620717411303424*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936)*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*Fp*cos(Q)-8499209376874803/18446744073709551616000*sin(Q+4851889062380631/72057594037927936)*Fp-125310456511243/147573952589676412928*sin(Q+4851889062380631/72057594037927936)^2*(-495/2*cos(100*pi*t)*pi^2-658845/2/(2116-121*sin(100*pi*t)^2)^(3/2)*sin(100*pi*t)^2*cos(100*pi*t)^2*pi^2-5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(100*pi*t)^2*pi^2+5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*sin(100*pi*t)^2*pi^2)-6140212369050907/737869762948382064640*sin(Q+4851889062380631/72057594037927936)^2+2579911819940018751/18446744073709551616000*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*(-495/2*cos(100*pi*t)*pi^2-658845/2/(2116-121*sin(100*pi*t)^2)^(3/2)*sin(100*pi*t)^2*cos(100*pi*t)^2*pi^2-5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(100*pi*t)^2*pi^2+5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*sin(100*pi*t)^2*pi^2)+126415679177060918799/92233720368547758080000*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)+4598773297575791/73786976294838206464*Fp*sin(Q)+50586506273333701/73786976294838206464*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*Fp*cos(Q))/(-1876299505410923/147573952589676412928+5174109172077131/9444732965739290427392*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936))


Qj =

-1/(-1876299505410923/147573952589676412928+5174109172077131/9444732965739290427392*sin(Q+4851889062380631/72057594037927936)*cos(Q+4851889062380631/72057594037927936))*(18625220548542048639/9007199254740992000000*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*(-495/2*cos(100*pi*t)*pi^2-658845/2/(2116-121*sin(100*pi*t)^2)^(3/2)*sin(100*pi*t)^2*cos(100*pi*t)^2*pi^2-5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(100*pi*t)^2*pi^2+5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*sin(100*pi*t)^2*pi^2)*cos(Q+4851889062380631/72057594037927936)+912635806878560383311/45035996273704960000000*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(Q+4851889062380631/72057594037927936)+33200036628417199/36028797018963968000*Fp*sin(Q)*cos(Q+4851889062380631/72057594037927936)+365200402912589189/36028797018963968000*sin(100*pi*t)/(2116-121*sin(100*pi*t)^2)^(1/2)*Fp*cos(Q)*cos(Q+4851889062380631/72057594037927936)+51/500000*Fp+54182459777576870783/288230376151711744000000*sin(Q+4851889062380631/72057594037927936)*(-495/2*cos(100*pi*t)*pi^2-658845/2/(2116-121*sin(100*pi*t)^2)^(3/2)*sin(100*pi*t)^2*cos(100*pi*t)^2*pi^2-5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*cos(100*pi*t)^2*pi^2+5445/2/(2116-121*sin(100*pi*t)^2)^(1/2)*sin(100*pi*t)^2*pi^2)+2654940529101266668367/1441151880758558720000000*sin(Q+4851889062380631/72057594037927936))
自由状态下的方程组，也是最简单的，其中Q为转动角度，Qs为转动速度 Qj为转动加速度，X为横向位移，Xs为横向速度，Xj为横向加速度，他们之间呈导数关系。

zk1985

菜鸟虚心向大家请教，不知道有无专业人士指点一二:lol ，难道要我自己去编个龙格库塔法，我看了ode45的愿程序，基本头是昏的，看不大懂，又长篇的。:'(

zk1985

搞忘了那个Fp，它是燃气爆发压力，我用高斯拟合出来的，如下所示
Fp=3800.*exp(-(2305843009213693952./1307182401923243.*t-30120074307853877248./816989001202026875).^2)+1957./5.*exp(-(9223372036854775808./2622204670077813.*t-618542387221573402624./8194389593993165625).^2)+2119.*exp(-(4611686018427387904./8794485237141029.*t-303938931651980034048./27482766366065715625).^2)+3311./5.*exp(-(9223372036854775808./3733621000518813.*t-627189298506124754944./11667565626621290625).^2)+881./2.*exp(-(288230376151711744./202568308359423.*t-19896903153722851328./633025963623196875).^2)+2156./5.*exp(-(1152921504606846976./3774665006082817.*t-20644500691866353664./2359165628801760625).^2)+2093.*exp(-(1152921504606846976./5034808210618101.*t-75732531333862260736./15733775658181565625).^2)-4077./10.*exp(-(2305843009213693952./5806112697200081.*t-126389019942525599744./18144102178750253125).^2)
可以将其代进去。

xjzuo

式子太长,可以用附件txt文档上传．
另外，有几点要注意：
1.希望你能了解本版发贴规则, 问问题一定要先将问题讲清楚，否则一般不予回答或直接删除；
2.你的这个问题虽然有点大，不过，原理就一个，给出"微分方程组及相应参数值",　即使非机械专业的人也可以提出建议，有时间也会帮你调试一下的．
3. 要有心理准备，别人不一定就能帮你解决问题，无论如何（即使没有解决问题）请尊重别人，因为别人都花了时间看你的问题，并且是无偿的．
4.问问题还是请虚心一点较好,就你这口气,我想很多人大概会觉得...

[ 本帖最后由 xjzuo 于 2007-4-6 17:31 编辑 ]

zk1985

我错了:@L

zk1985

可能用词不当，斑竹你误会了，我不是那个意思，只是相互探讨下。