rungekuta.m

内燃机转子仿真

clear all
clc
close all
%m=[2.00000 0.00000;0.00000 1.00000];
%c=[0.00000 0.00000;0.00000 0.00000];
%k=[6.00000 -2.00000;-2.00000 4.00000];
%F=[0.00000;10.00000];
m=2;
c=0.0;
k=8;
F=4;
n=1;
%n=length(F);
%w=0;
x=0;
y=0;
z0=inv(m)*F;

w=1;
%dt=60/w/360*2;
h=0.01;
t=10;
%t1=1;
dn=rem(t,h);
dn1=t-dn;
n1=dn1/h;

%k1=0.5;
%dn1=rem(t1,dt);
%dn2=t1-dn;
%n2=dn2/dt;
f=zeros(n,1);
df=zeros(n,1);
dx=zeros(n,1);
dy=zeros(n,1);
d2xx=zeros(n,1);
d2yy=zeros(n,1);
d3x=zeros(n,1);
d3y=zeros(n,1);
xyz1=zeros(2*n+1,n1);
xyz1(2:n+1,1)=x;
xyz1(4:2*n+1,1)=y;
for i=2:n1
    f=F*cos(w*i*h);
    df=0;
    d2f=0;
    dx=y;
    dy=inv(m)*(f-c*y-k*x);
    d2x=dy;
    d2y=inv(m)*(df-c*dy-k*y);
    d3x=d2y;
    d3y=inv(m)*(d2f-c*d2y-k*dy);
    x=x+h*(dx+0.5*h*(d2x+h*d2x/3));
    y=y+h*(dy+0.5*h*(d2y+h*d3y/3));
    %if (0.4<abs(y(1,1))<1.0|0.4<abs(y(2,1))<1.0)
    if 0.4<abs(y(1,1))<1.0
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h1=h/20;
        %i=i-1;
        for j=1:20
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h1*(dx+0.5*h1*(d2x+h1*d2x/3));
            y=y+h1*(dy+0.5*h1*(d2y+h1*d3y/3));
        end
    %elseif (abs(y(1,1))<=0.4|abs(y(2,1))<=0.4)
    elseif abs(y(1,1))<=0.4
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h1=h/40;
        for j=1:40
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h1*(dx+0.5*h1*(d2x+h1*d2x/3));
            y=y+h1*(dy+0.5*h1*(d2y+h1*d3y/3));
        end
    end
    xyz1(1,i)=i;
    xyz1((2:n+1),i)=x;
    xyz1((n+2):(2*n+1),i)=y;
end
save('sjjl1','xyz1');
figure(1);
plot(xyz1(1,:)*h,xyz1(2,:));
%figure(2);
%plot(xyz1(1,:)*h,xyz1(3,:));