pim1df.m

内燃机转子仿真

% Precise Integration Method
clear;
A=zeros(2);%完型矩阵
C=A;
D=[0.5,0;0,1];
B=[-6,2;2,-4];
f0=[0;0;0;10];
f1=zeros(size(f0));%完型矩阵
H=[A,D;B,C];%系数矩阵
I=eye(size(H));%单位矩阵，为雅可比矩阵泰勒展开之主项
iH=inv(H);

tf=20;
step=[2,0.5,0.1]; % different step size
N=20;
figure;
hold;
str=['o','x','b-'];
for jj=1:3%由计算机字长决定极值，泰勒截断阶次
  %PIM begin
  dt=step(jj)/2^N;%将h进一步细分
  Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12)/2;%泰勒展开高阶项
  for iter=1:N%求总Ta
    Ta=2*Ta+Ta*Ta;
  end
  T=I+Ta;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  vk=[0;0;0;0];
  for iter=1:tf/step(jj)
    iter;
    t(:,iter)=step(jj)*(iter-1);%求步长内精细积分点
    v(:,iter)=vk(1);%给定初值
    vk=T*(vk+iH*(f0+iH*f1))-iH*(f0+iH*f1+f1*step(jj));%计算精细积分点处的位移量
    %vk=T*vk+T*iH*f0+T*iH^2*f1-iH*f0-iH^2*f1-iH*f1*step(jj)
    %T*vk~vk  (T-I)*iH*f0     (T-I)*iH^2*f1  -iH*f1*step(jj)
    %(T-I)=Ta
    %
  end
  % PIM end
  % figure(jj);
  plot(t(1:tf/step(jj)),v(1,:),str(jj));
end