求助！wilson-θ法的MATLAB方法实现

sally1982

希望能帮忙，我的毕业设计里用到这个方法，可是到f(t+θ*dt)这项不知道怎么迭代，还希望高手给帮忙指点，

sally1982

自己顶一下
我参考了基本书，看到书里面对于f(t+θ*dt)得问题，都是把f(t)时刻的值代入了，而f(t+dt)-f(t)这一项好像忽略了，不知道怎么回事？
请大家帮忙

pengweicai

function [ua,ua1,ua2]=wilson(m,k,c,pt,u,u1,u2,n,theta,dertat)
% by swell2005
%---------------
a0=6/((theta*dertat)^2);
a1=3/(theta*dertat);
a2=2*a1;
a3=theta*dertat/2;
a4=a0/theta;
a5=-a2/theta;
a6=1-3/theta;
a7=dertat/2;
a8=dertat^2/6;
k=k+a0*m+a1*c;
[r,s]=ldl(k);
%----------------
for i=0:n,
  pta=pt+m*(a0*u+a2*u1+2*u2)+c*(a1*u+2*u1+a3*u2);
   %---pta=pta+theta*(pta-pt0)+m*(a0*u+a2*u1+2*u2)+c*(a1*u0+2*u1+a3*u2);
  uaa=r\pta;
  uaa=s\uaa;
  uaa=r'\uaa;
  ua2(:,i+1)=a4*(uaa-u)+a5*u1+a6*u2;
  ua1(:,i+1)=u1+a7*(ua2(:,i+1)+u2);
  ua(:,i+1)=u+dertat*u1+a8*(ua2(:,i+1)+2*u2);
  u=ua(:,i+1);
  u1=ua1(:,i+1);
  u2=ua2(:,i+1);
end

仅供参考(我没有运行过，请原谅)，希望楼主自己好好研究，这样才能提高啊。

sally1982

先感谢回帖的高手，我就是不太懂pta=pta+theta*(pta-pt0)+m*(a0*u+a2*u1+2*u2)+c*(a1*u0+2*u1+a3*u2);
中的为什么就能直接可以写成 pta=pt+m*(a0*u+a2*u1+2*u2)+c*(a1*u+2*u1+a3*u2); 那theta*(pta-pt0)这一项怎么为零的，
因为我查过相关的书籍，可能数学上还是不好，对于这项为零不懂，还请在赐教，感激不尽！！！