数值积分问题 simukink 和ode45的不同

ljq2008

我最近在进行动力学系统分析的时候，建立了一个一个实际物理对象的动力学方程，我用了两种仿真建模方法，一种是用simulink的S函数，搭建了基本的模块，一种是直接利用ode45进行数值积分，在同样的初始条件情况下，建模的基本部分完全一样，代码基本都是复制来去的，但是对同一个变量的结果却相差较大，我下面将我的一些代码贴出来，希望高手给我指正一下，
plantmodel0831.m这个文件是用S函数写的建模对象的动力学方程，初始条件也已经包含
basicmodel0831.mdl这个是基于S函数plantmodel0831.m搭建的simulink模型，只是为了将状态变量输出而已
plotbasicmoel0808.m这个函数将我关心的变量变化值图形显示
这三个文件都属于simulink的建模方法。
modelequation。m这个函数是我直接写的对象的模型，将来为ode45函数的调用
qiujiefangc.m这个函数是进行初始条件的设定和ode45函数的调用
为了方便大家直接运行并找出问题，我把我写的文件都放在附件里，同时也贴出具体代码
%%%%%%%%%%%%%%%%%%
function y=modelequation(t,x)
Msitasita = 0.256447668405926;
Mqq =[ 0.370100442297874   0.000008792825265; 0.000008792825265   0.303218779223962];
Msitaq =[ 0.259306820984583  -0.010542738101024];
Msitaata=0.256422679228154;
MsitaqT =[0.259306820984583; -0.010542738101024];
MataqT =[0.259306820984583; -0.010542738101024];
MsitaataT =0.256422679228154;
Mataq =[ 0.259306820984583  -0.010542738101024];
Mataata =0.256415026391973;
K =1.0e+003*[      0                   0                   0                     0
                   0         0.020844744563957        0.000271728078454          0
                   0         0.000271728078454        2.121545658107192          0
                   0                   0                   0                 1.829240205779346];
            
sita=x(1);  sitadot=x(5);  q=[x(2);x(3)];  qT=[x(2) x(3)];  qdot=[x(6);x(7)];
qdotT=[x(6) x(7)]; ata=x(4); atadot=x(8);
M=[ Msitasita+qT*Mqq*q           Msitaq        Msitaata+qT*Mqq*q
       MsitaqT                   Mqq                 MataqT
MsitaataT+qT*Mqq*q      Mataq        Mataata+qT*Mqq*q];

C=[   qdotT*Mqq*q                           sitadot*qT*Mqq                               0  
     -Mqq*q*sitadot                            zeros(2)                     -atadot*Mqq*q-2*sitadot*Mqq*q
   qdotT*Mqq*q                       sitadot*qT*Mqq+atadot*qT*Mqq               qdotT*Mqq*q ];      
X=[sita;q;ata];
Xdot=[sitadot;qdot;atadot];  
y=[zeros(4) eye(4);-M\K -M\C]*[X;Xdot];

同一个变量的两种方法的输出结果

说明：我只贴出了自己写的用于ode45函数运算的函数体，其余部分由于空间有限，没有贴出来，都在附件里，
我也试过ode113 ode15s 函数，好像效果都不好 希望高手能够帮我解答下，万分感谢！

ljq2008

哎，希望高手指点啊！

dfg25aw

可能算法不稳定 可以使用自己做的算法  比较好用有对角隐式RK方法