srungeandkutta2.m

内燃机转子仿真

clear all
clc
close all
%mm=[5.98000 1.08000 1.04000 2.91300 2.91300 2.91300 2.91300 2.91300 51.46300];
%cc=[0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000];
%kk=[0.00000 8.20000 39.20000 15.00000 11.78000 11.78000 11.78000 11.78000 11.78000 16.66000 0.00000];
%F=[10.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000];
%n=1;
m=[2 0;0 1];
c=[0 0;0 0];
k=[6 -2;-2 4];
F=[0;10];
n=length(F);
%m=zeros(n);
%c=zeros(n);
%k=zeros(n);
%for i=1:n
    %m(i,i)=mm(i);
%end
%kk=kk*1.0e5
%for i=1:n
   % k(i,i)=kk(i)+kk(i+1);
    %c(i,i)=cc(i)+cc(i+1);
%end
%for i=1:n-1
          %  k(i,i+1)=-kk(i+1);
           % c(i,i+1)=-cc(i+1);
           % k(i+1,i)=-kk(i+1);
           % c(i+1,i)=-cc(i+1);
%end
S=inv(m)*k;
[px pn]=eig(S);
p=diag(pn);
ff=max(p);
T=2*pi/ff^0.5;
h=T/10000;
t=20;
f=F;
i=1;
dn=rem(t,h);
dn1=t-dn;
n1=dn1/h;
x0=zeros(n,1);
y0=zeros(n,1);
z0=inv(m)*(f-c*y0-k*x0);
xyz(1,i)=i;
xyz(2:n+1,i)=x0;
xyz(n+2:2*n+1,i)=x0;
for i=2:n1
    
%    if i<=n2
        %f=F*(i-1)*dt;
   % else
       % f=F;
    %end
    
   % dk=k1*sin(w*(i-1)*h);
    %k(2,2)=k(2,2)+dk;
    %f=F;
    xz=x0+y0*h/2;
    yz=y0+z0*h/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5
    xz1=x0+yz*h/2;
    yz1=y0+zz*h/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h;
    y11=y0+zz1*h;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h*(y0+2*yz+2*yz1+y11)/6;
    y21=y0+h*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y11*h;
    y1=y0+z11*h;
    if (0.8<abs(y1(1,1))<1.0|0.4<abs(y1(2,1))<1.0)
        x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h1=h/2;
        %i=i-1;
        for j=1:2
            xz=x0+y0*h1/2;%预测中值
    yz=y0+z0*h1/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h1/2;
    yz1=y0+zz*h1/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
  
    end
    x11=x0+yz1*h1;
    y11=y0+zz1*h1;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h1*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h1*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h1;
    y1=y0+z21*h1;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
          elseif (0.6<abs(y1(1,1))<=0.8|0.6<abs(y1(2,1))<=0.8)
              x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h2=h/4;
        for jj=1:4
        xz=x0+y0*h2/2;%预测中值
    yz=y0+z0*h2/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h2/2;
    yz1=y0+zz*h2/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h2;
    y11=y0+zz1*h2;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h2*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h2*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h2;
    y1=y0+z21*h2;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
    elseif (0.4<abs(y1(1,1))<=0.6|0.4<abs(y1(2,1))<=0.6)
        x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h3=h/8;
        for jj=1:8
        xz=x0+y0*h3/2;%预测中值
    yz=y0+z0*h3/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h3/2;
    yz1=y0+zz*h3/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h3;
    y11=y0+zz1*h3;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h3*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h3*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h2;
    y1=y0+z21*h2;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
        elseif (0.2<abs(y1(1,1))<=0.4|0.2<abs(y1(2,1))<=0.4)
            x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h4=h/16;
        for jj=1:16
        xz=x0+y0*h4/2;%预测中值
    yz=y0+z0*h4/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h4/2;
    yz1=y0+zz*h4/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h4;
    y11=y0+zz1*h4;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h4*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h4*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h4;
    y1=y0+z21*h4;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
        elseif (0.1<abs(y1(1,1))<=0.2|0.1<abs(y1(2,1))<=0.2)
             x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h5=h/32;
        for jj=1:32
        xz=x0+y0*h5/2;%预测中值
    yz=y0+z0*h5/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h5/2;
    yz1=y0+zz*h5/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h5;
    y11=y0+zz1*h5;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h5*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h5*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h5;
    y1=y0+z21*h5;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
        elseif (abs(y1(1,1))<=0.1|abs(y1(2,1))<=0.1)
            x0=xyz(2:n+1,i-1);
        y0=xyz(n+2:2*n+1,i-1);
        h6=h/64;
        for jj=1:64
        xz=x0+y0*h6/2;%预测中值
    yz=y0+z0*h6/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*h6/2;
    yz1=y0+zz*h6/2;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*h6;
    y11=y0+zz1*h6;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+h6*(y0+2*yz+2*yz1+y11)/6;        %
    y21=y0+h6*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    x1=x0+y21*h6;
    y1=y0+z21*h6;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
        end
    end
    x0=x1;
    y0=y1;
    z0=z1;
    xyz(1,i)=i;
    xyz(2:n+1,i)=x0;
    xyz(n+2:2*n+1,i)=x0;
    %xyz(2*n+2:3*n+1,i)=z0;
    %zzl(1,i-1)=abs(xyz(2,i)-xyz(2,i-1))/h;
    %zzl(2,i-1)=abs(xyz(3,i)-xyz(3,i-1))/h;
end
figure(1);
plot(xyz(1,:),xyz(2,:));
figure(2);
plot(xyz(1,:),xyz(3,:));
%figure(3);
%plot(xyz(1,2:n1),zzl(1,:));