wilson-θ法中f(t+θ*dt)的问题

sally1982

各位好
       我的毕业设计里面也设计到利用wilson-θ法编程的问题，可是编到现在我开始怀疑自己的程序是否有问题，我把程序上穿一下，希望各位给指点一下，我现在怀疑是我的f(t+θ*dt)用的不当，希望各位帮看看，我实在是想不出哪里有问题
其中的MM，CC，KK都是已知的

1. d2=inv(MM)*(QF-CC*d1-KK*d0);                % d2是加速度的初始值
   
2. 
   
3. %%wilson循环
   
4. for i=1:TTT;
   
5. t(i)=dt*(i-1);
   
6. sit=1.4;
   
7. a0=6/(sit*dt)^2;
   
8. a1=3/(sit*dt);
   
9. a2=2*a1;
   
10. a3=sit*dt/2;
    
11. a4=a0/sit;
    
12. a5=-a2/sit;
    
13. a6=1-3/sit;
    
14. a7=dt/2;
    
15. a8=dt^2/6;
    
16. %%结束
    
17. nd=60*ii.*abs(d1)/Dt/pi;                       %电机在任意时刻的转速      r/min
    
18. 
    
19. NM=Mmax*2*(1+Sm)*Sm*n0.*(n0-nd)./((n0-nd).^2+(n0*Sm)^2+2*Sm^2.*(n0-nd)*n0);
    
20. 
    
21. QF=ii*2.*NM*yt/Dt;
    
22. QF=QF.*RR;
    
23. QDF(i,:)=QF;
    
24. 
    
25. d00=0.0;         %%d00初始位移
    
26. ff=(-AA1.*sign(d1)-AA2).*sign(d1);
    
27. pk=KK+a0*MM+a1*(CC+diag(f.*sign(d1)));
    
28. pk=inv(pk);
    
29. ff=ff+MM*(a0*d0+a2*d1+2*d2)+(CC+diag(f.*sign(d1)))*(a1*d0+2*d1+a3*d2);
    
30. pf1=QF+ff;
    
31. d00=d00+pk*pf1;                                          
    
32. d22=a4*(d00-d0)+a5*d1+a6*d2;                                         %%d22加速度
    
33. d11=d1+a7*(d22+d2);                                                    %%d11速度
    
34. d00=d0+dt*d1+a8*(d22+2*d2);
    
35. 
    
36. d0=d00;                                                             %%d0位移
    
37. d1=d11;                                                     %%d1速度，d2 加速度
    
38. d2=d22;
    
39. S(i,:)=d00;                                                  %%S位移，V速度，A加速度
    
40. V(i,:)=d11;
    
41. A(i,:)=d22/10;
    
42. F(i,:)=ff;
    
43. end
    
44. S
    
45. V
    
46. A
    
47. F

复制代码

sally1982

自己顶一下，
怎么没有人回应呢？
谢谢