【谷速软件】matlab源码-精细积分法.docVIP

function =precise(m,k,g,r,x,x1,N,step) % Precise Integration Method A=-inv(m)*g/2; C=g*inv(m)/2; D=inv(m); B=g*inv(m)*g/4-k; f0=[zeros(size(r));r]; f1=zeros(size(f0)); H=[A,D;B,C]; I=eye(size(H)); iH=inv(H); p=m*x1+g*x/2; vk=[x;p]; %--------------------------------- tf=N; dt=step/2^N; Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12)/2; for iter=1:N ???Ta=2*Ta+Ta*Ta; end T=I+Ta; for iter=1:tf/step ???vk=T*(vk+iH*(f0+iH*f1))-iH*(f0+iH*f1+f1*step); ???u(:,iter)=vk;%(1:2,; end %u=u(1:2,1:N+1);